Land Surface Temperature Retrieval from Fengyun-3D Medium Resolution Spectral Imager II (FY-3D MERSI-II) Data with the Improved Two-Factor Split-Window Algorithm

: Land surface temperature (LST) is an essential parameter widely used in environmental studies. The Medium Resolution Spectral Imager II (MERSI-II) boarded on the second generation Chinese polar-orbiting meteorological satellite, Fengyun-3D (FY-3D), provides a new opportunity for LST retrieval at a spatial resolution of 250 m that is higher than that of the already widely used Moderate Resolution Imaging Spectrometer (MODIS) LST data of 1000 m. However, there is no operational LST product from FY-3D MERSI-II data available for free access. Therefore, in this study, we developed an improved two-factor split-window algorithm (TFSWA) of LST retrieval from this data source as it has two thermal-infrared (TIR) bands. The essential coefﬁcients of the TFSWA algorithm have been carefully and precisely estimated for the FY-3D MERSI-II TIR thermal bands. A new approach for estimating land surface emissivity has been developed using the ASTER Global Emissivity Database (ASTER GED) and the International Geosphere-Biosphere Program (IGBP) data. A model to estimate the atmospheric water vapor content (AWVC) from the three atmospheric water vapor absorption bands (bands 16, 17, and 18) has been developed as AWVC has been recognized as the most important factor determining the variation of AT. Using MODTRAN 5.2, the equations for the AT estimate from the retrieved AWVC were established. In addition, the AT of the pixels at the far edge of FY-3D MERSI-II data may be strongly affected by the increase of the optical path. Viewing zenith angle (VZA) correction equations were proposed in the study to correct this effect on AT estimation. Field data from four stations were applied to validate the improved TFSWA in the study. Cross-validation with MODIS LST (MYD11) was also conducted to evaluate the improved TFSWA. The cross-validation result indicates that the FY-3D MERSI-II LST from the improved TFSWA are comparable with MODIS LST while the correlation coefﬁcients between FY-3D MERSI-II LST and MODIS LST over the Mid-East China region are in the range of 0.84~0.98 for different seasons and land cover types. Validation with 318 ﬁeld LST samples indicates that the average MAE and R 2 of the scenes at the four stations are about 1.97 K and 0.98, respectively. Thus, it could be concluded that the improved TFSWA developed in the study can be a good algorithm for LST retrieval from FY-3D MERSI-II data with acceptable accuracy.


Introduction
Land surface temperature (LST) is an important variable determining the process of land-atmosphere interactions over the Earth's surface [1,2]. Accurate LST retrieved from free Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) data from 2000 to 2008. Use of the multi-year average bare soil emissivity map from ASTER GED to replace the bare soil emissivity might be a good way to improve the LSE estimation for LST retrieval [43]. In these years, many studies have been reported to improve LSE estimation for various satellite data, such as MODIS [43], FY-3B/VIRR [30], Landsat8 OLI [44,45], VIIRS [46], and Sentinel-3A SLSTR [47]. To minimize the effect of four percent invalid pixels in the ASTER GED dataset, Meng [44,45] developed a gap-filling approach to use the database to estimate LSE for the Landsat8 OLI data in 30 m resolution using GLC Land cover data. This encouraged us to develop a similar approach to improve LSE estimation for the FY-3D MERSI-II data at 1000 m resolution.
Atmospheric transmittance (AT) is another necessary input parameter for LST estimation. In the TIR wavelength, water vapor is the main effective factor for AT variation. In the previous research [13,37], AT was generally estimated with water vapor content with the assumption that VZA = 0 for all the pixels. However, this method mostly neglects that because of the vast swath of polar-orbiting images, the optical path of the pixels at the far edge increases significantly, compared with those close to the nadir. Moreover, this will lead to the overestimation of AT for the pixels with relatively large SZAs, especially for the complex atmospheric situation, e.g., a tropical atmosphere. For FY-3D MERSI-II data, the vast swath (2400 km) could also significantly enlarge the optical path of the pixels at the far edge. A viewing zenith angle (VZA) correction model was necessary to remove the effects of the optical path on the atmospheric transmittance estimation.
In this study the aim was to develop a practical algorithm with input parameters improved, in support of generating LST products from FY-3D MERSI-II data, which has two TIR bands. Due to the perfect mathematical derivation from the thermal radiance transfer equation, the two-factor split-window algorithm (TFSWA) developed by Qin et al. [22] was selected as the LST retrieval algorithm in the study. Using the ASTER Global Emissivity Database (ASTER GED) and the International Geosphere-Biosphere Programme (IGBP) data, we developed a new approach to estimate the required key parameter, i.e., land surface emissivity (LSE), for the FY-3D MERSI-II data. In order to determine the atmospheric transmittance for LST retrieval from the data, we developed a model to estimate atmospheric water vapor content (WVC) from the three atmospheric water bands (bands 16, 17, and 18) of the FY-3D MERSI-II data. Using the Thermo-Dynamical Initial Guess Retrieval (TIGR) profile database [48,49], we developed the atmospheric transmittance estimation. The vast swath (2400 km) of FY-3D MERSI-II data may significantly increase the optical path to the pixels at the far edge compared to those close to the nadir. A viewing zenith angle (VZA) correction model was developed to remove the effects of the optical path on the atmospheric transmittance estimation. Cross-validation with the MODIS LST product (MYD11A1) and validation with field measurements was made in the study in order to demonstrate the applicability of the improved TFSWA for LST retrieval from FY-3D MERSI-II data.
The remaining parts of this paper are presented as follows. Materials used and the study area for validation of the proposed LST retrieval algorithm are firstly presented in Section 2, which is then followed by Section 3 to outline the methodology for LST retrieval, including the ASTER GED-based LSE scheme, AWVC, atmospheric transmittance estimation, and the LST determination. The results and discussion are presented in Sections 4 and 5, respectively. Finally, the conclusion is drawn in Section 6.

