Introduction

The land surface broadband emissivity (BBE) is a key variable for estimating land surface longwave net radiation, which is a surface radiation budget component and an important parameter in climate, weather, and hydrological models1,2,3,4,5,6. It is the most effective and accurate method to calculate the emissivity based on the obtained spectral data of the broadband surface emissivity.

Broadband emissivity spectra can be accurately measured by using Fourier transform infrared spectrometer (FTIR)7,8,9,10,11. Liu Yongqiang et al.12 and Liu et al.13 used FTIR to study the surface emissivity in the hinterland of Taklimakan Desert for the first time. Additional, thermal infrared remote sensing method is an effective method for getting the emissivity of the area. But thermal infrared remote sensing method can only obtain specific narrow-band emissivity and cannot represent the broadband emissivity14. However, due to the lack of reliable observations, a constant emissivity value or very simple parameterized models of the emissivity are currently used in Land Surface Models (LSM) and General Circulation Models (GCM)3,15. For example, the land model version 2 by the National Center for Atmospheric Research, USA calculates the canopy emissivity from the leaf area index and sets the soil and snow emissivities as 0.96 and 0.97, respectively16.

A sensitivity study of the simulated energy balance to the changes in emissivity in North Africa and the Arabian Peninsula indicated that a decrease of the soil BBE by 0.1 will cause the increases of the ground temperature and air temperature by an average of approximately 1.1 °C and 0.8 °C, respectively, and the decreases in the net and upward longwave radiation by approximately 6.6 W·m−2 and 8.1 W·m−2, respectively15. Jin and Liang3 also revealed the contributions of the BBE in improving the climate models simulation results.

Therefore, when estimating the broadband surface emissivity, most of the researchers would first select the emissivity remotely sensed data of the thermal infrared band to establish the broadband emissivity estimation model. Ogawa et al.17 determined the optimal broadband window for estimating longwave net radiation in the range from 8 μm to 13.5 μm, using the narrowband emissivity from MODIS and ASTER (Advanced Space borne Thermal Emission and Reflection Radiometer) to build the broadband emissivity model. They estimated the broadband land surface emissivity in the Sahara Desert and a cross-check was conducted.

Wang et al.18 established an emissivity estimation model for the band range from 8 μm to 13.5 μm using five thermal infrared bands linear combination from ASTER to estimate the longwave net radiation and analyzed the correlation between the emissivity and soil moisture, but did not verify the estimated emissivity accuracy. Wang et al.19 used the MODIS thermal infrared band data to estimate the surface emissivity of the window from 3 μm to 14 μm and verified the results using the MODIS emissivity library.

Tang et al.20 selected two broadband windows from 3 μm to14 μm and from 3 μm to ∞ μm respectively, then established broadband radiance estimation by using MODIS UCSB (University of California, Santa Barbara) emissivity library and ASTER spectral library (JHU & JPL spectral library) data. The model still used the spectral library data to validate the model and its estimation accuracy about soil, vegetation, water bodies, ice, and snow. Ogawa et al.21 also used MODIS data from the 7th band, but without using the ground observation data the applicability of these estimation equations in local areas is limited.

In view of the shortcomings of the above methods in estimating the surface BBE, this study used the FTIR measured data to optimize the coefficients in the GLASS BBE equation to calculate the BBE data in the Taklimakan Desert. Through the comparison and analysis between FTIR data with MODIS and GLASS BBE products data. In addition, the consistency of the GLASS BBE products on the spatial scale was evaluated. The GLASS and MODIS BBE distribution in the Taklimakan Desert was analyzed and FTIR data was used directly to verify and evaluate the GLASS BBE accuracy.

Results

Estimation equation

Li et al.22 pointed out that MODIS emissivity product data accuracy is preferably when FTIR measurement data is used and the coefficients of the formula (9) are fitted by using 10 measurement values as follows (Eq. 1):

$$y=0.0675{\varepsilon }_{29}+0.1326{\varepsilon }_{31}+0.7842{\varepsilon }_{32}-0.1206{\rho }_{7}+0.0071$$
(1)

The rest 15 values were used to verify the results (Fig. 1). It can be shown that The fitting result of the MODIS BBE estimation equation obtained by using FTIR data is better than the original, and R2 reaches 0.95. Therefore, the value of MODIS fitted by Eq. (1) can be used to verify GLASS BBE.

Figure 1
figure 1

