**1. Introduction**

Land surface temperature (LST) and land surface emissivity (LSE) are two key physical parameters characterizing the state of the land surface, which are applied in various areas such as mineral mapping [1–3], gas plume detection [4], soil moisture inversion [5] and oil-film thicknesses measurement [6]. The LST is sensitive to the external environment and reflects the energy budget of an object during a period; the LSE, determined by the object itself, reflects its physical and chemical properties. The LST and LSE of an object can be simultaneously derived from thermal infrared (generally 8–14 µm) data through a temperature-emissivity separation (TES) algorithm. TES is a crucial issue in thermal infrared remote sensing, which can be seen as an underdetermined equation problem, i.e., obtaining *N*+1 unknowns (*N* emissivities at each band and one temperature) by solving *N* equations (one sensor output radiance with *N* bands). Consequently, it needs at least one additional constraint to make the underdetermined equations solvable.

For thermal infrared data with low or medium spectral resolution, LST is the principal retrieval goal compared with LSE. In this case, the primary task of TES algorithms for multiband thermal infrared data is to determine LSE. Many algorithms were developed for the inversion of LST from multispectral thermal infrared data, such as the split-window method [7], grey body emissivity [8], reference channel method [9], normalized emissivity method (NEM) [10,11], temperature independent spectral indices (TISI) [12], spectral ratio method [13], NDVI-based emissivity method (NBEM) [14–17] and advanced spaceborne thermal emission and reflection radiometer temperature-emissivity separation (ASTER TES) [18]. Among them, NBEM, TISI and ASTER TES are widely used methods [19]. Based on the relationship between NDVI and LSE, NBEM determines the LSE at specific bands with the help of visible near-infrared data. The ASTER TES method iteratively adjusts the average emissivity according to the empirical relationship between the change range of emissivity and the average emissivity to obtain the optimal LST. TISI transforms ground-leaving radiance under some approximation of the Planck function to make it reflect the spectral shape of emissivity. The day and night algorithm, based on TISI, has obtained good results on LST inversion of multi-band thermal infrared remote-sensing data [17]. TISI is theoretically suitable for thermal infrared remote sensing data of any spectral resolution; however, it cannot obtain absolute values of the retrieval emissivity spectrum. Therefore, TISI requires other data or theoretical constraints to solve the problem. The improvement of hyperspectral thermal infrared (HTIR) remote sensing technology brings new opportunities for the development of thermal infrared remote sensing. Unlike multi-spectral remote sensing, HTIR remote sensing is applied to retrieve not only LST but also LSE more accurately, which gives full play to the role of LSE on the fields of classification, identification and parameter inversion of the land surface (especially urban land surfaces of various types). As a result, the application scenarios of HTIR remote sensing are different to those of multispectral thermal infrared remote sensing. To adapt to new application scenarios, the goal of TES of HTIR remote sensing is to obtain the LST and LSE with the highest accuracy under the smallest number of universal and reasonable additional constraints. At present, scholars have proposed or improved a number of TES algorithms for HTIR data [20–35]. Unlike multispectral thermal infrared inversion methods, most TES algorithms for HTIR data have far fewer constraints. One well-known assumption of HTIR TES algorithms is that the LSE of most objects [36] is much smoother than the spectral curve of atmospheric downwelling radiance.

Iterative spectrally smooth temperature-emissivity separation (ISSTES) [21] is a representative TES algorithm based on the concept of spectral smoothness for HTIR data. ISSTES uses a smoothness index as a cost function measuring the atmospheric residue in the estimated LSE, and iteratively optimizes the estimated LST until the cost function is minimized and the corresponding LSE is considered optimally smoothed. The automatic retrieval of temperature and emissivity using spectral smoothness (ARTEMISS) improved ISSTES [22] by updating its cost function. ARTEMISS uses the standard deviation of the estimated at-sensor radiance and true at-sensor radiance as the cost function, instead of the absolute deviation of the estimated LSE and true LSE. Based on ARTEMISS, the quick temperature-emissivity separation (QTES) algorithm [23] was modified to retrieve LST and LSE from HTIR data measured near the ground. QTES uses a narrow spectral range (9.5–10 µm) instead of the entire spectral range of data when searching for the optimal temperature. In addition to the three algorithms proposed by Borel mentioned above, researchers have successively proposed many other TES algorithms based on the concept of spectral smoothness. The spectral smoothness method (SpSm) [24] uses the first derivative of estimated LSE spectrum as cost function. The downwelling radiance residual index (DRRI) [25] algorithm uses downwelling radiance residual index to measure roughness of LSE curve on several selected three-channel groups. The correlation-based temperature and emissivity separation (CBTES) method [26] constructs a cost function to measure the correlation between the estimated LSE and the atmospheric downwelling radiation, based on the concept that LSE is not directly related to atmospheric downwelling radiation, but when estimated LST is not accurate, estimated LSE will include residual atmospheric absorption information. The stepwise refining temperature and emissivity separation (SRTES) [27] considers self-emitting radiation of the land surface a function of wavenumber in a narrow spectral range, where LSE can be seen as a constant and the Planck function can be expressed in a linear form. However, atmospheric downwelling radiation still contains unignorable atmospheric absorption features. SRTES utilizes the residual atmospheric downwelling radiation of calculated self-emitting radiation as a criterion and obtains both the LSE at specified bands and the LST by the stepwise refining method. The correlation-wavelet method [28] considers the land surface self-emission curve as a fundamental-frequency signal and atmospheric downwelling radiance as a high-frequency signal, and obtains the land surface self-emission from ground-leaving radiance by filtering out the atmospheric characteristics based on the spectral smoothness concept. The correlation-wavelet method calculates the correlation between the emissivity curves calculated by taking into the atmospheric downwelling radiance and by direct wavelet filtering without considering atmospheric downwelling radiance. At the maximum correlation, the corresponding temperature is seen as the optimal retrieval temperature. Meanwhile, the emissivity curve is synthesized by wavelet signals of different scales whose proportion is calculated based on the correlation. A TES algorithm for low-emissivity materials [29] is based on the bias characteristic of atmospheric downwelling radiance, namely when retrieval LSE deviates from its true value, atmospheric downwelling radiance also deviates from its true value and the deviations are the same at the bands of the atmospheric absorption peak and valley. Some researchers have tried to reduce the unknowns of equations by reducing the dimensions of the emissivity, making the underdetermined equations solvable. One way to reduce the dimensionality is to divide the emissivity curve into several segments, each of which is the linear function of the wavelength. Based on the segmented linear constraints, three TES methods have been proposed, namely, the linear spectral emissivity constraint (LSEC) [30], the improved LSEC (I-LSEC) [31] and the pre-estimate shape LSEC (PES-LSEC) [32], which all take the sum of squared residuals of the estimated and true ground-leaving radiance as cost function. Another way to reduce the dimensionality is using the wavelet transform. The wavelet transform method for separating temperature and emissivity (WTTES) [33], which expresses LSE as a function of low-frequency wavelet coefficients, reduces the number of unknowns from N + 1 to N/2 + 1, and iteratively searches for the optimal wavelet coefficient and LST. The multi-scale wavelet-based TES algorithm (MSWTES) [34] is based on the fact that both high frequencies of ground-leaving radiance and the LSE calculated according to inaccurate LST are closely correlated with the atmospheric downwelling radiance. From the perspective of the spectral curve, the above TES algorithms follow similar basic concepts. For eliminating the influence of the "thorns" of a curve due to atmospheric absorption lines, the TES algorithms first aim to obtain a degraded approximate spectral curve close to the true spectral curve. For example, ISSTES estimates the LSE from equations directly; LSEC and WTTES estimate LSE after modifying its expression. Then, they search for the optimal LST, which minimizes the cost function and calculates its corresponding LSE as the final retrieval LSE. In addition to the idea of spectral smoothness, there is a new idea from the perspective of statistics, which assumes that the LSE spectra of natural and man-made materials can be well represented in a given subspace of the original data space. Based on the new idea, the dictionary subspace based temperature and emissivity separation (D-SBTES) [35,36] uses a singular value decomposition to extract the basis matrix of the subspace from the emissivity spectra dictionary to obtain the retrieval emissivity. However, D-SBTES suffers from several factors, such as noise, the rank of basis matrix adopted to address the emissivity subspace and the true land surface temperature, and it is more suitable for high-emissivity objects.