The FY-3D MERSI-II Data
The FY-3D MERSI-II data have a total of 25 bands, but in the study, only 9 bands (Table 1) for the required parameter estimation towards the LST retrieval were used. The TIR bands 24 (centered at 10.8 µm) and 25 (centered at 12.0 µm) were the key bands for LST retrieval. LST retrieval with TFSWA from the two TIR bands required the knowledge of both AT and LSE. Therefore, other bands were also necessary, so that the estimation of both AT and LSE could be conducted. Specifically, bands 3 and 4 were used for LSE estimation and bands 15, 16, 17, 18, and 19 were used for AT estimation. These bands are in different pixel sizes and all of them were resampled to 1000 m spatial resolution with the nearest neighbor resampling. Figure 1 gives the spectral response functions of FY-3D MERSI-II and ASTER TIR bands.

The TIGR Profile Database
The TIGR database [48,49] was a widely-used atmospheric profile database established by the French Weather Dynamics Laboratory for various atmospheric applications [50]. The TIGR database contains 2311 atmospheric profiles selected from global atmospheric samples from 1968 to 1989. Generally, the atmospheric profiles with a relative humidity greater than 90% or relative humidity between two successive layers greater than 85% can be regarded as cloudy profiles [49,50]. Thus, those profiles that meet the abovementioned clear-sky conditions from the database were selected as the proposed clear-sky LST algorithm. Finally, 946 cloud-free profiles from the database were used, of which 236 ones belonged to the tropical atmosphere, 258 ones to the mid-latitude atmosphere, and 452 ones belonged to the polar atmosphere. The water vapor content of the cloud-free atmosphere profile varied from 0.06 to 6.27 g/cm 2 and the atmospheric temperature at the ground layer ranged from 231 K to 314 K. They represent various atmospheric conditions at different seasons around the world.

The ASTER GED Dataset
The ASTER Global Emissivity Database (ASTER GED) [43], created by the National Aeronautics and Space Administration's (NASA) Jet Propulsion Laboratory (JPL), California Institute of Technology, was used in the study for the LSE estimation.
Eight parameters are available in the database, including the mean and standard deviation of global emissivity, land surface temperature, NDVI, and ancillary products containing land/water masks, latitude, longitude, and ASTER GDEM. Referring to ASTER data, the parameters for five ASTER TIR bands centering at 8.3 µm, 8.6 µm, 9.1 µm, 10.6 µm, and 11.3 µm are available. The dataset is distributed as 1 • × 1 • tiles at~100 m and~1 km spatial resolution in HDF5 and binary formats [24,43]. The data were generated using all cloud-free ASTER data from 2000 to 2008 with the TES algorithm [23] in conjunction with a water vapor scaling atmospheric correction method [51,52]. Validation over a set of global sites indicated an average band error of~1% [52]. About 4% of pixels were found invalid probably due to incompletely expelling cloud contamination in high latitude regions [24,43].

The SURFRAD Measurements
Field measurements from four SURFRAD stations were applied in the study to evaluate the accuracy of the retrieved LST from FY-3D MERSI-II data. Established by the NOAA Office of Global Programs in 1993 to support satellite retrieval validation, modeling, climate, hydrology, and weather research for the United States [53], the SURFRAD network provides accurate, continuous, and long-term field measurements on the surface radiation budget every 1 min [54]. Seven stations are included and they are in a diverse climatological region ranging from the moist subtropical environment, the cool and dry northern plains to the hot and arid desert southwest, representing seven typical land cover types [54]. In this study, four SURFRAD stations, namely Bondvile (BND), Goodwin Creek (GWM), Penn State (PSU), and Sioux Fall (SXF), were selected for validation. Among them, the land cover type of the PSU station is crop, and the others are grass. Field LSTs were calculated from up-welling and down-welling long-wave radiation measurements of the SURFRAD stations through the Stefan-Boltzmann law as follows [55,56].

The Study Areas for Validation
Validation was necessary in this study. Two study areas are used to correspond two approaches of validation. The first validation was to use ground truth data from field measurements. Due to availability, we used the field measurements from 4 available SURFRAD stations in the USA to validate the improved TFSWA. Thus, for the validation, we needed to obtain FY-3D MERSI-II data covering the selected 4 SURFRAD stations in the USA and then compute the LST from the FY-3D MERSI-II data using the improved TFSWA. Comparison of the LST from field measurements with that from the retrieval would help us to understand the accuracy of the improved TFSWA for LST retrieval. The second approach of validation was cross-validation, i.e., to compare the retrieved LST from FY-3D MERSI-II data with the available satellite LST product. For this, we selected the MODIS LST product as a reference for the cross-validation. This is because the MODIS LST product has been widely applied to many studies and appreciated as one of the best satellite LST products in the world. Since FY-3D is a Chinese satellite for Earth observation, we selected Mid-East China as the region for the cross-validation with the MODIS LST product. The Mid-East China region has a latitude of 30 • 00 N to 45 • 00 N and longitude of 108 • 00 E to 123 • 00 E, including Henan, Shanxi, Hebei, and Shandong Province. Figure 2a shows the geographical location of the two study areas. Figure 2b gives the IGBP land cover types of the region in the USA and the four SURFRAD stations. Figure 2c shows the IGBP land cover types of the study area in China.

