A refined downscaling method of land surface temperature combining temporal and spatial lags

Abstract

The present application relates to a kind of surface temperature fine downscaling method of spatio-temporal hysteresis, belong to remote sensing image processing technical field.The method is first to Landsat8TIRS data is preprocessed and LST inversion, and calculate MODIS LST and DEM data.Subsequently, construct time series dataset, and establish the nonlinear relationship between MODIS LST and Landsat8NDVI, NDBI, MNDWI and SRTM DEM data based on spatio-temporal weighted lag random forest model and the spatio-temporal lag relationship of two consecutive time periods LST.Last, ordinary kriging interpolation method is used to downscale LST, obtain high-resolution ground temperature distribution.The method effectively improves the accuracy and reliability of LST downscaling result, provides scientific basis for mountain city climate regulation, ecological protection, thermal environment management etc..

Description

Technical Field

[0001] The present invention belongs to the technical field of remote sensing image processing and relates to a surface temperature fine downscaling method combining temporal and spatial hysteresis. Background Art

[0002] Land surface temperature (LST) is a key physical parameter that describes surface-atmosphere interactions and reflects the Earth's surface. It integrates the interaction between the Earth's surface and the atmosphere, as well as the energy exchange between the atmosphere and the land. As a key parameter in numerous basic disciplines and applied fields, LST is widely used in research areas such as urban thermal environment, agricultural drought, soil moisture estimation, ground-air water-heat exchange, and global climate change. In mountainous cities, the spatiotemporal distribution and evolution of LST is an extremely complex and dynamic issue, influenced by a variety of complex factors. These factors are the result of the intertwined interaction of multiple spatiotemporal characteristics, including the complexity of atmospheric motion, the direct impact of human activities, and spatiotemporal hysteresis and heterogeneity. Spatiotemporal hysteresis refers to the time delay or spatial dislocation between a geographical phenomenon or process and its preceding or subsequent moments. This is particularly evident in the distribution of LST in mountainous cities. For example, the impact of land use changes during urbanization on surface temperature often has a lagged effect. Furthermore, the spatiotemporal heterogeneity of LST in mountainous cities is more pronounced. This is due to the unique topographic and climatic conditions of mountainous cities, such as highly variable topographic changes, localized climate effects, and land use variations during urbanization. These factors lead to highly uneven spatial distribution of LST. For example, temperature differences between urban cores and surrounding rural areas can be quite pronounced, reflecting regional variations in heat island effects and vegetation cover. Furthermore, in mountainous cities, a more pronounced nonlinear relationship often exists between scale factors and LST. This suggests that the influence of scale factors on LST is not a simple linear superposition, but rather involves complex interactions and feedback mechanisms. For example, the influence of vegetation cover on surface temperature exhibits a distinct nonlinear characteristic. Therefore, when studying LST in mountainous cities, it is important to comprehensively analyze these complex factors to predict LST in mountainous cities with fine resolution. This will have important practical applications in climate regulation, ecological protection, heat wave response, and thermal environment management in mountainous cities.

[0003] Based on their principles, LST downscaling methods can be primarily categorized as those based on image fusion and statistical models based on scale-invariant relationships. Image fusion methods integrate multi-source data using methods such as similar pixel weighting or mixed pixel decomposition to generate dense, temporally dense, high-spatial-resolution images. Spatiotemporal fusion models, such as the flexible spatiotemporal data fusion method and the improved spatiotemporal data fusion method, were initially developed for reflectance data, but they can still be used to enhance LST. An LST spatiotemporal fusion model based on an improved U-STFM enhanced the spatial resolution of MODIS LST. The unbiased ESTARFM spatiotemporal fusion model highlights the potential of using LST data acquired from geostationary satellites, such as the Himawari-8. Spatiotemporal image fusion methods have gained widespread use due to their limited number of auxiliary parameters. However, these fusion methods rely on the quality of high-resolution LST and ignore the explicit physical context of satellite remote sensing data, making them inadequate for quantitative thermal infrared remote sensing.

[0004] LST downscaling models based on scale-invariant relationships have clear physical meaning, maintain consistency in thermal radiation information before and after downscaling, and are widely used. Their basic concept is to establish statistical regression relationships between LST and scaling factors such as vegetation index, building index, and topographic factors at low spatial resolution. These relationships are then applied to high-resolution scaling factors to achieve LST downscaling. Classic downscaling models include the DisTrad model and the TsHARP model, which use the Normalized Difference Vegetation Index (NDVI) and vegetation cover as scaling factors, respectively, and have achieved good LST downscaling results. However, these downscaling methods are unsuitable for urban and arid regions due to the limitations of a single variable and the simplified regression relationship. Surface albedo data is introduced into the TsHARP model to obtain high-resolution urban LST. TsHARP is coupled with the STARFM model to generate a robust LST downscaling framework. Based on the random forest method, a complex nonlinear relationship between multiple spectral indices and LST is established to achieve better downscaling results. A weighted LST downscaling method using kernel functions and fusion methods can achieve high-spatial-resolution daily LST data. The kernel-driven method is integrated into the fusion method, and a simple and effective downscaling (SED) algorithm is proposed to obtain LST with high spatiotemporal resolution.