Overall, spectral smoothness has been a reasonable strategy widely used in the field of TES. Among the TES algorithms mentioned above, the most commonly used cost function is the standard deviation of the simulated at-sensor radiance and true at-sensor radiance, and a classic and common method to approximately process the LSE curve is the boxcar average. Therefore, for this study we researched ARTEMISS, which is the representative algorithm of the spectral smoothness TES family and has first-rate performance.

There are many impact factors to the performance of a TES algorithm for HTIR data, mainly including (1) sensor-related parameters (e.g., sensor altitude, spatial resolution, spectral range, spectral resolution and instrument noise), (2) data pre-processing (before TES) related residual errors (e.g., radiation calibration, spectral calibration and atmospheric correction) and (3) the algorithm's limitations (e.g., the adopted assumptions or constraints) [37–48]. The sensor altitude affects the amount of atmospheric radiation entering the sensor. The higher the sensor altitude, the greater the inversion error [37,44]. If the spatial resolution is low, non-isothermal mixed pixels will appear in the thermal infrared hyperspectral image, which makes it difficult to separate the temperature emissivity accurately [39,41,42]. The spectral range of 7.5–8 µm, containing dense atmospheric absorption lines, is usually used to retrieve atmospheric parameters; however, the inversion requires a sufficient spectral resolution. The spectral range of 8–12.5 µm is usually used to retrieve the LST and LSE; the higher the spectral resolution, the smaller the retrieval errors [42,43]. A sufficient signal-to-noise ratio is necessary to obtain a unique solution for TES, and the retrieval error increases linearly with instrument noise [37,40,42]. The radiation calibration needs to have sufficient accuracy in order to make the simulated and measured radiance match [42]. The spectral calibration error greatly affects the temperature and emissivity errors even if the atmospheric correction has no error [40,45–48]. The influence of atmospheric correction errors (atmospheric upwelling radiance, atmospheric downwelling radiance and atmospheric transmittance) on the retrieval results varies by the TES algorithm. Among the three atmospheric parameters, atmospheric downwelling radiance has most of the influence on retrieval errors, especially for low-emissivity objects [41,43,45,46]. The adopted cost function is a mathematical expression of the adopted assumption for a TES algorithm, which certainly affects algorithm performance, regardless of whether the algorithm solves the underdetermined problem by increasing the constraints (e.g., spectral smoothness assumption) or reducing unknowns (e.g., piecewise linear constraint) [37,38,44,47,48].

The above factors, excluding the algorithm's limitations, can be classified as systematic errors and random errors. The systematic errors can be reduced or removed by some correction methods based on the study of error patterns. For example, besides laboratory calibration, scene-based spectral calibration is able to improve further the spectral calibration accuracy of the field-measured data, which selects some atmospheric absorption channels as the reference channel to perform spectral calibration on the measured data, and the minimum error of the center wavelength shift is within 1 nm [49]. The sensor altitude, spectral range and resolution depend on the equipment development level and the external conditions during data acquisition. They mainly affect errors of the atmospheric parameters (input parameters for separating temperature and emissivity), and the accuracy of atmospheric correction increases with the improvement of technology [40,50,51]. The random errors, however, are much more difficult to deal with. Like noise, it is uncontrollable and limited by the detector technology, and it is difficult to remove the noise non-destructively after data acquisition. Therefore, the influence of random errors on the TES algorithms is also widely researched [25–28,30–32,37,40–44,47,48]. However, most research on this issue has been about how much noise causes the temperature and emissivity inversion errors, and only a few studies have focused on the influence mechanism of noise on the TES algorithms.

Here, we selected the ARTEMISS algorithm—a representative of the iterative spectral smoothness TES algorithm family with a good application effect—as the research object. The influence mechanism of noise on the retrieval errors of ARTEMISS is derived from the cost function. The relationships between the instrument spectral resolution, noise level, the ARTEMISS parameter setting and the retrieval errors are investigated through simulation experiment. On the basis of the mechanism and law of the influence of noise on the retrieval errors obtained from the experiment, the resolution-degrade-based spectral smoothness (RDSS) algorithm, an improved method, is proposed. Corresponding suggestions on the instrument design are also provided. Section 1 reviews TES algorithms and impact factors. Section 2 theoretically analyzes how noise affects the retrieval results of ARTEMISS. Section 3 describes the data simulation experiment, implementation of TES and the evaluation metrics for the retrieval

results (Section 3.1) and the results of a sensitivity analysis of ARTEMISS regarding noise and smoothing window size (Section 3.2). Section 4 proposes the RDSS algorithm based on ARTEMISS, and presents the validation results. Sections 5 and 6 present the discussion and conclusions, respectively.

#### **2. Background**

The temperature-emissivity separation is based on radiative transfer theory. Ignoring multiple scattering, the radiative transfer equation (RTE) [9,52,53] is as follows:

$$L(\lambda\_{i\prime}T) = L\_{\mathcal{S}}(\lambda\_{i\prime}T)\tau(\lambda\_{i\prime}) + L\_u(\lambda\_i) \tag{1}$$

$$L\_{\mathcal{S}}(\lambda\_i, T) = \varepsilon(\lambda\_i) B(\lambda\_i, T) + [1 - \varepsilon(\lambda\_i)] L\_d(\lambda\_i) \tag{2}$$

$$B(\lambda\_{l\prime}T) = \frac{c\_1}{\lambda\_l^5 \cdot \exp\left(\frac{c\_2}{\lambda\_l T} - 1\right)}\tag{3}$$

where *L*(λ*<sup>i</sup>* , *T*) is the at-sensor radiance at *i*-th band and LST of *T*, *i* ∈ [1, *N*]; *N* is the number of bands; *Lg*(λ*<sup>i</sup>* , *T*) is the ground-leaving radiance, including the object's self-emitting radiance and reflected atmospheric downwelling radiance; *Lu*(λ*i*) is the atmospheric upwelling radiance; τ(λ*i*) is the atmospheric transmittance; ε(λ*i*) is the land surface emissivity; *B*(λ*<sup>i</sup>* , *T*) is the blackbody radiance at *T* and *L<sup>d</sup>* (λ*i*) is the atmospheric downwelling radiance.

The original ARTEMISS algorithm includes two parts: atmospheric correction and temperatureemissivity separation. This study focuses on the influence of instrument noise on the retrieval LST and LSE of ARTEMISS. Therefore, this study only discusses the temperature-emissivity separation process of ARTEMISS, assuming the three atmospheric parameters (τ(λ*i*), *L<sup>d</sup>* (λ*i*) and *Lu*(λ*i*)) known and without errors. Then, we know *Lg*(λ*<sup>i</sup>* , *T*) according to Equation (1). Given temperature estimation *T*ˆ, we can obtain the LSE according to Equation (2). The emissivity is calculated by,

$$\mathcal{E}(\lambda\_i) = \frac{L\_{\mathcal{g}}(\lambda\_{i\prime}T) - L\_d(\lambda\_i)}{B(\lambda\_{i\prime}\hat{T}) - L\_d(\lambda\_i)}\tag{4}$$

where εˆ(λ*i*) is the LSE estimation and *B* λ*i* , *T*ˆ is the blackbody radiance at *T*ˆ.

According to the concept of spectral smoothness, the goal of ARTEMISS is to find an optimal LST estimation where the corresponding LSE is the smoothest. ARTEMISS uses the standard deviation of estimated at-sensor radiance and true at-sensor radiance as the cost function. When the cost function reaches the minimum, the corresponding LST estimation is the optimal LST. The process of ARTEMISS temperature-emissivity separation is as follows:

Given an LST estimation, the LSE estimation can be calculated using Equation (4); then, a boxcar average is performed on it and one can obtain the smoothed LSE estimation εˆ(λ*i*):

$$\overline{\mathfrak{E}}(\lambda\_i) = \frac{1}{\mathfrak{Z}} \sum\_{j=i-1}^{i+1} \mathfrak{E}(\lambda\_j) \tag{5}$$

Then, the estimated at-sensor radiance, *Lfit*(λ, *T*ˆ, εˆ) is calculated according to the RTE:

$$L\_{\text{fit}}(\lambda, \Upsilon, \overline{\varepsilon}) = \varepsilon(\lambda) \mathcal{B}(\lambda, T) \tau(\lambda) + [1 - \varepsilon(\lambda)] L\_d(\lambda) \tau(\lambda) + L\_u(\lambda) \tag{6}$$

Then, the standard deviation (cost function) of the estimated at-sensor radiance and the true at-sensor radiance is calculated:

$$\sigma(\hat{T}, \overline{\hat{\varepsilon}}) = \sigma(L\_{\text{fil}}(\lambda, \hat{T}, \overline{\hat{\varepsilon}}) - L(\lambda, \hat{T}, \overline{\hat{\varepsilon}})) \tag{7}$$