The Scheme for LST Retrieval
The scheme for LST retrieval is shown in Figure 3. Among three components, the left part is to estimate AT required in the LST retrieval, the right part is to estimate the LSE, and the central part is finally to retrieve LST. As shown in Figure 3, five FY-3D MERSI-II bands (15, 16, 17, 18, and 19) were used for the AT estimation, including the MODTRAN simulation, AWVC estimation, AT estimation, and the angle correction. Two FY-3D MERSI-II bands (3 and 4) were used in associated with the ASTER GED data to conduct the LSE estimation. In addition to these two essential parameters, the brightness temperatures of two TIR bands were also the required inputs for LST retrieval.

TFSWA for LST Retrieval
Qin et al. [22] developed a TFSWA for LST retrieval from NOAA-AVHRR bands 4 and 5, completely based on the mathematical derivation from the thermal radiance transfer equation. The TFSWA for FY-3D MERSI-II has the functions as follows: C 24 = ε 24 τ 24n (7) where T s is the retrieved FY-3D MERSI-II LST (K); A 0 , A 1 and A 2 are the internal parameters; T i is the brightness temperature (BT) of FY-3D MERSI-II TIR band i (i = 24 or 25) (K); τ i is the AT for band i; ε i is LSE for the band i; a i, and b i are the constants for the band i. The constants for the two FY-3D MERSI-II TIR bands are determined as follows: where L i is the derivative of the Planck function B i (T) with temperature T for the FY-3D MERSI-II TIR band i, which could be approximated as [B i (T + ∆T) − B i (T))]/∆T for determination of the constants as follows.
where a i and b i are the constants for the band i. It assumes L i with T from −50 • C to 50 • C for the two FY-3D MERSI-II TIR bands and then the constants were determined finally as follows: a 24 = −53.477, b 24 = 0.3951, 24; a 25 = −57.087, and b 25 = 0.4292. With these constants, it is able to retrieve LST from the two FY-3D MERSI-II TIR bands under the conditions that the two required parameters AT and LSE are estimated. This algorithm is just for the cloud-free pixels' LST estimation. Thus, the cloud detection product was used as the cloud mask layer to filter clear-sky pixels.

Estimation of AWVC
As seen in the previous section, AT is an essential parameter for LST retrieval with the improved TFSWA. AWVC is the most important factor affecting the variation of AT. Thus, the estimation of AT was generally converted from AWVC [4,13]. Three bands (bands 16, 17, and 18) of FY-3D MERSI-II were within the spectral range of atmospheric water vapor absorption (AWVA) so that it was able to estimate the WVC from these AWVA bands. It was found that the surface reflectance differed almost linearly with the wavelength for various ground surface types around the spectral range of AWVA [59,60]. Using this relationship, Kaufman et al. [60] proposed a ratio approach for the AWVC estimation for MODIS data. Following this approach, a similar method in the study was used to estimate AWVC.
Together with the other two non-absorption (TNA) bands (bands 15 and 19, close to the AWVA bands in spectral), the computation of the band ratio R i for AWVC estimation is as follows: where R i is the ratio of the AWVC band i to the TNA bands; ρ i is top-of-atmosphere (TOA) reflectance for the AWVC band i; ρ 15 and ρ 19 are TOA reflectance of the TNA bands 15 and 19, respectively; λ i represents the central wavelength of the AWVA band i, which can be obtained in Table 1; ω i1 and ω i2 are the weight of the AWVC band i for TNA bands 19 and 15, respectively. To establish this relationship, the atmospheric simulation program MODTRAN 5.2 was used in the study to simulate the TOA reflectance for the 3 AWVA bands with various atmospheric profiles having different AWVCs (0 g/cm 2 to 7.5 g/cm 2 ). To represent the actual situations in the real world, the atmospheric simulations for various conditions were performed, including different ground surface types (snow cover, forest, farm, desert, ocean, burnt grass, maple leaf, decayed grass, cloud deck, and old grass), rural aerosol model, 23 km visibility, and six standard atmospheres (US-1976, mid-latitude winter/summer, sub arctic winter/summer, and tropical atmosphere). Then, using the results from the simulation, a look-up table (LUT) showing the corresponding relationship between the ratio R i and the AWVC was developed. This LUT or the regression functions (16)-(18) from this LUT were then used to estimate the AWVC (W i ) for each of the 3 AWVA bands.
Since the estimated AWVC for each of the 3 AWVA bands might be slightly different as a result of its spectral limitation in accurately reflecting the actual AWVC, a weighted average approach was developed in the study to determine the AWVC for the pixel as follows [59,60].
W* = f 16 where η i are computed numerically from simulated curves of transmittance versus AWVC as the sensitivity of the tranmittance (τ i ) to the AWVC (W i ) in each band i (16,17,18), which can be calculated from |∆τ i /∆W i |. The estimated W* from Equation (19) is a combined two-way water vapor amount with effects of the optical path from viewing and sun-shining. For AT estimation, the vertical AWVC is needed. Based on the solar zenith angle (SZA) and the viewing zenith angle (VZA), the derived two-way water vapor amount W* is converted to the vertical column water vapor amount W using the following equation [59,60].
where W is the estimated AWVC of the pixel over the vertical column for AT estimation; VZA and SZA are the viewing zenith angle and solar zenith angle of the pixel, respectively.

