Next Article in Journal
Morphodynamic Controls for Growth and Evolution of a Rubble Coral Island
Next Article in Special Issue
A Review of Reconstructing Remotely Sensed Land Surface Temperature under Cloudy Conditions
Previous Article in Journal
Long-Term and Emergency Monitoring of Zhongbao Landslide Using Space-Borne and Ground-Based InSAR
Previous Article in Special Issue
A Two-Step Integrated MLP-GTWR Method to Estimate 1 km Land Surface Temperature with Complete Spatial Coverage in Humid, Cloudy Regions
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Downscaling Land Surface Temperature Based on Non-Linear Geographically Weighted Regressive Model over Urban Areas

1
College of Global Change and Earth System Science, Beijing Normal University, Beijing 100875, China
2
State Key Laboratory of Remote Sensing Science, Beijing Normal University, Beijing 100875, China
3
Institute of Computer Science and Technology, Chongqing University of Posts and Telecommunications, Chongqing 400065, China
*
Author to whom correspondence should be addressed.
Submission received: 5 March 2021 / Revised: 15 April 2021 / Accepted: 16 April 2021 / Published: 19 April 2021
(This article belongs to the Special Issue Land Surface Temperature Estimation Using Remote Sensing)

Abstract

:
Land surface temperature (LST) is a vital physical parameter in geoscience research and plays a prominent role in surface and atmosphere interaction. Due to technical restrictions, the spatiotemporal resolution of satellite remote sensing LST data is relatively low, which limits the potential applications of these data. An LST downscaling algorithm can effectively alleviate this problem and endow the LST data with more spatial details. Considering the spatial nonstationarity, downscaling algorithms have been gradually developed from least square models to geographical models. The current geographical LST downscaling models only consider the linear relationship between LST and auxiliary parameters, whereas non-linear relationships are neglected. Our study addressed this issue by proposing an LST downscaling algorithm based on a non-linear geographically weighted regressive (NL-GWR) model and selected the optimal combination of parameters to downscale the spatial resolution of a moderate resolution imaging spectroradiometer (MODIS) LST from 1000 m to 100 m. We selected Jinan city in north China and Wuhan city in south China from different seasons as study areas and used Landsat 8 images as reference data to verify the downscaling LST. The results indicated that the NL-GWR model performed well in all the study areas with lower root mean square error (RMSE) and mean absolute error (MAE), rather than the linear model.

1. Introduction

Land surface temperature (LST) is one of the most important parameters for climate research and has been widely used in many research fields [1,2], including urban heat island monitoring [3], hydrological cycle [4], and climate change assessment [5], etc. In recent years, the spatial heterogeneity of land cover in urban areas has been drawing more and more research interest, which requires satellite remote sensing LST data with a higher spatiotemporal resolution [6]; however, there are many difficulties in acquiring high spatiotemporal data. According to the spatial and temporal resolution, the current remote sensing data can be divided into two categories. One is the data with finer temporal resolution and coarser spatial resolution. For instance, the temporal resolution of the FY-4/AGRI satellite is 15 min/full replay, but its spatial resolution of the thermal infrared band is 4 km; the temporal resolution of the Himawari-8/AHI is the same as FY-4/AGRI and its spatial resolution of the thermal infrared band is 2 km [7]; both the advanced very high-resolution radiometer (AVHRR) and the moderate resolution imaging spectroradiometer (MODIS) have a temporal resolution of twice a day, and the spatial resolution of thermal infrared bands are 1.1 and 1 km, respectively [8,9]. The other data category has finer spatial resolution and coarser temporal resolution. The Landsat satellite thematic mapper (TM), enhanced thematic mapper plus (ETM+), and thermal infrared sensor (TIRS) have a spatial resolution in thermal infrared bands of 120, 60 and 100 m, respectively, whereas their temporal resolutions are about 16 days; the imaging frequency of advanced spaceborne thermal emission and reflection radiometer (ASTER) is 16 days and the spatial resolution of its thermal infrared band is 90 m [10,11,12]. In general, a spatial resolution of 100 m and the temporal resolution of once or twice per day is desirable for research on urban heat source distribution; however, the spatiotemporal resolution of satellite remote sensing data cannot achieve this requirement [13]. There are two solutions capable of improving the spatial resolution of thermal infrared images. One solution is to improve the hardware, e.g., improving detector fabrication to increase the spatiotemporal resolution of the satellite remote sensors; however, this would have high production costs and take a significant amount of time. The second solution would be to improve image processing, which uses the visible-near infrared band of higher resolution data and selects an appropriate algorithm to improve the spatial resolution of the thermal infrared band [14]. LST disaggregation algorithms are widely applied to improve spatial resolution by unmixing the thermal infrared pixel. LST disaggregation is divided into two categories: temperature unmixing and thermal sharpening. The temperature unmixing means unmixing the LST of components while maintaining a constant spatial resolution of the pixels. Thermal sharpening refers to improving the spatial resolution by enhancing the spatial details of LST. Among these solutions, LST downscaling belongs to the category of thermal sharpening [15]. During the development of these algorithms, LST downscaling has branched out the following aspects: image fusion model, statistical linear regression model, modulation model, and hybrid model.
The downscaling algorithms based on the image fusion model obtain higher spatial resolution data with continuous time series by fusing the remote sensing data with different temporal and spatial resolutions. Fasbender et al. [16] proposed a Bayesian fusion method, which was proved to have a high accuracy and practicality when fusing visible and near-infrared bands of ASTER together. The spatial and temporal adaptive reflectance fusion model (STARFM) was proposed by Gao et al. [17], and achieved good results in reflectivity image fusion; however, when applied to LST data, the STARFM can only increase the spatial resolution of the effective LST data and is not continuous in the time series. Zhu et al. [18] put forward the flexible spatiotemporal data fusion (FSDAF) model, which required fewer auxiliary parameters to achieve a high prediction accuracy when there is little change in land covers. In the past few years, with the progress and development of artificial intelligence, machine learning methods were applied to image fusion by scholars. Bindhu et al. [19] proposed the backpropagation (BP) neural network to fit the non-linear downscaling non-linear DisTrad (NL-DisTrad) algorithm, which achieved excellent downscaling results when applied to estimating evapotranspiration. Dong et al. applied the neural network to image super-resolution reconstruction, then improved the super-resolution convolution neural network (SRCNN) model and proposed the accelerating the fast super-resolution convolutional neural network (FSRCNN) model, which can input low-resolution images into the network directly [20,21]. Ledig et al. [22] proposed the generative adversarial network (GAN), which has been widely applied to research. Recently, some researchers improved this model, including the dilated convolution generative adversarial network (DCGAN), Wasserstein generative adversarial network (WGAN), and boundary equilibrium generative adversarial network (BEGAN) [23,24,25]. Shao et al. [26] verified the GAN in thermal infrared image reconstruction. The results showed that this method could improve the image’s subjective visual effect while maintaining objective quality evaluation.
The algorithms based on the modulation model are used to distribute LST with a coarse resolution to each sub-pixel proportionally. Guo et al. [27] summarized past studies and proposed the pixel block intensity modulation (PBIM) algorithm, which selects the panchromatic band as the scale factor. The PBIM algorithm can retain the original thermal spectrum information and integrate detailed information into the thermal infrared band. Stathopoulou et al. [28], based on the PBIM method, used emissivity and two LST data types to simulate LST with high resolution to downscale the LST and achieved better results.
The algorithm based on a hybrid model utilizes many types of models. Zhukov et al. [29] proposed an multisensor multiresolution (MMT) method based on spectral decomposition to downscale the TM data and used the ASTER data to verify the downscaling results. Deng et al. [30] applied the spectral unmixing and thermal mixing (SUTM) algorithm to downscale the LST in urban areas.
Due to easier performance and better efficiency, the LST downscaling algorithms based on statistical models have been widely used in recent research. These algorithms can be divided into three types. The first one builds a statistical model between LST and NDVI to achieve the LST downscaling. The second one considers the complexity of the land covers and use variety auxiliary parameters for LST downscaling. The third one considers the spatial non-stationarity between LST and auxiliary parameters and applies the geographical model instead of the global model to downscale LST. These three types of models are under the assumption of scale invariability.
Kustas et al. [31] proposed disaggregation procedures for the radiometric (DisTrad) algorithm, which constructs a relationship between LST and NDVI to downscale the MODIS LST spatial resolution from 1000 to 250 m. Agam et al. [32] improved the DisTrad algorithm and proposed an algorithm for sharpening thermal imagery (TsHARP). This algorithm considered the non-linear relationship between LST and NDVI by adding the NDVI2 into the regression relationship. Yang et al. [33] pointed out that a single factor cannot reflect the difference of LST in different land covers and considered building a regressive relationship between LST and multiple factors. The downscaling results showed that this algorithm was better than DisTrad and indicated that selected multiple parameters would better explain LST. Zhu et al. [34] proposed an improved hierarchical regression method and compared it with DisTrad and TsHARP algorithms, and results showed that this method had obvious advantages and obtained higher accuracy than the other algorithms. Eswar improved the DisTrad algorithm to select NDVI, fractional vegetation cover (FVC), normalized difference water index (NDWI), and soil-adjusted vegetation index (SAVI) as the auxiliary parameters, and the results showed that the FVC and NDWI had higher accuracy for humid areas, whereas the NDVI was more suitable for dry areas [35]. Stefania et al. [36] considered the NDVI, SAVI, normalized difference built-up index (NDBI), and urban index (UI) parameters and selected different parameter combinations for LST downscaling, obtaining the optimal parameter combination for LST downscaling. Qi et al. proposed a new method to combine multiple variable and machine learning algorithms to downscale LST in urban areas [37]. Wang et al. compared and analyzed the downscaled results based on multiple linear regression (MLR), TsHARP, and random forest (RF) methods and indicated that the RF model is applicable to downscaling research in heterogeneous regions [38]. The above downscaling algorithms for LST are examples of global models that build relationships with a global scope and assume stable relationships between LST and auxiliary parameters; however, the non-stationarity relationship between independent variables and dependent variables was not considered. Researchers have paid attention to the non-stationarity between LST and auxiliary parameters and proposed geographical models to downscale LST in recent years. Duan et al. [39] proposed LST downscaling algorithm based on the geographically weighted regressive (GWR) model, which was the first time a geographical model was proposed to downscale LST and achieved better downscaling results than the global TsHARP algorithm. After that, scholars have made improvements to the GWR model. Pereira et al. [40] proposed the geographically weighted regressive kriging (GWRK) model, which downscales the LST of ASTER; Peng et al. [41] considered both the spatial and temporal non-stationarity and proposed the geographically temporally weighted regressive (GTWR) model, which was compared with the GWR and TsHARP models and obtained better downscaling results; Wang et al. [42] considered spatial non-stationarity and autocorrelation simultaneously and proposed the geographically weighted autoregressive (GWAR) model to downscale the spatial resolution of MODIS LST from 1000 to 100 m and achieved good downscaling results.
Based on the above, the statistical downscaling models evolved from a linear model to a geographical non-linear model, added different auxiliary parameters, and considered the non-stationarity relationship between LST and auxiliary parameters. However, the research into GWR and its improved models only considered the geographical linear relationship between LST and auxiliary parameters; the geographical non-linear relationship was overlooked. In this study, we addressed the complexity of the land covers in urban areas, and proposed a non-linear geographically weighted regressive (NL-GWR) model to downscale the MODIS LST.

