Surface temperature refined downscaling method combined with space-time hysteresis

Abstract

The invention relates to a land surface temperature refined downscaling method combined with space-time hysteresis, and belongs to the technical field of remote sensing image processing. The method comprises the following steps: firstly, carrying out preprocessing and LST inversion on Landsat8TIRS data, and calculating MODIS LST and DEM data; then, constructing a time sequence data set, and establishing a non-linear relationship between MODIS LST and Landsat8NDVI, NDBI, MNDWI and SRTM DEM data and a time-space lag relationship between two continuous time periods of LST based on a time-space weighted lag random forest model; and finally, downscaling the LST by using a common Kriging interpolation method to obtain high-resolution surface temperature distribution. According to the method, the precision and reliability of the LST downscaling result are effectively improved, and a scientific basis is provided for climate regulation, ecological protection, thermal environment management and the like of mountainous cities.

Description

Technical Field

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

[0002] Land surface temperature (LST) is an important physical parameter that describes the surface-atmosphere interaction process and reflects the surface. It is the result of the interaction between the surface and the atmosphere and the energy exchange between the atmosphere and the land. As a key parameter in many basic disciplines and application fields, LST is widely used in research fields 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 characteristics of LST are extremely complex and changeable issues. They are affected by a variety of complex factors. They are the result of the complexity of atmospheric movement, the direct impact of human activities, and the interweaving and joint action of multiple spatiotemporal characteristics such as spatiotemporal hysteresis and spatiotemporal heterogeneity. Spatiotemporal hysteresis means that there is a certain time delay or spatial dislocation between a geographical phenomenon or process and the moment before or after it. This is particularly evident in the distribution of LST in mountainous cities. For example, the impact of land use changes in the urbanization process on surface temperature often has a hysteresis effect. At the same time, the temporal and spatial heterogeneity of LST in mountainous cities is more prominent, which is due to the unique topographic and climatic conditions of mountainous cities, such as highly differentiated topographic changes, local climate effects, and differences in land use during urbanization. These factors lead to a highly uneven distribution of LST in space. For example, the temperature difference between the urban core and the surrounding rural areas may be quite obvious, reflecting the differences in the heat island effect and vegetation coverage in different regions. In addition, there is often a more obvious nonlinear relationship between scale factors and LST in mountainous cities, which means that the impact of scale factors on LST is not a simple linear superposition, but there are complex interactions and feedback mechanisms. For example, the impact of vegetation coverage on surface temperature shows obvious nonlinear characteristics. Therefore, when studying LST in mountainous cities, a comprehensive analysis of the above complex factors is conducted to predict LST in mountainous cities with fine resolution. This will have important practical application value in climate regulation, ecological protection, response to heat waves, and thermal environment management in mountainous cities.

[0003] According to different principles, LST downscaling methods can be mainly divided into models based on image fusion and statistical models based on relational scale invariance. The image fusion-based method integrates multi-source data according to similar pixel weighting or mixed pixel decomposition to generate time-intensive high-spatial resolution images. The spatiotemporal fusion model method was originally used for reflectance data, such as flexible spatiotemporal data fusion method and improved spatiotemporal data fusion method, but they can still be used to enhance LST. The LST spatiotemporal fusion model based on the improved U-STFM enhances the spatial resolution of MODIS LST. The unbiased ESTARFM spatiotemporal fusion model emphasizes the potential of using LST data obtained from geostationary satellites (such as Himawari-8) for applications. The spatiotemporal image fusion method has been widely used because it requires fewer auxiliary parameters. However, the fusion method relies on the quality of high-resolution LST and ignores the clear physical background in satellite remote sensing data, which makes it difficult to meet the quantitative requirements of thermal infrared remote sensing.