Emissivity comparison between FTIR observed values and MODIS estimated values. The graph was generated with OriginPro Portable 9.0.0. (https://www.originlab.com/).

Full size image

Comparative verification

First, extract the image value from the same position as the FTIR measured point from GLASS, were compared with MODIS emissivity products and the FTIR values, respectively. The results (Fig. 2) show that the GLASS products correlation with either the FTIR or MODIS is relatively low before the GLASS products are revised. The R2 with the FTIR was 0.421, and RMSE and Bias were 0.02899 and −0.02769, respectively. The correlation with MODIS was better than the correlation with FTIR, but its RMSE and Bias were both generally high. Therefore, GLASS product gives high emissivity values and its accuracy in describing the desert needs to be improved.

Figure 2
figure 2

Comparison between FTIR observation (left) and MODIS estimation (right) with GLASS BBE. The graph was generated with OriginPro Portable 9.0.0. (https://www.originlab.com/).

Full size image

We fitted the GLASS data by randomly selecting 10 figures from the measurement data and using the multiple regression method to obtain the equation coefficient (10) as follows (Eq. 2):

$${\varepsilon }_{8-13.5}=0.235{\alpha }_{1}-0.724{\alpha }_{2}-0.325{\alpha }_{3}+0.231{\alpha }_{4}+0.313{\alpha }_{5}+0.757{\alpha }_{6}-0.7126{\alpha }_{7}+0.964$$
(2)

The rest 15 figures were taken into the fitted Eq. (2) and the result was compared with the FTIR measurement data and MODIS data respectively as shown in Fig. 3.

Figure 3
figure 3

Comparison between FTIR observation (left) and MODIS estimation (right) with modified GLASS BBE estimations. The graph was generated with OriginPro Portable 9.0.0. (https://www.originlab.com/).

Full size image

It is shown in Fig. 3 that for the revised GLASS data, its R2 with FTIR data has changed from 0.421 to 0.948, and its R2 with MODIS has changed from 0.691 to 0.907. The corresponding RMSE was reduced by a magnitude order and the Bias was also significantly reduced.

Since MOD43A3 is the albedo data, the value in the desert area and the vegetation area is opposite to the emissivity. Therefore, the equation value in the oasis zone around the Taklimakan Desert is not accurate enough. Therefore, in combination with the GLASS broadband emissivity algorithm and the value of the MODIS NDVI product, the algorithm for the entire region is re-fitted. So combined with re-fitting broadband emissivity of GLASS and the product of MODIS NDVI. Out of the equation for the entire Taklimakan Desert as follows (Eq. 3):

$$\begin{array}{c}{\varepsilon }_{8-13.5}=0.036\cdot NDVI+0.235{\alpha }_{1}-0.724{\alpha }_{2}-0.325{\alpha }_{3}+0.231{\alpha }_{4}+0.313{\alpha }_{5}\\ \,+0.757{\alpha }_{6}-0.7126{\alpha }_{7}+0.964\end{array}$$
(3)

Emissivity distribution characteristics in the Taklimakan Desert

The emissivity distribution in the Taklimakan Desert was conducted by selecting the bands from 1 to 7 of the black-air albedo from MODIS MOD43A2 product on September 27, 2014. The surface broadband emissivity was calculated using Eq. (3). The area is 1300 km × 500 km, including the entire Taklimakan Desert and the oases around the Tarim Basin and the southern Tianshan Mountains and the northern Kunlun Mountains (Fig. 4). It can be seen that the broadband emissivity value from the Taklimakan Desert region was centered on the range between 0.89 and 0.92, covering the majority of the area distribution. The values in the regions with sparse vegetation were between 0.92 and 0.95, and the values for the oases near the edge of the desert, such as Alar, Qiemo, Luntai, Kashgar oasis, were between 0.95 and 0.98.

Figure 4
figure 4

Distribution for the broadband emissivity over the Taklimakan desert (Lines are the boundaries between deserts and oases). The map was generated with ArcMap 10.0. (http://www.esri.com/en/arcgis/arcgis-for-desktop/).

Full size image

Conclusion and Discussion

We used the FTIR measured data to establish a multivariate linear regression equation for GLASS BBE products and to establish the equation with MODIS bands 7, 29, 31, and 32. The accuracy of the modified GLASS BBE equation is significantly improved in the Taklimakan Desert and some results are as follows:

First, the GLASS emissivity products correlation analysis was carried out by using the data measured by FTIR, the R2 of the measured value is 0.421, RMSE and Bias are 0.029 and −0.02769, respectively. For MODIS BBE, the R2 is 0.691, the RMSE and Bias are 0.0268 and −0.026, respectively, and the GLASS data were higher than FTIR and MODIS data.

Second, using the FTIR data to perform multivariate linear regression on the calculation formula of GLASS emissivity products, the correlation between the calculated new values and the values not participating in the regression was obtained, and the obtained R2 is significantly improved, which was increased by 0.528, RMSE and Bias were reduced by 0.01713 and 0.02766, respectively. For MODIS data correlation analysis, the R2 was also increased by 0.217, RMSE and Bias were decreased by 0.024 and 0.02505, respectively, and the results were significantly improved.

Finally, using the calculated regression equation to perform the band calculation on the MODIS43A3 image, it can be seen that the surface emissivity value reached the highest of 0.95 near the oasis, and was in the range between 0.91 and 0.95 in the oasis and desert transition zone, in the range between 0.88 and 0.91 in the rest of the desert. The value in the desert hinterland area was less than 0.88.

Although the estimation method proposed in this paper achieved a reasonable broad-band surface emissivity, it was only validated and analyzed in the Taklimakan Desert area, and it requires more extensive verification if it is applied in other desert areas, or in other vegetation coverage areas. Its applicability in wet areas also remains to be further studied.

Study area and data

Study area

The Taklimakan Desert is located in the middle of the Tarim Basin in Xinjiang of China. It is adjacent to the Tomur Mountains and the Pamirs in the west, the Kunlun Mountains and the Altun Mountains in the south, the Tianshan Mountains in the north and the Lop Nur in the east. It is about 1070 km long from east to west, 420 km wide from north to south, with an area of 33.76 × 104 km2 and an altitude in the range from 800 to 1300 m. The terrain is high in the west and low in the east. It is the second-largest mobile desert in the world23,24,25,26. Desert soil type is single, mainly being silt, of which the most is fine sand, accounting for 43.8% to 75.5% of sediment transport, and a small amount of feldspar and muscovite12. According to the historical observation data of the Tazhong Meteorological Station in Xinjiang (39°00′N, 83°40′E, altitude 1099.3 m) from 1996 to 2014, the highest temperature in the region was up to 45.6 °C, the lowest temperature was −32.7 °C. The annual average wind speed was 2.5 m·s−1, and the instantaneous maximum wind speed was 24.0 m·s−1 (Fig. 5)27.

Figure 5
figure 5

FTIR observed sites along two desert highways. The map is generated with ArcMap Version 10.0(Image download from Google Earth) (http://www.esri.com/en/arcgis/arcgis-for-desktop/).

Full size image

Data

FTIR observation data

The field observation of the broadband emissivity in the Taklimakan Desert was conducted by using a portable Fourier transform thermal infrared spectroscopy (FTIR). The observed sites were picked up along two desert highways which cross through the desert from north to south (one is from Luntai to Niya, another from Alar to Hetian (Table 1). The soil composition and humidity change are small due to the single type of desert surface; thus, the suitable site was chosen every 50–100 km and the observed sites were increased in the oasis transition zone at the edge of the desert. In order to have effective observation data, the measurement was carried out on sunny days22. The field data collection was from Oct 16 to Oct 19, 2013, for the sites in the Luntai-Niya Desert Highway and from Sep 25 to Sep 27, 2014, for the sites in the Alar-Hetian Desert Highway. The data for a total of 25 sites were collected.

Table 1 Site place latitude and longitude coordinates.
Full size table

MODIS data

MODIS data provides 36 discrete-band images between 0.4 to 14.4 μm, with high spectral resolution and time resolution. Global observations are available every 1 to 2 days for large scales. Significant advantages in large-scale monitoring of global environmental parameters such as temperature28. The MODIS data has 36 channels, of which the channels 1–19 and 26 are distributed in the visible and near-infrared bands, and the other 16 channels are distributed in the thermal infrared band from 3–14 μm. The thermal infrared band atmospheric window that obtains the most suitable surface emissivity is in the range from 8 to14 μm, which correspond to the bands from 29 to 32. Because the 30th band has strong ozone absorption, it cannot be utilized. Therefore, the thermal infrared wave points were chosen band 29 (8.400–8.700 μm), band 31 (10.780–11.280 μm) and band 32 (11.770–12.270 μm). In addition, the 7th band (2.105–2.155 μm) in the near-infrared spectra is the reflection of surface and clouds29 which can represent the surface soil properties30, and thus the band 7 has the highest correlation with the broadband emissivity in all MODIS reflection channels27. The soils and rocks surface reflectivity depends on their mineral composition and emissivity. In general, the quartz-rich (SiO2) sand has higher reflectivity and lower emissivity, but the mafic minerals with lower SiO2 content generally have lower reflectance21. The Taklimakan Desert is rich in quartz sand and has a high surface reflectivity23,30, therefore the 7th band was added to the study.

We used the bands 31and 32 from MOD11A1 which is a MODIS product about the temperature versus emissivity, the band 29 from MOD11B1 (with a resolution of 1000 m) and the band 7 reflectance data (with a resolution of 500 m) from MOD09GA which is the land secondary standard product. The acquired image data covered the entire Taklimakan Desert in the time period from Oct 16 to Oct 18, 2013, and from Sep 25 to Sep 27, 2014. The MODIS data were downloaded from the National Aeronautics and Space Administration (NASA) website at http://modis.gsfc.nasa.gov.

GLASS data

The GLASS emissivity is a BBE (covering wavelength range from 8 μm to 13.5 μm) product derived from AVHRR VNIR data and MODIS albedos by using newly developed algorithms31,32. GLASS emissivity has two parts: (i) the global eight-day 1 km land surface BBE retrieved from MODIS albedos from 2000 through 2010 and (ii) the global eight-day 5 km land surface BBE retrieved from AVHRR VNIR reflectance from 1981 through 1999.

In the algorithm used to generate GLASS BBE from MODIS albedos, the land surface was classified into five types according to Normalized Difference Vegetation Index (NDVI) threshold values, such as water, snow or ice, bare soil, vegetation, and transition zone. The overlapping areas of bare soil and transition zone or of transition zone and vegetation were noted. The BBEs of water and snow or ice were set as 0.985 based on a combination of BBE calculated from the emissivity spectrum in the ASTER spectral library and the MODIS USCB emissivity library, and BBE simulated through radiative transfer models33.

BBEs of bare soil, vegetated areas, and transition zones were formulated as the linear function of seven MODIS narrowband black-sky albedos. When NDVI was lower than 0.1 or higher than 0.2, the formula for bare soil or vegetated areas was used to calculate individual BBE values. In overlapping bare soil and transition zone areas (0.1 < NDVI ≤ 0.156), BBE was taken as the average of values calculated by using formulas for bare soil and transition zones. By contrast, BBE for areas of the overlapping transition zone and vegetated area (0.156 <NDVI< 0.2) was taken as the average values calculated by using formulas for transition zones and vegetated areas. BBE derived from MODIS albedos was validated by field measurements in the US and China, and the absolute difference was 0.0231,34,35.

Methods

Broadband emissivity measured by FTIR

The observation of surface-specific emissivity spectral data was collected by using a portable FTIR spectrometer (Model 102 F, Designs and Prototypes Ltd.) with the following parameters: a luminous flux of 0.016 cm−2·sr, an operating spectral range from 2 μm to16 μm, and a spectral resolution of 2–24 cm−1. The measurement results can reach a standard deviation of less than 1%10. The formula for calculating the surface-specific emissivity spectrum is as follows (Eq. 4):

$${\varepsilon }_{\lambda }=\frac{{L}_{\lambda }-L{D}_{\lambda }}{{B}_{\lambda }(T)-L{D}_{\lambda }}$$
(4)

where εi is the spectrum of the surface specific emissivity when the wavelength is λ, Lλ (cm2·sr) is the surface radiance of the wavelength λ, Bλ (T)(cm2·sr) is the black body radiance when the wavelength is λ and the surface temperature is T(K), downward radiation Dλ (cm2·sr) is measured using a calibrated diffuse gold plate3,28. Liu et al. and Liu Yongqiang et al. gave the operation process and observation results of the surface specific emissivity spectrum observation in the Taklimakan Desert12,13.

In the land-surface patterns and numerical prediction models, the surface-specific emissivity uses the broadband emissivity spectra average. The calculation method of the surface emissivity between the wavelength range from λ1 to λ2 and is shown as follows (Eq. 5)36.

$${\varepsilon }_{\lambda 1-\lambda 2}=\frac{{\int }_{{\lambda }_{2}}^{{\lambda }_{1}}{\varepsilon }_{\lambda }{B}_{\lambda }(T)d\lambda }{{\int }_{{\lambda }_{2}}^{{\lambda }_{1}}{B}_{\lambda }(T)d\lambda }$$
(5)

where λ1 and λ2 is the upper and lower wavelength range limits. Rather than what Ogawa et al.17 pointed out that the wavelength ranges from 8 μm to13.5 μm is the best window for estimating long-wave radiation when using FTIR to observe the surface emissivity spectrum, higher precision can be obtained in the range from 8 μm to 14 μm13. Therefore, the broadband emissivity is calculated using a wavelength range from 8 μm to 14 μm. Since ελ and Bλ(T) are continuous functions, the integral equation is discretized as follows (Eq. 6):

$${\varepsilon }_{\lambda 1-\lambda 2}=\frac{{\sum }_{\lambda ={\lambda }_{1}}^{{\lambda }_{2}}{\varepsilon }_{\lambda }{B}_{\lambda }(T)\varDelta \lambda }{{\sum }_{\lambda ={\lambda }_{1}}^{{\lambda }_{2}}{B}_{\lambda }(T)\varDelta \lambda }$$
(6)

Broadband emissivity measured by MODIS

The MODIS thermal infrared single-band surface emissivity is calculated as follows (Eq. 7)37.

$${\varepsilon }_{i}=\frac{{\int }_{{\lambda }_{i1}}^{{\lambda }_{i2}}{f}_{i}(\lambda ){\varepsilon }_{\lambda }{B}_{\lambda }(T)\varDelta d\lambda }{{\int }_{{\lambda }_{i1}}^{{\lambda }_{i2}}{f}_{i}(\lambda ){B}_{\lambda }(T)\varDelta d\lambda }$$
(7)

where εi is the MODIS thermal infrared single band emissivity, fi(λ) is the spectral response function of the band i. According to the Eqs. (2) and (4), the broadband emissivity can be calculated by a linear combination of the emissivity in different bands (Eq. 8)21.

$${\varepsilon }_{{\lambda }_{1}-{\lambda }_{2}}=\frac{{\sum }_{i=1}^{n}{\int }_{{\lambda }_{1}}^{{\lambda }_{(i+1)}}{\varepsilon }_{\lambda }{B}_{\lambda }(T)d\lambda }{{\int }_{{\lambda }_{1}}^{{\lambda }_{2}}B(T)d\lambda }\approx {\sum }_{i=1}^{n}{g}_{i}{\varepsilon }_{i}$$
(8)

where gi is the combination coefficient, which is related to the black body thermal radiance when the temperature is T, and is independent of the single band-specific radiance and n is the number of bands. Thus, using the bands 29, 31 and 32 from MODIS data, the broadband emissivity estimation equation can be obtained by fitting the multiple linear regression equation. Since the band 7 of MODIS is highly correlated with the emissivity, it is added to the linear combination equation. Therefore, to calculate the broadband emissivity ε8–13.5, we have the following formula22 (Eq. 9).

$${\varepsilon }_{8-13.5}={\sum }_{i}{a}_{i}{\varepsilon }_{i}+b{\rho }_{7}+c$$
(9)

where, i = 29, 31, 32, ai, b, and c are the coefficients of the multiple regression equation, εi is the emissivity value for the bands 29, 31 and 32, ρ7 is the reflectance value of the band 7.

Broadband emissivity estimated by GLASS

The global land surface is divided into five types as water body, ice/snow, bare soil, vegetation cover, and transition area. The emissivity estimation method is given respectively based on the type. (i) For water bodies, ice/snow pixels, the value is 0.985; (ii) For other surface types, using ASTER emissivity products, MODIS reflectance and albedo products as data sources, the empirical relationships were established between the narrowband black sky albedo and the BBEs of the bare soil, transition areas, and vegetation-covered areas. In addition, the broadband emissivity spectral range is from 8 μm to 13.5 μm and the algorithm for calculating the bare soil broadband emissivity is shown as follows (Eq. 10):

$$\bar{\varepsilon }=e+{\sum }_{i=1}^{7}{d}_{i}BS{A}_{i}$$
(10)

where, \(\bar{\varepsilon }\) is broadband emissivity, BSAi is a narrowband black sky albedo, e is a constant, di are bands.

For bare soil: NDVI <= 0.1, pixel emissivity from bare soil algorithm; vegetation cover: NDVI > = 0.2, two cases fully vegetated (NDVI > 0.461) and Mixed Scenario (0.2 < = NDVI < 0.461). The emissivity is derived from the vegetation cover algorithm; the transition area is subdivided into two parts: the bare soil transition zone and the vegetation transition zone. The soil transition zone: NDVI is 0.1–0.156, and the pixel emissivity is derived from the transition zone and the bare soil zone. Mean of the algorithm; vegetation transition zone: NDVI is 0.156–0.2, and the pixel emissivity is derived from the mean of the vegetation zone and transition zone algorithm, and the broadband surface emissivity of the vegetation area is expressed as (Eq. 11):

$$\bar{\varepsilon }={A}_{0}\cdot NDVI+{\sum }_{i=1}^{7}{f}_{i}\cdot BS{A}_{i}+h$$
(11)

where, A0 and h are constants, NDVI is the normalized vegetation index, BSAi which are narrow-band black sky albedo, and fi are the coefficient of albedo.