We iteratively cycle through Equations (5)–(7) until the cost function reaches the minimum, and the corresponding temperature is seen as the optimal LST. Then we obtain the retrieval emissivity with the optimal temperature through Equation (4).

The above processes do not consider instrument noise. In fact, the output radiance of the thermal infrared hyperspectral sensor inevitably includes instrument noise. Instrument noise is a kind of random noise. In this study, the instrument noise is represented by additive Gaussian white noise, and the bands are uncorrelated. That is, to the right side of in Equation (1), add a noise term <sup>η</sup>(λ*i*). Accordingly, to the right side of in Equation (2), also add a noise term <sup>η</sup>(λ*<sup>i</sup>* ) τ(λ*<sup>i</sup>* ) . Then εˆ(λ*i*) in Equation (4) becomes,

$$\begin{split} \mathcal{E}(\lambda\_{i}) &= \frac{[\mathcal{B}(\lambda\_{i}, T) - \mathcal{L}\_{d}(\lambda\_{i})] \varepsilon(\lambda\_{i}) + \frac{\eta(\lambda\_{i})}{\tau(\lambda\_{i})}}{\mathcal{B}\left(\lambda\_{i}, \mathcal{T}\right) - \mathcal{L}\_{d}(\lambda\_{i})} \\ &= \frac{B(\lambda\_{i}, T) - \mathcal{L}\_{d}(\lambda\_{i})}{B(\lambda\_{i}, \mathcal{T}) - \mathcal{L}\_{d}(\lambda\_{i})} \varepsilon(\lambda\_{i}) + \frac{\frac{\eta(\lambda\_{i})}{\tau(\lambda\_{i})}}{\mathcal{B}\left(\lambda\_{i}, \mathcal{T}\right) - \mathcal{L}\_{d}(\lambda\_{i})} \end{split} \tag{8}$$

where εˆ(λ*i*) and ε(λ*i*) are the estimated emissivity and its true value respectively and *T*ˆ and *T* are the estimated temperature and its true value respectively.

For simplicity, εˆ(λ*i*), ε(λ*i*), *B*(λ*<sup>i</sup>* , *T*ˆ), *B*(λ*<sup>i</sup>* , *T*), τ(λ*i*), *Lu*(λ*i*), *L<sup>d</sup>* (λ*i*) and η(λ*i*) will be abbreviated in the following equations as εˆ*<sup>i</sup>* , ε*<sup>i</sup>* , *B*ˆ *i* , *B<sup>i</sup>* , τ*<sup>i</sup>* , *L* ↑ *i* , *L* ↓ *i* and η*<sup>i</sup>* respectively. We bring Equations (5) and (8) into Equation (7), and the cost function of ARTEMISS, i.e., Equation (7), becomes

σ *T*ˆ, εˆ = σ h ε*i*−1+ε*i*+ε*i*+<sup>1</sup> 3 *B*ˆ *<sup>i</sup>* − *L d i* τ*<sup>i</sup>* + *L d i* τ*i* i − h*B<sup>i</sup>* <sup>−</sup> *<sup>L</sup> d i* ε*i*τ*<sup>i</sup>* + *L d i* τ*<sup>i</sup>* + η*<sup>i</sup>* i = σ *Bi*−1−*L d i*−1 *B*ˆ *<sup>i</sup>*−1−*L d i*−1 ε*i*−<sup>1</sup> <sup>3</sup> + *Bi*−*L d i B*ˆ *<sup>i</sup>*−*L d i* ε*i* <sup>3</sup> + *Bi*+1−*L d i*+1 *B*ˆ *<sup>i</sup>*+1−*L d i*+1 ε*i*+<sup>1</sup> 3 +<sup>1</sup> 3 η*i*−1 τ*i*−1 *B*ˆ *<sup>i</sup>*−1−*L d i*−1 + <sup>1</sup> 3 η*i* τ*i B*ˆ *<sup>i</sup>*−*L d i* + <sup>1</sup> 3 η*i*+1 τ*i*+1 *B*ˆ *<sup>i</sup>*+1−*L d i*+1 *B*ˆ *<sup>i</sup>* − *L d i* τ*<sup>i</sup>* − *B<sup>i</sup>* − *L d i* ε*i*τ*<sup>i</sup>* − η*<sup>i</sup>* = σ 1 3 *Bi*−1−*L d i*−1 *B*ˆ *<sup>i</sup>*−1−*L d i*−1 *B*ˆ *<sup>i</sup>* − *L d i* ε*i*−1τ*<sup>i</sup>* + <sup>1</sup> 3 *Bi*+1−*L d i*+1 *B*ˆ *<sup>i</sup>*+1−*L d i*+1 *B*ˆ *<sup>i</sup>* − *L d i* ε*i*+1τ*<sup>i</sup>* − 2 3 *B<sup>i</sup>* − *L d i* ε*i*τ*<sup>i</sup>* +<sup>1</sup> 3 *B*ˆ *<sup>i</sup>*−*L d i B*ˆ *<sup>i</sup>*−1−*L d i*−1 η*i*−1τ*<sup>i</sup>* τ*i*−<sup>1</sup> + <sup>1</sup> 3 *B*ˆ *<sup>i</sup>*−*L d i B*ˆ *<sup>i</sup>*+1−*L d i*+1 η*i*+1τ*<sup>i</sup>* τ*i*+<sup>1</sup> − 2 3 η*i* (9)

where <sup>σ</sup>(·) in the right side of Equation (9) refers to <sup>q</sup> 1 *N* P*N i* = 1 [(·) 2 ].

In Equation (9), the first three items are not related to noise and the last three items are noise-related. The items not related to noise are fixed and determined by the temperature. In general, the cost function is smallest when *T*ˆ = *T*. However, after noise appears, the cost function may not be the minimum value when *T*ˆ = *T*. That is to say, after the coupling of noise, temperature and transmittance, different noise levels may cause different optimal *T*ˆ appearances for the same emissivity.

#### **3. Sensitivity Analysis of ARTEMISS**

This section analyzed the relationship between instrument spectral resolution, the noise level, the ARTEMISS parameter setting and the retrieval errors of ARTEMISS through simulation data. This section described the simulation experiment in Section 3.1 and the results in Section 3.2.

#### *3.1. Simulation Experiment*

A thermal infrared hyperspectral imager is the primary tool used to acquire HTIR remote sensing data. The instrument's settings and performance directly affect data quality, which in turn affects the accuracy of subsequent TES process. On the contrary, the TES algorithms also influence the development of HTIR instruments over time. Based on this background and recent research [48], this study mainly uses the instrument's key parameters, such as spectral range, noise level, spectral resolution and ARTEMISS parameter setting, to explore their influence on the accuracy of the TES algorithm. Note that the scenes simulated in this study were only isothermal homogeneous pixels, the complex radiation transfer process among non-isothermal heterogeneous pixels involves other

complex problems, such as pixel decomposition, which may not be conducive to focusing on evaluating the relationship between the accuracy of the TES algorithm and the selected impact factors.

During the simulation of at-sensor radiance, this study considered three impact factors unrelated to the algorithm itself (i.e., spectral range, noise level and spectral resolution) and four input parameters (i.e., atmospheric model, sensor altitude, surface temperature and surface emissivity), as shown in Table 1. They were set as follows:



(1) Spectral range: 7.5–13.5 µm is generally regarded as the atmospheric window of thermal infrared remote sensing. Considering the low detector response efficiency at 12.5–13.5µm and the atmospheric absorption at 7.5–8.0 µm, this study used two spectral ranges: 7.5–12.5 µm and 8–12.5 µm.

(2) Noise level: considering that a HTIR imager is a system, the noise in calibrated data is an accumulation of noise from the radiance entering the imager to have been calibrated, including detector noise, circuit noise, calibration source noise, etc. The calibration errors include both systematic and random errors. This study only considered the random calibration error, which is more difficult to be removed than systematic error. Therefore, the noise here refers to a comprehensive noise equivalent temperature difference, namely the random error in calibrated data, which is simulated by a Gaussian function [37,44]. Taking the current noise equivalent temperature difference (NEDT) of the instrument, which is generally below 0.3 K [54], and other sources of random errors into account, in this study we tested 11 noise levels, and the specific value was from 0 to 0.5 K with an incremental increase of 0.05 K.

(3) Spectral resolution: the HTIR imager splits the thermal infrared radiation spectrum into hundreds and thousands, approximately continuous and narrow bands, and the spectral response of each channel conforms to the Gaussian distribution [55,56]. When dealing with multi-spectral thermal infrared data, the spectral response function is usually obtained by direct measurement [57,58]. However, for the HTIR imager, the spectral response performance at each channel is usually fitted by a Gaussian function after spectral calibration by devices such as a monochromator, and described by the center wavelength and full width at half height (FWHM) obtained by fitting. Practically, the spectral response function of the HTIR imager is obtained by putting the center wavelength and FWHM at each band into a Gaussian function. In general, the spectral resolution of HTIR imagers is higher than 100 nm (one hundredth of a wavelength). The spectral resolution of most current airborne instruments is above 50 nm [54]. Four spectral resolutions (50 nm, 35.2 nm, 10 nm and 5 nm) were tested to analyze the influence of the spectral range and spectral resolution on the retrieval results (Table 1). The two spectral resolutions of 50 nm and 35.2 nm were from the airborne thermal-infrared hyperspectral imager system

(ATHIS) [55] and hyperspectral thermal emission spectrometer (HyTES) [59] respectively. The other two spectral resolutions of 10 nm and 5 nm were set to provide theoretical analysis for the design and development of HTIR instruments with higher spectral resolution in the future.

(4) Atmospheric model: to include as many atmospheric conditions as possible, we chose five atmospheric models from moderate resolution atmospheric transmission (MODTRAN) software developed by Spectral Sciences, Inc., Burlington, United States [60], including tropical atmospheric model, mid-latitude summer (winter) atmospheric model and sub-arctic summer (winter) atmospheric model.

(5) Sensor altitude: the development of HTIR imagers is currently in its infancy. There are no spaceborne HTIR instruments in orbit yet, and airborne HTIR imagers are still the mainstream. Thus, airborne HTIR data are the main source currently available. The sensor altitude directly affects atmospheric transmittance and path radiance. With due consideration of the airborne HTIR remote sensing scene, this study sets three sensor altitudes of 2 km, 10 km and 750 km for low aerial, high aerial and spaceborne instruments, respectively.

(6) Surface temperature: the heat between the land surface and atmosphere is frequently exchanged. LST generally fluctuates more than atmospheric bottom temperature. Therefore, the surface temperature is set around the atmospheric temperature. That is, the surface temperature is equal to the atmospheric temperature plus an offset. In the simulation, when the atmospheric temperature was higher than 280 K, the offset was set from −5 to 20 K with an incremental step of 5 K; when the atmospheric temperature was lower than 280 K, the offset was set from −10 to 15 K with an incremental step of 5 K.

(7) Surface emissivity: the theoretical basis of ARTEMISS is that the smoothness of LSE is higher than the downwelling radiance spectrum of the atmosphere. It is necessary to choose a dataset containing a large number of spectra to test the influence of the fluctuation of spectral curves of different ground objects on inversion accuracy. We chose the advanced spaceborne thermal emission and reflection radiometer (ASTER) 2.0 spectral library [61], which is often adopted in a remote sensing numerical simulation experiment. There are 2445 spectra of different ground objects in the library. Among them, 1524 emissivity curves covering thermal infrared wavelengths were selected, simulating the majority of scenarios in practical applications. The selected dataset involves common ground object types including stony mineral (909), rock (388), manmade (84), soil (75), stony meteorite (60), vegetation (4), water (4), snow (4) and ice (1) [48].

The influence of ARTEMISS's smooth window size on the inversion accuracy was studied during the process of TES in this study. The method used in [48] was adopted in the terms of searching for the optimal temperature; the specific settings were as follows:

(1) Smooth window size setting: here, the algorithm parameter refers to the best window size for smoothing, where the retrieval error is the smallest among the several window sizes for a group of data at a certain noise level. We used the ARTEMISS algorithm with 17 smoothing window sizes (from 3 to 35 with an incremental step of 2) to retrieve the simulated at-sensor radiances to analyze the relationship between the ARTEMISS spectral smoothing window sizes and the retrieval errors.

(2) The strategy of optimal temperature solution: to reduce the computational complexity and avoid the local minimum of the cost function, the strategy of the optimal temperature solution adopted in this study was the same as in [48]. First, a search range of temperature was set, as in [48]. Then, the cost function was calculated for each temperature in the temperature search range. Finally, the temperature corresponding to the minimum cost function value was seen as the optimal temperature. According to the literature [48], the optimal temperature obtained with a search range of true LST ± 100 K is basically same with that obtained with a search range of true LST ± 20 K. Therefore, in order to save calculation time, the true surface temperature ± 20 K was set as the temperature search range in this study.

We used the root mean square error (RMSE) metrics to evaluate the retrieval temperature error. The RMSE of a group of retrieval LSTs with a certain spectral parameter and noise level is defined as:

$$\text{RMSE}\_{T} = \sqrt{\frac{1}{N\_s} \sum\_{j=1}^{N\_s} \left(\hat{T}\_j - T\_j\right)^2} \tag{10}$$

where *j* ∈ [1, *N<sup>s</sup>* ], *N<sup>s</sup>* is the number of the samples in the group of data to be evaluated; and *T*ˆ *<sup>j</sup>* and *T<sup>j</sup>* are the retrieval LST and the real LST of a sample in the group of data to be evaluated respectively.

The RMSE is sensitive to the amplitude of error at some bands (especially large errors), but it is not sensitive to the overall deviation tendency of retrieval LSE from its true LSE. The overall deviation means that the amplitudes of error at all bands are about the same. Therefore, we used both the RMSE and median absolute deviation (MAD) to evaluate the retrieval LSE. This strategy can simultaneously measure the emissivity error in terms of both details and overall morphologies.

The RMSE and MAD of a group of retrieval LSEs with a certain spectral parameter and noise level are defined as:

$$\text{RMSE}\_{\varepsilon} = \frac{1}{N\_s} \sum\_{j=1}^{N\_s} \sqrt{\frac{1}{N} \sum\_{i=1}^{N} \left( \varepsilon\_j(\lambda\_i) - \varepsilon\_j(\lambda\_i) \right)^2} \tag{11}$$

$$\text{MAD}\_{\mathcal{E}} = \frac{1}{N\_s} \sum\_{j=1}^{N\_s} \text{median} \{ \hat{\varepsilon}\_j(\lambda\_i) - \varepsilon\_j(\lambda\_i) \} \tag{12}$$

where *j* ∈ [1, *N<sup>s</sup>* ], *N<sup>s</sup>* is the number of samples in a group of data to be evaluated; *i* ∈ [1, *N*], *N* is the number of bands of a retrieval LSE spectrum; <sup>q</sup> 1 *N* P*N i* = 1 εˆ*j*(λ*i*) − ε*j*(λ*i*) 2 is RMSE of retrieval LSE spectrum; *median* εˆ*j*(λ*i*) − ε*j*(λ*i*) is MAD of retrieval LSE spectrum and εˆ*j*(λ*i*) and ε*j*(λ*i*) is the retrieval and true LSE of sample j at the i-th band, respectively.

The data simulation included two processes: the ultra-high-resolution (1 nm) at-sensor radiance simulation and sensor output radiance simulation, as detailed in article [48]. The simulation experiment was based on the interface data language (IDL), Exelis Visual Information Solutions, Inc. MODTRAN was called by an IDL program to simulate the atmospheric parameters. All other experimental processes were implemented by the IDL code. In order to facilitate the drawing and explanation of the results, the simulation results were classified into eight groups according to spectral range and spectral resolution in this article, as shown in Table 2.


**Table 2.** Groups of simulated hyperspectral thermal infrared (HTIR) data.

#### *3.2. Results Analysis*

3.2.1. The Relationship of the Retrieval Errors of ARTEMISS vs. the Noise Level and Spectral Parameters

This study retrieved eight groups of simulated HTIR data with two spectral ranges (7.5–12.5 µm and 8–12.5 µm) and four spectral resolution (50 nm, 35.2 nm, 10 nm and 5 nm) by the ARTEMISS algorithm and measured the retrieved LST and LSE errors. Figure 1 shows the retrieved LST errors of the eight groups of simulated HTIR data as a function of the noise level.

μ

μ

**Figure 1.** Plot of automatic retrieval of temperature and emissivity using spectral smoothness (ARTEMISS) land surface temperature (LST) inversion errors varying with noise equivalent temperature difference (NEDT), spectral range and spectral resolution.

μ μ μ μ μ μ μ μ μ μ μ For all eight groups of simulated HTIR data shown in Figure 1, the retrieved LST errors of ARTEMISS increased linearly, roughly with the noise level. For the four groups of simulated data with high spectral resolution of 10 nm and 5 nm, after the NEDT of 0.2 K, the retrieved LST errors increased with the noise level at a slower growth rate. For the eight groups of simulated data except the two with the spectral resolution of 5 nm, the retrieved LST errors of ARTEMISS decreased as the spectral resolution increased. However, the two groups of data with a spectral resolution of 5 nm show differences to each other. For the spectral range of 8–12.5 µm, the retrieved LST errors of 5 nm at 8–12.5 µm (refers to Group8 with a spectral resolution of 5 nm and in a spectral range of 8–12.5 µm) were smaller than those with a lower spectral resolution (50 nm and 35.2 nm) and were consistent with those of 10 nm and a spectral range of 8–12.5 µm in terms of both amplitude and trend. For the spectral range of 7.5–12.5 µm, the retrieved LST errors of 5 nm at 7.5–12.5 µm were larger than those with a lower spectral resolution (35.2 nm and 10 nm) and were numerically similar with those of 50 nm at 7.5–12.5 µm. For the spectral resolutions of 50 nm and 32.5 nm, the retrieved LST errors with a spectral range of 7.5–12.5 µm were smaller than those in the spectral range of 8–12.5 µm. Among the eight groups of results, the group with 10 nm at 8–12.5 µm had the smallest retrieved LST error (e.g., LST RMSE = 0.11 K, 0.55 K and 1.10 K when NEDT = 0.05 K, 0.20 K and 0.50 K, respectively), 50 nm at 8–12.5 µm had the largest retrieval LST error (e.g., LST RMSE = 0.46 K, 1.20 K and 2.56 K when NEDT = 0.05 K, 0.20 K and 0.50 K, respectively). The results suggest that for temperature inversion, a higher spectral resolution and wider spectral range do not necessarily lead to smaller LST inversion errors. The results of this study show that, considering noise, 10 nm at 8–12.5 µm was the most suitable setting for thermal infrared hyperspectral temperature inversion among the eight spectral settings investigated.

Figure 2a,b shows the RMSE and MAD of ARTEMISS retrieved LSE for eight groups of simulated HTIR data. Figure 2a shows that the retrieved LSE RMSEs of ARTEMISS increased approximately linearly with the noise level for all eight groups of simulated HTIR data, which is consistent with the retrieved LST RMSEs in Figure 1. However, the retrieved LSE RMSEs show a significant difference between the two spectral ranges. The retrieved LSE RMSEs of all four groups of data in the spectral range of 8–12.5 µm were smaller than the retrieved LSE RMSEs of four groups of data with in spectral range of 7.5–12.5 µm. Moreover, the relationship between the retrieved LSE RMSEs and the spectral resolution were different in the two spectral ranges. For the four groups of data in spectral range of

7.5–12.5 µm, the retrieved LSE RMSEs of ARTEMISS increased as the spectral resolution increased; in contrast, for the four groups of data in the spectral range of 8–12.5 µm, the retrieved LSE RMSEs of ARTEMISS decreased as the spectral resolution increased (except the group of 5 nm at 8–12.5 µm). Among the eight groups, 10 nm at 8–12.5 µm had the smallest LSE RMSEs (0.0021, 0.0082 and 0.0203 when NEDT = 0.05 K, 0.20 K and 0.5 K, respectively), while the 5 nm at 7.5–12.5 µm had the largest LSE RMSEs (0.0173, 0.0418 and 0.0796 when NEDT = 0.05 K, 0.20 K and 0.5 K, respectively). Figure 2b shows that the trend of the retrieved LSE MADs of ARTEMISS was close to that of the retrieved LST RMSEs for the eight groups. For the eight groups of data with two spectral ranges, the retrieved LSE MADs of ARTEMISS roughly increased linearly with the noise level. Moreover, the relationship between the retrieved LSE MADs and the spectral resolution were the same for the two spectral ranges. The retrieved LSE MADs decreased as the spectral resolution increased, except for 5 nm at 7.5–12.5 µm. Among the eight groups, 5 nm at 8–12.5 µm had the smallest retrieved LSE MADs and 50 nm at 8–12.5 µm had the largest retrieved LSE MADs. μ μ μ μ μ μ μ μ

μ

μ

μ

**Figure 2.** Plots of ARTEMISS retrieved land surface emissivity (LSE) errors varying with NEDT, spectral range and spectral resolution. (**a**) The retrieved LSE root mean square error (RMSE) and (**b**) the retrieved LSE median absolute deviation (MAD).

μ μ The LSE RMSE represents the degree to which the retrieved emissivity spectrum deviates from the true value at each band numerically, while LSE MAD represents the degree to which the retrieved emissivity spectrum deviates from the true spectral curve in the overall shape. Figures 1 and 2b show that the smaller the temperature inversion error, the smaller the emissivity MAD, indicating that the more accurate the temperature inversion, the closer the overall trend of the inversion emissivity curve is to the true spectrum curve. Both the retrieved LST errors and the retrieved LSE MADs of the four groups of data with a spectral range of 7.5–12.5 µm are between those of 10 nm at 8–12.5 µm and 50 nm at 8–12.5 µm with a spectral range of 8–12.5 µm. However, the retrieved LSE RMSEs of all four groups of data with a spectral range of 7.5–12.5 µm were larger than those with a spectral range of 8–12.5 µm. The reason for this seemingly contradictory phenomenon is that, compared to 8–12.5 µm, the retrieved LSE spectra with a spectral range of 7.5–12.5 µm (especially 7.5–8.0 µm is where dense atmospheric absorption lines exist) have large outliers in some bands, as shown in Figure 3. The existence of the huge outliers causes the cases that although the retrieved LSE has a large RMSE, the overall shape of the retrieved emissivity is closer to the true curve than some retrieved LSE with a small RMSE. The finding that in terms of the spectral range of 7.5–12.5 µm, the higher the spectral resolution is, the larger the retrieved LSE RMSE, suggests that for the cases with a spectral range of 7.5–12.5 µm,

μ

μ μ

the higher the spectral resolution was, the more outliers of the retrieved LSE RMSE occurred. This also suggests that for thermal infrared hyperspectral remote sensing, the spectral range of 8–12.5 µm is more suitable for temperature inversion, which is similar to the statement in [42]. μ

μ

μ μ

μ μ

μ

**Figure 3.** A case of outliers occurring in the retrieved LSE spectrum.

3.2.2. The Relationship of the Retrieval Errors vs. the Optimal Window Size Spectral Smoothing

To investigate the effect of the window size for the spectral smoothing index in ARTEMISS on the retrieved LST and LSE errors, we tested 17 spectral smoothing window sizes for ARTEMISS, going from 3 to 35 in increments of 2 for eight groups of simulated data with two spectral ranges and four spectral resolutions. In the process of a simulated radiance spectrum inversion, we looped the window size for spectral smoothing, calculated the cost functions and found the smallest one among all of the cost functions under 17 window sizes. The temperature at the minimum cost function was considered as the retrieved LST, the corresponding emissivity was the retrieved LSE, and the corresponding window size was the optimal spectral smoothing window size because the retrieved LST error was close to zero in the majority of cases. Figure 4 has eight subplots for the eight groups of simulated data, each of which shows the distribution of the optimal spectral smoothing window size vs. the noise level under a certain spectral setting. The distributions of the optimal spectral smoothing window size vs. the noise level show a consistent trend for the eight spectral settings. For each noise level, the window size with the largest proportion is 3. Generally, the larger the window size, the smaller the proportion (except 35, at the end of the window size range). In the ideal situation without noise, a window size of 3 has the most advantageous proportion in the optimal spectral smoothing window size distribution, especially when the spectral resolution is high (e.g., 10 nm and 5 nm); as the noise level increased, the proportion of the window size 3 gradually decreased for a certain spectral setting. The results suggest that when the hyperspectral thermal infrared data is noisy, it may be beneficial to improve the accuracy of temperature inversion by increasing the spectral smoothing window size in some cases.

**Figure 4.** The distribution of the optimal spectral smoothing window size varying with the noise level.

̂ߝ ,ܶ൫ߪ

൯ ݂ఎ൫ܤ

ܤఌ൫݂

ߟ ൯

ܤఎ൫݂

ߟ

ሜ൯ = ߪ ൬ቀ݂ఌ൫ܤ

൯ = <sup>ܤ</sup>ିଵ − ܮିଵ

ܤ

൯ = <sup>1</sup> 3

൯

ܤఌ൫݂

ܤఎ൫݂

൯ + ݂ఎ൫ܤ