2. Study Area and Data Preparation

2.1. Study Area

In this study, we selected two Chinese provincial capital cities, Jinan and Wuhan, as the study areas. Figure 1 shows the false-color images generated from Landsat 8 data of the study areas.
Jinan is a typical city in the North China Plain. The selected study area covers a total area of 7998 km2. Jinan has a warm temperate continental monsoon climate with a cold winter and hot summer. The annual average temperature is approximately 13.8 °C, where the lowest and highest temperatures are −19.7 and 42.5 °C, respectively. The average annual precipitation is approximately 685 mm. The study area selected in this paper is the Jinan urban area, and the main types of land cover include building, vegetation, water, and bare soil.
Wuhan covers a total area of 8569 km2. Wuhan represents cities in the mid-south of China where rivers and lakes are widely distributed. It has a subtropical monsoon humid climate, which has abundant rainfall, sufficient heat, and four distinct seasons. The weather is hotter and wetter than Jinan, with an annual average temperature of about 15.8 to 17.5 °C, and the annual precipitation is 1150–1450 mm. The land covers are similar to Jinan.

2.2. Datasets and Preprocessing

The remote sensing data used in this study include Landsat 8 raw reflectance data and MODIS LST data. Three pairs of images were collected for each study area to test the algorithm’s performance in different seasons. Table 1 shows the main characteristics of these images in Jinan and Wuhan. Table 1 indicates that the Landsat and MODIS images were acquired within half an hour time. Due to the short time difference, the solar geometries, orbital parameters, and the viewing (near-nadir) of the MODIS Terra platform are highly consistent those of the corresponding Landsat 8.

2.2.1. MODIS Data