[0005] Based on a downscaling model that maintains local scale invariance, the mapping relationship between LST and scale factors varies with local location. Duan et al., taking into account spatial heterogeneity, first applied the geographically weighted regression (GWR) model to the study of LST downscaling. This model investigated the local relationship between LST and vegetation indices and DEMs. The GWR model was combined with kriging interpolation to generate geographically weighted kriging models and multiscale geographically weighted kriging models. Some researchers, taking into account spatiotemporal constraints such as LST spatial autocorrelation and spatiotemporal nonstationarity, as well as local nonlinear mapping relationships, have constructed downscaling models based on spatiotemporal geographic weighting, geographically weighted autoregression, spatiotemporal geographic weighted autoregression, and nonlinear geographical weighting, achieving excellent performance. Wu et al. proposed a weighted autoregressive model based on geographic and temporal neural networks to improve the overall downscaling accuracy. Mu et al. proposed a geographically weighted random forest downscaling model by coupling the geographically weighted model with the random forest algorithm. This model, based on the geographically weighted regression model, incorporates the random forest regression model to account for the nonstationarity and nonlinearity of the relationship between surface temperature and scale factors.

[0006] Overall, while surface temperature downscaling models have yielded numerous valuable research results, they still face numerous challenges, primarily in the following areas: First, the consistent consideration of spatiotemporal constraints is insufficient. Existing downscaling methods often fail to account for the spatiotemporal lags of LST. Consequently, many downscaling models fail to effectively correct for the temporal variations of LST, struggling to accurately capture its dynamics across time periods. This can lead to biased results in practical applications. Second, the complex nonlinear relationship between LST and scale factors in mountainous cities remains under-considered. Existing models still struggle to adequately address this nonlinear relationship, limiting their adaptability to complex environments. Finally, LST exhibits significant nonlinear characteristics in both spatial and temporal dimensions. Accurately capturing these variations during downscaling remains a pressing challenge. Summary of the Invention

[0007] In view of this, the purpose of the present invention is to provide a refined surface temperature downscaling method that combines temporal and spatial hysteresis. This method incorporates the impact of temporal and spatial hysteresis on LST, comprehensively considers the temporal and spatial hysteresis of LST, temporal and spatial heterogeneity, the nonlinear relationship between LST and scale factors, and the nonlinear variation of LST in time and space dimensions, providing a new approach for downscaling complex and changeable mountainous urban LST from coarse spatial resolution to fine spatial resolution.

[0008] In order to achieve the above object, the present invention provides the following technical solutions:

[0009] A surface temperature fine-grained downscaling method incorporating temporal and spatial hysteresis includes the following steps:

[0010] The acquired 30m spatial resolution Landsat8 TIRS data of the corresponding area were subjected to radiometric correction, atmospheric correction, and geometric correction, and NDVI, NDBI, and MNDWI were calculated. The atmospheric correction method was used to invert the LST of the Landsat8 TIRS data.

[0011] The MODIS LST is calculated using the 1000m spatial resolution MOD11A1 product data provided by MODIS;

[0012] The terrain correction and terrain factor data extraction were completed by using the acquired 90m spatial resolution DEM data;

[0013] The Landsat8 NDVI, NDBI, and MNDWI images with a spatial resolution of 30 m and the SRTM DEM images with a spatial resolution of 90 m were resampled to 100 m and 1000 m spatial resolutions by using the cubic convolution interpolation method.

[0014] A time series dataset of the study area was constructed, and the MODIS LST of two consecutive time periods were arranged in chronological order and used as the dependent variable in the LST downscaling process;

[0015] The aggregated Landsat8 NDVI, NDBI, and MNDWI data, 1000-meter spatial resolution SRTM DEM data, and land cover type data were arranged in chronological order to construct a time series dataset. This dataset was used to establish the nonlinear relationship between MODIS LST and the coarse spatial resolution Landsat8 NDVI, NDBI, MNDWI data and SRTM DEM data, as well as the spatiotemporal lag relationship between LST in two consecutive time periods.

[0016] The aggregated Landsat8 NDVI, NDBI, MNDWI data and 100-meter spatial resolution SRTM DEM data were arranged in chronological order to construct a time series dataset of independent variables for high-resolution LST downscaling.

[0017] Using the spatiotemporal series dataset of explanatory variables and MODIS LST, the temporal and spatial heterogeneity relationship between the current MODIS LST and Landsat8 NDVI, NDBI, MNDWI and SRTM DEM, the spatiotemporal lag relationship with the previous MODIS LST, and the spatiotemporal nonlinear relationship with all explanatory variables were established based on the GTWLRF model at a coarse spatial resolution.

[0018] The LST was downscaled and the regression residuals at 1000m spatial resolution were interpolated to 100m spatial resolution using ordinary kriging interpolation;

[0019] Assuming that the spatiotemporal lag, spatiotemporal heterogeneity, and spatiotemporal nonlinear relationship between LST and explanatory variables are scale invariant, the regression relationship between the current spatiotemporal LST and the current spatiotemporal NDVI, NDBI, MNDWI, DEM, and the previous spatiotemporal LST at a spatial resolution of 100 m after downscaling is expressed as follows:

[0020]

