Land Surface Temperature Retrieval Using Airborne Hyperspectral Scanner Daytime Mid-Infrared Data

1 University of Chinese Academy of Sciences, Beijing 100049, China; E-Mails: zhaoenyusdu@163.com (E.Z.); hongyuanh@gmail.com (H.H.) 2 Key Laboratory of Quantitative Remote Sensing Information Technology, Academy of Opto-Electronics, Chinese Academy of Sciences, Beijing 100094, China; E-Mail: caixiagao2010@hotmail.com 3 College of Geography and Planning, Ludong University, Yantai 264025, China; E-Mail: emails305@163.com 4 College of Geomatics and Geoinformation, Guilin University of Technology, Guilin 541004, China


Introduction
Land surface temperature (LST) is an important indicator for monitoring the changing of earth resources and one of the most critical parameters in the physical process of surface energy and water balance at local and global scales [1][2][3][4][5][6].The knowledge of LST plays a valuable role in urban climate, evapotranspiration, vegetation assessment, heat flux estimation, hydrological cycle, and environment studies [7][8][9][10][11][12][13][14][15].Therefore, it is necessary to find a reliable way to acquire LST in regional and global scales.With the development of remote sensing, infrared sensors, such as Moderate Resolution Imaging Spectroradiometer (MODIS), Advanced Very High Resolution Radiometer (AVHRR), and Spinning Enhanced Visible and Infrared Imager (SEVIRI), provide a valuable way for measuring LST over the entire globe.
To date, different methods have been proposed to retrieve LST from Thermal Infrared (TIR) remotely sensed data, such as the single channel algorithm and the split-window algorithm [16][17][18][19][20].By comparison, the study on LST retrieval from Mid-Infrared (MIR) data is under-developed because the radiance measured during daytime in the MIR spectrum contains both the surface emitted thermal radiance and the reflected solar radiance, which are equal in magnitude [21,22].Moreover, it is difficult to eliminate solar effects in the MIR spectrum during daytime because the separation of solar radiance from the total energy requires the accurate atmospheric information and the knowledge of the surface bidirectional reflectivity.However, the MIR data has its own advantages, such as higher detection sensitivity for high temperatures, higher atmospheric transmittance in the atmospheric window, and less sensitive to water vapor content (WVC) [22,23].Therefore, Sun and Pinker proposed a split algorithm, with three TIR channels and one MIR channel, to retrieve LST from SEVIRI data [24].The Visible Infrared Imaging Radiometer Suite (VIIRS) workgroup developed a dual split window day/night LST algorithm for 17 IGBP surface types by using two TIR bands (10.8 μm and 12 μm) and two MIR bands (3.75 μm and 4.005 μm), with a solar zenith angle cosine correction during the daytime [25].Still, LST retrieval only from the Airborne Hyperspectral Scanner (AHS) daytime MIR data has not been studied.In this paper, a new method is proposed to retrieve LST from AHS daytime data in two MIR channels (channel 66: 3.746~4.084μm and channel 68: 4.418~4.785μm) after reducing/eliminating the effect of direct solar radiance.
Section 2 describes the theory associated with the LST retrieval and the procedures for AHS data simulation.The details of direct solar radiance estimation and LST retrieval, as well as the sensitivity analyses, are presented in Sections 3 and 4, respectively.In Section 5, the proposed method is applied to AHS data, and the results are validated with in situ measurements.The conclusion is drawn in Section 6.

Basic Theory
Under local thermodynamic equilibrium during a clear-sky day, the radiative transfer equation (RTE) in the MIR region (3~5 μm) can be written as [15,26]: where Bi is the Planck function.Bi(Ti) is the radiance measured at the top of the atmosphere (TOA) in channel i, and Ti is the brightness temperature.εi and Ts are the surface emissivity and surface temperature, respectively.τi is the transmittance of the atmosphere from the ground to the TOA along the viewing angle.are the upward and downward solar diffusion radiances, respectively, which result from atmospheric scattering of the solar radiance.ρbi is the surface bidirectional reflectivity.
is solar radiance at ground level.In addition, , ⁄ , where Ei is the solar irradiance at TOA, , is the transmittance of the atmosphere from TOA to the ground along the solar angle, and and are the solar zenith and azimuth angle, respectively.Equation ( 1) can be rewritten as follows: where is the radiance after extracting the reflected solar direct radiance, and is the equivalent brightness temperature.

Data
To develop the LST retrieval method, the at-sensor radiances should be simulated.For this purpose, the MODTRAN 4.0 radiative transfer code has been used to predict the radiances for the AHS MIR channels (CH66 and CH68) in terms of the channel filter functions [27].In total, 705 atmospheric profiles, with the atmospheric bottom temperature (Ta) of 250~310 K and the WVC of 0.06~5.39g/cm 2 , extracted from the TOVS Initial Guess Retrieval (TIGR) database [28,29], are used to analyse atmospheric effects.The attenuation of the surface radiance has been considered by adding the uniformly mixed gases (CO2, N2O, CO and CH4) and ozone, included in the standard atmospheres of the MODTRAN 4.0 code, to the water vapor taken from profiles in the TIGR radiosoundings [30].To accomplish this, a previous classification of the surface temperature is made with the rule that the surface temperatures are from Ta − 5 K to Ta + 15 K with a step of 5 K. Furthermore, the VZAs are set to be 0°, 33.56°, 44.42°, 51.32°, 56.25°, and 60° (corresponding values of 1/cos(VZAs) are 1, 1.2, 1.4, 1.6, 1.8, and 2.0), respectively, so that 1/cos(VZAs) could be sampled with a step of 0.2.The SZAs are set as 0°, 25.84°, 36.87°,45.57°, 53.13°, and 60° (cos(SZAs) are 1, 0.9, 0.8, 0.7, 0.6, and 0.5), respectively, so that the cos(SZAs) could be sampled with a step of 0.1.Also, 70 different emissivities obtained from the Johns Hopkins University (JHU) Spectral library (soils, vegetation, and water, etc.) are considered.Once the simulations are made, TOA radiance could be determined according to Equation (1).In total, for the TIGR database and the JHU Spectral library, 8,883,000 different situations are simulated for retrieval.

LST Retrieval Method from Two AHS MIR Channels
The daytime MIR radiance contains not only the radiance emitted by the land surface and atmosphere but also the solar radiance reflected by the land surface and scatted by the atmosphere.To retrieve LST from MIR data, the direct solar radiance should be estimated first to eliminate the effect of solar radiance, and then a split window method should be developed to estimate the LST using the radiance ( ) after extracting the reflected solar direct radiance.

Estimation of Direct Solar Radiance
For the daytime MIR data, estimation of direct solar radiance is premise for LST retrieval because the radiance measured during daytime is strongly affected by the reflected direct solar radiance, which is coupled by the bi-directional reflectivity of the surface (ρbi), the solar radiance at ground level ( ), and the transmittance from ground to sensor ( ).As we all know, is related to WVC and VZA, while is related to WVC and SZA.Therefore, the relationship between direct solar radiance ( ) and WVC, VZA and SZA is investigated to estimate the direct solar radiance by assuming that the surface is Lambertian and the LSE is known.

Relationship between Direct Solar Radiance and WVC
To investigate the relationship between Di and WVC with the aid of simulated data, a scatter plot between Di and ln(WVC) is shown in Figure 1.The data is shown for CH66 and CH68 with different VZAs and SZAs at the LST conditions of 250~310 K and WVC of 0~5.5 g/cm 2 .Figure 1a,b shows the relationships at the six different VZAs when SZA = 0°, while Figure 1c,d shows those at six different SZAs when VZA = 0°.It is noted that Di and ln(WVC) can be fitted using a quadratic polynomial with a formula as Equation (3) with a correlation coefficient of 0.985.Similar results also can be obtained for other combinations of SZAs and VZAs.

(d) Direct Solar Radiance at Different VZAs
To analyse the effect of VZA on Di, Figure 2a,b express the relationships between coefficients a, b, c and 1/cos(VZA) at SZA = 0° and SZA = 60°, respectively.It is found that the coefficients a, b, and c can be fitted using the formulations of a = a1/cos(VZA) + a2, b = b1/cos(VZA) + b2, and c = c1/cos(VZA) + c2.The direct solar radiance can be described as a function of WVC and VZA as Equation ( 4) with a correlation coefficient of 0.992.