ௗ

ଵିߟ ߬ିଵ

ିଵ − ܮିଵ

ଵିߝ 3 +

> ௗ + 1 3

ିଵ − ܮିଵ ௗ

ܤ

൯ቁ ൫ܤ

ܮ − ௗ

ܮ − ܤ ௗ

> ߟ ߬ ܤ ܮ − ௗ + 1 3

ܤ ܮ − ௗ ߝ 3 +

ܮ − ܤ൫ − ߬൯

ௗ

ାଵܮ − ାଵܤ ௗ

> ାଵ − ܮାଵ ௗ

> > ାଵߟ ߬ାଵ

ାଵ − ܮାଵ ௗ

> ܤఌ൫݂ ൯

ܤ

ܤ

൰ߟ − ߬ߝ൯

ାଵߝ 3

> ܤఎ൫݂ ൯

ܤఌ൫݂

൯ ݂ఎ൫ܤ

൯

#### **4. The Proposal and Validation of the RDSS Algorithm**

#### *4.1. RDSS—An Improved TES Algorithm Based on ARTEMISS*

The coupling of temperature and emissivity, as shown in Equation (9), suggests that the temperature is easier to determine than the emissivity in the TES process. Therefore, most algorithms adopt the strategy of first determining the temperature and then calculating the emissivity. Hence, accurately determining the temperature has also become the key to TES. We rewrite Equation (9) by combining the items determined by the temperature estimation and the items determined by both the temperature estimation and noise as follows:

$$\sigma(\hat{\mathcal{T}}, \overline{\hat{\varepsilon}}) = \sigma(\left(f\_{\hat{\varepsilon}}(\nexists\_{l}) + f\_{\eta}(\nexists\_{l})\right)(\hat{\mathcal{B}}\_{l} - L\_{i}^{d})\tau\_{i} - \left(B\_{l} - L\_{i}^{d}\right)\varepsilon\_{i}\tau\_{l} - \eta\_{i}\right) \tag{13}$$

where *f*<sup>ε</sup> *B*ˆ *i* and *f*<sup>η</sup> *B*ˆ *i* are the emissivity estimation after box-average smoothing and the residual containing noise respectively. For the algorithms that use the standard deviation of the simulated at-sensor radiance and the actual at-sensor radiance as the cost function, *f*<sup>ε</sup> *B*ˆ *i* and *f*<sup>η</sup> *B*ˆ *i* can also be considered as the approximately processed emissivity and residual noise, respectively. They can be expressed as

$$f\_{\varepsilon}(\mathfrak{B}\_{i}) = \frac{\mathcal{B}\_{i-1} - L\_{i-1}^{d}}{\mathfrak{B}\_{i-1} - L\_{i-1}^{d}} \frac{\mathcal{E}\_{i-1}}{\mathfrak{B}} + \frac{\mathcal{B}\_{i} - L\_{i}^{d}}{\mathfrak{B}\_{i} - L\_{i}^{d}} \frac{\mathcal{E}\_{i}}{\mathfrak{B}} + \frac{\mathcal{B}\_{i+1} - L\_{i+1}^{d}}{\mathfrak{B}\_{i+1} - L\_{i+1}^{d}} \frac{\mathcal{E}\_{i+1}}{\mathfrak{B}} \tag{14}$$

$$f\_{\eta}(\boldsymbol{\hat{\mathcal{B}}\_{i}}) = \frac{1}{3} \frac{\frac{\eta\_{i-1}}{\tau\_{i-1}}}{\boldsymbol{\hat{\mathcal{B}}\_{i-1}} - \boldsymbol{L}\_{i-1}^{d}} + \frac{1}{3} \frac{\frac{\eta\_{i}}{\tau\_{i}}}{\boldsymbol{\hat{\mathcal{B}}\_{i}} - \boldsymbol{L}\_{i}^{d}} + \frac{1}{3} \frac{\frac{\eta\_{i+1}}{\tau\_{i+1}}}{\boldsymbol{\hat{\mathcal{B}}\_{i+1}} - \boldsymbol{L}\_{i+1}^{d}} \tag{15}$$

When the emissivity approximation processing can remove the "thorns" in a curve caused by atmospheric background radiation and obtain relatively accurate *f*<sup>ε</sup> *B*ˆ *i* , and the noise residuals *f*<sup>η</sup> *B*ˆ *i* with a noise item η*<sup>i</sup>* of 0 or close to 0, the searched optimal temperature is likely to be equal or much closer to the true temperature. As the magnitude of noise residuals *f*<sup>η</sup> *B*ˆ *i* and noise term η*<sup>i</sup>* increase, it is likely that the signs of the items will cancel each other out, resulting in the searched optimal temperature deviating from the true surface temperature. The results of ARTEMISS sensitivity on the spectral smoothing window size in Section 3.2.2 suggest that increasing the spectral smoothing window size did not effectively reduce the inversion error; different window sizes are needed to find the optimal temperature as close to the true one possible when different noises, emissivity and transmittances are coupled together. However, it will add computation complexity and cost much more computation time. That is to say, it is difficult to remove the influence of noise on the optimal temperature search simply by emissivity approximation processing. Therefore, in this study, from the perspective of reducing the noise of the data, we filtered the variables participating in the calculation of the emissivity estimation with a unified filter. The unified filter ensured the variables were at the same spectral resolution while eliminating the noise in the at-sensor radiances. We named the TES algorithm proposed in this study the resolution-degrade-based spectral smoothness (RDSS) algorithm. The calculation process of RDSS is:

(1) Unified filtering: before calculating the emissivity with Equation (4), we filtered each variable participating in the calculation of the emissivity estimation with mean filtering. The calculation method was as follows:

$$
\overline{V} = F(V) \tag{16}
$$

where *F* is the filter function, *V* is the variable to be filtered and *V*e is the filtered variable. Here, the variables refer to the ground-leaving radiance, the atmospheric downwelling radiance and the blackbody radiance.

*Remote Sens.* **2020**, *12*, 2295

Specifically, the calculation method for filtering the variables in Equation (4) with a mean filter was as follows:

$$\left[\widetilde{L}\_3(\lambda\_i, T), \widetilde{L}\_d(\lambda\_i), \widetilde{\mathcal{B}}(\lambda\_i, \mathcal{T})\right] \\ = \frac{1}{2N\_r + 1} \sum\_{j = -N\_r}^{N\_r} \left[L\_3(\lambda\_{i+j}, T), L\_d(\lambda\_{i+j}), \mathcal{B}(\lambda\_{i+j}, \mathcal{T})\right] \\ \tag{17}$$

where e*Lg*(λ*<sup>i</sup>* , *T*), e*L<sup>d</sup>* (λ*i*) and e*B* λ*i* , *T*ˆ are the filtered *L<sup>g</sup>* λ*i*+*<sup>j</sup>* , *T* , *L<sup>d</sup>* λ*i*+*<sup>j</sup>* and *B* λ*i*+*<sup>j</sup>* , *T*ˆ , respectively, with the mean filter; and 2*N<sup>r</sup>* + 1 is the spectral smoothing window size.

(2) Emissivity estimation and approximate processing: to avoid the local minimum value issue, taking advantage of the fact that the object temperature has a finite value range as a physical parameter, this study adopted the strategy of an exhaustive search to find the optimal temperature. We determined a temperature search range *Tsupin f* and the incremental step δ*T* (also named as the temperature resolution), then calculated the emissivity estimations for all of the temperatures in the search range and approximately processed the emissivity estimations (i.e., boxcar average):

$$\mathcal{E}\left(\lambda\_{i\prime}, T\_{inf}() \frac{\widetilde{L}\_{\mathcal{S}}(\lambda\_{i\prime}T) - \widetilde{L}\_{d}(\lambda\_{i})}{\widetilde{B}(\lambda\_{i\prime}T - \widetilde{L}\_{d\inf})}\right) \tag{18}$$

$$\left\{\lambda\_{i\nu}, T\_{inf}(\cdot)\frac{1}{3}\sum\_{j=-1}^{1\sum}\varepsilon\{\lambda\_{i+j\nu}, T\_{inf}(\cdot)\}\right\}\tag{19}$$

where εˆ λ*i* , *<sup>T</sup>in f*() is the emissivity estimation at the temperature of *Tin f* , *k* is the steps, *k* = 0, 1, 2, . . . , *<sup>T</sup>*sup <sup>−</sup> *<sup>T</sup>*inf /δ*<sup>T</sup>* and <sup>e</sup><sup>ε</sup> λ*i* , *<sup>T</sup>in f*() is the corresponding emissivity estimation of εˆ λ*i* , *<sup>T</sup>in f*() after approximate processing.