[0021] In the formula represents the Landsat8LST with a spatial resolution of 100m at point i at time t-1, represents the k points adjacent to point i at time t-1, with a Landsat8 LST spatial resolution of 100m. W represents the set of Landsat8 NDVI, NDBI, MNDWI and SRTM DEM with a spatial resolution of 100m at point i at time t, respectively. 100 (u i , v i , τ i ) represents the spatiotemporal weight matrix of point i with a spatial resolution of 100m, represents the regression residual of point i at time t from 1000m spatial resolution to 100m through ordinary kriging interpolation. It is the downscaled LST result of point i at time t100m spatial resolution.

[0022] Furthermore, the GTWLRF model is shown as follows:

[0023]

[0024] Where (u i , v i , τ i ) is the space-time position of point i, W(u i , v i , τ i ) it is the spatiotemporal adjacency matrix, X it is the independent variable matrix at point i at time t, y i(t-1) is the time lag term, is the spatial lag term, y it is the target variable at time t, ∑ p f p is a random forest regression model, p represents the number of decision trees, ε it is the error term; where the spatiotemporal weight matrix W(u i , v i , τ i )=diag(wi1 , w i2 ,…,w in ) is a square matrix; the diagonal elements w ij Represents the geographical and temporal weights of target point i to surrounding observation points j.

[0025] Furthermore, the spatiotemporal weight matrix W(u i , v i , τ i ) by using space-time distance and is constructed based on the Gaussian distance decay function; the spatiotemporal weight matrix W(u i , v i , τ i )'s diagonal elements w ij Obtained by the following formula:

[0026]

[0027] Where h sT It is the space-time bandwidth, which is used to control the decay rate of the weight and determine the degree of influence of distance on the weight. The calculation of space-time distance includes the weight combination of space and time, as shown in the following formula:

[0028]

[0029] In the formula is the spatial distance, is the time distance, λ and are weight coefficients for balancing the spatial and temporal scales; λ is set to 1, and μ is determined by the coefficient of determination R 2 Statistics and modified Akaike information criterion are determined and optimized by cross validation; spatiotemporal weight matrix W(u i , v i , τ i )'s diagonal elements w ij Rewritten as:

[0030]

[0031] where h S is the time bandwidth, h T is the spatial bandwidth; take the spatiotemporal distance from the target point i to the kth neighboring point as the adaptive kernel bandwidth h ST ,Right now Ensure that the kernel bandwidth of each target point i is dynamically adjusted based on the distance of its neighborhood, and the weight calculation automatically adapts to different point densities.

[0032] Furthermore, the LST downscaling process is shown as follows:

[0033]

[0034] In the formula represents the MODIS LST with a spatial resolution of 1000m at point i at time t, represents the MODIS LST with a spatial resolution of 1000m at point i at time t-1, represents the MODIS LST of k points adjacent to point i at time t-1 with a spatial resolution of 1000m, W represents the set of Landsat8 NDVI, NDBI, MNDWI and SRTM DEM with a spatial resolution of 1000m at point i at time t, respectively. 1000 (u i , v i , τ i ) represents the spatiotemporal weight matrix of point i with a spatial resolution of 1000m, represents the regression residual at point i at time t with a spatial resolution of 1000m.

[0035] Furthermore, the LST downscaling process is shown as follows:

[0036]

[0037] In the formula represents the Landsat8LST with a spatial resolution of 100m at point i at time t-1, represents the k points adjacent to point i at time t-1, with a Landsat8 LST spatial resolution of 100m. W represents the set of Landsat8 NDVI, NDBI, MNDWI and SRTM DEM with a spatial resolution of 100m at point i at time t, respectively. 100 (u i ,v i ,τ i ) represents the spatiotemporal weight matrix of point i with a spatial resolution of 100m, represents the regression residual of point i at time t from 1000m spatial resolution to 100m through ordinary kriging interpolation. It is the downscaled LST result of point i at time t100m spatial resolution.

[0038] The beneficial effects of the present invention are:

[0039] First, this paper proposes a refined surface temperature downscaling method that incorporates spatiotemporal hysteresis, innovatively incorporating the impact of spatiotemporal hysteresis on LST. The spatiotemporal hysteresis effect of LST refers to its delayed response to various influencing factors in both time and space. This effect is of great significance in ecological and environmental research. First, LST temporally lags behind changes in solar radiation. If LST was high at a previous point in time, the surface will also store relatively more heat, which may result in a relatively high LST at the current point in time, even if external conditions such as solar radiation have not changed significantly. For example, although solar radiation peaks at noon, the surface temperature typically peaks after noon because the surface takes time to absorb and conduct heat. Similar delays also occur in seasonal variations. For example, after an increase in solar radiation in spring, the surface temperature rises with a lag due to the influence of winter snow cover and soil moisture. Spatially, the conduction of surface heat causes temperature changes to gradually diffuse. For example, the urban heat island effect causes temperature changes to be transmitted to surrounding areas, resulting in a spatial delay in temperature changes. Topography also influences this effect. Valleys and low-lying areas experience slower temperature drops at night due to the accumulation of cold air, while higher altitudes experience faster drops. Different land cover types also respond differently to temperature changes. For example, evapotranspiration in vegetation-covered areas slows temperature changes, causing them to lag behind those of bare land. Furthermore, atmospheric conditions such as cloud cover and precipitation interact with the physical properties of the Earth's surface, further exacerbating the complexity of the LST's spatiotemporal lag. This spatiotemporal lag helps us better understand the complex interactions between surface temperature and the climate system, and has broad application in fields such as climate prediction, urban planning, and agricultural management. Combining temporal and spatial lags in analysis can more comprehensively capture the changing characteristics of LST and improve its prediction accuracy.

[0040] Second, the present invention proposes a land surface temperature (LST) downscaling model suitable for mountainous cities - the temporal and spatial weighted lag random forest model (GTWLRF). During the model design process, the temporal and spatial lag, temporal and spatial heterogeneity of LST, and the complex nonlinear relationship between LST and scale factors are fully considered. Traditional LST downscaling methods are usually difficult to accurately reflect the dynamic change characteristics of LST in the temporal and spatial dimensions, especially in complex terrain such as mountainous cities. By introducing a temporal and spatial weighting mechanism, the GTWLRF model effectively captures the temporal and spatial nonlinear relationship between LST and scale factors, thereby improving the adaptability of the model at different time nodes and spatial locations. At the same time, the lag consideration in the model can utilize LST data at past moments to improve the prediction accuracy of the current LST. In addition, the random forest algorithm, as the core of the model, not only has powerful nonlinear fitting capabilities, but can also effectively reduce the noise impact brought by multi-dimensional feature input and enhance the robustness of the model. The advantage of this approach is that the GTWLRF model, taking advantage of the rugged terrain and high environmental diversity found in mountainous cities, can more accurately capture local trends in LST and achieve refined predictions of its spatial distribution. This approach has important applications in climate change research, urban thermal environment assessment, and ecological and environmental monitoring, and also provides a scientific basis for refined environmental management in mountainous cities. This model can better understand the spatiotemporal variations of LST in mountainous cities, thereby supporting urban planning, ecological and environmental protection, and sustainable development.

[0041] Other advantages, objects, and features of the present invention will be described in part in the following description and, in part, will be apparent to those skilled in the art upon examination of the following description or may be learned from practice of the present invention. The objects and other advantages of the present invention may be realized and obtained through the following description. BRIEF DESCRIPTION OF THE DRAWINGS

[0042] In order to make the purpose, technical solutions and advantages of the present invention more clear, the present invention will be described in detail below with reference to the accompanying drawings, in which:

[0043] Figure 1 This is a flow chart of the surface temperature refined downscaling method combining temporal and spatial hysteresis proposed in the present invention;

[0044] Figure 2 Schematic diagram of the space-time lag effect;

[0045] Figure 3 This is a schematic diagram of the downscaling effect of the method of the present invention;

[0046] Figure 4 This is the absolute error map of downscaling using the method of the present invention. DETAILED DESCRIPTION

[0047] The following describes the embodiments of the present invention by means of specific examples, and those skilled in the art can easily understand other advantages and effects of the present invention from the contents disclosed in this specification. The present invention can also be implemented or applied through other different specific embodiments, and the details in this specification can also be modified or changed in various ways based on different viewpoints and applications without departing from the spirit of the present invention. It should be noted that the illustrations provided in the following embodiments are only schematic illustrations of the basic concept of the present invention, and the following embodiments and features in the embodiments can be combined with each other without conflict.

[0048] Among them, the accompanying drawings are only for illustrative purposes and represent only schematic diagrams rather than actual pictures, and should not be understood as limiting the present invention. In order to better illustrate the embodiments of the present invention, some parts of the accompanying drawings may be omitted, enlarged or reduced, and do not represent the dimensions of actual products. For those skilled in the art, it is understandable that some well-known structures and their descriptions may be omitted in the accompanying drawings.

[0049] The same or similar numbers in the drawings of the embodiments of the present invention correspond to the same or similar parts; in the description of the present invention, it should be understood that if there are terms such as "upper", "lower", "left", "right", "front", "back", etc. indicating directions or positional relationships, they are based on the directions or positional relationships shown in the drawings. They are only for the convenience of describing the present invention and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific direction, be constructed and operate in a specific direction. Therefore, the terms describing the positional relationship in the drawings are only used for illustrative purposes and cannot be understood as limiting the present invention. For ordinary technicians in this field, the specific meanings of the above terms can be understood according to specific circumstances.

[0050] First, this paper proposes a refined surface temperature downscaling method that incorporates spatiotemporal hysteresis, innovatively incorporating the impact of spatiotemporal hysteresis on LST. The spatiotemporal hysteresis effect of LST refers to its delayed response to various influencing factors in both time and space. This effect is of great significance in ecological and environmental research. First, LST temporally lags behind changes in solar radiation. If LST was high at a previous point in time, the surface will also store relatively more heat, which may result in a relatively high LST at the current point in time, even if external conditions such as solar radiation have not changed significantly. For example, although solar radiation peaks at noon, the surface temperature typically peaks after noon because the surface takes time to absorb and conduct heat. Similar delays also occur in seasonal variations. For example, after an increase in solar radiation in spring, the surface temperature rises with a lag due to the influence of winter snow cover and soil moisture. Spatially, the conduction of surface heat causes temperature changes to gradually diffuse. For example, the urban heat island effect causes temperature changes to be transmitted to surrounding areas, resulting in a spatial delay in temperature changes. Topography also influences this effect. Valleys and low-lying areas experience slower temperature drops at night due to the accumulation of cold air, while higher altitudes experience faster drops. Different land cover types also respond differently to temperature changes. For example, evapotranspiration in vegetation-covered areas slows temperature changes, causing them to lag behind those of bare land. Furthermore, atmospheric conditions such as cloud cover and precipitation interact with the physical properties of the Earth's surface, further exacerbating the complexity of the LST's spatiotemporal lag. This spatiotemporal lag helps us better understand the complex interactions between surface temperature and the climate system, and has broad application in fields such as climate prediction, urban planning, and agricultural management. Combining temporal and spatial lags in analysis can more comprehensively capture the changing characteristics of LST and improve its prediction accuracy.

[0051] Second, the present invention proposes a land surface temperature (LST) downscaling model suitable for mountainous cities - the temporal and spatial weighted lag random forest model (GTWLRF). During the model design process, the temporal and spatial lag, temporal and spatial heterogeneity of LST, and the complex nonlinear relationship between LST and scale factors are fully considered. Traditional LST downscaling methods are usually difficult to accurately reflect the dynamic change characteristics of LST in the temporal and spatial dimensions, especially in complex terrain such as mountainous cities. By introducing a temporal and spatial weighting mechanism, the GTWLRF model effectively captures the temporal and spatial nonlinear relationship between LST and scale factors, thereby improving the adaptability of the model at different time nodes and spatial locations. At the same time, the lag consideration in the model can utilize LST data at past moments to improve the prediction accuracy of the current LST. In addition, the random forest algorithm, as the core of the model, not only has powerful nonlinear fitting capabilities, but can also effectively reduce the noise impact brought by multi-dimensional feature input and enhance the robustness of the model. The advantage of this approach is that the GTWLRF model, taking advantage of the rugged terrain and high environmental diversity found in mountainous cities, can more accurately capture local trends in LST and achieve refined predictions of its spatial distribution. This approach has important applications in climate change research, urban thermal environment assessment, and ecological and environmental monitoring, and also provides a scientific basis for refined environmental management in mountainous cities. This model can better understand the spatiotemporal variations of LST in mountainous cities, thereby supporting urban planning, ecological and environmental protection, and sustainable development.

[0052] The embodiment of the present invention is: downscaling the surface temperature of a certain area (106°28′-106°39′ east longitude, 29°22′-29°40′ north latitude).

[0053] Figure 1 The figure shows a flow chart of the surface temperature fine-grained downscaling method proposed in this invention that combines temporal and spatial hysteresis. Specifically, the acquired 30m spatial resolution Landsat8 TIRS data for the corresponding area is subjected to radiation correction, atmospheric correction, and geometric correction, and NDVI, NDBI, and MNDWI are calculated. Finally, the atmospheric correction method is used to invert the Landsat8 TIRS data for LST.

[0054] The NDVI calculation formula is shown in formula (1):

[0055]

[0056] Where b nir Indicates the near-infrared band band5 of Landsat8 TIRS, b red Indicates the red band 4 of Landsat8 TIRS.

[0057] The calculation formula of NDBI is shown in formula (2):

[0058]

[0059] Where b nir Indicates the near-infrared band 5 of Landsat8 TIRS. Indicates the first shortwave infrared band band 6 of Landsat8 TIRS.

[0060] The calculation formula of MNDWI is shown in formula (3):

[0061]

[0062] Where b nir Indicates the near-infrared band band5 of Landsat8 TIRS, b green Indicates the green band 2 of Landsat8 TIRS.

[0063] The inversion process of the surface temperature of the 10th band of Landsat8 is as follows: Get the atmospheric radiation brightness B of the 10th band 10 (T 10 ), the calculation formula is shown in formula (4):

[0064] B 10 (T 10 )=M 10 ×Q al +A 10 (4)

[0065] Among them, Q al is the DN value of the pixel, M 10 is the gain value of the 10th band, A 10 is the bias value of the 10th band. Calculate the surface radiation brightness B according to the radiation transfer equation 10 (T s ), the calculation formula is shown in formula (5):

[0066]

[0067] In this formula, B 10 (T 10 ) represents the radiation brightness received by band 10, ∈ 10 is the surface emissivity in band 10, τ 10 is the atmospheric transmittance. represents the atmospheric downward radiance, represents the upward radiation brightness of the atmosphere. The surface temperature is calculated by inverse operation of the Planck function, and the calculation formula is shown in formula (6):

[0068]

[0069] Where K1 and K2 are constant parameters in the inverse operation of the Planck function

[0070] The MODIS LST is calculated using the 1000m spatial resolution MOD11A1 product data provided by MODIS. The temperature value of this product is Kelvin temperature, with a scaling factor of 0.02 and a range of 7500-65535K. The formula for converting MOD11A1 data to degrees Celsius is shown in Equation (7):

[0071] T(℃)=T(K)×0.02-273.15 (7)

[0072] Where T(℃) represents Celsius temperature, and T(K) represents the temperature value of MOD11A1 product.

[0073] The terrain correction and terrain factor data extraction are completed by acquiring the 90m spatial resolution DEM data, as shown in formula (8):

[0074] Shown: DEM = SRTM elevation ×SF (8)

[0075] Among them, DEM represents the extracted elevation value, SRTM elevation It represents the original value of SRTM data in GEE, and SF represents the scaling factor, which is usually set to 1.

[0076] The 30m spatial resolution Landsat8 NDVI, NDBI, and MNDWI images and the 90m spatial resolution SRTM DEM images were resampled to 100m and 1000m spatial resolutions to match the spatial resolutions of Landsat8 LST and MODIS LST images. The expression is shown in Equation (9):

[0077]

[0078] f(i, j) = NDVI, NDBI, MNDWI, DEM (9) where f(x, y) is the value of the interpolation point, (i, j) is the integer grid point where the interpolation point (x, y) is located, f(i+m, j+n) is the value of the surrounding known grid points, w(m, xi) and w(n, yj) are one-dimensional cubic convolution weight functions, defined as:

[0079]

[0080] Where d represents the distance between the interpolation point and the integer grid point, and m and n range from -1 to 2, indicating the 4x4 grid involved in the cubic convolution interpolation.

[0081] A time series dataset for the study area was constructed, chronologically arranging MODIS LST data from two consecutive time periods to serve as the dependent variable in the LST downscaling process. Aggregated Landsat8 NDVI, NDBI, and MNDWI data, 1000-meter spatial resolution SRTM DEM data, and land cover type data were also chronologically arranged to construct a time series dataset. This dataset was used to establish the nonlinear relationship between MODIS LST and the coarse spatial resolution Landsat8 NDVI, NDBI, MNDWI, and SRTM DEM data, as well as the spatiotemporal lag relationship between LST data from two consecutive time periods. Furthermore, aggregated Landsat8 NDVI, NDBI, and MNDWI data, along with 100-meter spatial resolution SRTM DEM data, were chronologically arranged to construct a time series dataset of the independent variables used for high-resolution LST downscaling.

[0082] The temporal and spatial weighted lag random forest model (GTWLRF) is a downscaling model that comprehensively considers the temporal and spatial lag of LST, the temporal and spatial heterogeneity, the nonlinear relationship between LST and scale factor, and the nonlinear change of LST in time and space, as shown in Equation (11).

[0083]

[0084] Where (u i , v i , τ i ) is the space-time position of point i, W(u i , v i , τ i ) it is the spatiotemporal adjacency matrix, X it is the independent variable matrix at point i at time t, y i(t-1) is the time lag term, is the spatial lag term, y it is the target variable at time t, ∑ p f p is a random forest regression model, p represents the number of decision trees, ε it Is the error term. Where the spatiotemporal weight matrix W(u i , v i , τ i )=diag(w i1 , w i2 ,…,w in ) is a square matrix. The diagonal elements w ij Represents the geographical and temporal weights of the target point i to the surrounding observation points j. The spatiotemporal weight matrix W(u i , v i , τ i ) can be achieved by using space-time distance and is constructed based on the Gaussian distance decay function. The spatiotemporal weight matrix W(u i , v i , τ i )'s diagonal elements w ij It can be obtained by the following formula:

[0085]

[0086] Where h ST is the spatiotemporal bandwidth, which controls the decay rate of the weight and determines the degree of influence of distance on the weight. The calculation of spatiotemporal distance includes the combination of spatial and temporal weights, as shown in Equation (13).

[0087]

[0088] In the formula is the spatial distance, is the time distance, λ and are weight coefficients that balance the spatial and temporal scales. λ can be set to 1, and μ can be adjusted based on the coefficient of determination (R 2 ) statistic and modified Akaike information criterion are determined and optimized by cross validation. According to formula (12) and formula (13), the spatiotemporal weight matrix W(u i , v i , τ i )'s diagonal elements w ij can be rewritten as:

[0089]

[0090]

[0091] where h S is the time bandwidth, h T This paper takes the spatial and temporal distance from the target point i to the kth neighboring point as the adaptive kernel bandwidth h ST ,Right now This ensures that the kernel bandwidth of each target point i is dynamically adjusted based on the distance to its neighborhood, and the weight calculation automatically adapts to different point densities.

[0092] Using the spatiotemporal series dataset of explanatory variables and MODIS LST, the MODIS LST and Landsat8 NDVI at the current moment are established based on the GTWLRF model at a coarse spatial resolution. 1000 ,NDBI 1000 ,MNDWl 1000 and SRTM DEM 1000The spatiotemporal heterogeneity relationship of , the spatiotemporal lag relationship with the previous moment MODIS LST, and the spatiotemporal nonlinear relationship with all explanatory variables are expressed as shown in formula (15):

[0093]

[0094] In the formula represents the MODIS LST with a spatial resolution of 1000m at point i at time t, represents the MODIS LST with a spatial resolution of 1000m at point i at time t-1, represents the MODIS LST of k points adjacent to point i at time t-1 with a spatial resolution of 1000m, W represents the set of Landsat8 NDVI, NDBI, MNDWI and SRTM DEM with a spatial resolution of 1000m at point i at time t, respectively. 1000 (u i , v i , τ i ) represents the spatiotemporal weight matrix of point i with a spatial resolution of 1000m, represents the regression residual at point i at time t with a spatial resolution of 1000m.

[0095] The LST was downscaled and the regression residuals at 1000m spatial resolution were interpolated using ordinary kriging. Interpolate to a spatial resolution of 100 m. Then, assuming that the spatiotemporal lag relationship, spatiotemporal heterogeneity, and spatiotemporal nonlinear relationship between LST and explanatory variables are scale-invariant, the regression relationship between the current spatiotemporal LST and the current spatiotemporal NDVI, NDBI, MNDWI, DEM, and the previous spatiotemporal LST at a spatial resolution of 100 m after downscaling can be expressed as shown in Equation (16):

[0096]

[0097] In the formula represents the Landsat8LST with a spatial resolution of 100m at point i at time t-1, represents the k points adjacent to point i at time t-1, with a Landsat8 LST spatial resolution of 100m. W represents the set of Landsat8 NDVI, NDBI, MNDWI and SRTM DEM with a spatial resolution of 100m at point i at time t, respectively. 100 (u i ,v i ,τ i) represents the spatiotemporal weight matrix of point i with a spatial resolution of 100m, represents the regression residual of point i at time t from 1000m spatial resolution to 100m through ordinary kriging interpolation. It is the downscaled LST result of point i at time t100m spatial resolution.

[0098] Figure 2 Shown is a schematic diagram of the spatiotemporal hysteresis effect. The principle is that the present invention enhances spatiotemporal constraints by introducing spatiotemporal hysteresis, thereby improving the prediction accuracy of LST. Specifically, the LST prediction at the current moment depends on the LST information at the previous moment, which includes the LST of the pixel at the previous moment and the LST of its surrounding pixels, providing time lag information and spatial lag information, respectively. Because the spatiotemporal variation characteristics of LST are nonlinear, spatiotemporal lagged temperature data can be input into the random forest model along with other predictive factors for training. This method not only focuses on the changes in the current pixel, but also considers its mutual relationship with surrounding pixels, thereby enhancing spatiotemporal constraints. It also effectively utilizes historical data to capture the spatiotemporal variation patterns of temperature, helping the model better learn the pattern of temperature changes over time and space, thereby improving the accuracy and stability of the prediction.

[0099] Figure 3 Figure 2 shows a schematic diagram of the downscaling effect of the method of the present invention. By considering the impact of the spatiotemporal hysteresis of LST, the method of the present invention more comprehensively incorporates spatiotemporal constraints. It also comprehensively considers the spatiotemporal hysteresis, spatiotemporal heterogeneity, spatiotemporal nonstationarity of LST, the nonlinear relationships between predictors, and the nonlinear variations of LST in both temporal and spatial dimensions. This comprehensive spatiotemporal constraint and thorough consideration of nonlinear variations enable the method of the present invention to more accurately capture the complex spatiotemporal characteristics of LST while reducing noise during the downscaling process, and to more comprehensively restore LST spatial details.

[0100] Figure 4Shown is the absolute error map of the downscaling method of the present invention. A significant advantage of the method of the present invention is that it effectively reduces the occurrence of local high-error areas and makes the error distribution more uniform throughout the study area. In traditional LST downscaling methods, due to insufficient capture of the temporal and spatial variation laws of LST, large errors are easily generated in certain areas with complex terrain or significant environmental characteristics. These error areas not only affect the accuracy of the overall downscaling results, but also weaken the reliability of the model in practical applications. By introducing the method of the present invention, this problem is effectively overcome. This uniform error distribution characteristic not only improves the adaptability of the model under various terrains, but also enhances the stability and reliability of the LST downscaling results, enabling it to achieve more accurate surface temperature simulation and prediction in complex environments such as mountainous cities. This is of great significance for improving the scientific nature of regional environmental monitoring, urban planning and ecological protection, and can provide more accurate temperature data support for practical applications.

[0101] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not limiting. Although the present invention has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that the technical solutions of the present invention can be modified or replaced by equivalents without departing from the purpose and scope of the technical solutions, which should all be included in the scope of the claims of the present invention.

Claims

1. A refined downscaling method for surface temperature incorporating temporal and spatial hysteresis, characterized by: The following steps are involved: The acquired 30m spatial resolution Landsat8 TIRS data of the corresponding area were subjected to radiometric, atmospheric, and geometric corrections, and NDVI, NDBI, and MNDWI were calculated. The atmospheric correction method was used to invert the LST of the Landsat8 TIRS data. The MODIS LST is calculated using the 1000m spatial resolution MOD11A1 product data provided by MODIS; The terrain correction and terrain factor data extraction were completed by using the acquired 90m spatial resolution DEM data; The Landsat8 NDVI, NDBI, and MNDWI images with a spatial resolution of 30 m and the SRTM DEM images with a spatial resolution of 90 m were resampled to 100 m and 1000 m spatial resolutions by using the cubic convolution interpolation method. A time series dataset of the study area was constructed, and the MODIS LST of two consecutive time periods were arranged in chronological order and used as the dependent variable in the LST downscaling process; The aggregated Landsat8 NDVI, NDBI, and MNDWI data, 1000-meter spatial resolution SRTM DEM data, and land cover type data were arranged in chronological order to construct a time series dataset. This dataset was used to establish the nonlinear relationship between MODIS LST and the coarse spatial resolution Landsat8 NDVI, NDBI, MNDWI data and SRTM DEM data, as well as the spatiotemporal lag relationship between LST in two consecutive time periods. The aggregated Landsat8 NDVI, NDBI, MNDWI data and 100-meter spatial resolution SRTM DEM data were arranged in chronological order to construct a time series dataset of independent variables for high-resolution LST downscaling. Using the spatiotemporal series dataset of explanatory variables and MODIS LST, the temporal and spatial heterogeneity relationship between the current MODIS LST and Landsat8 NDVI, NDBI, MNDWI and SRTM DEM, the spatiotemporal lag relationship with the previous MODIS LST, and the spatiotemporal nonlinear relationship with all explanatory variables were established based on the GTWLRF model at a coarse spatial resolution. The LST was downscaled and the regression residuals at 1000m spatial resolution were interpolated to 100m spatial resolution using ordinary kriging interpolation; Assuming that the spatiotemporal lag, spatiotemporal heterogeneity, and spatiotemporal nonlinear relationship between LST and explanatory variables are scale invariant, the regression relationship between the current spatiotemporal LST and the current spatiotemporal NDVI, NDBI, MNDWI, DEM, and the previous spatiotemporal LST at a spatial resolution of 100 m after downscaling is expressed as follows: In the formula represents the Landsat8LST with a spatial resolution of 100m at point i at time t-1, represents the k points adjacent to point i at time t-1, with a Landsat8 LST spatial resolution of 100m. W represents the set of Landsat8 NDVI, NDBI, MNDWI and SRTM DEM with a spatial resolution of 100m at point i at time t, respectively. 100 (u i ,v i ,τ i ) represents the spatiotemporal weight matrix of point i with a spatial resolution of 100m, represents the regression residual of point i at time t from 1000m spatial resolution to 100m through ordinary kriging interpolation. It is the downscaled LST result of point i at time t100m spatial resolution.

2. The surface temperature refined downscaling method incorporating spatiotemporal hysteresis according to claim 1 is characterized by: The GTWLRF model is shown as follows: Where (u i ,v i ,τ i ) is the space-time position of point i, W(u i ,v i ,τ i ) it is the spatiotemporal adjacency matrix, X it is the independent variable matrix at point i at time t, y i(t-1) is the time lag term, is the spatial lag term, y it is the target variable at time t, ∑ p f p is a random forest regression model, p represents the number of decision trees, ε it is the error term; where the spatiotemporal weight matrix W(u i ,v i ,τ i )=diag(w i1 ,w i2 ,…,w in ) is a square matrix; the diagonal elements w ij Represents the geographical and temporal weights of target point i to surrounding observation points j.

3. The surface temperature refined downscaling method incorporating spatiotemporal hysteresis according to claim 2 is characterized by: The spatiotemporal weight matrix W(u i ,v i ,τ i ) by using space-time distance and is constructed based on the Gaussian distance decay function; the spatiotemporal weight matrix W(u i ,v i ,τ i )'s diagonal elements w ij Obtained by the following formula: Where h ST It is the space-time bandwidth, which is used to control the decay rate of the weight and determine the degree of influence of distance on the weight. The calculation of space-time distance includes the weight combination of space and time, as shown in the following formula: In the formula is the spatial distance, is the time distance, λ and are weight coefficients for balancing the spatial and temporal scales; λ is set to 1, and μ is determined by the coefficient of determination R 2 Statistics and modified Akaike information criterion are determined and optimized by cross validation; spatiotemporal weight matrix W(u i ,v i ,τ i )'s diagonal elements w ij Rewritten as: where h S is the time bandwidth, h T is the spatial bandwidth; Take the spatiotemporal distance from target point i to the kth neighboring point as the adaptive kernel bandwidth h ST ,Right now Ensure that the kernel bandwidth of each target point i is dynamically adjusted based on the distance of its neighborhood, and the weight calculation automatically adapts to different point densities.

4. The surface temperature refined downscaling method incorporating spatiotemporal hysteresis according to claim 1 is characterized by: The LST downscaling process is as follows: In the formula represents the MODIS LST with a spatial resolution of 1000m at point i at time t, represents the MODIS LST with a spatial resolution of 1000m at point i at time t-1, represents the MODIS LST of k points adjacent to point i at time t-1 with a spatial resolution of 1000m, W represents the set of Landsat8 NDVI, NDBI, MNDWI and SRTM DEM with a spatial resolution of 1000m at point i at time t, respectively. 1000 (u i ,v i ,τ i ) represents the spatiotemporal weight matrix of point i with a spatial resolution of 1000m, represents the regression residual at point i at time t with a spatial resolution of 1000m.

5. The surface temperature refined downscaling method incorporating spatiotemporal hysteresis according to claim 1 is characterized by: The LST downscaling process is as follows: In the formula represents the Landsat8 LST with a spatial resolution of 100m at point i at time t-1, represents the k points adjacent to point i at time t-1, with a Landsat8 LST spatial resolution of 100m. W represents the set of Landsat8 NDVI, NDBI, MNDWI and SRTM DEM with a spatial resolution of 100m at point i at time t, respectively. 100 (u i ,v i ,τ i ) represents the spatiotemporal weight matrix of point i with a spatial resolution of 100m, represents the regression residual of point i at time t from 1000m spatial resolution to 100m through ordinary kriging interpolation. It is the downscaled LST result of point i at time t100m spatial resolution.

Images