Estimation of LST
Based on the differential absorption (especially for WVC) in two TIR channels in 10~12.5 μm, a split-window method was improved for LST retrieval from TIR data by expressing LST as a linear function of the brightness temperatures Ti and Tj measured in the two adjacent TIR channels [14,15].In consideration of the similar RTEs in MIR and TIR without the influence of solar direct radiance, this paper extends the split-window method to the MIR spectral region for LST retrieval after eliminating the effect of direct solar radiance.The new method is expressed as follows: where ε = (εi + εj)/2, Δε = εi − εj, and k0, k1, k2, k3, k4, k5, and k6 are unknown coefficients, which can be derived from simulated AHS data.and are the TOA equivalent brightness temperatures in two MIR channels.εi and εj are the LSEs in channel i and j, respectively.ε is the averaged emissivity, and Δε is the emissivity difference between the two MIR channels.

Estimated Result of Direct Solar Radiance
The above analyses show that Di can be expressed as a function of WVC, SZA and VZA, and the fitting coefficients can be obtained using the simulated data (see Table 1).To evaluate the accuracy of

Coefficients of LST Retrieval Method
After eliminating the effect of direct solar radiance, WVC and LST are divided into several tractable sub-ranges to improve the accuracy of LST retrieval.WVCs are divided into five sub-ranges: [0-1.5],[1-2.5],[2-3.5],[3-4.5], and [4-5.5]g/cm 2 , and LSTs are divided into three sub-ranges: 265 K ≤ LST ≤ 295 K, 290 K ≤ LST ≤ 310 K, and 305 K ≤ LST ≤ 325 K [31].Then, the coefficients in Equation ( 6) can be obtained through a statistical regression method for each sub-range under different VZAs.As an example, Figure 5 displays the coefficients as functions of the secant of VZAs at the sub-ranges of LSTs, which vary from 305 K to 325 K, for the two WVC groups.The coefficients k0~k6 for other VZAs can be linearly interpolated as function of the secant of VZA.Similar results are obtained for the other sub-ranges.

Result of LST Retrieval
In practice, three steps are needed to retrieve LST.First, direct solar radiance is estimated with Equation ( 5).Second, approximate LSTs are estimated according to Equation ( 6) with the coefficients derived for the whole range of LST (provided that the sub-ranges of emissivity and WVC are known).Finally, more accurate LSTs are estimated using Equation ( 6) but with the coefficients k0~k6 of the LST sub-range that is determined according to the approximate LST. Figure 6 gives the RMSEs between the actual and estimated LST as functions of the secant of VZA for different sub-ranges.The RMSEs are shown to increase with the increase of VZAs.The RMSEs are less than 1 K for all sub-ranges; the minimum value is 0.16 K (LST = 305~325 K, WVC = 4~5.5 g/cm 2 , and VZA = 0°).

Sensitivity Analysis
LST retrieval accuracy is affected by several factors such as the instrumental noises, uncertainties of LSEs and atmospheric properties, and the uncertainty of the method.In this study, the instrument noises (Noise Equivalent difference Temperature, NEΔT), the uncertainties of LSEs and WVCs are taken into account.

Sensitivity Analysis to Instrumental Noises
The expected NEΔT in AHS channels 66 and 68 is approximately 0.33 K [32].To analyse the effect of NEΔT on the LST retrieval using Equation ( 6), a Gaussian random distribution error of 0.33 K is added to the TOA brightness temperatures in CH66 and CH68 in Equation (1).Also, the errors between the LSTs retrieved from the noise-added brightness temperatures and those determined from the noise-free brightness temperatures for the sub-ranges of LST varying from 305 K to 325 K in the five WVC sub-ranges are shown in Figure 7.The results show that the NEΔT of 0.33 K produces approximately an error within 0.63 K when LST is varying from 305 K to 325 K.

