Practical retrieval of land surface emissivity spectra in 8-14μm from hyperspectral thermal infrared data

Hua Wu; Ning Wang; Li Ni; Bo-Hui Tang; Zhao-Liang Li

doi:10.1364/OE.20.024761

1. Introduction

Land surface emissivity (LSE) is an intrinsic property of natural materials and is of great importance in mineral mapping, resource exploration, climate research, and short/medium range forecast applications among others [1–3]. Many attempts have been made to retrieve regional or global LSE from remotely sensed thermal infrared data. However, many LSE products, such as MOD11/MYD11 from Moderate Resolution Imaging Spectroradiometer (MODIS) and AST05 from Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER), only provide the LSE at several fixed wavelength [4, 5], which means the details on the wavelength-variation of the whole emissivity spectra are still lost at most spectral locations. Consequently, a constant emissivity instead of the emissivity spectra was often used in energy balance studies and general circulation models. This substitute may greatly influence the estimation of the surface radiation budget. It was reported that a 10% error on the emissivity approximately corresponds to a 10% error in the energy emitted from the surface [6].

The appearance of the hyperspectral sounders, such as atmospheric infrared sounder (AIRS) on Earth observing system (EOS) satellites, interferometer atmospheric sounding instrument (IASI) on European meteorological operational satellite programme (METOP-A), and cross-track infrared sounder (CrIS) on the next-generation National polar-orbiting operational environmental satellite system (NPOESS), improves the spectral resolution greatly and provides a new opportunity to retrieve LSE with high spectral resolution. Several methods from semi-empirical to physics-based have been developed to retrieve LSE from those sounders [2, 7–10]. Generally, the semi-empirical regression methods are often of high retrieval efficiency and easy to implement. On the contrary, physical methods are time consuming and complicated in essence. The retrieval results of physical methods are often more accurate than those of semi-empirical methods because semi-empirical methods are more dependent on the regression training database. However, when the training database is optimal in the sense of representing various situation, the semi-empirical methods should be accurate as those from physical methods [9].

To meet the demand of fast, accurate LSE spectra retrieval from hyperspectral infrared data, this paper focuses on estimating the thermal infrared windows (8-14μm) emissivity spectra under clear sky conditions from the simulated IASI data. A practical physics-based regression method will be proposed with improved training and coefficients. Section 2 presents the methodology for LSE spectra retrieval from space. Section 3 describes the data used for this study. Section 4 gives some preliminary results. Finally, the conclusions are summarized in the last section.

2. Methodology

2.1 Radiative transfer equation

For a cloud-free atmosphere under local thermodynamic equilibrium, neglecting the scattering effects of atmosphere, the radiative transfer equation (RTE) in the thermal infrared region can be written as:

\begin{matrix} \begin{array}{l} R (v, θ) = ε (v, θ) B (v, T_{s}) τ (v, θ, p_{s}) - \int_{0}^{p_{s}} B (v, T_{p}) \frac{\partial τ (v, θ, p)}{\partial p} d p + \\ (1 - ε (v, θ)) \int_{0}^{p_{s}} B (v, T_{p}) \frac{\partial τ^{*} (v, θ, p)}{\partial p} d p \end{array} & , \end{matrix}

where R is the spectral radiance measured at the top of the atmosphere (TOA) with wavenumber ν and viewing zenith angle θ. For simplicity, the wavenumber and viewing zenith angle are ignored in the following expressions. ε is the effective LSE. B is the Planck function of the temperature. T_s and T_p are land surface temperature and atmospheric temperature at pressure level p, respectively. τ is the transmittance from the pressure level p to the TOA along the viewing angle. τ*(ν,θ,p) is the transmittance from surface p_s to the pressure level p. The subscript s denotes surface values. Equation (1) can thus be rewritten as follows:

\begin{matrix} \begin{array}{l} R (v, θ) = [B (v, T_{s}) τ (v, θ, p_{s}) - \int_{0}^{p_{s}} B (v, T_{p}) \frac{\partial τ^{*} (v, θ, θ_{s}, p)}{\partial p} d p] ε (v, θ) - \\ [\int_{0}^{p_{s}} B (v, T_{p}) \frac{\partial τ (v, θ, p)}{\partial p} d p - \int_{0}^{p_{s}} B (v, T_{p}) \frac{\partial τ^{*} (v, θ, p)}{\partial p} d p] \\ = M ε (v, θ) + N \end{array} & , \end{matrix}

where both M and N are regression coefficients and are regarded as the function of land surface temperature and atmospheric temperature and moisture profiles. Consequently, the spectral radiance at TOA could be thought as a linear function of the emissivity. Introducing a background mean condition and inverting Eq. (2), we can get:

δ ε = {(M^{T} M)}^{- 1} M^{T} δ T_{B} = K δ T_{B},

where superscript T denotes the transpose of matrix. Symbol δ is the perturbation with respect to an a priori estimated or mean condition. T_B is a brightness temperature (BT) vector corresponding to the radiance R at TOA. Obviously, K is the regression coefficients from BT of TOA to emissivity spectra.

2.2 Solution in the eigenvector domain

To reduce the uncertainties in the solution and take advantage of spectral correlations, the principal component analysis (PCA) technique has been suggested and widely used to represent the interested parameters [2]. Consequently, both the BT at TOA and the emissivity spectra are represented by several components, making only a few unknowns, i.e. coefficients of eigenvectors, need to be determined in the retrieval:

\begin{matrix} \begin{array}{l} δ T_{B} = \sum_{i = 1}^{N_{B T}} f_{i}^{B T} q_{i}^{B T} = Q_{B T} f_{B T} \\ δ ε = \sum_{i = 1}^{N_{ε}} f_{i}^{ε} q_{i}^{ε} = Q_{ε} f_{ε} \end{array} & , \end{matrix}

where q_i is the ith eigenvector, f_i is the associated coefficient, and N is the number of eigenvectors used. Q_BT and f_BT are the corresponding eigenvector matrix and coefficient vector for the BT of TOA, while, Q_ε and f_ε for the emissivity spectra. Combining Eqs. (3) and (4), one gets:

δ ε = Q_{ε} f_{ε} = K δ T_{B} = K Q_{B T} f_{B T} .

Because of the high correlation within the equations, the observation errors may have a great impact on the retrieval accuracy. To further reduce this effect, a ridge regression is introduced to stabilize the solution and to obtain reliable regression coefficients at the expense of losing small accuracies. Therefore, the solution of transformation coefficients K in Eq. (3) can be given by:

K = Q_{ε} f_{ε} {(Q_{B T} f_{B T})}^{T} {(Q_{B T} f_{B T} {(Q_{B T} f_{B T})}^{T} + r I)}^{- 1},

where r is the ridge parameter and I is the identity matrix. Small positive values of r will improve the conditioning of the problem. The coefficients establish the relationship between the perturbation of BT and that of the emissivity spectra with respect to the background mean condition. Once regression coefficients are known, the emissivity spectra departure from the a priori mean condition can be estimated from Eq. (3).

3. Data

Because there are few campaigns that measure the emissivity spectra and atmospheric profiles simultaneously in coincidence with satellite overpass, one possible way to develop a practical regression method is to use numerical simulation with emissivity spectra, atmospheric profiles, as well as the hyperspectral atmospheric radiative transfer models. Taking into account the feasibility, the ASTER Spectral Library, the MODIS UCSB (University of California, Santa Barbara) emissivity library, the Thermodynamic Initial Guess Retrieval (TIGR) database TIGR2000, and the atmospheric radiative transfer model 4A/OP (Operational Release for Automatized Atmospheric Absorption Atlas) [11] for hyperspectral infrared data were used to generate the simulated IASI data, to get the regression coefficients, and to evaluate the regression accuracies.

The ASTER and UCSB spectral libraries include spectra of various samples. The spectra of main materials (soils, vegetation, water, and snow/ice) of the terrestrial ecosystem covering from 8 to 14μm were chosen in this study. In total 173 spectral samples including 122 soil types, 32 vegetation types, 10 water types, and 9 snow/ice types from ASTER and UCSB spectral libraries were selected to develop this regression method.

The TIGR atmospheric profiles constructed by the Laboratoire de Meteorologie Dynamique (LMD) represent a worldwide set of atmospheric situations from polar to tropical atmosphere with total precipitable water (TPW) ranging from 0.1 to 8.0 g/cm². For LSE retrieval, only atmospheric profiles in clear sky are taken into account. Consequently, a profile selection procedure was performed to choose the cloud-free atmospheric situations. In this procedure, the profiles with relative humidity at any layer greater than 90% or at two consecutive layers greater than 85% were considered as cloudy [12]. Moreover, a subset of the remained profiles was further elaborately selected to insure that there was a nearly uniform probability distribution for TPW. At last, 348 atmospheric profiles were extracted from TIGR2000.

To well establish and evaluate the regression model, both the spectral samples and the atmospheric profiles were randomly divided into two data sets. One is the training data set used to establish the regression, and the other is the evaluating data set used to validate the accuracies. The training data set contains 86 emissivity spectra and 174 atmospheric profiles, while the evaluating data set contains 87 emissivity spectra and 174 atmospheric profiles. The general information about these data sets is displayed in Fig. 1 . Figure 1(a) shows the spectral variation of the emissivity spectra. Meanwhile, Fig. 1(b) and Fig. 1(c) display the histograms of the bottom atmospheric temperature (T_a) and TPW of the selected profiles. It is obvious that both training and evaluating data sets almost behave the similar properties and distributions for T_a and TPW. However, the vertical distribution of temperature and moisture for those atmospheric situations may be different.

Fig. 1 Characterization of variability and distribution in selected emissivity spectra and atmospheric profiles. The character ‘T’ denotes the training data set and ‘E’ the evaluating data set. (a) the surface emissivity spectra: mean emissivity (solid), mean emissivity ± 1 standard deviation (dashed); (b) the bottom atmospheric temperature, and (c) the TPW of profiles.

Download Full Size | PDF

4. Results

4.1 Determination of regression coefficients

To determine the regression coefficients, nine different viewing zenith angles (VZAs) (0°, 24.62°, 33.56°, 39.72°, 44.42°, 48.19°, 51.32°, 56.25°, 60°) were used in simulations by taking the angular dependence into account. Then, the VZAs and the atmospheric profiles were inputted into 4A/OP to estimate the spectral atmospheric parameters (transmittance, upwelling radiance, downwelling radiance). In the simulation, the sampling interval was set to 0.25cm⁻¹ and the spectral resolution was 0.5cm⁻¹ according to the configuration of IASI. Subsequently, the BT at TOA could be easily calculated with the atmospheric parameters, LSE and land surface temperature (LST). To make the simulation more representatives, the reasonable variations of LST were thus varied in a wide range according to T_a of the atmospheric profiles used. That is, LST varied from T_a−5K to T_a + 15K in steps of 5 K for T_a≥290K and from T_a−5K to T_a + 5K in steps of 5 K for T_a< 290K. During the simulation, Gaussian white noises were also added with the noise equivalent temperature differences (NEΔT) of IASI [13].

With the PCA technique, the first 16 components for emissivity spectra and 10 components for the spectral BT at the TOA are enough to describe the variations. The reconstruction errors of emissivity spectra were less than 0.01 and those of the BT at the TOA were less than 1K. Finally, the regression coefficients K can be obtained directly through Eq. (6).

4.2 Experiment results

To analyze the results obtained with regression coefficients, the retrieved emissivity spectra were compared with the selected emissivity spectra. Three parameters were calculated for each sample during the simulation. The first parameter is the mean spectral difference between the retrieved and selected emissivity in simulation. This parameter is expressed as:

b i a s = \frac{1}{N_{s}} \sum_{j = 1}^{N_{s}} (ε_{v} [j] - ε'_{v} [j]),

where N_S is the total number of sample during the simulation. ε_i[j] and ε'_i[j] are the selected and retrieved emissivity for sample j at wavenumber v. The second parameter is the root mean square error (RMSE) of the retrieved emissivity at wavenumber v, deðned as follows:

R M S E = \sqrt{\frac{1}{N_{s}} \sum_{j = 1}^{N_{s}} {(ε_{v} [j] - ε'_{v} [j])}^{2}} .

The third parameter is the overall RMSE_all for the selected spectral domain:

R M S E_{a l l} = \sqrt{\frac{1}{N_{D} N_{S} - 1} \sum_{v = 1}^{N_{D}} \sum_{j = 1}^{N_{S}} {(ε_{v} [j] - ε'_{v} [j])}^{2}},

where N_D is the total number of wavenumber in the spectral domain from 714 to 1250cm⁻¹.

Figure 2 shows the bias and the RMSE at each wavenumber for VZA 0°. Figure 2(a) gives the statistics for training data set, while Fig. 2(b) for evaluating data set. Obviously, the retrieved emissivity spectra have no bias at each wavenumber for both the training and evaluating data sets. However, the RMSE behaves wavelength-dependence. The RMSE is generally lower than 0.015 for spectral range 714-970cm⁻¹, while may be reach up to 0.03 for spectral range 970-1250cm⁻¹. The main reason may come from the fact that the variation in emissivity spectra is larger in that spectral range, which can be found in Fig. 1. It is worth noting the accuracy of retrieved emissivity spectra in 714-970cm⁻¹ is more accurate than that in 970-1250cm⁻¹, which means the retrieved emissivity spectra in 714-970cm⁻¹ may meet the required needs of LST retrieval and may be used as inputs to split-window algorithms to retrieval LST.

Fig. 2 Comparison between the selected and retrieved emissivity at each wavenumber for spectral region (8-14μm) at VZA 0°. The color solid is the mean spectral bias and the color dashed is the RMSE. (a) statistics for training data set, and (b) statistics for evaluating data set.

Download Full Size | PDF

Figure 3 gives the overall histogram of emissivity spectra difference between the retrieved and selected emissivity at VZA 0°. The RMSE_all of the retrieved emissivity spectra is about 0.02 for both training data set and evaluating data set. The similar performance on the training and evaluating data sets indicates the proposed method is robust and applicable. Because the data sets were generated by the world-wide library, the proposed method is promising at various atmospheric conditions.

Fig. 3 The overall histogram of emissivity spectra difference between the retrieved and selected emissivity.

Download Full Size | PDF

The experiment results at other VZAs are nearly similar to those at VZA 0°. Table 1 gives the statistical mean bias, RMSE, and RMSE_all between the retrieved and selected emissivity spectra for both training and evaluating data sets. Due to the small intervals of VZAs used in simulations, the regression coefficients for other VZAs might be linearly interpolated in function of the secant VZA.

Table 1. Statistical mean bias, RMSE, and RMSEall for both training and evaluating data sets.

View Table

5. Conclusion

In this work, a practical regression method to retrieve the hyperspectral emissivity spectra within the spectral range of 8-14μm has been proposed to balance the tradeoff between the retrieval speed and accuracy. ASTER Spectral Library, UCSB emissivity library, and TIGR database were used to generate simulated data that covered all the possible situations. The combination of PCA and ridge regression technique insured the robustness of the regression. Therefore, the proposed model may be suitable for recovering the emissivity spectra from world-wide hyperspectral thermal infrared data.

To validate the proposed method, the evaluating data set that were independent to the training data set were used. The results showed that the proposed method might be unbiased for various land surface and atmospheric conditions and be of robustness and speediness. The emissivity spectra retrieval accuracies under clear sky conditions were better than 0.015 and 0.03 for 714-970cm⁻¹ and 970-1250cm⁻¹, respectively.

All those preliminary work and results showed the proposed method is promising. Due to its simplicity in reality, this method might be easily used for all hyperspectral infrared sensors. Since this method was only evaluated in the simulated data, further studies are required to apply this method with the actual remotely sensed IASI data. Furthermore, how to improve the retrieval accuracy in spectral range 970-1250cm⁻¹ and how to use the retrieved emissivity spectra as a priori knowledge to direct the retrieval of LST and atmospheric profiles will be studied in the future.

Acknowledgments

This work was jointly supported by the “Strategic Priority Research Program” of the Chinese Academy of Sciences, Climate Change: Carbon Budget and Relevant Issues, Grant No. XDA05040204, by the National High Technology Research and Development Program of China under Grant No. 2009AA122102 and No. 2012AA121103, and by the National Natural Science Foundation of China under Grant No. 41101332 and No. 41231170. The authors would also like to sincerely thank NOVELTIS Inc. for providing the 4A/OP model, Jet Propulsion Laboratory for providing the ASTER spectral library, MODIS LST group for providing the UCSB spectral library, and the Laboratoire de Meteorologie Dynamique for providing the TIGR database.

References and links

1. A. R. Gillespie, S. Rokugawa, T. Matsunaga, J. S. Cothern, S. Hook, and A. B. Kahle, “A temperature and emissivity separation algorithm for advanced spaceborne thermal emission and reflection radiometer (ASTER) images,” IEEE Trans. Geosci. Rem. Sens. 36(4), 1113–1126 (1998). [CrossRef]

2. J. Li, E. Weisz, and D. K. Zhou, “Physical retrieval of surface emissivity spectrum from hyperspectral infrared radiances,” Geophys. Res. Lett. 34(16), L16812 (2007). [CrossRef]

3. B.-H. Tang, H. Wu, C. Li, and Z.-L. Li, “Estimation of broadband surface emissivity from narrowband emissivities,” Opt. Express 19(1), 185–192 (2011). [CrossRef] [PubMed]

4. Z. Wan, “MODIS land-surface temperature algorithm theoretical basis document, version 3.3” (NASA/GSFC, 1999). http://modis.gsfc.nasa.gov/data/atbd/atbd_mod11.pdf.

5. A. R. Gillespie, S. Rokugawa, S. J. Hook, T. Matsunaga, and A. B. Kahle, “Temperature/emissivity separation algorithm theoretical basis document, Version 2.4” (NASA/GSFC, 1996). http://eospso.gsfc.nasa.gov/eos_homepage/for_scientists/atbd/docs/ASTER/atbd-ast-05-08.pdf.

6. K. Ogawa, T. Schmugge, F. Jacob, and A. French, “Estimation of land surface window (8-12μm) emissivity from multi-spectral thermal infrared remote sensing - a case study in a part of Sahara Desert,” Geophys. Res. Lett. 30(2), 1067 (2003). [CrossRef]

7. J. Li and J. L. Li, “Derivation of global hyperspectral resolution surface emissivity spectra from advanced infrared sounder radiance measurements,” Geophys. Res. Lett. 35(15), L15807 (2008). [CrossRef]

8. L. Zhou, M. Goldberg, C. Barnet, Z. Cheng, F. Sun, W. Wolf, T. King, X. Liu, H. Sun, and M. Divakarla, “Regression of surface spectral emissivity from hyperspectral instruments,” IEEE Trans. Geosci. Rem. Sens. 46(2), 328–333 (2008). [CrossRef]

9. D. K. Zhou, A. M. Larar, X. Liu, W. L. Smith, L. L. Strow, P. Yang, P. Schlüssel, and X. Calbet, “Global land surface emissivity retrieved from satellite ultraspectral IR measurements,” IEEE Trans. Geosci. Rem. Sens. 49(4), 1277–1290 (2011). [CrossRef]

10. E. Péquignot, A. Chédin, and N. A. Scott, “Infrared continental surface emissivity spectra retrieved from AIRS hyperspectral sensor,” J. Appl. Meteorol. Climatol. 47(6), 1619–1633 (2008). [CrossRef]

11. L. Chaumat, C. Standfuss, B. Tournier, R. Armante, and N. A. Scott, “4A/OP reference documentation” (NOVELTIS, 2012). http://4aop.noveltis.com/sites/4aop.noveltis.loc/files/NOV-3049-NT-1178v4.3.pdf.

12. J. M. Galve, C. Coll, V. Caselles, and E. Valor, “An atmospheric radiosounding database for generating land surface temperature algorithms,” IEEE Trans. Geosci. Rem. Sens. 46(5), 1547–1557 (2008). [CrossRef]

13. F. Aires, W. B. Rossow, N. A. Scott, and A. Chedin, “Remote sensing from the infrared atmospheric sounding interferometer instrument. 1. compression, denoising, and first-guess retrieval algorithms,” J. Geophys. Res. 107(D22), 4619 (2002). [CrossRef]

VZAs		0°	24.62°	33.56°	39.72°	44.42°	48.19°	51.32°	56.25°	60°
Training data set	bias	0.000	0.000	0.000	0.000	0.000	0.000	0.000	0.000	0.000
	RMSE	0.019	0.020	0.020	0.021	0.021	0.021	0.022	0.022	0.023
	RMSE_all	0.019	0.020	0.020	0.021	0.021	0.021	0.022	0.022	0.023
Evaluating data set	bias	0.001	0.001	0.002	0.002	0.002	0.002	0.002	0.002	0.002
	RMSE	0.020	0.020	0.021	0.021	0.021	0.022	0.022	0.023	0.023
	RMSE_all	0.020	0.020	0.021	0.021	0.021	0.022	0.022	0.023	0.023

Practical retrieval of land surface emissivity spectra in 8-14μm from hyperspectral thermal infrared data

Abstract

1. Introduction

2. Methodology

2.1 Radiative transfer equation

2.2 Solution in the eigenvector domain

3. Data

4. Results

4.1 Determination of regression coefficients

4.2 Experiment results

5. Conclusion

Acknowledgments

References and links

Cited By

Figures (3)

Tables (1)

Equations (9)

Optics Express