(3) Cost function calculation and the optimal temperature determination. We calculated the cost functions for all temperatures in the temperature search range according to the approximate processed emissivity. The calculation method was as follows:

$$\sigma(T\_{\text{inf}} + \delta T \cdot k) = \sigma(\left(\widetilde{\mathcal{B}}(\lambda\_{i\prime} T\_{\text{inf}} + \delta T \cdot k) - \widetilde{L}\_{\text{d}}(\lambda\_{i})\right) \widetilde{\varepsilon}(\lambda\_{i\prime} T\_{\text{inf}} + \delta T \cdot k) + \widetilde{L}\_{\text{d}}(\lambda\_{i}) - \widetilde{L}\_{\text{g}}(\lambda\_{i\prime} T) \tag{20}$$

$$\inf\_{k} \underset{k}{\text{arg min}} \int\_{k} \sigma(T\_{\text{inf}}())() \underset{\text{ $\mathcal{T}\_{\text{opt}}}}{\text{$ \mathcal{T}\_{\text{opt}}}} = \text{ $\mathcal{T}$ } \tag{21}$$

where *T*ˆ *opt* is the optimal temperature estimation, i.e., the retrieved temperature.

(4) Final emissivity calculation: we substituted the retrieved temperature from (3) and the other three variables before filtering them into Equation (4) to obtain the final emissivity. The calculation method was:

$$\mathcal{E}\_{opt}(\lambda\_i) = \frac{L\_{\mathcal{S}}(\lambda\_{i\prime}T) - L\_d(\lambda\_i)}{B(\lambda\_{i\prime}\hat{T}\_{opt}) - L\_d(\lambda\_i)}\tag{22}$$

where εˆ*opt*(λ*i*) is the final emissivity estimation, i.e., the retrieved emissivity.

#### *4.2. Validation Results Analysis*

#### 4.2.1. The Retrieval Errors of the RDSS Algorithm

We retrieved the eight groups of simulated HTIR data with two spectral ranges (7.5–12.5 µm and 8–12.5 µm) and four spectral resolutions (50 nm, 35.2 nm, 10 nm and 5 nm) using the RDSS algorithm and measured the retrieved LST and LSE errors. Figures 5 and 6 show the retrieved LST errors and the retrieved LSE errors of the eight groups of simulated HTIR data as a function of the noise level, respectively. The change trends of the retrieved LST and LSE errors of the RDSS algorithm with noise, spectral resolution and spectral range were basically the same as those of the ARTEMISS algorithm, μ

μ

namely both the LST and LSE errors increase with noise linearly. The difference is that in terms of the magnitude of errors, the LST and LSE errors of RDSS were smaller than those of ARTEMISS except for some LSE RMSEs.

μ

μ

**Figure 5.** The plot of resolution-degrade-based spectral smoothness (RDSS) algorithm LST inversion errors varying with NEDT, spectral range and spectral resolution.

**Figure 6.** Plots of RDSS retrieved LSE errors varying with NEDT, spectral range and spectral resolution. (**a**) The retrieved LSE RMSE and (**b**) the retrieved LSE MAD.

Figure 7 shows a comparison chart of RDSS and ARTEMISS on the retrieved LST error, and Table 3 shows the corresponding error list. They show that under most noise and spectral settings, the retrieved LST error of RDSS are reduced by varying degrees, with a maximum reduction of 0.75 K (Group4, NDET = 0.5 K). The reduced values of the retrieved LST errors vary for different noise levels and spectral settings. Specifically, the reduced degrees of RDSS LST errors increase with the noise level; however, LST errors of RDSS are larger than ARTEMISS when there are no noise or less noise. The initial noise level at which LST errors of RDSS became smaller than ARTEMISS decreased as the spectral resolution increased. For example, the initial noise level decreased from 0.15 (Group1) and 0.2 K

(Group5) with a spectral resolution of 50 nm to 0.05 K and 0.05 K (Group4 and Group8) with the spectral resolution of 5 nm, respectively. The above trends indicate that in terms of temperature inversion, RDSS was more effective for data with high noise levels and a high spectral resolution. The same temperature inversion accuracy can be obtained by RDSS, while the instrument design specifications or data correction accuracy were reduced.


**Table 3.** The retrieved LST error list of the RDSS and ARTEMISS and their differences.

<sup>1</sup> Group 1. <sup>2</sup> ARTEMISS. <sup>3</sup> RDSS. <sup>4</sup> Difference between RDSS and ARTEMISS (i.e., RDSS – ARTEMISS). <sup>5</sup> Relative difference between RDSS and ARTEMISS (RDSS – ARTEMISS)/ARTEMISS). <sup>6</sup> The gray values show how much RDSS performs better than ARTEMISS and the bold values are the maximum reduction in error. Figure 8 shows a comparison of RDSS and ARTEMISS on the retrieved LSE RMSEs, and Table 4 is a list of the RMSE corresponding to Figure 8. Figure 9 shows a comparison chart of RDSS and ARTEMISS on the retrieval LSE MADs, and Table 5 is a list of the MAD corresponding to Figure 9. The results of LSE RMSE show that unlike the retrieval LST RMSE, for a few of the eleven noise levels and eight spectral settings (31 of 88 cases), the retrieved LSE RMSEs from the RDSS are smaller than those of the ARTEMISS algorithm. From the perspective of LSE RMSE, it seems that the RDSS had limited improvement on the emissivity inversion accuracy, and was not as good as ARTEMISS in terms of the overall effect. However, the results for LSE MAD tell another story.


**Table 4.** The retrieved LSE RMSE for RDSS and ARTEMISS and their differences.

**Figure 7.** The comparison of the retrieved LST errors of RDSS and ARTEMISS.

**\**


**Table 5.** The retrieved LSE MAD for the RDSS and ARTEMISS algorithms and their differences.

**Figure 8.** The comparison of the retrieved LSE RMSEs of RDSS and ARTEMISS.

**Figure 9.** A comparison of the retrieved LSE MADs using the RDSS and ARTEMISS algorithms.

Figure 8 shows a comparison of RDSS and ARTEMISS on the retrieved LSE RMSEs, and Table 4 is a list of the RMSE corresponding to Figure 8. Figure 9 shows a comparison chart of RDSS and ARTEMISS on the retrieval LSE MADs, and Table 5 is a list of the MAD corresponding to Figure 9. The results of LSE RMSE show that unlike the retrieval LST RMSE, for a few of the eleven noise levels and eight spectral settings (31 of 88 cases), the retrieved LSE RMSEs from the RDSS are smaller than those of the ARTEMISS algorithm. From the perspective of LSE RMSE, it seems that the RDSS had limited improvement on the emissivity inversion accuracy, and was not as good as ARTEMISS in terms of the overall effect. However, the results for LSE MAD tell another story.

Figure 9 and Table 5 show that the retrieved LSE MADs of RDSS have similar distribution with the retrieved LST RMSE for the eight groups of data. For the majority of cases (79 out of 88 cases), the retrieved LSE MADs of RDSS were smaller than those of ARTEMISS. The results of LSE RMSE and LSE MAD behaved differently, probably due to the RMSE itself. The RMSE of an emissivity will be much larger because of some large abnormal values, although the retrieved temperature was accurate. As shown in Figure 10, there were some cases where the retrieved LSE of RDSS with a small LST error and LSE MAD was much closer to the true value from the perspective of the trend of the emissivity curve; however, the RMSE of the emissivity calculated using the accurate temperature was relatively large due to the influence of individual outliers. The phenomenon that the RMSE of an emissivity did not decrease but increased when the retrieved temperature was more accurate was also described in another paper [34].

**Figure 10.** A case with a small LST error and LSE MAD but large LSE RMSE.

This study used MAD to evaluate the retrieved LSE error to measure the quality of the retrieved LSE curve from the perspective of the overall degree of deviation of emissivity from the true value. However, MAD cannot be compared with RMSE in terms of values, because MAD focuses on the deviation degree of the overall trend, which is different from the commonly used RMSE in terms of magnitude. We used the percentage of MAD to measure the reduction of the retrieved LSE error, as an indicator reflecting how close the trend of the emissivity curve retrieved by RDSS was to the true emissivity curve. Table 5 shows that the maximum relative decline in the percentage of LSE MAD can reach 31% (G5, 0.45 K). As the noise increased, the decline in the percentage of the LSE MAD of RDSS increased accordingly. When there was no noise or less noise, however, the MAD of RDSS instead increased. This is similar to the retrieved LST error in terms of change characteristics. However, the decline in the percentage of MAD of RDSS did not increase with the spectral resolution. In the range of 7.5–12.5 µm (G1–G4), the decline in the percentage of retrieval LSE MAD firstly decreased, and then increased along with the spectral resolution (it reached a minimum when the spectral resolution was 10 nm); in the range of 8.5–12.5 µm (G5−G8), the decline in the percentage of retrieved LSE MAD gradually decreased as spectral resolution increased. Overall, when the noise level was greater than 0.1 K, the decline in the percentage of the retrieved LSE MAD of RDSS was greater than 10%, which proved the effectiveness of the RDSS algorithm.

### 4.2.2. The Relationship of RDSS and the Window Setting

During spectral degradation, the window setting (2*N<sup>r</sup>* + 1 in Equation (17)) is a key parameter. Different window settings may lead to different results, and temperature inversion is the key of the iterative spectral smoothness TES algorithms. The results of this study, at least the LSE MAD results, suggest that the improvement of temperature inversion accuracy is beneficial to the improvement of emissivity inversion accuracy. Therefore, we set 17 spectral smoothing window sizes, ranging from 3

increasing incrementally by 2–35, and separated the temperature and emissivity of the eight sets of simulated data to study the effect of window size on the performance of RDSS. First, we performed a mean filter cyclically with a window size on a radiance curve to obtain spectral degraded radiance curves under different window settings; second, we retrieved the temperatures under different window settings by RDSS, and looped the two steps for the eight groups of data; finally, we measured the retrieved temperature errors under different window settings.

Figure 11 shows the retrieved LST errors for the RDSS algorithm varying with 11 noise levels under 17 window settings for the eight groups of data. For the four groups of data with lower spectral resolutions (G1–2 and G5–6), 3 was the optimal window setting using RDSS, where the retrieved LST errors were the lowest among the 17 window settings under all or most noise levels. For the four groups of data with higher spectral resolutions (G3–4 and G7–8), the optimal window setting increased rapidly with the noise level. For the two groups of data with a spectral range of 7.5–12.5 µm (G3 and G4), the optimal window setting tended to be stable when the noise level was above 0.15 K, i.e., 27 for G3 and 31 for G4. For the group of data in the spectral range of 8–12.5 µm and the spectral resolution of 10 nm (G7), the optimal window setting tended to be stable at a noise level above 0.25 K—5 for G7. For the group of data in a spectral range of 8–12.5 µm and with a spectral resolution of 5 nm (G8), the optimal window setting was not stable but fluctuated in the range 5–21 when the noise level was above 0.05 K. It is worth noting that for the cases in G8 (5 nm at 8–12.5 µm), as the noise increased, the retrieved LST errors corresponding to the optimal window setting and the second optimal one, or even some other window settings, were very close to each other, and their difference was within 0.01 K. Therefore, it can be considered that when separating the LST and LSE with RDSS for HTIR data, we obtained the minimum LST error with a window setting of 3 for HTIR data with a spectral resolution in the tens of nanometers, and need to increase the window setting according to the spectral range and spectral resolution for HTIR data with a spectral resolution of 10 nm.

**Figure 11.** The retrieved LST errors of RDSS under different window size settings.

#### **5. Discussion**

This study focused on the sensitivity of an iterative spectral smoothness TES algorithm to the spectral range, spectral resolution, noise and smooth window size. We analyzed the characteristics of the influence of the four variables on temperature and emissivity inversion accuracy, and proposed an improved algorithm called RDSS. Compared with previous studies, this study analyzed the influence of the key spectrometer indices on the inversion error more comprehensively, especially the spectral resolution and the smooth window size.

The improved method, RDSS, performs spectral degradation on the input data, which can eliminate the effect of random noise on the search for the optimal temperature and improve the accuracy of temperature and emissivity inversion. In terms of method validation, the effect of random noise on RDSS was mainly studied in this study, while the effect of atmospheric error on RDSS needs further study. However, RDSS can also further improve the inversion accuracy by effectively removing random noise in atmospheric correction errors. For the systematic errors in atmospheric correction, the performance of the RDSS algorithm is similar to the ARTEMISS algorithm.

The results of this study also provide a valuable reference for the performance improvement of hyperspectral thermal infrared instruments and data processing. The improvement in spectral resolution is beneficial to reducing the error of temperature inversion. However, for a HTIR imager, the spectral resolution is dependent on the spatial resolution. In other words, the spatial resolution may decrease while the spectral resolution is improved, which will lead to the increase of mixed pixels. Without a good pixel decomposition algorithm, the accuracy of TES may decrease. Therefore, it is not a case of "the higher the spectral resolution, the better". The choice of spectral resolution needs to consider the application requirements, inversion algorithm, instrument design and other factors, and seek the optimal solution. This may be a very meaningful research topic.

From the perspective of the coupling of temperature and emissivity, accurately determining the temperature is an extremely critical step in the TES process. It is naturally believed that the improvement on the accuracy of temperature inversion will bring corresponding improvements in the accuracy of emissivity inversion. However, there are some anomalies that when the retrieval LST error is reduced to a certain value, the RMSE of retrieved LSE increases instead of decreasing. The emissivity of an object is a spectral curve composed of a set of digits, and the evaluation of the retrieval LSE error is equivalent to the measurement of the relationship between the two sets of digits, i.e., the retrieved LSE and the true value. The temperature of an object, however, is a unique value. Therefore, there is a certain difference between the retrieved LSE errors and the retrieved LST errors on the evaluation method. Moreover, the emissivity retrieved from the HTIR data is mostly used for target identification or surface parameters inversion. We adopted MAD alongside RMSE by combining the evaluation index and the application target in the evaluation of the retrieved LSE. In addition, the retrieved LSEs were not denoised before accuracy evaluation in the experiment, resulting in the retention of many outliers. In the practical application of HTIR data, however, the retrieved LSE curves can be denoised to eliminate outliers; at this moment, the evaluation of the trend of LSE becomes particularly important. Furthermore, the highest spectral resolution in the simulation experiment of the study was set to 5 nm, which is mainly considering that the spectral resolution of the emissivity of the ground objects in the ASTER spectral library is 10 nm or lower. Therefore, cases with a higher spectral resolution than 5 nm need further research when there are emissivity spectra of ground objects with sufficiently high spectral resolution.

### **6. Conclusions**

We focused on the study of instrument noise, a widely existing factor that is difficult to eliminate losslessly in hyperspectral thermal infrared data, and focused on optimizing instrument design. Taking the ARTEMISS algorithm—the representative of the iterative spectral smoothness TES algorithm family with a good application effect—as the object of study. Supplemented by atmospheric radiation transfer simulation software and emissivity spectral library, we carried out simulation and inversion experiments; studied the relationship between the spectral resolution of the instrument, noise level, the ARTEMISS parameter setting and the inversion error; and proposed an improved method—RDSS—based on the mechanism and law of the influence of noise on the inversion error.

The law of the spectral response range, spectral resolution and noise on the inversion error obtained in this study provide a reference for the future development of airborne and spaceborne thermal infrared hyperspectral imagers and optimize the instrument design. The improved TES algorithm can be used for temperature and emissivity inversion of various thermal infrared hyperspectral remote sensing data. The resistance of this method also provides a certain surplus space for the design of the instrument in NEDT.

**Author Contributions:** Conceptualization, H.S., C.L. (Chengyu Liu), F.X. and J.W.; methodology, H.S. and C.L. (Chengyu Liu); software, H.S.; validation, H.S. and C.L. (Chengyu Liu); formal analysis, H.S.; investigation, H.S.; resources, H.S.; data curation, H.S. and C.L. (Chengyu Liu); writing—original draft preparation, H.S.; writing—review and editing, H.S. and C.L. (Chengyu Liu); visualization, H.S.; supervision, C.L. (Chengyu Liu), F.X., C.L. (Chunlai Liu) and J.W.; project administration, F.X., C.L. (Chunlai Li) and J.W.; funding acquisition, F.X., C.L. (Chunlai Li) and J.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Project of Science and Technology Commission of Shanghai Municipality (grant number 18511102202) and the Project of civil space technology pre-research of the 13th five-year plan (grant number D040104).

**Acknowledgments:** The authors would like to sincerely thank Jet Propulsion Laboratory for the ASTER emissivity library. The authors would also like to sincerely thank the editors and anonymous reviewers for their careful reading and constructive comments.

**Conflicts of Interest:** The authors declare no conflict of interest.