Sensitivity Analysis to LSEs
To analyse the effect of the LSE uncertainty, a Gaussian random distribution error, where the mean of the distribution is 0 and the standard deviation is 0.01, is added to emissivities εi and εj in Equation ( 6).Table 2 shows the effect of emissivity on the accuracy of LST retrieval at the condition of VZA = 0°.It is worth noting the LST retrieval errors (the LSTs retrieved from LSE-uncertainty-added conditions minus those determined from no-LSE-uncertainty conditions) vary from 0.4 K to 2.83 K by assuming that the uncertainty of the emissivity is 0.01; errors increase with the decrease of LSTs.The reason may be that the proportion of direct solar radiance is larger in the total radiance when the LST is lower, and the same emissivity error may produce a larger effect on TOA radiance with lower LSTs.Meanwhile, the retrieval accuracy will be increasing with the increase of WVC.The possible reason is that the impact of the atmosphere on LST retrieval is more significant when the WVC increases; thus, the same emissivity error would produce smaller errors when the WVC is larger.

Sensitivity Analysis to WVC
To investigate how significant is the effect of the uncertainty of the WVC on the retrieval of LST, a Gaussian random distribution error, where the mean of the distribution is 0 and the standard deviation is 10%, is added to WVC.The WVC error will affect the accuracy of direct solar radiance estimation and cause the wrong selection of sub-range retrieval coefficients.Figure 8a,b shows the histograms of LST errors caused by a WVC uncertainty of 10% under the conditions of VZA = 0°, LST = 305~325 K, and atmospheric WVC = 0~1.5 g/cm 2 (dry atmosphere) or 4~5.5 g/cm 2 (wet atmosphere).It can be seen from Figure 8 that RMSE and Bias are 0.21 K and 0.02 K under dry atmosphere, respectively, and 0.21 K and 0.047 K under wet atmosphere, respectively.
Considering the instrument noises ( sin cos • cos • cos sin • sin ), the uncertainties of LSEs ( ) and WVC ( ), and the accuracy of the algorithm ( ), the overall error on the LST ( ) can be described by the following: Inserting the values illustrated in Figure 6 for sub-ranges of LST = 305~325 K and LST = 265~295 K at all combinations of VZAs and WVCs, the overall error is shown in Figure 9. Figure 9 shows that the biggest RMSE is approximately 3.39 K for the sub-range of LST = 265~295 K and WVC = 0~1.5 g/cm 2 at VZA = 60°.The smallest RMSE is approximately 0.72 K for the sub-range of LST = 305~325 K and WVC = 4~5.5 g/cm 2 at VZA = 0°.

Data Processing
The AHS instrument has 80 spectral bands covering the visible and near infrared (VNIR), short wave infrared (SWIR), mid-infrared (MIR), and thermal infrared (TIR) spectral ranges.The instrument is operated by Instituto Nacional de Técnica Aerospacial (INTA) and it has been involved in several field campaigns since 2004.The arrangement of the AHS mid-infrared bands covered by 3.3~5.4µm has seven bands (i.e., CH 64 ~CH 70), and its bandwidth is 35 nm with λ/∆λ (minimum) of 110 [33].In this study, the proposed method is applied to AHS data collected on 4 July 2008 over a study area in the city and suburbs of Madrid covered by vegetation, water, buildings, roads, bare soil, etc.  Due to the influence of the attitude and speed of the spacecraft and the Earth's rotation, the AHS image would emerge a certain degree of geometric distortions, such as compression, stretching and offset compared to the actual position of the ground target pixels.Hence, to obtain the accurate location of the image pixels, the image geometric correction is carried out using the module of IGM in ENVI 4.7 with the aid of the geographic coordinate of each pixel.Then, the SZA for every pixel can be calculated using Equation (8) with the information of latitude and longitude, while the VZA of each pixel can be acquired according to the field of view (FOV) of the AHS sensor (90°) and the pixel numbers.sin( ) cos( ) cos( ) cos( ) sin( ) sin( ) ) where θ is the solar zenith angle.h is the hour angle at the local sidereal system.δ is the solar declination.φ is the local latitude.
Because it is difficult to obtain the ground measured LSEs in practice, the LSEs are determined using the supervised classification methods with the VNIR hyperspectral data equipped on the AHS sensor.Five categories, i.e., grass, water, roads, houses and bare soil, have been classified.The LSEs in AHS CH 66 and CH 68 are obtained by integrating emissivity spectral curves extracted from the JHU library with the spectral response function.In addition, WVC is extracted from the radiosonde atmospheric profiles acquired at the imaging time (approximately 0.76 g/cm 2 ).

Results and Validation
After data processing, the LSTs can be derived using Equation ( 6) from AHS MIR data on 4 July 2008 (N-S data 11:32 and W-E data at 11:16, both are at 2497 m with scan-rate 18.7 rpt, see Figure 10a,b).The selected surfaces for validation were (see Figure 11): (i) Water body, as cold target, in the lake at "El Retiro" park (40°25'1.65"N,3°41'2.65"W),(ii) Bare soil, as hot target, in the soccer field at UAM (40°32'52.44"N,3°41'48.45"W),(iii) Green grass, as cold target, in the Rugby field at UAM (40°32'51.71"N,3°41'54.33"W).The in situ measurements were carried out making transects at regular steps in representative and vast surfaces during the flight overpasses.Various instruments in the thermal infrared domain were used, including the CIMEL models CE312-1 and CE312-2, Heitronics KT 19.85 and NEC TH9100, whose technical specifications are shown in Table 3 [33].(a) (b) (c) For each validation area, the LSTs of the 2 × 2 pixels nearest to the geographic coordinate of the in situ measurement are selected.The estimated average LST of the four pixels is compared with the average of in situ measurements acquired within 5 minutes.The differences between estimated and measured LST for the three areas are approximately 0.7 K, 0.9 K, and 2.3 K (see Table 4).The vegetation area presents a larger LST retrieval error, which would be caused by its larger bias of measured LST (the bias is approximately 2 K, see Figure 12) due to the nonuniformity of the grass area covered by grass and soil; whereas, water and bare soil received better results because the compositions of water and bare soil are relatively uniform.(a) (b) (c)

Conclusions
In this study, a new algorithm for LST retrieval from daytime AHS data in two MIR channels, CH66 and CH68, is proposed.In this method, to retrieve LST more accurately, the atmospheric WVC and the LST are divided into several tractable sub-ranges.The coefficients of each sub-range in the algorithm are calculated using a statistical regression method from the numerical values simulated with the atmospheric radiative transfer model MODTRAN 4.0 under different atmospheric and surface conditions.The simulation analysis shows that the LST could be estimated by this algorithm with the RMSE less than 1 K for all of the sub-ranges with the knowledge of the LSEs.
In addition, the sensitivity analysis is performed in terms of the instrumental noise, the uncertainty of the LSE and WVC.The results show that, for the sub-range with VZA = 0° and LST = 305~325 K, an instrumental error of 0.33 K would produce a LST retrieval error of 0.6 K; the LST error caused by an LSE uncertainty of 0.01 varies from 0.4 K to 2.8 K.A LST retrieval error of approximately 0.2 K is caused by the WVC uncertainty of 10%.The largest total uncertainty is approximately 3.39 K for the sub-range of LST = 265~295 K and WVC = 0~1.5 g/cm 2 at VZA = 60°.The smallest one is approximately 0.72 K for the sub-range of LST = 305~325 K and WVC = 4~5.5 g/cm 2 at VZA = 0°.
Finally, the proposed LST retrieval algorithm is applied to the AHS data and is validated over three land surface types with in situ measurements.The validation results show that the differences between estimated and measured LST for water, bare soil and grass areas are approximately 0.7 K, 0.9 K, and 2.3 K, respectively.
where a, b and c are the fitting coefficients.Di is the direct solar radiance.
retrieval, Figure4a,b shows the histograms of the difference between the estimated and actual Di in CH66 and CH68, respectively.The root mean square errors (RMSEs) are 0.0123 W/(m

Figure 4 .
Figure 4. Histogram of the difference between the estimated and actual direct solar radiance (Di) for CH66 (a) and CH68 (b).

Figure 7 .
Figure 7. LST retrieval error caused by NEΔT when LST is varying from 305 K to 325 K.

Table 4 .
Validation results over three areas.