Estimation of AT
An AT estimation from AWVC is required to know the relationship between AT and AWVC [22,35,37]. To establish this relationship, MODTRAN 5.2 again was used to simulate the change of AT with AWVC for the two FY-3D MERSI-II TIR bands under various conditions including the atmospheric profiles, viewing angles, and surface types. Generally, AT decreases with AWVC under various atmospheric conditions. The simulation was performed for each of the six available standard atmospheres with AWVC changing from 0.1 g/cm 2 to 7.0 g/cm 2 for various viewing angles. Results of the simulation were used to construct an LUT for AT estimation from AWVC at the pixel level.
Due to the obvious VZA difference among the pixels, a two-step approach was established in the study for AT estimation. Firstly, τ i (0), the AT at the nadir viewing with the AWVC (W) over the vertical column of the atmosphere for the two FY-3D MERSI-II bands 24 and 25, was estimated. This was done with the established LUT for AT estimation. The second step was to consider the effect of the pixel's viewing zenith angle on AT. As is well-known, the pixels apart from the nadir had a VZA greater than 0 compared with the pixels at the nadir. This would increase the viewing path through the atmosphere. Consequently, the AT of the pixel would decrease [35]. The effect of VZA on AT had to be considered for the pixels apart from nadir. Assuming the AT with the VZA = 0 was demonstrated as τ i (0) and the AT with VZA = θ was demonstrated as τ i (θ) for the two FY-3D MERSI-II TIR bands, to minimize this angular effect, the atmospheric simulation using MODTRAN 5.2 to establish the relationship between τ i (0) and τ i (θ) was performed. Since the FY-3D data had a scanning swath of 2400 km, and the pixels at the far edge part might have had a VZA of~55 • , the simulation was conducted with the VZA varying from 0 to 65 • . A total of 946 clear-sky atmospheric profiles from the TIGR database were used for the simulation. Figure 4 shows the results of the simulation for FY-3D MERSI-II TIR bands 24 and 25. The obvious difference as seen in Figure 4 was that the AT decreased remarkably with VZA. This was especially true when AWVC was very high or AT was quite low. An angular correction model was established in the study to involve the effect of VZA on AT estimation as follows [35]: where a 1i , a 2i , a 3i , b 1i , b 2i , b 3i , c 1i , c 2i , and c 3i are the coefficients of the angular correction models for band i, as shown in Table 2. Determination of the coefficients was done over the results from 26,656 simulations in MODTRAN 5.2 with 946 TIGR atmospheric profiles, six standard atmospheres, and 14 VZAs (the range is from 0 • to 65 • , the step is 5 • ). Table 2 clearly shows that correlation coefficients' square R 2 value were high up to 0.9999 and 1, implying that the fitting models are feasible for angular effect correction.

Estimation of LSE
Variation of LSE at the pixel scale is related to the surface types and their structures [39,61]. Under the condition that the LSE variations caused by internal reflection and distribution of natural surface [10,43,45] were not considered, LSE estimation could be done with a linear combination of the pixel's components, mainly including soil and vegetation [39,41].
where ε i is the emissivity for band i, ε vi and ε si are the vegetation and bare soil components for band i, and P v is the vegetation cover fraction calculated from NDVI [40]. Usually, such as spectral databases as the ASTER spectral dataset [62] and University of California at Santa Barbara (UCSB) emissivity library [63] were used to compute the ε vi and ε si in Equation (27) for LSE estimation [39]. Moreover, the computed ε vi and ε si from the databases were set as the same for all images. However, soil emissivity among different soil types may have great differences in magnitude. Consequently, this relatively wide range of emissivity variation among different soil types may cause significant effects on the accuracy of LST retrieval. To avoid these possible effects, an approach of estimating the soil emissivity from the ASTER GED product was proposed. The approach was composed of four steps: soil emissivity separation from ASTER GED product, IGBP based gap-filling for soil component emissivity, emissivity adjustment to FY-3D MERSI-II data, and the computation of LSE for the pixels of FY-3D MERSI-II data. Details are given as follows.

Estimation of Soil Emissivity from ASTER GED Product
ASTER GED product provides LSE for the entire land surface of the Earth. Viewing the grid of ASTER GED products as a mixture of vegetation and soil, it was able to separate the soil emissivity from the product. The separation can be done as follows [43]: where P v a is the vegetation cover fraction computed from ASTER GED mean NDVI data [64]; NDVI max and NDVI min are the maximum and minimum NDVI values in the ASTER GED product, which was obtained from 95% and 5% histogram intervals; ε vi a , ε si a , and ε i a are the emissivity of vegetation component, bare soil component, and mean emissivity of the mixed grid for ASTER TIR band i, respectively. In the study, the emissivity of the vegetation component was computed as an average emissivity by convolving with the ASTER spectral response function (SRF) with vegetation emissivity spectrum from the ASTER spectral database [62] and UCSB spectral database [63]. The vegetation emissivity of 0.8981 for the ASTER TIR band 13 and 0.983 for band 14 were adopted in this study.

Gap Filling of ASTER Soil Emissivity
As mentioned previously in Section 2.3, the ASTER GED product includes some pixels (~4%) that involve invalid data, probably due to the effects of long-prone cloud cover [43,64]. These pixels with invalid data should be filled with proper data for LSE estimation. The MODIS land cover product (MCD12Q1) was used in the study to fill the gaps of LSE for the pixels with invalid data. Among the five schemes of the MODIS MCD12Q1 product, the International Geosphere Biosphere Program (IGBP) global vegetation classification was selected to determine the soil component emissivity of mixed pixels. The gap was filled according to the soil LSE look-up table (LUT) of the IGBP classes shown in Table 3 where soil types correspond to specific IGBP land surface types [45]. Furthermore, soil emissivity of ASTER TIR bands 13 and 14 were obtained from the spectral response functions (SRFs), the ASTER spectral database [62] and the UCSB emissivity [63]. Mean and standard deviation values were calculated for soil LSEs of ASTER TIR bands 13 and 14. It was found that the soil LSE ranges from 0.954 to 0.973 for different soil types. This relatively broad dynamic range of vegetation LSEs indicates that it is not reasonable to use the same soil emissivity for LSE for each pixel of FY-3D MERSI-II.

