Evaluation of MODIS LST Products Using Longwave Radiation Ground Measurements in the Northern Arid Region of China

This study presents preliminary results of the validation of the Moderate Resolution Imaging Spectroradiometer (MODIS) daily LST products (MOD/MYD11A1, Version 5) using longwave radiation ground measurements obtained at 12 stations in the North Arid and Semi-Arid Area Cooperative Experimental Observation Integrated Research program. In this evaluation process, the broadband emissivity at each station was obtained from the ASTER Spectral Library or estimated from the MODIS narrowband emissivity Collection 5. A comparison of the validation results based on those two methods shows that no significant differences occur in the short-term validation, and a sensitivity analysis of the broadband emissivity demonstrates that it has a limited effect on the evaluation results. In general, the results at the 12 stations indicate that the LST products have a lower accuracy in China’s arid and semi-arid areas than in other areas, with a mean absolute error of 2–3 K. Compared with the temporal mismatch, the spatial mismatch has a stronger effect on the validation results in this study, and the stations with homogeneous land cover have more comparable MODIS LST accuracies. Comparisons between the stations indicate that the spatial mismatch can increase the influence of the temporal mismatch. OPEN ACCESS Remote Sens. 2014, 6 11495


Introduction
Land surface temperature (LST) is a key parameter in climatological and environmental studies [1,2].Obtaining LST data is important considering its common use in environmental studies and resource management.Measuring LST from ground-based instruments at the regional and global scales is practically impossible; however, the use of satellites in the thermal infrared (TIR) region is a viable option [3].Evaluating satellite LST retrievals is critical to their application, and the feedback from validation activities also helps improve the generation of these products [4,5].LSTs are very difficult to validate because of their large spatial and temporal variation, particularly during the daytime.Therefore, the careful selection of validation sites is essential [5][6][7].
Multiple daily LST products (on the global scale) are generated by the science team of the Moderate Resolution Imaging Spectroradiometer (MODIS) onboard the NASA Terra and Aqua Earth Observation System (EOS) satellites.MODIS LST products have been used in various studies since their release.However, the errors caused by LSTs have been mostly disregarded in these studies.Although MODIS LSTs provide a potentially inexpensive means of validating and improving existing land surface and climate models, these products were often ignored by the modeling community until recently.The major concern regarding the use of these LST products is that their accuracy has not been adequately assessed on a global scale, although some validation work has been conducted in previous studies [6,8,9].
Surface longwave radiation is related to LST and emissivity [10,11].Recently, high-quality long-term longwave radiation measurements have become globally available, such as measurements from FLUXNET [12,13], Surface Radiation Budget Monitoring (SURFRAD) [14,15], Baseline Surface Radiation Network (BSRN) [16,17], and Atmospheric Radiation Measurement (ARM) [18][19][20].It is helpful to investigate the efficacy of these measurements when evaluating satellite LST products.MODIS Collection 4 LST/emissivity products were evaluated using longwave radiation measurements collected at FLUXNET sites [21].It was concluded that the MODIS Collection 4 LST/emissivity has an obvious negative bias as high as -3 K. Substantial improvements have been made to MODIS Collection 5 LST/emissivity products [22].In this study, the MODIS Collection was evaluated using the longwave radiation measurement data from observation sites in the northern arid region of China.Additionally, the effects of the parameter emissivity and the spatial and temporal mismatches are discussed.

Study Area
China has large arid regions.In these regions, the LST is a sensitive factor to the local ecology and environment.Therefore, the study area is selected in the northern arid region of China, and the validation stations (shown in Table 1) are mainly located in the eastern part of the region (see Figure 1).The arid region in northern China covers approximately 4.0 × 10 6 km 2 , which is over 30 percent of the total land area in China.The region includes most of the Xinjiang Uygur Autonomous Region, the Inner Mongolia Autonomous Region, and the Ningxia Hui Autonomous Region, as well as parts of Gansu, Shanxi, and Hebei provinces.The study area spans 35°N to 50°N and a 50° longitudinal swath.The temporal change across the region is over two hours.In this study, the arid and semi-arid regions were defined according to the study of Song et al. [23], as shown in Figure 1.The data are presented for the months of July, August, and September each year listed in the last column.
Figure 1 shows the location of each site.

Method
The methodology of this study consisted of two parts: evaluation of the MODIS LST MOD/MYD11A1 using the ground-based LSTs obtained from the station's longwave radiation data, and evaluation of the LST spatial heterogeneity of the validating stations.The study flowchart is shown in Figure 2. The detailed methodology is as follows.

LST Retrieval from the Ground Measurements
Longwave radiation ground measurements were obtained from the stations (Figure 1).Based on thermal radiative transfer theory, the upwelling longwave radiation at the surface level depends on the LST, emissivity, and downwelling longwave radiation [10]: Where L↑ is the surface upwelling longwave radiation, L↓ is the surface downwelling longwave radiation, εb is the surface broadband emissivity, δ is the Stefan-Boltzmann constant (5.67 and TS is the land surface temperature.Therefore, the ground-measured surface temperature can be obtained according to Equation (2): ( ) In Equation ( 2), εb is a key parameter because L↑ and L↓ are observational data parameters, and δ is a constant.In this study, 3-14-µm broadband emissivity values are assumed to be equal to the emissivity (εb) of the entire longwave range of Earth's natural surface material, which peaks at 9.7 µm according to Wien's displacement law.The error associated with this assumption can be ignored [21].The 3-14-µm broadband emissivity for each site was derived in two ways.One method involved retrieving the information from the spectra of the ASTER Spectral Library [24], and the other method involved estimating the information from the MODIS narrow emissivity Collection 5 in the thermal infrared region [11].

Evaluation of Results
In the first method, the emissivity at each site was derived using the spectra from the ASTER Spectral Library, based on its land cover type.In support of the Advanced Spaceborne Thermal Emission Reflection Radiometer (ASTER) on NASA's Terra platform, the ASTER Spectral Library was compiled and made available from the website [25].The library is a collection of contributions in a standard format with ancillary data from the Jet Propulsion Laboratory (JPL), the Johns Hopkins University (JHU), and the United States Geological Survey (USGS).The latest version of the library, Vision 2.0, is available online or via CD, and it includes the spectra of minerals, rocks, lunar and terrestrial soils, manmade materials, meteorites, vegetation, snow, and ice; in total, more than 2300 types of spectrums are available [26].The library is one of the most comprehensive collections of spectra.
In the second method, the MODIS narrowband emissivity derived from the MODIS day/night LST algorithm was used to estimate the broadband emissivity parameter.Following Ogawa et al. [27], a multiple regression was obtained, as shown in Equation (3): where W ε is the broadband emissivity and ε29, ε31, and ε32 are the narrow emissivities of MODIS bands 29, 31, and 32, respectively.a, b and c are coefficients and statistical parameters that demonstrate the suitability of the regression.
Wang et al. [11] demonstrated that there are no significant differences in the coefficients of different types of materials, except for ice, water, and snow, although different materials influence the coefficients.Therefore, the use of one set of coefficients to calculate the broadband emissivity will not result in a significant error.Considering the diversity of land covers of the sites, a set of coefficients that are suitable for all the data was used to calculate the broadband emissivity [24], as demonstrated in Equation (4): 0.2122 0.3859 0.4029 The satellite-estimated broadband emissivity was compared with three-year ground-based emissivity measurements at Gaize (32.30°N, 84.06°E, 4420 m) in the western Tibetan Plateau; the result shows that the broadband emissivity calculation from MODIS narrowband emissivities agrees with the ground measurements quite well, with a standard deviation of 0.0085 and a bias of 0.0015 [11].To match the 1-km spatial resolution LST products, the LST/Emissivity products with a 5-km spatial resolution at each site were resampled to 1-km resolution.Then, the broadband emissivity at each station was calculated following Equation ( 4) by resampling the multiple narrow emissivity data.

Spatial Heterogeneity Analysis
Because of the stable performance and larger footprints of the pyrgeometers equipped at the validation stations, the longwave radiation ground measurements are more suitable for validating the MODIS LSTs than the other direct ground-measured LST data, such as the infrared radiometer SI-100 series.However, the pixels of MODIS MOD11A1 LSTs, with a 1-km resolution, are large scale compared with those of the four-component radiation sensor, and the differences caused by the spatial mismatch effect at different sites warrant further discussion.In this study, the semi-variance analysis method was used to analyze the effects on the validation results based on the higher resolution land surface temperature data retrieved from Landsat Thematic Mapper (TM) images.
Matheron [28] synthesized a geographical statistical theory, namely, the regionalized variable theory, from some of the scattered statistical application results.The conventional variogram of the geostatistics, defined formally in Equation ( 5), can distinguish the spatial variation in regions within a magnitude of one or two spatial scales simultaneously and can describe the randomness and structure of the regionalized variable.We assume that Z is a realization of a two-dimensional random process, Z(x), using the following variogram: where Z(x) and Z(x + h) are the random variables at positions x and x + h separated by vector h, with a spatial lag in both distance and direction.In this context, we consider the land surface temperature, which is scattered from the surface grid of the TM land 2 × 2-km image, as regionalized variable Z(x), in which lag h is treated as the distance only.Then, we can estimate the variogram as follows: Based on Equation ( 6), the relevant γ°(h) was calculated according to the distance variable h.Taking h as the x-axis and γ°(h) as the y-axis, the variogram can be drawn.The diagram directly shows the characteristics of the spatial variation of the regionalized variable Z(xi), which is the land surface temperature in this study.Therefore, we can analyze the properties of the spatial variable and the structure of the land surface temperature variable, which indicates the land surface temperature heterogeneity at the validation sites.
Four important parameters are included in the variogram diagram: the nugget (C0), sill (C + C0), structural variance ratio (C0/(C + C0), and range (A0).The γ°(h) function increases from a nonzero value to a relatively stable constant value that varies with h.When h = 0, the γ°(h) value is a nonzero value C0, namely, γ°(0) = C0 and C0 is the nugget of the variogram.The stable constant value is the sill (C + C0) of the variogram.When γ°(h) reaches the sill, the value of the variable h is A0, namely, the range of the variogram.To analyze the spatial heterogeneity of the regionalized variable LST, C0 was employed to indicate the level of Z(x) randomness, which may be caused by the internal variation in Z(x) over a smaller distance h than the sampling distance, or it may be derived from the sampling error.C + C0 can represent the largest extent of the regionalized variation.C0/(C + C0) represents the ratio of the nugget to the total spatial variation, and it can be used to weigh the randomness of the regionalized variable.A0 is the mean distance over which the spatial autocorrelation of Z(x) exists.

MODIS Data
As a part of the NASA EOS project, two MODIS instruments were placed on the Terra and Aqua satellite platforms to provide information for global studies of atmosphere, land, and ocean processes [29].Aqua passes overhead at approximately 1:30 p.m. (ascending mode) and 1:30 a.m.(descending mode) local time, whereas Terra passes at 10:30 a.m.(descending mode) and 10:30 p.m. (ascending mode) local time.The MODIS instruments have distinct advantages in terms of their global coverage, high radiometric resolution and dynamic ranges, and accurate calibration in the thermal infrared (TIR) bands [30].MODIS land surface temperature products provide global temperatures and emissivities with daily or quasi-daily temporal resolutions and either a 1-km spatial resolution, retrieved using the generalized split-window algorithm [31], or a 5-km spatial resolution, retrieved using the MODIS day/night LST algorithm [32].In this study, MODIS Collection-5 LSTs at a 1-km spatial resolution (MOD/MYD11A1) were evaluated, while the synchronous Collection-5 LSTs at a 5-km spatial resolution (MOD/MYD11B1) were downloaded to calculate the land surface broadband emissivity of each station based on Equation (4).
The MOD/MYD11A1 product pixels were determined based on the location of the validation stations.Then, the validated pixels were further selected through the LST product quality control field (QC = 0) and the sensor view zenith angle field (SVZ < 40°).To calculate the broadband emissivities of the stations, emissivity data in band 29, band 31, and band 32 were selected from the MOD/MYD11B1 corresponding to the MOD/MYD11A1 pixels that had been previously selected.

TM Data
The Landsat Thematic Mapper (TM) images have been extensively studied for various purposes [33,34].The remote sensor has a thermal band (TM6) with a wavelength range of 10.45-12.45µm and a nominal ground resolution of 120 m × 120 m.This spatial resolution of Landsat TM6 is high enough for analyzing the detailed spatial patterns of thermal variations across the Earth's surface.Therefore, at each station, a TM image with a clear sky and a pass time close to MODIS was selected to retrieve the land surface temperature at each site for the scale-effect analysis.
The TM land surface temperature was obtained from the following expression of the radiative transfer equation (RTE) applied to the thermal infrared region: where Lsensor,λ is the at-sensor radiance or top of atmospheric (TOA) radiance, ελ is the land surface emissivity, Bλ(TS) is the blackbody radiance given by Planck's law, TS is the LST, , atm L λ ↓ is the downwelling atmospheric radiance, τλ is the total atmospheric transmissivity between the surface and sensor, and L λ ↑ can be calculated using MODTRAN [35].Therefore, from Equation ( 7), TS is calculated by inverting Planck's law.The inversion of Equation ( 7) can be interpreted as a correction of the atmospheric and emissivity effects on the data obtained by the sensor.
A convenient atmospheric correction tool, which was developed on a public access website [36] for the thermal bands of the Landsat-5 and Landsat-7 sensors, was used to calculate the atmospheric parameters.The Atmospheric Parameter Calculator uses the National Centers for Environmental Prediction (NCEP) modeled atmospheric global profiles [37], which were interpolated to a particular date, time, and location, as input.The upwelling and downwelling radiances were retrieved using the MODTRAN radiative transfer code and a suite of integration algorithms for the site-specific atmospheric transmission [38,39].The atmospheric parameter results are shown in Table 2.Because some sites are in close proximity, 10 TM image files can completely cover the 12 sites.To calculate the ελ in Equation ( 7), a modification of the NDVI Thresholds Method-NDVI THM was used, which depends on the NDVI data obtained from TM3 and TM4 bands.This method works well compared with reference methods, such as the method based on the TISI indices [40] noted by Sobrino et al. [41].
The emissivity retrieved from the TM3 and TM4 bands has a resolution of 30 m.To match the emissivities retrieved from the TM3 and TM4 bands, the TM6 thermal infrared band images were resampled at a 30-m spatial resolution.Next, the TM land surface temperature was calculated using the method described above.To analyze the site's heterogeneity, a 2 × 2-km region centered on the site was selected from the TM retrieval land surface temperature images as the analysis region for each site.In this context, the semi-variance was introduced to compare the diversity of the scale effect between the sites.The GS+ software was applied to analyze the variogram for each site.Prior to this analysis, the TM land surface image grids were discretized to scatter points with a 30-m distance.

Ground Measurements
The China Northern Arid and Semi-Arid Region Cooperative Experimental Observation Integrated Research program was implemented to obtain a series of field observations focused on China's arid and semi-arid areas.The program is an allied observation experiment by several research institutions, including the Institute of Atmospheric Physics-Chinese Academy of Sciences (CAS), the Cold and Arid Regions Environmental and Engineering Research Institute-CAS, Lanzhou University, Beijing Normal University, the Gansu Meteorological Bureau, the Institute of Soil and Water Conservation-CAS, the Institute of Botany-CAS, and the Xinjiang Meteorological Bureau.The continuous observational dataset covers hydrology, soil, atmosphere, and ecology in the northern arid and semi-arid regions of China.
The longwave radiation is measured by pyrgeometers.The field-of-view of the pyrgeometers mounted on a 10-meter tower is approximately 74.6 m in diameter.For the different observation objects, the sensors were not all installed at the same height; the heights ranged from 10 to 35 m.The diameters of the sensor footprints ranged from 74.6 m to 130.6 m.High-quality surface longwave radiation measurements acquired with the pyrgeometers at the 12 stations (in Figure 1) were obtained during three months in 2008 and 2009, namely, July, August, and September, which compose the growing season.All the ground measurements at the 12 stations are 30-min averages.For consistency with the MODIS LSTs, the ground measurements were selected based on the viewing time field associated with the MODIS LST products.

Presentation of Validation Results
Based on the scatter plot for each site in Figure 3, the validation results derived from the two different emissivity retrieval methods for obtaining the land surface temperature from the ground-measured longwave radiation data agree well.The two fit lines in each figure are very close.In general, the scatter is mostly near the 1:1 line (the black line in the scatter plot).However, there are a few outliers positioned far from the y = x line that are mostly located below the y = x line; thus, the MODIS LSTs that have large errors correspond to these underestimated data.One explanation for the outliers is that the cloud effect was not completely eliminated; as a result, the pixels in the MODIS land surface product images appear to represent a lower temperature than the ground-measured temperature.Additionally, the representativeness of the ground measurements and the scale effect may be questionable.
Table 3 shows the validation results; there is a high level of agreement between the MODIS 1-km land surface temperature products (MOD/MYD11A1) and the ground-measured land surface temperatures for the AR, NM, JZ, TYG, TYF, and HZZ stations, with a bias of <1 K.However, at the DS, MQ, and YK stations, the biases are within 1-2 K.In general, most of the biased values are positive, which means that the MODIS LST products are underestimated compared with the ground-measured data.The poorest validation result occurred at the SPT station, with a 3.8 K error.With the exception of the SPT station, the mean absolute errors (MAEs) of the stations are within 2-3 K, and the root-mean-square errors (RMSEs) are all less than 4 K.The validation results showed that MODIS LST is obviously underestimated at the SPT station.The surface cover at the SPT station is desert, which is different than that at the other stations.Therefore, a possible explanation for this result is that the emissivities of the MODIS V5 land surface temperature products are universally overestimated in desert areas, as noted by Hulley [42], which leads to an underestimated LST by MODIS V5 in deserts.In contrast, at the stations with other land surface covers, the MODIS LSTs had a high level of accuracy.Measurements-MODIS LSTs|), and RMSE for the land surface temperature ground measurements obtained using the emissivities from the ASTER Spectral Library (T ASTER ) and MODIS emissivity products, which are both compared with the MODIS LSTs (T MODIS ).

Temporal Mismatch
In this study, the ground-measured data are half-hour time aggregates (30-min temporal average), whereas the Terra and Aqua imagery are effectively temporally instantaneous (acquired in 5-min temporal swath/scenes).From the WATER (Watershed Allied Telemetry Experimental Research) experiment [43] conducted in the Heihe River Basin, which is China's second largest inland river basin, high-frequency longwave radiation data are available for the AR, HZZ, and YK stations.The data are 10-min averages that have been used to validate the MODIS 1-km LSTs (MOD/MYD11A1) [44].A comparison of the validation results of the 10-min averages of the ground-measured data with the results in this study shows similarities (except for the results from the YK station): at the AR station, the biases are 0.30 K and −0.49K, respectively; at the HZZ station, the biases are 0.51 K and 0.11, respectively.At the YK station, the bias between the GMS (ground-measured data) and MOD (MODIS land surface temperature) is 0.06 K based on the 10-min averages, whereas in this study, the bias is −1.74 K.The scatter diagram shows that large biases mainly exist in the high-temperature regions.Comparisons between the different stations indicate that a difference occurs in the temporal mismatch effect.The only explanation for the discrepancies between the three stations is that the fragmented land cover at the YK station exacerbates the mismatch effect.

Spatial Mismatch
Figure 4 provides 2 × 2-km land surface temperature images retrieved from the TM images, which are centered on each station.The analysis focuses on the central 1 × 1-km region because the spatial resolution of the MODIS LSTs is 1 km.The 2 × 2-km land surface temperatures at AR can be divided into two main subareas, namely, the low-temperature area in the upper-left corner and bottom-left corner, and the high-temperature zone in the center.The land surface temperature in the central 1 × 1-km region is approximately 24.7-25.2°C.In the 2 × 2-km region at the HZZ station, with the exception of the upper-right low-temperature region, the land surface temperature ranges within 47.6-48.8°C.The YK station is surrounded by fragmented cropland, so the land surface temperature at this station is composed of small sub-regions and is greatly heterogeneous.Similarly, we found that at the JZ, MQ, TYF, and TYG stations, the central 1 × 1-km land surface temperatures are more uniform than those at the MY, NM, SPT, YZ, and DS stations.From the images shown in Figure 4, the YK station's land surface temperature is more heterogeneous than that at the other stations, and the land surface status of the YK station leads to a larger temporal mismatch.Therefore, these findings indicate that strong land surface spatial heterogeneity can exacerbate the temporal mismatch, particularly in high-temperature regions.
Figure 5 shows the variogram and the optimal model curve at each station.The optimal models of the 12 stations are different.For the HZZ, JZ, MY, TYG, TYF, and YZ stations, the optimal models are exponential models.For the YK, DS, and SPT stations, the optimal models are Gaussian models.For the NM, MQ, and AR stations, the optimal models are spherical models.The optimal models are determined by the r 2 and residual sum of squares (RSS) of the models.Table 4 provides the important parameters of the variogram.The C0 (nugget) and C0/C + C0 (structure variance ratio) provided in Table 4 demonstrate that a small amount of spatial heterogeneity is caused by a random factor, mostly less than 5 percent.For the HZZ, TYG, TYF, and MQ stations, C + C0 < 1 indicates that the variations are all small.A0 indicates the autocorrelation scale of the regionalized variable.The A0 values of the AR, DS, HZZ, MQ, and NM stations are all close to 1 km, which means the autocorrelation scales of these stations are approximately 1 km and the land surface temperature units at these stations are close to the resolution of the MODIS LSTs.Therefore, the station validation results that were obtained using dozens of meters of high-resolution ground-based data to validate the 1-km resolution MODIS LSTs are reasonable.The validation biases of these stations are primarily less than 1 K.The validation at the SPT station is inferior to that at the other stations (see Section 5.1 and Table 3), with a bias of 3.8 K (GMS-MOD), an A0 of approximately 244 m, and a C + C0 value of 3.041 (Table 4).The A0 and C + C0 values at the SPT station are not the smallest values among all the validation stations, and the effect of the heterogeneity of the land surface temperature at the SPT station is not expected to be the strongest.The analyses suggest that the MODIS LST at the SPT station is less accurate than that at the other stations.The SPT station is the only station located in the desert region.Thus, the land cover type is responsible for the underestimation of the MODIS LST products at the SPT station; this finding supports the explanation in section 5.1, which agrees with Hulley's study [42].

Sensitivity Analysis of Emissivity
To quantify the effect of the broadband emissivity estimation error, we established a series of biases Δεb (±0.001, ±0.002, ±0.003, ±0.004, ±0.005, and 0.01) based on the εb value obtained from the ASTER Spectral Library (shown in Figure 6).The results provided in Table 2 show that for the two methods, the obtained emissivities closely correspond.The emissivities obtained from the ASTER Spectral Library are fixed, whereas the emissivities from the narrowband MODIS emissivity products vary with time.When the data used to validate the results only include three months in one or two years (2008 and 2009), the changes in the emissivities are too small to affect the validation results.Thus, the error of the ground-based land surface temperature (ΔTS) varies with time, as shown in Equation ( 8): where εb is the broadband emissivity of the station obtained from the ASTER Spectral Library, Δεb is the setting bias, and TS(εb) and TS(εb + Δεb) are the corresponding land surface temperatures when the land emissivities of the station are εb and (εb + Δεb), according to Equation (2). Figure 6 shows the curve of the results for ΔTS at the 12 stations, and Table 5 provides a summary of the results of the emissivity.Based on Figure 6 and Table 5, the following are concluded: (1) ΔTS varies directly with Δεb in parallel form.Based on Equations ( 8) and ( 2), Equation ( 9) can be formulated as: where F↑ − F↓ is the value of the ground-measured longwave upwelling radiation minus the value of the ground-measured longwave downwelling radiation.In the equation, the values Ts, σ, and F↑ − F↓ are fixed when considering different εb values for the same ground-measured longwave radiation. 2b ε is very close to 1 and can be considered 1 when its value is much smaller than F↑ − F↓ and 3 S T .
Therefore, ΔTS primarily yields a value proportional to Δεb, and the curved lines for different Δεb values vary parallel to each other in every image.There are two primary sections for every station's curve because the ground-measured data include both the day-measured longwave data and the night-measured data; the beginning section of the figure includes the daytime results, and the second section includes the nighttime results.This result also indicates that the same emissivity estimation error yields a greater land surface temperature error during the day than during the night, as calculated using Equation (2) from the ground-measured longwave radiation data.
(2) The maximum MAE is <0.3 when the estimation emissivity |error| is ≤0.01 (see Table 5), whereas the MAE of the validation results ranges within 2-3 K.The results show that the broadband emissivity calculated from the MODIS narrowband emissivity using Equation (4) agrees well with the ground-based measurements of emissivity at Gaize (32.30°N, 84.06°E, with an elevation of 4420 m) in the western Tibetan Plateau, with a standard deviation of 0.085 and a bias of 0.0015, according to Wang et al. [11].Therefore, the validation error of all the stations that results from that method of obtaining the broadband emissivity is less than 0.037 K, which is the maximum average absolute error value of all the stations, when the emissivity bias = 0.002.Compared with the validation results with a MAE of 3 K, the effect of the emissivity estimation with a MAE of 0.037 K can be considered insignificant.
(3) In summary, the validation results from the two emissivity methods are nearly the same.However, the results between the various stations differ.An extreme value appears at the HZZ station, with a discrepancy of 0.2 K between the two methods.According to Equation (2), the only parameter that can cause this difference is the broadband emissivity.The broadband emissivity value that is calculated using the ASTER Spectral Library is a constant value based on the cover types of the stations.In contrast, the broadband emissivity value calculated from the narrowband emissivities in the thermal regions is a value that varies with time and space.Therefore, because the emissivity at the HZZ station changes more over time compared to the other stations for the three months studied, a large difference occurs.The two methods for obtaining broadband emissivity are both suitable for short-term validation, and significant emissivity changes do not occur regionally.In contrast, for land areas with emissivities that change quickly over time and space or over a long period, estimating the stations' broadband emissivities from the MODIS narrowband emissivity Collection 5 in the thermal infrared region is the optimal choice.

Summary
The MODIS LST emissivity products provide global temperature and narrowband emissivity data in thermal infrared bands at daily, eight-day, and monthly time intervals.This study presents the validation results of the MODIS 1-km-resolution daily LST products (MOD/MYD11A1) using ground-measured longwave radiation data from 12 stations located near the semi-arid and arid regions of North China, based on thermal radiative transfer theory.The important factors that impact the validation results are discussed.
The results show that there are general underestimations of the MODIS LST product MOD/MYD11A1, which corroborates Wan's finding [6], in which the values of the MODIS emissivity products are usually overestimated in semi-arid and arid regions; as a result, the MODIS LST products are underestimated.The validation results indicate that the accuracy of the fine-quality MODIS LST products (the QC field is 0) is less than 1 K at most of the stations (7 out of 12), with the MAE within 2-3 K at all stations, except at the SPT station.The results at the SPT station exhibit lower accuracy MODIS LST products in the desert region.The discussion of the temporal and spatial mismatches indicates that spatial heterogeneity can affect the results obtained using this validation method.The use of ground-measured data from stations with the most spatial heterogeneity can yield more reasonable validation results.Moreover, strong spatial heterogeneity can aggravate the effect of the temporal mismatch on the validation results.At the stations with strong spatial heterogeneity, the radiance-based method (R-based) is an advanced alternative method for validation, which does not rely on ground-measured LST value but does require both LSE spectra and measured atmospheric profiles over the validation site at the time of the satellite overpass [45].The comparisons that used the two methods to obtain the broadband emissivity for the stations yielded nearly the same validation results.Based on the different results obtained from the two methods at the different stations, an optimal method for obtaining the broadband emissivity is proposed, although the assumption is limited because the emissivity of the material does not vary with its LST when this method is used for estimating broadband emissivity during daytime and nighttime.Based on a comparison of the error attributable to the broadband emissivity estimation in the daytime with that in the nighttime, the effect of the validation results in the daytime is significantly larger than that in the nighttime.However, the overall validation error caused by the emissivity estimation is minor, with a MAE of less than 0.037 K, and it can even be ignored when compared with the MAE of the validation results (2-3 K).

Figure 1 .
Figure 1.Study area and locations of the validation stations.

Figure 2 .
Figure 2. The flowchart of the study."MODIS NB_Emissivity" is for the MODIS narrowband emissivity products retrieved by the day/night algorithm."GB_LW Data"is the ground-based longwave radiation data."ASTER SL Data" is the ASTER spectral library data."NB_AW Data" is the ground-based LST obtained from "GB_LW Data" and the broadband emissivity retrieved based on "MODIS NB_Emissivity"."ASTER_AW Data" is the ground-based LST retrieved from "GB_LW Data" and broadband emissivity from the "ASTER SL Data"."BE_MODIS Results" is a comparison of results between the "NB_AW Data" and "MODIS LST MOD/MYD11A1"."ASTER_MODIS Results" is a comparison of results between the "ASTER_AW Data" and "MODIS LST MOD/MYD11A1".

Figure 3 .
Figure 3. Plots of the ground-measured LSTs vs. MODIS LSTs at the 12 sites.The black diamonds (ASTER_MODIS(C)) indicate the evaluation results for the ground-station LSTs based on the broadband emissivity data from the ASTER Spectral Library compared with the MODIS LSTs.The pink circles (BE_MODIS(C)) indicate the evaluation results for the ground station LSTs based on the broadband emissivity data from the MODIS narrowband emissivity products compared with the MODIS LSTs.MS_LST (x-axis) is the LST obtained from a station's measured longwave radiation.The linear fit of ASTER_MODIS corresponds to the ASTER_MODIS(C) plots, and the linear fit of BE_MODIS(C) corresponds to BE_MODIS(C).(a) The plot of AR station; (b) the plot of DS station; (c) the plot of HZZ station; (d) the plot of JZ station; (e) the plot of MQ station; (f) the plot of MY station; (g) the plot of NM station; (h) the plot of SPT station; (i) the plot of TYF station; (j) the plot of YK station; (k) the plot of TYG station; (l) the plot of YZ station.

Figure 6 .
Figure 6.Results of the sensitivity analyses of the broadband emissivities for each station.The y-axis indicates the ΔTS in Equation (8).The x-axis indicates the sequence numbers of longwave radiation data used to evaluate the MODIS LSTs during the daytime or nighttime.(a) AR station; (b) DS station; (c) HZZ station; (d) JZ station; (e) MQ station; (f) MY station; (g) NM station; (h) SPT station; (i) TYF station; (j) YK station; (k) TYG station; (l) YZ station.

Table 1 .
Information regarding the stations in this study.

Table 2 .
Atmospheric parameter results.The first column in the table lists the TM image names used to analyze the spatial heterogeneity, and the last eight digits of the image file names indicate the acquisition data.The second and third columns show the center location of each image.The last three columns are the atmospheric parameters from the atmospheric correction tool. a

Table 3 .
Summary of validation sites.

Table 4 .
Summary of the variogram model.

Table 5 .
Summary of the emissivity sensitivity analyses for all stations a .In this table, the second column shows the average of ΔT S at each station, based on Equation (2), by adding a ±0.01 difference to the ASTER estimation broadband emissivity value.Similarly, the last five columns are the average of ΔT S obtained by adding ±0.005, ±0.004, ±0.003, ±0.002, and ±0.001 to the broadband emissivities, respectively.