[0004] The LST downscaling model based on the invariant relationship scale has a clear physical meaning, maintains the consistency of thermal radiation information before and after downscaling, and is widely used. Its basic idea is to establish a statistical regression relationship between LST and scale factors such as vegetation index, building index, and terrain factor at low spatial resolution, and apply the statistical relationship to high-resolution scale factors to achieve LST downscaling. Classical downscaling models include the DisTrad model and the TsHARP model, which use the normalized difference vegetation index (NDVI) and vegetation coverage as scale factors, respectively, and obtain good LST downscaling results. The above downscaling methods are not suitable for cities and arid areas due to the limitation of a single variable and the simplification of the 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 is established between multiple spectral indices and LST to obtain better downscaling effects. The LST downscaling method weighted by kernel function and fusion method can achieve daily LST data with high spatial resolution. The kernel-driven approach is integrated into the fusion method, and a simple but effective downscaling (SED) algorithm is proposed to obtain LST with high spatiotemporal resolution.

[0005] Based on the downscaling model with invariant local relationship scale, the mapping relationship between LST and scale factor varies with local location. Duan et al. considered spatial heterogeneity and applied the geographically weighted regression (GWR) model to the research field of LST downscaling for the first time. The model studied the local relationship between LST and vegetation index and DEM. The GWR model was combined with Kriging interpolation to generate geographically weighted Kriging model and multi-scale geographically weighted Kriging model. Some scholars comprehensively considered the spatial autocorrelation of LST, spatiotemporal non-stationarity and other spatiotemporal constraints, as well as local nonlinear mapping relationships, and constructed downscaling models based on spatiotemporal geographic weighting, geographically weighted autoregression, spatiotemporal geographic weighted autoregression and nonlinear geographical weighting, achieving good performance. Wu et al. proposed a weighted autoregression model of geographic and temporal neural network to improve the overall accuracy of downscaling. Mu et al. proposed a geographically weighted random forest downscaling model by coupling the geographically weighted model with the random forest algorithm. On the basis of the geographically weighted regression model, the random forest regression model was introduced to express the non-stationarity and nonlinearity of surface temperature and scale factor.

[0006] In general, the surface temperature downscaling model has achieved a lot of valuable research results, but it still faces many challenges, which are mainly reflected in the following aspects: First, the unified consideration of spatiotemporal constraints is insufficient. Existing downscaling methods often fail to consider the spatiotemporal lag of LST. Therefore, many downscaling models fail to effectively correct the temporal variation characteristics of LST and find it difficult to accurately capture its dynamic changes in different time periods, which leads to deviations in the results in practical applications. Secondly, for mountainous urban LST, the complex nonlinear relationship between it and the scale factor has not been fully considered. The existing models still do not adequately handle this nonlinear relationship, which limits their adaptability in complex environments. Finally, the changes in LST in spatial and temporal dimensions show significant nonlinear characteristics. How to accurately capture these changing laws in the downscaling process is still a difficult problem that needs to be solved urgently. Summary of the invention

[0007] In view of this, the purpose of the present invention is to provide a refined downscaling method for surface temperature combined with spatiotemporal hysteresis. The method incorporates the influence of spatiotemporal hysteresis on LST, comprehensively considers the spatiotemporal hysteresis of LST, spatiotemporal heterogeneity, the nonlinear relationship between LST and scale factor, and the nonlinear change of LST in time and space dimensions, providing a new way to downscale the complex and changeable LST of mountainous cities 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 refined downscaling method for surface temperature combined with temporal and spatial hysteresis includes the following steps:

[0010] Perform radiation correction, atmospheric correction and geometric correction on the acquired 30m spatial resolution Landsat8TIRS data of the corresponding area, calculate NDVI, NDBI, MNDWI, and use the atmospheric correction method to perform LST inversion on 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 acquiring the 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 dependent variables in the LST downscaling process;

[0015] The aggregated Landsat8 NDVI, NDBI, 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, which 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-m spatial resolution SRTM DEM data were arranged in chronological order to construct an independent variable time series dataset for high-resolution LST downscaling;

[0017] Using the spatiotemporal series dataset of explanatory variables and MODIS LST, the spatiotemporal heterogeneous 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 are 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 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 is expressed as:

[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 spatial resolution of 100m Landsat8 LST, W represents the collection of Landsat8 NDVI, NDBI, MNDWI and SRTM DEM with a spatial resolution of 100 m at point i at time t. 100 (u i , v i , τ i ) represents the spatiotemporal weight matrix of point i with a spatial resolution of 100m, It 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 the spatial resolution of 100m at point i at time t.

[0022] Further, 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 speed of the weight and determine the 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 according to the coefficient of determination R 2 Statistics and modified Akaike information criterion are determined and optimized by cross-validation; the spatiotemporal weight matrix W(u i , v i , τ i )’s diagonal elements w ij Rewritten as:

[0030]

[0031] where h S is the temporal bandwidth, h T is the spatial bandwidth; the spatiotemporal distance from the target point i to the kth neighboring point is taken 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 in the following formula:

[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 collection of Landsat8 NDVI, NDBI, MNDWI and SRTM DEM with a spatial resolution of 1000m at point i at time t. 1000 (u i , v i , τ i ) represents the spatiotemporal weight matrix of point i with a spatial resolution of 1000m, represents the regression residual with a spatial resolution of 1000m at point i at time t.

[0035] Furthermore, the LST downscaling process is shown in the following formula:

[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 spatial resolution of 100m Landsat8 LST, W represents the collection of Landsat8 NDVI, NDBI, MNDWI and SRTM DEM with a spatial resolution of 100 m at point i at time t. 100 (u i ,v i ,τ i ) represents the spatiotemporal weight matrix of point i with a spatial resolution of 100m, It 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 the spatial resolution of 100m at point i at time t.

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

[0039] First, the present invention proposes a method for fine-grained downscaling of surface temperature combined with spatiotemporal hysteresis, which innovatively incorporates the impact of spatiotemporal hysteresis on LST. The spatiotemporal hysteresis effect of LST refers to the delayed response of LST to various influencing factors in time and space. This effect is of great significance in ecological environment and climate research. First, LST will lag behind the change of solar radiation in time. If the LST at the previous time point is high, the heat stored in the surface will also be relatively more, which may cause the LST at the current time point to be relatively high, even if the external conditions such as solar radiation have not changed significantly. For example, in a day, although the peak of solar radiation occurs at noon, the peak of surface temperature usually occurs after noon because the surface needs time to absorb and conduct heat. In addition, there are similar delays in seasonal changes. For example, after the increase of solar radiation in spring, the surface temperature will increase with hysteresis due to the influence of winter snow and soil moisture. In space, the conduction of surface heat causes the temperature change to gradually spread. For example, the urban heat island effect will be conducted to the surrounding areas, resulting in a spatial delay in temperature change. Topography also affects this effect. Valleys and low-lying areas will cool down more slowly at night due to the accumulation of cold air, while high-altitude areas will cool down relatively quickly. Different land cover types respond to temperature changes at different speeds. For example, vegetation-covered areas will delay temperature changes through evaporation, causing their temperature changes to lag behind bare land. In addition, atmospheric conditions such as cloud cover and precipitation interact with the physical properties of the surface, further exacerbating the complexity of the LST spatiotemporal lag effect. This spatiotemporal lag effect helps us better understand the complex interaction between surface temperature and the climate system, and has a wide range of application value in climate prediction, urban planning, and agricultural management. Combining time and space lags for analysis can more comprehensively capture the changing characteristics of LST and improve the prediction accuracy of LST.

[0040] Second, the present invention proposes a land surface temperature (LST) downscaling model suitable for mountainous cities, namely the spatiotemporal weighted lag random forest model (GTWLRF). In the process of model design, the spatiotemporal hysteresis, spatiotemporal heterogeneity, 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 time and space dimensions, especially in complex terrains such as mountainous cities. The GTWLRF model effectively captures the spatiotemporal nonlinear relationship between LST and scale factors by introducing a spatiotemporal weighting mechanism, thereby improving the adaptability of the model at different time nodes and spatial locations. At the same time, the hysteresis considerations in the model can utilize LST data at past moments to improve the prediction accuracy of the current LST. In addition, as the core of the model, the random forest algorithm not only has a powerful nonlinear fitting ability, but also can effectively reduce the noise impact brought by multi-dimensional feature input and enhance the robustness of the model. The advantage of this method is that the GTWLRF model can more accurately capture the local change trend of LST and realize the refined prediction of LST spatial distribution in view of the characteristics of large terrain fluctuations and high environmental diversity in mountainous cities. This has important application value in the fields of climate change research, urban thermal environment assessment, ecological environment monitoring, etc., and also provides a scientific basis for the refined environmental management of mountainous cities. Through this model, we can better understand the temporal and spatial variation of LST in mountainous cities, thereby providing support for urban planning, ecological environment protection and sustainable development.

[0041] Other advantages, objectives and features of the present invention will be described in the following description to some extent, and to some extent, will be obvious to those skilled in the art based on the following examination and study, or can be taught from the practice of the present invention. The objectives and other advantages of the present invention can 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 in conjunction with the accompanying drawings, wherein:

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

[0044] Figure 2 This is a schematic diagram of the space-time lag effect;

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

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

[0047] The following describes the embodiments of the present invention by 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 only illustrate the basic concept of the present invention in a schematic manner, and the following embodiments and features in the embodiments can be combined with each other without conflict.

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

[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 the terms "upper", "lower", "left", "right", "front", "rear", etc. indicate the orientation or position relationship, they are based on the orientation or position relationship shown in the drawings, which is only for the convenience of describing the present invention and simplifying the description, rather than indicating or implying that the device or element referred to must have a specific orientation, be constructed and operate in a specific orientation. Therefore, the terms describing the position 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, the present invention proposes a method for fine-grained downscaling of surface temperature combined with spatiotemporal hysteresis, which innovatively incorporates the impact of spatiotemporal hysteresis on LST. The spatiotemporal hysteresis effect of LST refers to the delayed response of LST to various influencing factors in time and space. This effect is of great significance in ecological environment and climate research. First, LST will lag behind the change of solar radiation in time. If the LST at the previous time point is high, the heat stored in the surface will also be relatively more, which may cause the LST at the current time point to be relatively high, even if the external conditions such as solar radiation have not changed significantly. For example, in a day, although the peak of solar radiation occurs at noon, the peak of surface temperature usually occurs after noon because the surface needs time to absorb and conduct heat. In addition, there are similar delays in seasonal changes. For example, after the increase of solar radiation in spring, the surface temperature will increase with hysteresis due to the influence of winter snow and soil moisture. In space, the conduction of surface heat causes the temperature change to gradually spread. For example, the urban heat island effect will be conducted to the surrounding areas, resulting in a spatial delay in temperature change. Topography also affects this effect. Valleys and low-lying areas will cool down more slowly at night due to the accumulation of cold air, while high-altitude areas will cool down relatively quickly. Different land cover types respond to temperature changes at different speeds. For example, vegetation-covered areas will delay temperature changes through evaporation, causing their temperature changes to lag behind bare land. In addition, atmospheric conditions such as cloud cover and precipitation interact with the physical properties of the surface, further exacerbating the complexity of the LST spatiotemporal lag effect. This spatiotemporal lag effect helps us better understand the complex interaction between surface temperature and the climate system, and has a wide range of application value in climate prediction, urban planning, and agricultural management. Combining time and space lags for analysis can more comprehensively capture the changing characteristics of LST and improve the prediction accuracy of LST.

[0051] Second, the present invention proposes a land surface temperature (LST) downscaling model suitable for mountainous cities, namely the spatiotemporal weighted lag random forest model (GTWLRF). In the process of model design, the spatiotemporal hysteresis, spatiotemporal heterogeneity, 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 time and space dimensions, especially in complex terrains such as mountainous cities. The GTWLRF model effectively captures the spatiotemporal nonlinear relationship between LST and scale factors by introducing a spatiotemporal weighting mechanism, thereby improving the adaptability of the model at different time nodes and spatial locations. At the same time, the hysteresis considerations in the model can utilize LST data at past moments to improve the prediction accuracy of the current LST. In addition, as the core of the model, the random forest algorithm not only has a powerful nonlinear fitting ability, but also can effectively reduce the noise impact brought by multi-dimensional feature input and enhance the robustness of the model. The advantage of this method is that the GTWLRF model can more accurately capture the local change trend of LST and realize the refined prediction of LST spatial distribution in view of the characteristics of large terrain fluctuations and high environmental diversity in mountainous cities. This has important application value in the fields of climate change research, urban thermal environment assessment, ecological environment monitoring, etc., and also provides a scientific basis for the refined environmental management of mountainous cities. Through this model, we can better understand the temporal and spatial variation of LST in mountainous cities, thereby providing support for urban planning, ecological environment protection and sustainable development.

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

[0053] Figure 1 What is shown is a flow chart of the refined downscaling method of surface temperature combined with time and space hysteresis proposed in the present invention, specifically: the acquired 30m spatial resolution Landsat8 TIRS data of the corresponding area is subjected to radiation correction, atmospheric correction and geometric correction, and NDVI, NDBI, MNDWI are calculated, and the atmospheric correction method is used to perform LST inversion on the Landsat8 TIRS data.

[0054] The calculation formula of NDVI 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 band6 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 band2 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 obtained by inverse calculation 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 Planck function

[0070] MODIS LST is calculated using the 1000m spatial resolution MOD11A1 product data provided by MODIS. The temperature value of this product is Kelvin temperature, the scaling ratio is 0.02, and the value range is 7500-65535K. The formula for converting MOD11A1 data to 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 Landsat8 NDVI, NDBI, MNDWI images with a spatial resolution of 30 m and the SRTM DEM images with a spatial resolution of 90 m are resampled to 100 m and 1000 m spatial resolutions to match the spatial resolutions of Landsat8 LST and MODIS LST images. The expression is shown in formula (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 of the study area was constructed, and the MODIS LST of two consecutive time periods was arranged in chronological order and used as the dependent variable in the LST downscaling process. At the same time, the aggregated Landsat8 NDVI, NDBI, MNDWI data and 1000-meter spatial resolution SRTM DEM data and land cover type data were arranged in chronological order to construct a time series dataset for establishing the nonlinear relationship between MODIS LST and the coarse spatial resolution Landsat8 NDVI, NDBI, MNDWI data and SRTM DEM data and the spatiotemporal lag relationship between LST in two consecutive time periods. In addition, the aggregated Landsat8 NDVI, NDBI, MNDWI data and 100-meter spatial resolution SRTM DEM data were arranged in chronological order to construct an independent variable time series dataset for high-resolution LST downscaling.

[0082] The spatiotemporal weighted lag random forest model (GTWLRF) is a downscaling model that comprehensively considers the spatiotemporal lag of LST, spatiotemporal 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. 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 geographic 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 the 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 is used to control the decay speed of the weight and determine the influence of distance on the weight. The calculation of spatiotemporal distance includes the weight combination of space and time, as shown in formula (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 set based on the determination coefficient (R 2 ) statistic and modified Akaike information criterion are determined and optimized by cross validation. According to equations (12) and (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 temporal bandwidth, h T is the spatial bandwidth. In this paper, the temporal and spatial distance from the target point i to the kth neighboring point is taken 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 of 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 spatial and temporal heterogeneity relationship of the MODIS LST at the previous moment, the spatial and temporal lag relationship with the MODIS LST at the previous moment, and the spatial and temporal 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 collection of Landsat8 NDVI, NDBI, MNDWI and SRTM DEM with a spatial resolution of 1000m at point i at time t. 1000 (u i , v i , τ i ) represents the spatiotemporal weight matrix of point i with a spatial resolution of 1000m, represents the regression residual with a spatial resolution of 1000m at point i at time t.

[0095] The LST was downscaled and the regression residuals at a spatial resolution of 1000 m were interpolated using ordinary kriging. Interpolate to 100m spatial resolution. 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 after downscaling at 100m spatial resolution can be expressed as shown in formula (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 spatial resolution of 100m Landsat8 LST, W represents the collection of Landsat8 NDVI, NDBI, MNDWI and SRTM DEM with a spatial resolution of 100 m at point i at time t. 100 (u i ,v i ,τ i) represents the spatiotemporal weight matrix of point i with a spatial resolution of 100m, It 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 the spatial resolution of 100m at point i at time t.

[0098] Figure 2 What is shown is a schematic diagram of the spatiotemporal hysteresis effect. The principle is: the present invention enhances the 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 adjacent pixels, which respectively provide time lag information and space lag information. Since the spatiotemporal variation characteristics of LST are nonlinear, the spatiotemporal lagged temperature data can be input into the random forest model together with other predictors for training. This method not only pays attention to the changes in the current pixel, but also considers the relationship between it and the surrounding pixels, thereby enhancing the spatiotemporal constraints. And effectively use historical data to capture the spatiotemporal variation of temperature, help the model better learn the pattern of temperature variation over time and space, and thus improve the accuracy and stability of the prediction.

[0099] Figure 3 What is shown is a schematic diagram of the downscaling effect of the method of the present invention. Since the method of the present invention takes into account the influence of the spatiotemporal lag of LST, the spatiotemporal constraints are more comprehensively integrated. And the spatiotemporal lag, spatiotemporal heterogeneity, spatiotemporal non-stationarity of LST, the nonlinear relationship between the prediction factors, and the nonlinear changes of LST in time and space dimensions are comprehensively considered. This comprehensive spatiotemporal constraint and full consideration of nonlinear changes 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 restore the spatial details of LST more comprehensively.

[0100] Figure 4What is shown 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 distribution of errors more uniform in the entire study area. In the traditional LST downscaling method, due to the insufficient capture of the temporal and spatial variation laws of LST, large errors are easily generated in some 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 solution of the present invention rather than to limit it. 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 solution of the present invention can be modified or replaced by equivalents without departing from the purpose and scope of the technical solution, which should be included in the scope of the claims of the present invention.

Claims

1. A refined downscaling method for surface temperature combined with temporal and spatial hysteresis, characterized by: The following steps are involved: Perform radiation correction, atmospheric correction and geometric correction on the acquired 30m spatial resolution Landsat8 TIRS data of the corresponding area, calculate NDVI, NDBI, MNDWI, and use the atmospheric correction method to perform LST inversion on 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 acquiring the 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 dependent variables in the LST downscaling process; The aggregated Landsat8 NDVI, NDBI, 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 for establishing 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-m spatial resolution SRTM DEM data were arranged in chronological order to construct an independent variable time series dataset for high-resolution LST downscaling; Using the spatiotemporal series dataset of explanatory variables and MODIS LST, the spatiotemporal heterogeneous 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 are 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 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 is expressed as: 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 spatial resolution of 100m Landsat8 LST, W represents the collection of Landsat8 NDVI, NDBI, MNDWI and SRTM DEM with a spatial resolution of 100 m at point i at time t. 100 (u i ,v i ,τ i ) represents the spatiotemporal weight matrix of point i with a spatial resolution of 100m, It 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 the spatial resolution of 100m at point i at time t.

2. The method for fine-grained downscaling of land surface temperature in combination with 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 method for fine-grained downscaling of land surface temperature in combination with 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 speed of the weight and determine the 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 according to the coefficient of determination R 2 Statistics and modified Akaike information criterion are determined and optimized by cross-validation; the spatiotemporal weight matrix W(u i ,v i ,τ i )’s diagonal elements w ij Rewritten as: where h S is the temporal 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.

4. The method for fine-grained downscaling of land surface temperature in combination with 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 collection of Landsat8 NDVI, NDBI, MNDWI and SRTM DEM with a spatial resolution of 1000m at point i at time t. 1000 (u i ,v i ,τ i ) represents the spatiotemporal weight matrix of point i with a spatial resolution of 1000m, represents the regression residual with a spatial resolution of 1000m at point i at time t.

5. The method for fine-grained downscaling of land surface temperature in combination with 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 spatial resolution of 100m Landsat8 LST, W represents the collection of Landsat8 NDVI, NDBI, MNDWI and SRTM DEM with a spatial resolution of 100 m at point i at time t. 100 (u i ,v i ,τ i ) represents the spatiotemporal weight matrix of point i with a spatial resolution of 100m, It 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 the spatial resolution of 100m at point i at time t.

Images