Adjustment of ASTER Soil Emissivity to FY-3D MERSI-II TIR Bands
Firstly, the emissivity for the target spectral bands of MRESI-II and ASTER were estimated with the UCSB and ASTER spectral databases. Then, the least-square fitting method was used to establish the conversion functions (Equations (30) and (31)) of bare soil emissivity, from the bare soil LSEs of ASTER channel 13 ε s13 a and channel 14 ε s14 a to those of MERSI-II bands 24 ε s24 m and 25 ε s25 m . Figure 5 compares the estimated emissivity with the actual one computed from UCSB and ASTER spectral databases. The linear Pearson correlation coefficient R values for FY-3D MERSI-II TIR bands 24 and 25 were 0.9950 and 0.8734, respectively, and the RMSEs were 0.0006 and 0.0017, respectively. This excellent fitting performance indicates that the improved TFSWA could be used to estimate the bare soil emissivity for LST retrieval from the FY-3D MERSI-II data. ε s24 m = 1.038 ε s13 a + 0.032 ε s14 a − 0.069 (30) ε s25 m = −0.360 ε s13 a + 0.978 ε s14 a + 0.375 (31)

Calculation of the LSE for FY-3D MERSI-II
After the emissivity of bare soil ε si m for FY-3D MERSI-II TIR bands 24 and 25 was calculated from the bare soil emissivity of ASTER TIR bands 13 ε s13 a and 14 ε s14 a , the land surface emissivity ε i m of each FY-3D MERSI-II pixel could be further estimated with the mixed pixel emissivity model, i.e., Equation (32), provided that FY-3D MERSI-II vegetation emissivity ε vi m and the corresponding vegetation cover fraction P v m had been calculated from FY-3D MERSI-II NDVI data. The emissivity for FY-3D MERSI-II data can be calculated as follows: where ε i m is the LSE of FY-3D MERSI-II data; ε vi m is the vegetation emissivity used to estimate the LSE for FY-3D MERSI-II TIR bands 24 and 25, which has the value as ε v24 m = 0.982 and ε v25 m = 0.984, respectively; ε si m is the bare soil emissivity for LSE estimation. The FY-3D MERSI-II multi-year average bare soil emissivity (ε si m ) image could be obtained from the ASTER GED database's mean emissivity to replace the bare soil emissivity constant usually used in the VCM method, which can reflect the variation of bare soil emissivity in space, thus improving the emissivity inversion accuracy.

Validation of the Improved TFSWA
Validation is necessary to evaluate the developed algorithm and to understand the possible errors caused by the major factors involved. Three methods were used in the study for the validation, including atmospheric radiation transfer (RTE) simulation by MODTRAN, T-based validation using SURFRAD field measurements, and cross-validation with the MODIS MYD11A1 product.

Validation with Atmospheric Radiation Transfer Simulation
Since the derivation of the TFSWA from the radiance transfer equation (RTE) was completed with several approximations to the relevant parameters of the equations, possible errors of the TFSWA were unavoidable as a result of these approximations. One way to understand this possible error was to compare the LST calculated from the atmospheric radiation transfer simulation and the retrieved MERSI-II LST by the TFSWA and the simulated dataset [65,66]. In the study, MODTRAN 5.2 was used to conduct this atmospheric radiation transfer simulation for the FY-3D MERSI-II TIR bands. The sensed TIR radiance of the FY-3D MERSI-II TIR band i could be computed as follows: where L i toa , L ↑ , L ↓ and τ i are the TOA radiance, up-welling radiance, down-welling radiance, and AT of the FY-3D MERSI-II TIR band i, respectively; B i (T s ) is the Planck function to describe the emitted thermal radiance of ground as a blackbody with LST of T s for the band i; ε i is the LSE of the band i. The value of all these parameters for the specific conditions (such as AWVC and LST) used for the simulation could be obtained from the output of the simulation with the MODTRAN 5.2 program. The simulation was conducted with the LSE = 0.978 for band 24 and 0.983 for band 25 of FY-3D MERSI-II. The emissivities used here were calculated from the UCSB emissivity dataset and ASTER spectral library, which can represent mixed land types according to the actual estimated emissivity values using the MERSI-II LSE algorithm. Therefore, RTE LST (T s ) used for validation can be obtained from the inverse Planck function derived from Equation (33).
Since the altitude of the FY-3D satellite was 750 km, which was much higher than the TOA height, the converted brightness temperature (BT) from L i toa with Planck function could be viewed as the corresponding at-sensor brightness temperature for FY-3D MERSI-II TIR bands 24 and 25. Accordingly, the corresponding FY-3D MERSI-II LST could then be retrieved from the brightness temperature with the improved TFSWA in combination with the assumed LSE used for the simulation and the simulated AT (τ i ).

Validation with Field Measurements
In the study, field measurements from the available four stations in the central USA ( Figure 2b) were obtained from SURFRAD websites to validate the accuracy of the improved TFSWA. Land cover types of the four stations mainly included cropland (PSU) and grassland (BND, GWM, and SXF). The whole-year field measurement datasets of each station in 2020 were downloaded and the FY-3D MERSI-II LST of the region covering the four stations were retrieved. Using the latitude and longitude of the four stations, the pixel values were extracted from the retrieved FY-3D MERSI-II LST of each station from all the clear-sky scenes in 2020 for the validation. A total of 318 clear-sky FY-3D MERSI-II LST samples were used covering the whole year of 2020. Then, the field measurements of the LST of each station at the time closest to FY-3D MERSI-II overpass time were obtained and used to compare with the retrieved FY-3D MERSI-II LST. It was assumed that the field measurements of LST were the true LST of the station at the FY-3D MERSI-II overpass time. Despite the incompatibility between the field measurements of LST as a point and the FY-3D MERSI-II LST as a measurement from a pixel with a size of 1000 m × 1000 m, which is far from the point measurement of each station, the comparison of the difference between the retrieved LST and the measured LST still provides some general information about the accuracy of the TFSWA for LST retrieval from the FY-3D MERSI-II data in the real world. The field measurement of LST was calculated with Equation (1) given in Section 2.4.

Cross-Validation with MODIS LST Products
The MODIS LST product MYD11A1 has been widely used for various applications. Since the overpass time of the FY-3D satellite is close to that of MODIS, a comparison of the retrieved FY-3D MERSI-II LST with the MYD11A1 LST can also provide useful information about the accuracy of the TFSWA for LST retrieval from the FY-3D MERSI-II data. Different seasons and land surface types were considered for the comparison. The daytime FY-3D MERSI-II data in 2020 in Mid-East China were obtained for LST retrieval with the TFSWA. The daily MODIS LST products in 2020 in the same region of Mid-East China were collected for comparison.
Generally, cross-validation is done over single-day images from two sources. In our study, we computed the FY-3D MERSI-II LST for each day of the entire year 2020 and downloaded the daily MODIS LST product covering the same Mid-East China region for the same year. The retrieval of the LST from satellite images was only possible for the clearsky condition when thermal infrared radiance could transfer from the ground through the atmosphere to the sensor at space for remote sensing of LST. Therefore, our daily cross-validation was completed for only the clear-sky pixels.
In addition to the accuracy analysis of single-day LST data, the monthly average LST between MERSI-II LST and MODIS LST (LST MERSI-II − LST MODIS ) were also further calculated and analyzed in the study. This is because the monthly average LST can provide some more information about the general trend of spatial variation of the LST over the same region in comparison to the single-day LST image. Thus, cross-validation with the monthly average LST between the two LST datasets was also necessary. Since there was less cloud contamination in July 2020, the monthly average LST results of MERSI-II and MODIS for Mid-East China (Figure 2c) in this month were calculated.
In this study, for evaluating the difference of FY-3D LST and MODIS LST in different seasons of the year, MERSI-II and MODIS LSTs in January, April, July, and October 2020 were selected to represent the four seasons of the year, respectively. In addition, a series of conditions for observation time and angle to select proper points was used, to ensure that the pixel pairs selected from the two kinds of datasets had the most similar observation conditions. The local solar time is used for the observing time of the two LST datasets, and the clear-sky pixels in the two LST datasets with a time gap within 5 min were applied [67][68][69]; based on these selected pixels of two datasets, pixels with the difference between the tangent ratio of the view zenith angle less than 0.01 were regarded as having a similar atmospheric optical path [68,[70][71][72]. Figure 6a shows the scatter plots of the simulated RTE LST from MODTRAN and the retrieved LST with the improved TFSWA for VZA = 0, which have 30 samples (six atmospheric models, five water vapor contents, and VZA = 0). In addition, Figure 6b shows the scatter plots for the entire five VZAs, which have 150 samples (six atmospheric models, five water vapor contents, and five VZAs). Since the LST in different atmospheric models has a great difference, for example, 256.58 K in the sub-Arctic winter atmosphere and 299.20 K in the tropical atmosphere, scatter plots of the samples shown in Figure 6 are not evenly distributed but in a cluster for each atmospheric condition. All the sample dots in both Figure 6a,b are distributed around the 1:1 line, indicating that the retrieved LSTs with the TFSWA are very close to the RTE LSTs in all cases for the six atmospheric models. Despite this very good closeness between the simulated RTE LST and TFSWA LST for the six atmospheric conditions, slight differences between them still exist for different VZAs for a specific atmosphere condition. The statistical result of the difference indicates that~70% of the difference ∆Ts are within 1 K. The mean average value (MAE) shown in Figure 6a for 30 samples of VZA = 0 is 0.33 K, while the MAE in Figure 6b for a total of 150 samples of five VZAs is 0.73 K. These statistical results prove that the improved TFSWA is very good for LST retrieval from FY-3D MERSI-II data.  Figure 7 presents the validation results with the field measurement LST datasets. Scatter plots of the field measurement LST with the retrieved LST are shown in Figure 7a for Bondville station. Though the RMSE in the Bondville station is 1.97 K, most points in Figure 7a are on the non-bias line, implying the high accuracy between the retrieval LST and the field measurement LST. A similar situation is also seen in Figure 7b for Goodwin Creek station even though the retrieved LST is slightly underestimated for the field measurement LST and the RMSE increases to 2.70 K. This is probably due to the mismatch of the station in spatial location over the FY-3D MERSI-II images. Another reason may be the relatively spatial heterogeneous ground surface around the station. Relatively high RMSE is also observed in Figure 7c for PSU station where ground heterogeneity is high. Figure 7d    day was also about 13:30 h, indicating they were comparable for the LST variation. According to the statistical analysis, LSTs in the study area range from 268.6 K to 309.30 K. In addition, from Northwest to Southeast, the LSTs increase evidently. The LSTs in the Southeast near the ocean are higher, mainly above 290 K. In comparison, those in Northwest China are mostly less than 290 K. Some relatively lower LSTs below 273.15 K are mainly distributed in Shanxi and Hebei province where Taihang Mountain, Lvliang Mountain, and Yanshan Mountain are located. Furthermore, there is quite a similar trend in the MERSI-II and MODIS LST images. In addition, as shown in Figure 8c, we analyze the difference (LST MERSI-II − LST MODIS ) between the two LST maps. According to the statistical analysis, for the LST difference LST image in Figure 8c, the RMSE, bias, STD, and R 2 are 1.74 K, 0.92 K, 1.47 K, and 0.96, respectively. For the frequency histogram graph shown in Figure 8d, the LST difference follows the normal distribution. In Figure 8, some pixels with invalid data result from cloud contamination distributed in both LST images. Because of the different transit times of the two sensors, the spatial distribution of invalid data in the two images are not the same for different cloud coverage.  Figure 9 shows the monthly average LST of FY-3D MERSI-II (a), MODIS (b), the LST difference between them (c), and the frequency histogram of LST difference (d). As shown in Figure 9a,b, the LST distribution in the MERSI-II and MODIS monthly average LST is highly consistent in the spatial distribution, and the difference also follows the normal distribution in Figure 9d. For the LST difference distribution in Figure 9a,b, relatively low LSTs are mainly distributed in Northern Hebei, Shanxi, and Eastern Henan where the mountains are located. In the meantime, it can be seen from the LST difference image in Figure 9c that the negative LST difference value mainly appears in the same location. According to the statistic result, for the difference of the monthly average LST in July 2020, the RMSE, bias, R 2 , and STD are 2.97 K, −1.58 K, 0.84, and 2.52 K, respectively.  Figure 10 shows the scatter plots between MERSI-II and MODIS LST in January, April, July, and October 2020. According to the statistical analysis in Figure 10, there are significant differences in LST for different seasons. The R 2 range from 0.84 (July) to 0.96 (January); the RMSE (2.52 K), MAE (1.96 K), and standard deviation (2.46 K) are relatively large in July, while in January, the RMSE (1.57 K), MAE (1.21 K) and standard deviation (1.57 K) are relatively low.

Cross-Validation over Different Land Surface Types
To clarify the different land cover types of MERSI-II, the variation of surface temperature accuracy of FY-3D MERSI-II with different types of surface coverage is analyzed. There are four types of land cover in the study area: deciduous broad-leaved forest (DBF), woody grassland (WS), grassland (GL), and farmland (CL). Figure 11 shows the correlation coefficient (R 2 ), root mean squared error (RMSE), mean absolute error (MAE), standard deviation (STD), and other precision evaluation indexes of FY-3D MERSI-II LST over different land cover types in the summer of 2020 (June to August). Figure 11 shows the scatter plots of FY-3D MERSI-II and MODIS LST in the summer of 2020 (June to August). According to the chart, the land cover type has a significant impact on LST. The correlation coefficients of the four types of land cover are all relatively high, from 0.93 of the woody grassland to 0.98 of the other land cover types. In addition, the RMSE (1.98 K) and STD (1.85 K) of deciduous broad-leaved forests are the smallest, while the RMSE (2.87 K) and STD (2.87 K) of woody grasslands are the largest. In addition, the deciduous broadleaf forest has a minor MAE (1.46 K) and a more significant correlation coefficient (0.98), while the woody grassland has a more significant MAE (1.74 K) and the smallest correlation coefficient (0.93).

Discussion
The accuracy of the TFSWA algorithm is mainly influenced by the AWVC, transmittance, and LSE. The BT, LSE, and transmittance changed with different VZAs, for the pixels with extreme VZA that will induce more error. In the research, transmittance angular correction was performed carefully to minimize part of the angle effects in MERSI-II LST retrieval. Some additional factors should also not be neglected. Some errors were induced in the validation using ground data. As described above, the error for the SURFRAD crop and grass sites during daytime was about −0.3-0 • C [58]. This is another non-negligible error source for validation. Cloud detection errors inevitably exist in the cloud mask layer, especially when thin clouds exist. Moreover, some pixels with lousy quality and strips are inside in MERSI-II data. In this study, all test points for in-situ validation were checked carefully to exclude undetected clouds or bad-quality pixels to minimize the effect of cloud contamination and image quality in the validation work. In addition, a series of potential relating factors that may lead to LST errors, such as soil water content, surface air, humidity, wind speed, air temperature, should be analyzed in the near future [32].
To validate the accuracy and adaptability of the TFSWA algorithm, the MERSI-II LST with MYD11A1 LST for different seasons and land surface types were compared. According to the single-day and single-month LST comparison in Section 4.3, the LST difference along with mountains was more significant. This may be because the spatialtemporal variation of surface humidity in the rugged areas caused more uncertainty of cross-validation results. [68,71,72] As described in the monthly cross-validation, one can find that the LSTs in July had relatively low accuracy, the maximum standard deviation (2.46 K) existed in July and the minimum one (1.57 K) is in January. These results are mainly consistent with the previous studies [68,70,72]. This result may be due to the considerable spatial variation of LST because the difference of LST between bare soil and vegetation is more significant in July [72]. The influence of land surface emissivity may cause the difference. Moreover, the error of LSE may lead to more considerable differences between the two LST data sources in July. Another reason for the significant standard deviation of July LST may be the increasing influence of atmospheric change. The monsoon climate in Central and Eastern Asia is warm and humid, and the atmospheric water vapor content increases, which leads to the increase of LST error to a certain extent.
For the cross-validation results of various land cover types, the LST accuracy of the deciduous broad-leaved forest was relatively high than that of other land surface types, while the accuracy of woody grassland was the lowest. The reason for these results is mainly that the spatial structure of the evergreen broad-leaved forest is relatively uniform. At the same time, the woody grassland usually has more significant spatial heterogeneity, and the LSTs observed are easily affected by observation and illumination conditions [72]. Therefore, for land cover types with substantial spatial homogeneity, such as evergreen broad-leaved forest and grassland, the difference of LST retrieved from the two sensors is relatively small. In contrast, for rugged land surfaces, especially for poor or sparse vegetation land cover, the difference is relatively considerable. Besides, about three pixels of geometric positioning difference can be found between the MODIS and FY-3D MERSI-II data; the geographical position error is another influencing factor that should not be negligible for cross-validation at the 1000 m scale.
A limitation of this LST estimation scheme that should be noticed is that, the transmittance estimation models were simulated according to the standard atmospheric models inside MODTRAN, which was decided by the latitude and observing date of the image. For example, for the LST retrieval in the USA, the 1976 USA atmospheric model was used. However, for the global LST estimation, how to properly use these standard atmospheric models is an issue to address. In the next step, we plan to explore transmittance models that can be used for global LST estimation from numerous representative profiles. In addition, in the process of LSE estimation shown in Section 3.5, we calculated the NDVI max and NDVI min values using 95% and 5% percentiles for 1 km × 1 km tiles, while for different tiles of ASTER GED datasets, the NDVI max and NDVI min values are different. This could lead to some strips in the mosaic edges of ASTER GED soil emissivity tiles, which will result in some uncertainties in the process of LSE estimation. In the meantime, to derive soil emissivities with Equations (28) and (29), we applied average vegetation emissivities from ASTER and USCB databases in channels ASTER 13 and ASTER 14. Thus, the unique calculated values can have a significant impact on the soil emissivities. Although the dynamic range of vegetation specific emissivity is relatively small compared with soil specific emissivity, the soil specific emissivity can be underestimated or overestimated in certain areas where the vegetation specific emissivity is higher or lower than the average. In the near future, it should also be considered to develop an algorithm to obtain the vegetation emissivity distribution map to substitute the average values to further minimize the influence of mean vegetation LSE.

Conclusions
FY-3D MERSI-II, the second generation of the Chinese meteorological satellite for Earth observation, has two TIR bands for LST retrieval, and hence may have great potential for many environmental applications such as agricultural drought monitoring. In this study, an improved TFSWA has been developed for LST retrieval from FY-3D MERSI-II data. Essential coefficients of the TFSWA were carefully and precisely estimated according to the spectral features of the FY-3D MERSI-II TIR bands for LST retrieval. The ASTER GED and IGBP products were used to estimate the key parameter LSE required for LST retrieval with the improved TFSWA. The AT estimation equations and band ratio LUT for AWVC estimation were established based on the simulation database obtained from MODTRAN 5.2 and TIGR atmospheric profiles. The swath of FY-3D MERSI-II data is 2400 km, which may significantly increase the optical path of the pixels apart from the nadir and such an increase in optical path may perform an obvious effect on the magnitude of AT for TIR remote sensing. A viewing zenith angle (VZA) correction equation was proposed in the study to consider this effect on AT estimation for precise LST retrieval from the FY-3D MERSI-II data.
Validation is very necessary for the improved TFSWA developed in the study. Three approaches were adopted to assess the applicability of the improved TFSWA: simulated RTE LSTs, field data, and MODIS LSTs. Compared with the RTE LSTs simulated by MODTRAN 5.2, the improved TFSWA had an LST estimation error MAE of <1 K, implying the applicability of the improved TFSWA for LST retrieval from the FY-3D MERSI-II data. Validation with 318 field samples from the available four stations of the SURFRAD network indicated that the average MAE and R 2 of the scenes at the four stations were about 1.79 K and 0.98. Therefore, it could be concluded that the improved TFSWA developed in the study can be a good algorithm for LST retrieval from FY-3D MERSI-II data with acceptable accuracy. Cross-validation with MODIS LST products MYD11A1 was also completed over the Mid-East China region for different land surface types and seasons. Results from the cross-validation showed that the LST accuracy varied considerably in different seasons and over various land cover types. According to the statistical analysis, the land cover type had a significant impact on LST. The LST accuracy of the deciduous broad-leaved forest was relatively higher than that of other land surface types, while the accuracy of woody grassland was the lowest. For the cross-validation results of monthly average LSTs, the standard deviation values range from 1.57 K (in January) to 2.46 K (in July). Correlation coefficients between FY-3D MERSI-II LST and MODIS LST over the Mid-East China region were high up to 0.84~0.98 in different seasons and land cover types, implying that the MERSI-II LST generally had a high consistency with the MYD11 LST. The high consistency proves that the improved TFSWA has an acceptable accuracy for LST inversion of FY-3D MERSI-II data.
Author Contributions: W.D. and Z.Q. wrote the manuscript and worked on the methodology. J.F. collected and preprocessed the FY-3D MERSI-II data. C.Z., Q.H. and K.C. participated in the processing of LSE, the ground measurements, and the cross-validation. K.C. and B.A. gave support on TIGR data processing. Z.Q., J.F. and Q.H. provided suggestions for the correction of the final manuscript. All authors have read and agreed to the published version of the manuscript.