The MODIS LST product selected in this study is the MOD11_L2, collection 6 data with a spatial resolution of 1000 m; the MODIS LST data were downloaded from the NASA website (http://reverb.echo.nasa.gov/ accessed on 31 October 2020). The MOD11_L2, collection 6 product contained segmented data with a spatial resolution of 1000 m and was generated by the split-window algorithm. Extended validation works indicated that the accuracy of the MOD11_L2 product is approximately 1 K [43,44,45] and is widely used in LST downscaling algorithm research.

2.2.2. Landsat 8 Data

The Landsat 8 data of this study were downloaded from the geospatial data cloud (http://www.gscloud.cn/ accessed on 31 October 2020). The downloaded data were systematically processed with radiometric and geometric correction. In this study, the Landsat 8 data were processed with calibration and atmospheric correction modules from the ENVI 5.3 software to derive surface reflectance. In addition, accurate geometric correction between MODIS and Landsat 8 data was necessary. In this study, we used the image-to-image module of ENVI 5.3 software for geometric correction. The Landsat 8 data were selected as the reference images, and feature points such as road and river intersections were chosen as control points to correct MODIS images. Thus, the conjugated pixels from different sensor were matched and stacked to ensure the geometric consistency between images.
In this study, the MODIS LST was vital to model the relationship in the coarse spatial resolution. The Landsat 8 data serve two purposes in this study. One is to provide auxiliary parameters for establishing model, the auxiliary parameters including NDVI, SAVI, NDBI, UI and NDWI. These auxiliary parameters were retrieved from Landsat 8 data and resampled to 100 m and 1000 m for establishing the fine resolution and coarse resolution downscaling model, respectively. The resample method of spatial aggregations is adopted in this study to upscale the auxiliary parameters to coarse resolution of 100 m and 1000 m, so the values of coarse pixels are simply spatial means of the fine pixels. Another purpose is to provide validation reference data. The Landsat 8 LST was retrieved by the Mono-window algorithm and served as the validation data to verify and analyze the LST downscaling results.
The remote sensing inversion methods for LST are maturing and, the accuracy is constantly improving. Scholars have proposed different LST inversion methods for different remote sensing data. The common methods include a radiative transfer equation algorithm, Mono-window algorithm, and split-window algorithm [46,47,48,49]. MODIS LST usually uses the day/night algorithm and the split-window algorithm, whereas the Landsat series satellites use the Mono-window algorithm with the support of in situ measured surface and atmosphere parameters. Compared with the radiative transfer equation, the Mono-window algorithm includes the influence of the surface and atmosphere in the calculation formula, making the application more convenient [50,51].
The United States Geological Survey (USGS) demonstrated that the calibration of the TIRS band 11 of Landsat 8 is temporarily unstable; therefore, the Mono-window algorithm, in combination with band 10 of Landsat 8, was used to retrieve LST data. The formula of the Mono-window algorithm is written as follows [52]:
T S = a 10 ( 1 C 10 D 10 ) + [ b 10 ( 1 C 10 D 10 ) + C 10 + D 10 ] T 10 D 10 T a C 10
where T S , T 10 and T a represent land surface temperature, the brightness temperature of TIRS band 10, and effective average temperature of the atmosphere, respectively. Both a 10 and b 10 are constant. For band 10 of Landsat 8, the a 10 and b 10 can be seen in Table 2 [53].
The coefficient C 10 and D 10 can be calculated by the following equations:
C 10 = ε 10 τ 10
D 10 = ( 1 τ 10 ) [ 1 + ( 1 ε 10 ) τ 10 ]
where ε 10 and τ 10 represent the land surface emissivity and atmospheric transmittance, respectively. It can be indicated that the effective atmospheric average temperature, atmospheric transmittance and surface emissivity are three significant parameters required when using the Mono-window algorithm to retrieve the land surface temperature. In general, the effective mean atmospheric temperature is obtained using the near-surface temperature by using a linear equation. Many factors influenced the atmospheric transmittance, including water vapor, aerosol, wavelength and ozone, although the atmosphere is the most important factor [54]. Two steps were used to estimate the surface emissivity. The first step used the land cover map to distinguish different land cover types in 30 m spatial resolution. The second step was the NDVI threshold method, used to estimate the surface emissivity [55].

3. Method

3.1. Introduction of Non-Linear Geographically Weighted Regressive Model

The relationships between LST and auxiliary parameters are complex and require in-depth exploration and research [56]. The development of the least square method starts from building linear relationship between LST and a single parameter to a linear relationship with multiple auxiliary parameters and further considering non-linear relationship. The traditional least square method ignores the non-stationarity in the relationship between the LST and auxiliary parameters. As an extension of the least square model, the GWR model is proposed as the typical model to consider the non-stationarity in downscaling LST; however, the non-linear relationship is ignored in the GWR model. In this study, we proposed a non-linear geographically weighted regressive (NL-GWR) model to address this issue. Compared with the GWR model, the NL-GWR added a non-linear auxiliary parameter in the process of building the model. The formula of NL-GWR is as follows:
y i = β 0 ( u i , v i ) + β 1 ( u i , v i ) x i 1 2 + β 2 ( u i , v i ) x i 1 + β 3 ( u i , v i ) x i 2 + + β n ( u i , v i ) x i n + ε i
where, y i means dependent variable; x i 1 , x i 2 , x i n represent different independent variables; x i 1 2 represents the quadratic of parameters, enabling the non-linear relationship between LST and auxiliary parameters; n is the count of independent variables; ( u i , v i ) represents the coordinate of the ith pixel. β 0 ( u i , v i ) refers to the intercept of the regression; the β 1 ( u i , v i ) , β 2 ( u i , v i ) and β n ( u i , v i ) are the 1st, 2nd, and nth regression coefficients of the pixel, respectively. ε i is the random error. Fotheringham et al. used Tobler’s first law of geography to determine the weight and proposed the weighted least square method to estimate the regression parameters [57]. The formula of regression parameters is as follows:
β ^ ( u i , v i ) = ( X T W ( u i , v i ) X ) 1 X T W ( u i , v i ) Y
where β ^ is the estimated value of β ; X , Y are the vectors of independent and dependent variables, respectively; W ( u i , v i ) represents the kernel function, which is used to ensure the weight of observation pixel. The weight of the observation pixel close to pixel i will be larger, as the contribution of the observation pixel far away from pixel i will be relatively small, the weight will also be relatively smaller. The kernel function selected in this study is obtained by the following formula:
W i j = exp ( d i j 2 b 2 )
where d i j is the Euclidean distance between the pixel i and j; b is the adaptive bandwidth, which can calculate by the cross-validation (CV) method of local regression analysis [58]. The relationship between CV and bandwidth can be shown as follows:
C V = 1 n n i = 1 [ y i y ^ i ( b ) ] 2
where y ^ i ( b ) means that the regression parameter estimation does not include the regression pixel itself and builds the relationship between the around pixels. When the relationship between the bandwidth and CV is built, the minimum CV corresponds to the optimal bandwidth.

3.2. LST Downscaling Algorithm Based on NL-GWR Model

In this study, we proposed the LST downscaling algorithm based on the NL-GWR model; a flow chart of the NL-GWR model downscaling LST algorithm is shown in Figure 2.
As shown in Figure 2, the NL-GWR model downscaling LST algorithm can be divided into three parts: data processing, LST downscaling model establishment, and verification and analysis of the downscaling results.
(1)
Data processing. Firstly, the Landsat 8 reflectance data needed preprocessing, including radiometric calibration, atmospheric correction and geometric correction, etc., and then calculated the auxiliary parameters, including NDVI, SAVI, NDBI, UI, and NDWI using the Landsat 8 data, and resampled these indices to 1000 and 100 m, respectively. The data with a spatial resolution of 1000 and 100 m are the input parameters for fitting relationship in the coarse resolution and the fine resolution, respectively. For the MODIS LST data (MOD11_L2, collection 6 data), we used the MODIS reprojection tool (MRT) to registered to a UTM WGS 1984 reference system. In addition, then, the MODIS LST data were used to establish model in 1000 m resolution.
(2)
LST downscaling model establishment. We used the coarse resolution auxiliary parameters and the MODIS LST to establish the NL-GWR model at a resolution of 1000 m, which is as follows:
L S T i C R = β 0 ( u i , v i ) + β 1 C R ( u i , v i ) i n d e x i 1 2 C R + β 2 C R ( u i , v i ) i n d e x i 1 C R + β 3 C R ( u i , v i ) i n d e x i 2 C R + + β n C R ( u i , v i ) i n d e x i ( n 1 ) C R + ε i C R
where the superscript CR denotes the data with a coarse resolution. L S T i C R is the MODIS LST at the pixel i; i n d e x i 1 2 C R , i n d e x i 1 C R , i n d e x i 2 C R and   i n d e x i ( n 1 ) C R represent the quadratic and linear components of the auxiliary parameters at the pixel i, respectively. β 1 C R , β 2 C R , β 3 C R , β n C R are the coefficient of quadratic component and coefficient of multiple one power parameters, respectively; ε i C R is the fitting error at a coarse resolution.
Secondly, the coefficients and error are interpolated to 100 m, which was used to build a regressive relationship at the fine resolution in the Kriging interpolation module in ArcMap 10.2 software.
Finally, we assumed that the fitting relationship between LST and auxiliary parameters is irrelevant to the spatial resolution in the LST downscaling algorithm, which means the fitting relationship established for the coarse spatial resolution can be directly used for fine resolution modeling [37,59,60]; therefore, the relationship between LST and auxiliary parameters for the fine resolution can be expressed as follows:
L S T j F R = β 0 F R ( u j , v j ) + β 1 F R ( u j , v j ) i n d e x j 1 2 F R + β 2 F R ( u j , v j ) i n d e x j 1 F R + β 3 F R ( u j , v j ) i n d e x j 2 F R + + β n F R ( u j , v j ) i n d e x j ( n 1 ) F R + ε j F R
where the superscript FR represents the data at a fine resolution in Formula (9); the L S T j F R is the downscaled LST based on the NL-GWR model.
  • Verification and analysis of the downscaling results. The LST from Landsat 8 retrieved by the Mono-window algorithm was used as reference data to verify the downscaled LST, and root mean square error (RMSE) and mean absolute error (MAE) were chosen as evaluating indicators. RMSE is the deviation between the observed value and its predicted value and illustrates the sample’s dispersion degree. MAE represents the average value of the absolute error between the predicted value and the observed value. The smaller of RMSE and the MAE, the better the downscaling results.

4. Selection of Optimal Index

4.1. Candidates of the Remote Sensing Indices

NDVI can reflect the influence of the vegetation and eliminate the interference of soil and atmosphere [61]. The purpose of proposing SAVI is to correct the sensitivity of NDVI to the soil background. Compared with NDVI, SAVI adds the soil adjustment coefficient L into the formula, which is determined according to actual conditions. The range of SAVI is [0,1], and L = 1 means the vegetation coverage is high and the influence of soil background is low; L = 0 means the vegetation coverage is small. Normally, L = 0.5 can eliminate the influence of soil background, so L = 0.5 was selected in this study [62]. Zha et al. used the NDBI to describe the intensity of urbanization and distinguish building information. The range of NDBI is [−1,1], and areas of NDBI greater than 0 are considered urban land cover [63]. Urban index (UI) was proposed by Kawamura et al. to extract urban areas [64]. NDWI was used to extract water information in the study areas. NDWI is a measurement of liquid water molecules interacting with solar radiation [65] and was adopted in our study as an indicator of water information in the study area. The calculation formulas of the remote sensing indices used in this study are shown in Table 3.
In Table 3, R N I R , R R E D , R S W I R 1 , R S W I R 2 and R G R E E N are the reflectance values of the near-infrared band, red band, the first shortwave infrared band, the second shortwave infrared band, and green band, respectively, and correspond to band 5, band 4, band 6, band 7, and band 3 of Landsat 8, respectively.

4.2. Optimal Index Combination of Research Areas

4.2.1. Downscaling with Single Remote Sensing Index

In this analysis, NDVI, SAVI, NDBI, UI, and NDWI were selected as auxiliary parameters to downscale the LST. The purpose of our first analysis was to choose the most sensitive parameter among the auxiliary parameters for building the downscaling relationship. For the selection of the optimal auxiliary parameters, we choose the stationary least square method that can explain the fitting accuracy of the auxiliary parameters and LST through the R2 of fitting. We used the forward approach by successively adding significant terms into the model to select the optimal indices. The NDVI can reflect the influence of the vegetation and eliminate the interference of soil and atmosphere and is an important index for urban area LST downscaling. Therefore, we selected the NDVI as the first index, establish regressive relationship between Landsat LST and NDVI and calculate the R2. Then, we added other index into regressive relationship and count R2. At last, according to the R2, we determined the optimal combination of auxiliary parameters. As all studied images have similar characteristics, one image for each city (Jinan: 11 July 2014 and Wuhan: 24 July 2016) was used to demonstrate the statistics of optimal index selection in Figure 3, Figure 4, Figure 5, Figure 6 and Figure 7 and Table 4, Table 5, Table 6 and Table 7.
As shown in Figure 3, the trend of the two areas is consistent. Compared with only using NDVI, when SAVI was added as index, the R2 barely increased at all. Then, when NDBI was added as index, the R2 has great improved, from 0.3224 to 0.6657 (Jinan, 11 July 2014) and from 0.2756 to 0.5857 (Wuhan 24 July 2016). Later on, When UI was added into the regressive relationship, the R2 has little improvements. Finally, added NDWI into the regressive relationship, there were even a slight decrease in R2. According to this analysis, SAVI, UI and the NDWI is not selected as the downscaling index. In total, we chose the NDBI, NDVI, and their quadratics as the auxiliary parameters in this study.

4.2.2. Least Square Regression with Combined Remote Sensing Indices

Having analyzed each auxiliary parameter’s sensitivity, we chose the optimal index combination of research areas. At first, we assumed that there are stationary relationships between the auxiliary parameters and LST and tested the performance of the combination of parameters in the sense of least square relationship. At a spatial resolution of 1000 m, we built a least square relationship between Landsat 8 LST and NDVI + NDBI, NDVI2 + NDBI, NDBI2 + NDVI, NDVI2 + NDVI + NDBI, and NDBI2 + NDBI + NDVI. The results of the two study areas are shown in Table 3 and Table 4.
As shown in Table 4 and Table 5, the two parameter combinations of NDVI2 + NDVI + NDBI and NDBI2 + NDBI + NDVI achieved better results in the global least square regression. The combination of NDVI2 + NDBI also obtained better results in the Jinan area. In the global model, the fit results are better, obtain a better R2, and lower RMSE when there are more parameters. Furthermore, the statistical results show that adding the quadratic of the auxiliary parameter can improve the least squares relationship’s accuracy; therefore, it is necessary to consider the non-linear relationship between LST and auxiliary parameters.

4.2.3. Geographical Weighted Regression with Combined Remote Sensing Indices

In fact, there is spatial non-stationarity between LST and the auxiliary parameters during the modeling process, which means the relationship will change when the geographic location changes [66]. In this section, we considered the non-stationary between auxiliary parameters to establish model. The GWR model takes the spatial location into account in the modeling process that can reflect the independent variable’s change rule on the dependent variable with the spatial location.
Table 6 and Table 7 show the GWR and NL-GWR models’ statistical results in the two study areas. As shown in Table 5 and Table 6, the geographical model’s statistical results are better than the least square model. Furthermore, by comparing Table 3, Table 4, Table 5, and Table 6, we can see that adding the quadratic of the parameters into the model has a great effect on the improvement of results; therefore, considering the non-stationarity of the parameters with a geographical model and the non-linear relationship between LST and auxiliary parameters can improve the results, especially in the Wuhan area, where the coefficient of the determination of R2 of the geographical model increased significantly.

5. Results and Discussion

5.1. Downscaling Results

Figure 4 and Figure 5 show the downscaling results at a spatial resolution of 100 m using different parameter combinations. Among them, Figure 4a and Figure 5a are the Landsat 8 LST retrieved with the Mono-window algorithm applied to Jinan (11 July 2014) and Wuhan (24 July 24, 2016), respectively; Figure 4 and Figure 5b are the downscaling results based on GWR and NL-GWR models with the parameter combination of NDVI + NDBI, NDVI2 + NDBI, NDVI + NDBI2, NDVI2 + NDVI + NDBI, and NDBI2 + NDBI + NDVI, respectively. Figure 4 and Figure 5b are the traditional GWR model, whereas Figure 4 and Figure 5c–f added the quadratic of the parameters considering the non-linear relationship between LST and auxiliary parameters.
As shown in Figure 4, the downscaling results with different parameter combinations are consistent with the reference Landsat 8 LST visually speaking. Compared to the NL-GWR model, the GWR model results only consider the linear relationship between the LST and NDVI + NDBI parameter combination, which is substandard in the Yellow River passing through in the north of the Jinan area. The LST of the water is higher than the reference LST because there is a negative linear correlation between the LST and NDVI in the area with high vegetation coverage; however, the NDVI and LST is not a linear relationship in the area with abundant soil moisture. In the previous studies of geographical models, the non-linear relationship of the heterogeneous urban areas and areas with abundant soil water content is neglected, which caused large downscaling errors. The difference is quite large in the north of the study area, where a reservoir is marked with a black rectangle in Figure 4, and different combinations of parameters have different results around the reservoir. As shown in the Landsat 8 LST (Figure 4a), the LST of the reservoir is the lowest, but the LST of downscaling results in Figure 4d,e is slightly higher. From the visible results, the NL-GWR model’s result with the parameter combination NDVI2 + NDBI is closest to the reference data and shows the best visible results.
As Figure 5 shows, Wuhan has a similar visual effect to Jinan. The downscaling image trend is consistent with the reference image, but the downscaled LST has a certain smoothing effect in some areas. Two possible reasons can explain this phenomenon. One reason could be that the processes of regression and interpolation are both based on the minimum mean square error (MMSE) method. There is a common problem for the MMSE method, i.e., that low values tend to be overestimated and high values tend to be underestimated. Another reason is that, in the process of spatial aggregation from 100 to 1000 m of NDBI and NDVI, values of NDVI and NDBI will be replaced by the average values of surrounding pixels, making the resampled image smoother than the original image; this may cause the smoothing effect of the downscaling results. The downscaling LST is the highest in the Yangtze River’s tributary when the NDVI + NDBI is used as the auxiliary parameter combination (marked with a black rectangle in the figure). This situation is not consistent with the reference image, which causes a relatively significant error. By considering the non-linear relationship between LST and auxiliary parameters, the LST of the tributaries has decreased. In the visual images, the downscaling results of the NL-GWR downscaling model with the NDVI2 + NDBI parameter combination are closest to the reference image; this result is consistent with the Jinan area.
In the Wuhan study area, the Yangtze River runs through it, and there are many lakes. Given that there is much water in and around Wuhan, it is intuitive to add the NDWI index into the downscaling model; therefore, the NDWI was added into the optimal parameter combination NDVI2 + NDBI to downscale the LST and verify the accuracy of the downscaling results. The downscaling results of NDVI2 + NDBI and NDVI2 + NDBI + NDWI downscaling were evaluated by analyzing their errors. Figure 6a,b show the error distributions between reference Landsat 8 LST and the downscaling results of the NL-GWR model with NDVI2 + NDBI and NDVI2 + NDBI + NDWI as parameter combinations, respectively. The error range of Figure 6a,b is –4~4 °C. In Figure 6a, the large errors are mainly in rivers, lakes, and high soil water content, and the intersecting areas of different land cover when using the NDVI2 + NDBI parameter combination. Though the errors of using the NDVI2 + NDBI + NDWI parameter combination decreased in the lake, river, and high soil water content areas, in the vegetation, building and other areas were increased, as shown in Figure 6b.
Figure 7 shows the histogram of the error image, and Figure 7a,b correspond to Figure 6a,b. The number of pixels in the interval of [−1,1] in Figure 7a is higher than that in Figure 7b. We selected the mean value and standard deviation (StdDev) as evaluation indicators. The mean value was calculated based on all pixels and is preferably as close to 0 as possible, which indicates unbiased prediction. By comparing Figure 7a,b, we found that the mean values are 0.1271 and −0.2289 °C and the StdDev values are 1.6390 and 1.7131, respectively. So, the mean value of Figure 7a is closer to 0, and the StdDev of Figure 7a is smaller than (b), which means that the dispersion is also smaller.
In summary, adding the NDWI into the downscaling model can improve the downscaling results of water areas; however, the errors of the other land covers have increased. The downscaling results with NDVI2 + NDBI have a larger error in the water areas, but the overall accuracy is higher than that of adding the NDWI and obtaining better downscaling results. The conclusion is that NDWI is not needed in LST downscaling applications.
The visual results show that the downscaling results are improved when considering the non-linear relationship between LST and auxiliary parameters. It is necessary to consider the non-linear relationship in the geographical downscaling model. In the two study areas, the parameter combination of NDVI2 + NDBI obtained the best visual results, and adding the NDWI into the parameter combination does not improve the error statistics of the NL-GWR model.

5.2. Landsat 8 LST as the Reference Data

It is necessary to quantitatively verify the downscaling effect of different parameter combinations instead of relying on simple visual inspection. Ideally, the land surface temperature detected in the ground is the most accurate and direct source for verification data; however, LST cannot accurately be measured from the ground due to its extreme variability as well as scale effect, not to mention the impractical workload to collect representative samples for the whole study area. As a compromise, the inversed high-resolution LST from Landsat 8 data is usually selected as the reference data in the LST downscaling studies.
Significant amounts of research have shown that if there is no ground monitoring data, the LST accuracy retrieved by the MODIS can be used as a reference [67,68,69]. Peng and wang et al. proved that the relationship between the MODIS LST and the LST retrieved by the Landsat 8 is preferable and that the RMSE is less than 2 °C [41,42]. Duan et al. used the ASTER LST as reference data to evaluate the downscaling LST based on the GWR model and demonstrated an error of 2.1 K for the ASTER LST and MODIS LST, which is bigger than the RMSE found between MODIS LST and Landsat 8 retrieved LST at a 1000 m spatial resolution [39]. For this reason, the Landsat 8 LST retrieved by the Mono-window algorithm can be used as reference data to verify the accuracy of downscaling results.
In this study, the Landsat 8 LST was used as reference data to verify the downscaling LST, with RMSE and MAE used as indicators for verification. The statistics of RMSE and MAE of Jinan and Wuhan are shown in Table 8 and Table 9, respectively. The results are based on the GWR and NL-GWR models with different parameter combinations.
As shown in Table 8 and Table 9, the following conclusion can be drawn from the six datasets: the quantitative analysis results are consistent with the visual results; the NL-GWR model with the parameter combination of NDVI2 + NDBI obtains the best downscaling LST. In all the datasets, RMSE and MAE are smaller when considering the non-linear relationship between LST and auxiliary parameters. In the perspective of seasonality, the summer data achieved better downscaling results than other seasons and gained the smallest RMSE (Jinan (11 July 2014): 1.3208 °C, Wuhan (24 July 2016): 1.4957 °C) and MAE (Jinan: 0.9208 °C, Wuhan: 1.0686 °C). The statistical results show that consideration of the non-linear relationship improves the accuracy of downscaling results.
Figure 8 and Figure 9 show the density scatter plots between the reference Landsat 8 LST and the downscaling LST-based NL-GWR model with the parameter combination of NDVI2 + NDBI in the six pairs of study images with a 100 m spatial resolution.
Figure 8 and Figure 9 show that all fitting coefficients of determination are above 0.90 and obtain satisfactory fitting results, indicating that the NL-GWR model can achieve better downscaling results than the GWR model and improve the accuracy of the downscaling algorithm. This conclusion provides new approaches to the study of downscaling algorithms.

6. Conclusions

Due to the spatial heterogeneity of the urban area’s underlying surface, it is necessary to gain LST data with high spatial and temporal resolutions; however, high spatial and temporal resolution can not be simultaneously achieved by current satellite sensors. To solve this problem, we proposed an NL-GWR model, which simultaneously considers the spatial non-stationarity and non-linearity in order to downscale the spatial resolution of MODIS LST from 1000 to 100 m and verify the downscaling results comparing with the Landsat 8 retrieved LST.
In this study, Jinan and Wuhan were selected as study areas. We selected the auxiliary parameters among the NDVI, SAVI, NDBI, UI, and NDWI, and chose the NDVI and NDIBI as auxiliary parameters; then, we selected the different combinations of NDVI + NDBI, NDVI2 + NDBI, NDVI + NDBI2, NDVI2 + NDVI + NDBI, and NDBI2 + NDBI + NDVI to downscale LST; lastly, by comparing and verifying the downscaling results from different parameter combinations, we concluded that the NL-GWR model with the NDVI2 + NDBI parameter combination obtained the best downscaling results. Due to Wuhan’s many water bodies, we tried to add NDWI into the downscaling parameter combination. The downscaling results were improved for areas containing water, but the error in other areas increased; therefore, the NDWI was not used as an auxiliary parameter in this study. The retrieved Landsat 8 LST by the Mono-window algorithm was used as reference data to verify the LST downscaling results. In six datasets, the coefficients of determination between downscaling LST and the Landsat 8 LST both reached above 0.90. The downscaling algorithm based on the NL-GWR model has achieved better downscaling results than the former GWR model, which only considers linear relationships.
Although the NL-GWR model proposed in this study has achieved good downscaling results, many problems still need to be improved in subsequent studies. In this study, only the auxiliary parameters’ quadratic parameters were added to the model the non-linear relationship between LST and auxiliary parameters, although more non-linear relationships, e.g., exponential or rational functions, could be considered. Another aspect that needs improvement is that more indices besides NDVI, SAVI, NDBI, UI, and NDWI should be considered to adapt to the urban land cover diversities. These problems need to be further explored in the future studies.
We consider our study a successful attempt to build a non-linear relationship in geographically weighted regressive algorithms where the downscaling result’s accuracy is significantly improved due to the introduction of non-linear terms. Two cities with different climate, as well as different seasonal data, were investigated, and the proposed method gets better downscaling results in all datasets. Therefore, we believe that the NL-GWR model can be applied to other urban area around the world. However, in extremely heterogeneous areas such as mountains, using nonlinear term in LST downscaling may introduce unstable result and may need extra constraints. It is also recommended that before applying this algorithm to scenarios other than urban areas, selection of more auxiliary parameters and optimization of their combination should be performed. For example, slope angle or incoming solar radiation may be more relevant in mountainous areas than NDBI. The overall procedure was repeatable, and the data are available for many different practices, thus providing a new urban heat island study tool.

Author Contributions

Conceptualization, Q.L. and X.L. (Xiaobo Luo); methodology, S.W.; soft-ware, Y.L. and K.Y.; validation, X.L. (Xia Li), Y.L. and X.L. (Xiuhong Li); writing—original draft preparation, Q.L. and S.W.; writing, review and editing, S.W., Q.L., Y.L., X.L. (Xia Li) and K.Y. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Natural Science Foundation of China (No. 20180913), the National Key Research and Development Program of China (No. 2020YFA0608703), and Research on comprehensive management and system innovation of water environment quality in Beijing Tianjin Hebei region (No. 2018ZX07111).

Acknowledgments

The authors would like to thank the contributions of the anonymous reviewers and the Institute of Remote Sensing and Digital Earth Chinese Academy of Sciences. The authors would also like to thank NASA for providing the satellite data and insightful comments that helped significantly improve this manuscript.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

References

  1. Ebrahimy, H.; Azadbakht, M. Downscaling MODIS land surface temperature over a heterogeneous area: An investigation of machine learning techniques, feature selection, and impacts of mixed pixels. Comput. Geosci. 2019, 124, 93–102. [Google Scholar] [CrossRef]
  2. Weng, Q.; Lu, D. A sub-pixel analysis of urbanization effect on land surface temperature and its interplay with impervious surface and vegetation coverage in Indianapolis, United States. Int. J. Appl. Earth Obs. 2008, 10, 68–83. [Google Scholar] [CrossRef]
  3. Peng, F.; Weng, Q. A time series analysis of urbanization induced land use and land cover change and its impact on land surface temperature with Landsat imagery. Remote Sens. Environ. 2016, 175, 205–214. [Google Scholar]
  4. Sun, Z.; Wang, Q.; Ouyang, Z.; Watanabe, M.; Matsushita, B.; Fukushima, T. Evaluation of MOD16 algorithm using MODIS and ground observational data in winter wheat field in North China Plain. Hydrol. Process. 2007, 21, 1196–1206. [Google Scholar] [CrossRef]
  5. Wu, P.; Shen, H.; Zhang, L. Integrated fusion of multi-scale polar-orbiting and geostationary satellite observations for the mapping of high spatial and temporal resolution land surface temperature. Remote Sens. Environ. 2015, 156, 169–181. [Google Scholar] [CrossRef]
  6. Rajasekar, U.; Weng, Q. Spatio-temporal modelling and analysis of urban heat islands by using Landsat TM and ETM+ imagery. Remote Sens. Environ. 2009, 30, 3531–3548. [Google Scholar] [CrossRef]
  7. Wu, M.; Zhang, X.; Huang, W.; Niu, Z.; Wang, C.; Li, W.; Hao, P. Reconstruction of Daily 30 m Data from HJ CCD, GF-1 WFV, Landsat, and MODIS Data for Crop Monitoring. Remote Sens. 2015, 7, 16293–16314. [Google Scholar] [CrossRef] [Green Version]
  8. Kandasamy, S.; Fernandes, R. An approach for evaluating the impact of gaps and measurement errors on satellite land surface phenology algorithms: Application to 20 year NOAA AVHRR data over Canada. Remote Sens. Environ. 2015, 164, 114–129. [Google Scholar] [CrossRef]
  9. Verbesselt, J.; Hyndman, R.; Zeileis, A.; Culvenor, D. Phenological change detection while accounting for abrupt and gradual trends in satellite image time series. Remote Sens. Environ. 2010, 114, 2970–2980. [Google Scholar] [CrossRef] [Green Version]
  10. Luo, J.; Li, X.; Ma, R.; Li, F.; Duan, H.; Hu, W.; Qin, B.; Huang, W. Applying remote sensing techniques to monitoring seasonal and interannual changes of aquatic vegetation in Taihu Lake, China. Ecol. Indic. 2016, 60, 503–513. [Google Scholar] [CrossRef]
  11. Mo, Y.; Momen, B.; Kearney, M.S. Quantifying moderate resolution remote sensing phenology of Louisiana coastal marshes. Ecol. Model. 2015, 312, 191–199. [Google Scholar] [CrossRef]
  12. Fisher, J.I.; Mustard, J.F. Cross-scalar satellite phenology from ground, Landsat, and MODIS data. Remote Sens. Environ. 2007, 109, 261–273. [Google Scholar] [CrossRef]
  13. Weng, Q.; Fu, P. Modeling diurnal land temperature cycles over Los Angeles using downscaled GOES imagery. ISPRS J. Photogramm. 2014, 97, 78–88. [Google Scholar] [CrossRef]
  14. Sun, L.; Anderson, M.C.; Gao, F. Investigating water use over the Choptank River Watershed using a multisatellite data fusion approach. Water Resour. Res. 2017, 53, 5298–5319. [Google Scholar] [CrossRef] [Green Version]
  15. Zhan, W.; Chen, Y.; Zhou, J.; Wang, J.; Liu, W.; Voogt, J.; Zhu, X.; Quan, J.; Li, J. Disaggregation of remotely sensed land surface temperature: Literature survey, taxonomy, issues, and caveats. Remote Sens. Environ. 2013, 131, 119–139. [Google Scholar] [CrossRef]
  16. Fasbender, D.; Tuia, D.; Bogaert, P. Support-Based Implementation of Bayesian Data Fusion for Spatial Enhancement: Applications to ASTER Thermal Images. IEEE Geosci. Remote Sens. Lett. 2008, 5, 598–602. [Google Scholar] [CrossRef]
  17. Gao, F.; Masek, J.G.; Schwaller, M.R. On the Blending of the Landsat and MODIS Surface Reflectance: Predicting Daily Landsat Surface Reflectance. IEEE Trans. Geosci. Remote Sens. 2006, 44, 2207–2218. [Google Scholar]
  18. Zhu, X.L.; Helmer, E.H.; Gao, F.; Liu, D.S.; Chen, J.; Lefsky, M.A. A flexible spatiotemporal method for fusing satellite images with different resolutions. Remote Sens. Environ. 2016, 172, 165–177. [Google Scholar] [CrossRef]
  19. Bindhu, V.M.; Narasimhan, B.; Sudheer, K.P. Development and verification of a non-linear disaggregation method (NL-DisTrad) to downscale MODIS land surface temperature to the spatial scale of Landsat thermal data to estimate evapotranspiration. Remote Sens. Environ. 2013, 135, 118–129. [Google Scholar] [CrossRef]
  20. Dong, C.; Loy, C.C.; He, K. Learning a Deep Convolutional Network for Image Super-Resolution; Springer: Cham, Switzerland, 2014; pp. 184–199. [Google Scholar]
  21. Dong, C.; Loy, C.C.; Tang, X. Accelerating the Super-Resolution Convolutional Neural Network. European Conference on Computer Vision; Springer: Cham, Switzerland, 2014; pp. 391–407. [Google Scholar]
  22. Ledig, C.; Theis, L.; Huszar, F. Photo-Realistic Single Image Super-Resolution Using a Generative Adversarial Network. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, Honolulu, HI, USA, 21–26 July 2017; pp. 4681–4690. [Google Scholar]
  23. Yu, Y.; Gong, Z.; Zhong, P.; Shan, J. Unsupervised Representation Learning with Deep Convolutional Neural Network for Remote Sensing Images. In International Conference on Image and Graphics; Springer: Cham, Switzerland, 2017; pp. 97–108. [Google Scholar]
  24. Arjovsky, M.; Chintala, S.; Bottou, L. Wasserstein gan. arXiv 2017, arXiv:1701.07875. [Google Scholar]
  25. Berthelot, D.; Schumm, T.; Metz, L. Began: Boundary equilibrium generative adversarial networks. arXiv 2017, arXiv:1703.10717. [Google Scholar]
  26. Shao, B.; Tang, X.; Jin, L.; Li, Z. Single frame infrared image super-resolution algorithm based on generative adversarial nets. J. Infrared Millim. Waves 2018, 37, 427–432. [Google Scholar]
  27. Guo, L.J.; Moore J, M. Pixel block intensity modulation: Adding spatial detail to TM band 6 thermal imagery. Int. J. Remote Sens. 1998, 19, 2477–2491. [Google Scholar] [CrossRef]
  28. Stathopoulou, M.; Cartalis, C. Downscaling AVHRR land surface temperatures for improved surface urban heat island intensity estimation. Remote Sens. Environ. 2009, 113, 2592–2605. [Google Scholar] [CrossRef]
  29. Zhukov, B.; Oertel, D.; Lanzl, F.; Rtinhackel, G. Unmixing-based multisensor multiresolution image fusion. IEEE Trans. Geosci. Remote Sens. 1999, 37, 1212–1226. [Google Scholar] [CrossRef]
  30. Deng, C.; Wu, C. Examining the impacts of urban biophysical compositions on surface urban heat island: A spectral unmixing and thermal mixing approach. Remote Sens. Environ. 2013, 131, 262–274. [Google Scholar] [CrossRef]
  31. Kustas, W.P.; Norman, J.M.; Anderson, M.C.; French, A.N. Estimating subpixel surface temperatures and energy fluxes from the vegetation index–radiometric temperature relationship. Remote Sens. Environ. 2003, 85, 429–440. [Google Scholar] [CrossRef]
  32. Agam, N.; Kustas, W.P.; Anderson, M.C. A vegetation index based technique for spatial sharpening of thermal imagery. Remote Sens. Environ. 2007, 107, 545–558. [Google Scholar] [CrossRef]
  33. Yang, G.; Pu, R.; Huang, W.; Wang, J.; Zhao, C.A. Novel Method to Estimate Subpixel Temperature by Fusing Solar-Reflective and Thermal-Infrared Remote-Sensing Data with an Artificial Neural Network. IEEE Trans. Geosci. Remote Sens. 2010, 48, 2170–2178. [Google Scholar] [CrossRef]
  34. Zhu, S.; Guan, H.; Millington, A.C.; Zhang, G. Disaggregation of land surface temperature over a heterogeneous urban and surrounding suburban area: A case study in Shanghai, China. Int. J. Remote Sens. 2013, 34, 1707–1723. [Google Scholar] [CrossRef]
  35. Eswar, R.; Sekhar, M.; Bhattacharya, B. Disaggregation of LST over India: Comparative analysis of different vegetation indices. Int. J. Remote Sens. 2016, 37, 1035–1054. [Google Scholar] [CrossRef]
  36. Bonafoni, S. Downscaling of Landsat and MODIS Land Surface Temperature Over the Heterogeneous Urban Area of Milan. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 2016, 9, 2019–2027. [Google Scholar] [CrossRef]
  37. Qi, P.; Cui, Y.; Zhang, H. Evaluating Multivariable Statistical Methods for Downscaling Nighttime Land Surface Temperature in Urban Areas. IEEE Access 2020, 8, 162085–162098. [Google Scholar] [CrossRef]
  38. Wang, R.; Gao, W.; Peng, W. Downscale MODIS Land Surface Temperature Based on Three Different Models to Analyze Surface Urban Heat Island: A Case Study of Hangzhou. Remote Sens. 2020, 12, 2134. [Google Scholar] [CrossRef]
  39. Duan, S.-B.; Li, Z.-L. Spatial Downscaling of MODIS Land Surface Temperatures Using Geographically Weighted Regression: Case Study in Northern China. IEEE Trans. Geosci. Remote Sens. 2016, 54, 6458–6469. [Google Scholar] [CrossRef]
  40. Osvaldo, P.; Adolpho, M.; Célia, M.; Lucas, Y. Downscaling of ASTER Thermal Images Based on Geographically Weighted Regression Kriging. Remote Sens. 2018, 10, 633. [Google Scholar]
  41. Peng, Y.; Li, W.; Luo, X. A Geographically and Temporally Weighted Regression Model for Spatial Downscaling of MODIS Land Surface Temperatures Over Urban Heterogeneous Regions. IEEE Trans. Geosci. Remote Sens. 2019, 57, 5012–5027. [Google Scholar] [CrossRef]
  42. Wang, S.; Luo, X.; Peng, Y. Spatial Downscaling of MODIS Land Surface Temperature Based on Geographically Weighted Autoregressive Model. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 2020, 54, 6458–6469. [Google Scholar]
  43. Wan, Z.; Dozier, J. A generalized split-window algorithm for retrieving land-surface temperature from space. IEEE Trans. Geosci. Remote Sens. 1996, 34, 892–905. [Google Scholar]
  44. Duan, S.B.; Li, Z.L.; Wu, H. Radiance-based validation of land surface temperature products derived from Collection 6 MODIS thermal infrared data. Int. J. Appl. Earth Obs. 2018, 70, 84–92. [Google Scholar] [CrossRef]
  45. Duan, S.B.; Li, Z.L.; Cheng, J.; Leng, P. Cross-satellite comparison of operational land surface temperature products derived from MODIS and ASTER data over bare soil surfaces. ISPRS J. Photogramm. Remote Sens. 2017, 126, 1–10. [Google Scholar] [CrossRef]
  46. Qinhuo, L.; Xiru, X.; Jiayi, C. The Retrieval of Land Surface Temperature and Emissivity byRemote Sensing Data: Theory and Digital Simulation. J. Remote Sens. 1998, 2, 1–9. [Google Scholar]
  47. McMillin Larry, M. Estimation of sea surface temperatures from two infrared window measurements with different absorption. Int. J. Remote Sens. 1975, 80, 5113–5117. [Google Scholar] [CrossRef]
  48. Becker, F. The impact of spectral emissivity on the measurement of land surface temperature from a satellite. Int. J. Remote Sens. 1987, 8, 1509–1522. [Google Scholar] [CrossRef]
  49. Qin, Z.; Karnieli, A.; Berliner, P. A mono-window algorithm for retrieving land surface temperature from Landsat TM data and its application to the Israel-Egypt border region. Int. J. Remote Sens. 2001, 22, 3719–3746. [Google Scholar] [CrossRef]
  50. Luo, J.H.; Zhang, J.C.; Huang, W.J.; Yang, G.; Gu, X.; Yang, H. The Analysis of Consistency between HJ-1B and Landsat 5 TM for Retrieving LST Based on the Single-Channel Algorithm. Spect. Anal. 2010, 30, 3285–3289. [Google Scholar]
  51. Hu, D.Y.; Qiao, K.; Wang, X.L. Land surface temperature retrieval from Landsat 8 thermal infrared data using mono-window algorithm. Int. J. Remote Sens. 2015, 19, 964–976. [Google Scholar]
  52. Balcik, F.B.; Ergene, E.M. Determining the impacts of land cover/use categories on land surface temperature using landsat8-oli. Remote Sens. Spat. Inf. Sci. 2016, 41, 251–256. [Google Scholar]
  53. Wang, F.; Qin, Z.; Song, C.; Tu, L.; Karnieli, A.; Zhao, S. An Improved Mono-Window Algorithm for Land Surface Temperature Retrieval from Landsat 8 Thermal Infrared Sensor Data. Remote Sens. 2015, 7, 4268–4289. [Google Scholar] [CrossRef] [Green Version]
  54. Sobrino, J.A.; Coll, C.; Caselles, V. Atmospheric correction for land surface temperature using NOAA-11 AVHRR channels 4 and 5. Remote Sens. Environ. 1991, 38, 19–34. [Google Scholar] [CrossRef]
  55. Jimenez-Munoz, J.C.; Sobrino, J.A. Split-Window Coefficients for Land Surface Temperature Retrieval From Low-Resolution Thermal Infrared Sensors. IEEE Geosci. Remote Sens. 2008, 5, 806–809. [Google Scholar] [CrossRef]
  56. Hengl, T.; Heuvelink, G.B.M.; Tadi, M.P.; Pebesma, E.J. Spatio-temporal prediction of daily temperatures using time-series of MODIS LST image. Theor. Appl. Climatol. 2012, 107, 265–277. [Google Scholar] [CrossRef] [Green Version]
  57. Fotheringham, A.S.; Brunsdon, M.C.; Brunsdon, C. The geography of parameter space: An investigation of spatial non-stationarity. Int. J. Geogr. Inf. Sci. 1996, 10, 605–627. [Google Scholar] [CrossRef]
  58. Calvo, E.; Escolar, M. The Local Voter: A Geographically Weighted Approach to Ecological Inference. Am. J. Political Sci. 2003, 47, 189–204. [Google Scholar] [CrossRef]
  59. Zhan, W.; Chen, Y.; Zhou, J. Sharpening Thermal Imageries: A Generalized Theoretical Framework from an Assimilation Perspective. IEEE Trans. Geosci. Remote Sens. 2011, 49, 773–789. [Google Scholar] [CrossRef]
  60. Zhan, W.; Chen, Y.; Wang, J.; Zhou, J.; Quan, J.; Liu, W.; Li, J. Downscaling land surface temperatures with multi-spectral and multi-resolution images. Int. J. Appl. Earth Obs. 2012, 18, 23–36. [Google Scholar] [CrossRef]
  61. Schnell, J.A. Monitoring the vernal advancement and retrogradation (greenwave effect) of natural vegetation. In Great Plains Corridor; NASA/GSFC Type III Final Report: Greenbelt, MD, USA, 1974. [Google Scholar]
  62. Huete, A.R. A soil-adjusted vegetation index (SAVI). Remote Sens. Environ. 1988, 25, 295–309. [Google Scholar] [CrossRef]
  63. Zha, Y.; Gao, J.; Ni, S. Use of normalized difference built-up index in automatically mapping urban areas from TM imagery. Int. J. Remote Sens. 2003, 24, 583–594. [Google Scholar] [CrossRef]
  64. Kawamura, M. Relation between social and environmental conditions in Colombo Sri Lanka and the urban index estimated by satellite remote sensing data roc. In Proceedings of the 51st Annual Conference of the Japan Society of Civil Engineers, Vienna, Austria, 9–19 July 1996. [Google Scholar]
  65. Gao, B.C. NDWI—A normalized difference water index for remote sensing of vegetation liquid water from space. Remote Sens. Environ. 1996, 58, 257–266. [Google Scholar] [CrossRef]
  66. McMillen, D.P. Geographically Weighted Regression: The Analysis of Spatially Varying Relationships. Am. J. Agric. Econ. 2004, 86, 554–556. [Google Scholar] [CrossRef]
  67. Mutiibwa, D.; Strachan, S.; Albright, T. Land Surface Temperature and Surface Air Temperature in Complex Terrain. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 2016, 8, 4762–4774. [Google Scholar] [CrossRef]
  68. Gemma, S.; Vicente, G.; Maria, J.; Martínez-Villagrasa, D.; Picos, R.; Caselles, V.; Cuxart, J. Landsat and Local Land Surface Temperatures in a Heterogeneous Terrain Compared to MODIS Values. Remote Sens. 2016, 8, 849. [Google Scholar]
  69. Zhu, X.; Wang, X.; Yan, D. Analysis of remotely-sensed ecological indexes’ influence on urban thermal environment dynamic using an integrated ecological index: A case study of Xi’an, China. Int. J. Remote Sens. 2019, 40, 3421–3447. [Google Scholar] [CrossRef]
Figure 1. Landsat 8 false-color images (R: band 5, G: band 4, B: band 3) in Jinan and Wuhan.
Figure 1. Landsat 8 false-color images (R: band 5, G: band 4, B: band 3) in Jinan and Wuhan.
Remotesensing 13 01580 g001
Figure 2. Flow chart of the non-linear geographically weighted regressive (NL-GWR) model downscaling LST algorithm.
Figure 2. Flow chart of the non-linear geographically weighted regressive (NL-GWR) model downscaling LST algorithm.
Remotesensing 13 01580 g002
Figure 3. R2 statistics of establishing regression relationship with different combination of indices. (a) and (b) are the R2 statistics of Jinan (11 July 2014) and Wuhan (24 July 2016), respectively.
Figure 3. R2 statistics of establishing regression relationship with different combination of indices. (a) and (b) are the R2 statistics of Jinan (11 July 2014) and Wuhan (24 July 2016), respectively.
Remotesensing 13 01580 g003
Figure 4. Downscaling results at a resolution of 100 m using different parameter combinations in Jinan (11 July 2014). (a) is the Landsat 8 LST retrieved by Mono-window algorithm; (bf) are the downscaling results based on GWR and NL-GWR models and with the parameter combinations of NDVI + NDBI, NDVI2 + NDBI, NDVI + NDBI2, NDVI2 + NDVI + NDBI, and NDBI2 + NDBI + NDVI, respectively.
Figure 4. Downscaling results at a resolution of 100 m using different parameter combinations in Jinan (11 July 2014). (a) is the Landsat 8 LST retrieved by Mono-window algorithm; (bf) are the downscaling results based on GWR and NL-GWR models and with the parameter combinations of NDVI + NDBI, NDVI2 + NDBI, NDVI + NDBI2, NDVI2 + NDVI + NDBI, and NDBI2 + NDBI + NDVI, respectively.
Remotesensing 13 01580 g004
Figure 5. Downscaling results at a resolution of 100 m using different parameter combinations in Wuhan (24 July 2016). (a) is the Landsat 8 LST retrieved by Mono-window algorithm; (bf) are the downscaling results based on GWR and NL-GWR models and with the parameter combinations of NDVI + NDBI, NDVI2 + NDBI, NDVI + NDBI2, NDVI2 + NDVI + NDBI and NDBI2 + NDBI + NDVI, respectively.
Figure 5. Downscaling results at a resolution of 100 m using different parameter combinations in Wuhan (24 July 2016). (a) is the Landsat 8 LST retrieved by Mono-window algorithm; (bf) are the downscaling results based on GWR and NL-GWR models and with the parameter combinations of NDVI + NDBI, NDVI2 + NDBI, NDVI + NDBI2, NDVI2 + NDVI + NDBI and NDBI2 + NDBI + NDVI, respectively.
Remotesensing 13 01580 g005
Figure 6. Error distribution between reference Landsat 8 LST and the downscaling results in Wuhan, 24 July 2016. (a,b) are the NL-GWR model results with NDVI2 + NDBI and NDVI2 + NDBI + NDWI parameter combinations, respectively.
Figure 6. Error distribution between reference Landsat 8 LST and the downscaling results in Wuhan, 24 July 2016. (a,b) are the NL-GWR model results with NDVI2 + NDBI and NDVI2 + NDBI + NDWI parameter combinations, respectively.
Remotesensing 13 01580 g006
Figure 7. The histogram of error image corresponding to Figure 6. (a,b) are the error statistics of the parameter combination of NDVI2 + NDBI and NDVI2 + NDBI + NDWI.
Figure 7. The histogram of error image corresponding to Figure 6. (a,b) are the error statistics of the parameter combination of NDVI2 + NDBI and NDVI2 + NDBI + NDWI.
Remotesensing 13 01580 g007
Figure 8. Density scatter plots between the reference Landsat 8 LST and the downscaling LST with a spatial resolution of 100 m in Jinan: (a) 11 July 2014; (b) 25 April 2015; (c) 10 December 2017.
Figure 8. Density scatter plots between the reference Landsat 8 LST and the downscaling LST with a spatial resolution of 100 m in Jinan: (a) 11 July 2014; (b) 25 April 2015; (c) 10 December 2017.
Remotesensing 13 01580 g008
Figure 9. Density scatter plots between the reference Landsat 8 LST and the downscaling LST with a spatial resolution of 100 m in Wuhan: (a) 23 January 2014; (b) 24 July 2015; (c) 30 October 2017.
Figure 9. Density scatter plots between the reference Landsat 8 LST and the downscaling LST with a spatial resolution of 100 m in Wuhan: (a) 23 January 2014; (b) 24 July 2015; (c) 30 October 2017.
Remotesensing 13 01580 g009
Table 1. The Landsat 8 and moderate resolution imaging spectroradiometer (MODIS) land surface temperature (LST) data collected in this study.
Table 1. The Landsat 8 and moderate resolution imaging spectroradiometer (MODIS) land surface temperature (LST) data collected in this study.
Study AreaAcquisition Time
(Landsat 8 Data)
Acquisition Time
(MODIS LST Data)
Jinan11 July 2014 02:48:2311 July 2014 03:00:00
25 April 2015 02:47:5125 April 2015 02:55:00
10 December 2017 02:48:4010 December 2017 03:00:00
Wuhan23 January 2014 02:57:2623 January 2014 04:00:00
24 July 2016 02:56:1724 July 2016 03:05:00
30 October 2017 02:56:3630 October 03:05:00
Table 2. The value of a 10 and b 10
Table 2. The value of a 10 and b 10
Range of LST (°C).a10b10R2
20–70 70.1775 0.4581 0.9997
0–50 62.7182 0.4339 0.9996
−20–30 55.4276 0.4086 0.9996
Table 3. Remote sensing indices used and their calculation formulas.
Table 3. Remote sensing indices used and their calculation formulas.
IndicesAbbreviationFormulation
Normalized Difference Vegetation IndexNDVI N D V I = R N I R R R E D R N I R + R R E D (10)
Soil-Adjusted Vegetation IndexSAVI S A V I = ( R N I R R R E D ) ( 1 + L ) R N I R + R R E D + L ( L = 0.5 ) (11)
Normalized Difference Built-up IndexNDBI N D B I = R S W I R 1 R N I R R S W I R 1 + R N I R (12)
Urban indexUI U I = R S W I R 2 R N I R R S W I R 2 + R N I R (13)
Normal Difference Water IndexNDWI N D W I = R G R E E N R N I R R G R E E N + R N I R (14)
Table 4. Least squares relationship results in Jinan (11 July 2014).
Table 4. Least squares relationship results in Jinan (11 July 2014).
RelationshipR2RMSE (°C)
L S T = a + b N D V I + c N D B I 0.661.38
L S T = a + b N D V I 2 + c N D B I 0.711.26
L S T = a + b N D V I + c N D B I 2 0.601.61
L S T = a + b N D V I 2 + c N D V I + d N D B I 0.751.15
L S T = a + b N D B I 2 + c N D B I + d N D V I 0.721.22
Table 5. Least squares relationship results in Wuhan (24 July 2016).
Table 5. Least squares relationship results in Wuhan (24 July 2016).
RelationshipR2RMSE (°C)
L S T = a + b N D V I + c N D B I 0.420.96
L S T = a + b N D V I 2 + c N D B I 0.540.72
L S T = a + b N D V I + c N D B I 2 0.600.63
L S T = a + b N D V I 2 + c N D V I + d N D B I 0.700.59
L S T = a + b N D B I 2 + c N D B I + d N D V I 0.710.54
Table 6. Geographical model results of Jinan (11 July 2014).
Table 6. Geographical model results of Jinan (11 July 2014).
RelationshipR2RMSE (°C)
L S T = a ( u i , v i ) + b ( u i , v i ) N D V I + c ( u i , v i ) N D B I 0.800.99
L S T = a ( u i , v i ) + b ( u i , v i ) N D V I 2 + c ( u i , v i ) N D B I 0.860.85
L S T = a ( u i , v i ) + b ( u i , v i ) N D V I + c ( u i , v i ) N D B I 2 0.810.98
L S T = a ( u i , v i ) + b ( u i , v i ) N D V I 2 + c ( u i , v i ) N D V I + d ( u i , v i ) N D B I 0.810.94
L S T = a ( u i , v i ) + b ( u i , v i ) N D B I 2 + c ( u i , v i ) N D B I + d ( u i , v i ) N D V I 0.830.89
Table 7. Geographical model results of Wuhan (24 July 2016).
Table 7. Geographical model results of Wuhan (24 July 2016).
RelationshipR2RMSE (°C)
L S T = a ( u i , v i ) + b ( u i , v i ) N D V I + c ( u i , v i ) N D B I 0.600.76
L S T = a ( u i , v i ) + b ( u i , v i ) N D V I 2 + c ( u i , v i ) N D B I 0.820.39
L S T = a ( u i , v i ) + b ( u i , v i ) N D V I + c ( u i , v i ) N D B I 2 0.770.60
L S T = a ( u i , v i ) + b ( u i , v i ) N D V I 2 + c ( u i , v i ) N D V I + d ( u i , v i ) N D B I 0.790.53
L S T = a ( u i , v i ) + b ( u i , v i ) N D B I 2 + c ( u i , v i ) N D B I + d ( u i , v i ) N D V I 0.760.57
Table 8. Accuracy statistics of downscaling results in Jinan.
Table 8. Accuracy statistics of downscaling results in Jinan.
Parameter
Combination
11 July 201425 April 201510 December 2017
RMSEMAERMSEMAERMSEMAE
NDVI + NDBI2.49971.40542.55161.46872.33281.3574
NDVI2 + NDBI1.32080.92081.48580.97211.52311.0063
NDVI + NDBI21.89261.24511.79761.15262.18101.2762
NDVI2 + NDVI + NDBI2.02241.37082.13641.38402.31801.4308
NDBI2 + NDBI + NDVI1.69211.17642.12211.27652.31001.2364
Table 9. Accuracy statistics of downscaling results in Wuhan.
Table 9. Accuracy statistics of downscaling results in Wuhan.
Parameter
Combination
23 January 201424 July 201630 October 2017
RMSEMAERMSEMAERMSEMAE
NDVI + NDBI2.91172.10172.81651.97572.03431.4837
NDVI2 + NDBI1.72351.23841.49571.06861.41681.1172
NDVI + NDBI22.56112.16093.02111.15241.68211.2407
NDVI2 + NDVI + NDBI2.50392.14362.02241.28441.85661.2280
NDBI2 + NDBI + NDVI2.54452.21402.52541.62981.79751.2451
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Wang, S.; Luo, Y.; Li, X.; Yang, K.; Liu, Q.; Luo, X.; Li, X. Downscaling Land Surface Temperature Based on Non-Linear Geographically Weighted Regressive Model over Urban Areas. Remote Sens. 2021, 13, 1580. https://0-doi-org.brum.beds.ac.uk/10.3390/rs13081580

AMA Style

Wang S, Luo Y, Li X, Yang K, Liu Q, Luo X, Li X. Downscaling Land Surface Temperature Based on Non-Linear Geographically Weighted Regressive Model over Urban Areas. Remote Sensing. 2021; 13(8):1580. https://0-doi-org.brum.beds.ac.uk/10.3390/rs13081580

Chicago/Turabian Style

Wang, Shumin, Youming Luo, Xia Li, Kaixiang Yang, Qiang Liu, Xiaobo Luo, and Xiuhong Li. 2021. "Downscaling Land Surface Temperature Based on Non-Linear Geographically Weighted Regressive Model over Urban Areas" Remote Sensing 13, no. 8: 1580. https://0-doi-org.brum.beds.ac.uk/10.3390/rs13081580

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop