A three-stage high-precision near-surface air temperature remote sensing estimation method

Abstract

The present invention discloses a three-stage high-precision near-surface temperature remote sensing estimation method, comprising: performing temporal normalization on remotely sensed surface temperature data to eliminate imaging time differences; clustering a study area into multiple sub-areas based on differences in natural conditions; and using four machine learning models to estimate the temperature in each sub-area based on the temporally normalized surface temperature data and spatial auxiliary variables to obtain spatially determined temperature estimates; and integrating the temperature estimation results of the four machine learning models using a generalized additive model to obtain a more accurate near-surface temperature. This method not only effectively eliminates the uncertainty caused by temporal differences in surface temperature but also considers the differences in model adaptability in different geographical environments. By integrating multiple single machine learning models based on an integrated approach, the accuracy of remote sensing temperature estimation is improved.

Description

Technical Field

[0001] The present invention belongs to the technical field of thermal infrared remote sensing, and in particular relates to a three-stage high-precision near-surface air temperature remote sensing estimation method. Background Art

[0002] Near-surface air temperature is a key meteorological parameter reflecting the thermal environment, playing an important role in many fields, including climate change monitoring, urban heat island effect research, public health, and agricultural production. High-temporal and spatial resolution near-surface air temperature data are crucial for better understanding and modeling complex land-atmosphere processes. While traditional ground-based observations can provide accurate near-surface air temperature data at a point scale, they suffer from limitations such as limited availability and uneven spatial distribution, making them inadequate for describing detailed temperature variations over large areas.

[0003] The development of satellite remote sensing has provided a new means for thermal environment monitoring. Satellite thermal infrared remote sensing, with its advantages of large-scale, long-term continuous observation, has been widely used in near-surface air temperature estimation. Four methods for estimating near-surface air temperature using remote sensing data can be categorized: the temperature-vegetation index method, the energy balance method, statistical methods, and machine learning. The temperature-vegetation index method estimates air temperature based on the assumption that the air temperature within the vegetation canopy is equal to the canopy surface temperature. This method requires few input parameters, all of which can be obtained through remote sensing, but is not suitable for areas with sparse vegetation. The energy balance method uses the surface energy balance equation to estimate the relationship between surface and air temperature and then calculates air temperature from remotely sensed surface temperatures. This method has clear physical principles and is highly portable and versatile. However, some variables in the surface energy balance equation, such as dynamic impedance and surface roughness, cannot be obtained through remote sensing, limiting its large-scale application. Statistical methods estimate air temperature by establishing regression equations between air temperature and environmental factors such as surface temperature. This method is easy to implement. For example, Patent Publication No. CN113255148A discloses a method for estimating all-weather temperature and its spatiotemporal distribution based on MODIS product data, and Patent Publication No. CN112016052B discloses a method, system, and terminal for estimating near-surface daily maximum temperature based on multi-source data. However, regression models often cannot fit complex relationships between variables. Machine learning methods estimate temperature based on nonlinear machine learning models. Compared with traditional regression methods, machine learning methods can effectively handle complex nonlinear relationships and achieve higher accuracy, making them the most widely used temperature estimation methods.

[0004] However, existing methods still overlook several issues: 1) Due to the inherent observational characteristics and wide field of view of polar-orbiting satellite sensors, surface temperature products obtained from these sensors can have time differences of up to two hours between different pixels on the same image, or between the same pixel on different days within a revisit period. This time difference introduces significant uncertainty when using remotely sensed surface temperature products to estimate near-surface air temperature in large-scale or long-term studies. 2) Most existing studies develop a unified model for the entire region to estimate near-surface air temperature. However, due to the diversity of natural conditions and the differences in the relationship between near-surface air temperature and environmental factors across large-scale regions, unified models fail to account for spatial heterogeneity in temperature estimates, leading to bias in near-surface air temperature estimates. 3) Existing studies use only a single machine learning algorithm to estimate near-surface air temperature, or compare multiple machine learning algorithms and select the best-performing model for temperature estimation. However, single models often have limitations and uncertainties, as the performance of different algorithms varies depending on geographic location and environmental conditions. Summary of the Invention

[0005] To address the above-mentioned problems, the present invention proposes a three-stage high-precision near-surface air temperature remote sensing estimation method, which not only effectively eliminates the uncertainty caused by the temporal differences in surface temperature, but also takes into account the differences in model adaptability in different geographical environments. Based on the idea of integration, it integrates multiple single machine learning models to improve the accuracy of remote sensing estimation of air temperature.

[0006] The present invention discloses a three-stage high-precision near-surface air temperature remote sensing estimation method, which comprises the following steps:

[0007] S1, time normalization of remotely sensed surface temperature data to eliminate imaging time differences;

[0008] S2: The study area is divided into multiple sub-areas based on differences in natural conditions. Based on the time-normalized surface temperature data and the spatial auxiliary variables corresponding to the meteorological stations in the sub-areas, four machine learning models are used to estimate the temperature in each sub-area to obtain spatial temperature estimation results.

[0009] S3 integrates the temperature estimation results of four machine learning models through a generalized additive model to obtain near-surface temperature with higher accuracy.

[0010] Step S1 further includes the following steps:

[0011] The surface temperature normalization formula is used to perform time normalization on the surface temperature data collected by the MODIS remote sensing sensor:

[0012] ;

[0013] Where, T norm is the surface temperature normalized to the reference time, T ac is the actual surface temperature during the solar period, a and b are coefficients, t ac is solar time, t norm is the base time;

[0014] Hourly surface temperatures were extracted based on the ERA5_land reanalysis meteorological data. The total cloud cover information of ERA5 was used for masking. The clear-sky phase surface temperature was taken for hourly monthly synthesis to obtain the spatiotemporally continuous monthly composite hourly ERA5_Land surface temperature data. The following formula was used to convert the UTC coordinated universal time corresponding to the ERA5_Land surface temperature data to solar time:

[0015] ;

[0016] Where, t solar_time is solar time, t UTC For Coordinated Universal Time, is the time difference between the standard time of the time zone where the study area is located and UTC, and Lon is the longitude value of the pixel;

[0017] Based on the monthly composite hourly ERA5_Land surface temperature data and the corresponding solar time, the coefficients a and b were obtained by least squares fitting based on the surface temperature normalization formula for each pixel, and the coefficients a and b were resampled to 1 km resolution to be consistent with the resolution of the MODIS surface temperature data. The MODIS surface temperature data were time-normalized using the surface temperature time normalization formula for each pixel to obtain surface temperature data normalized to 11 o'clock.

[0018] Furthermore, in step S2, the process of clustering the study area into multiple sub-areas based on differences in natural conditions includes the following steps:

[0019] The MODIS reflectance product MOD09A1, MODIS albedo product MCD43A3, MODIS albedo product MCD43A3, ERA5-Land, and SRTM / DEM data were resampled to 1 km resolution to match the resolution of MODIS land surface temperature data. Reflectances of the red, near-infrared, green, and shortwave infrared 1 bands were extracted from the preprocessed MODIS reflectance product MOD09A1. Albedo data and longitude and latitude data were extracted from the preprocessed MODIS albedo product MCD43A3. Downward solar radiation (DSR) was extracted from the preprocessed ERA5-Land data. Elevation and terrain index (CTI) were extracted from the preprocessed SRTM / DEM data.

[0020] The vegetation index NDVI and water body index MNDWI at a resolution of 500 m are calculated using the following formulas:

[0021] ;

[0022] ;

[0023] Where, ρ NIR and ρ RED Represent the reflectance of the near-infrared band and the red band respectively; ρ Green is the green band reflectivity, ρ SWIR1 is the reflectivity of the shortwave infrared band 1.

[0024] Furthermore, the albedo, elevation, downward solar radiation (DSR), vegetation index (NDVI), and water body index (MNDWI) were annually averaged and normalized to eliminate the influence of value range differences during the classification process. Longitude and latitude were combined as classification features, and an iterative self-organizing data analysis model was used for clustering and partitioning. The study area was divided into several sub-areas through iteration.

[0025] Furthermore, in step S2, the process of estimating the temperature in each sub-area using the four machine learning models includes the following steps:

[0026] By performing spatiotemporal matching of longitude and latitude and time, all spatial variables corresponding to the meteorological stations are extracted. Spatial variables include normalized surface temperature and spatial auxiliary variables including vegetation index NDVI, water body index MNDWI, albedo, elevation, downward solar radiation DSR, topographic index CTI, annual cumulative day DOY, longitude and latitude.

[0027] With all spatial variables as independent variables and the measured temperature values at meteorological stations as dependent variables, four machine learning methods, random forest, extreme gradient boosting, classification boosting, and Cubist, were used to construct temperature estimation models in each sub-area:

[0028] ;

[0029] Where, represents the estimated temperature at station i; represents the spatial coordinates at site i; is the intercept of the equation; is the regression coefficient of the corresponding independent variable; are the normalized surface temperature, vegetation index NDVI, water body index MNDWI, albedo, elevation, downward solar radiation DSR, topographic index CTI, annual cumulative day DOY, longitude and latitude at site i; is the random error term;

[0030] The four constructed machine learning models were applied to the time-normalized surface temperature and spatial auxiliary variables including vegetation index NDVI, water body index MNDWI, albedo, elevation, downward solar radiation DSR, topographic index CTI, annual cumulative days DOY, longitude and latitude to obtain spatial temperature estimation results.

[0031] Furthermore, in step S3, the temperature estimation results obtained by each machine learning model in step S2 are used as independent variables, and the measured temperature at the meteorological station is used as the dependent variable. The generalized additive model (GAM) is used to construct an integrated model:

[0032] ;

[0033] Where, Ta To observe the near-surface air temperature; Ta model 1 、Ta model 2 、Ta model 3 、Ta model 4 are the temperature estimates of the four machine learning models: random forest, extreme gradient boosting, classification enhancement, and Cubist; S() represents the smooth function of the corresponding parameters; e i,j represents the residual error term for each observation;

[0034] The constructed ensemble model is applied to the estimated temperatures obtained by the four machine learning models to obtain the final integrated temperature estimate.

[0035] The advantages and beneficial effects of the present invention are:

[0036] First, the three-stage high-precision near-surface air temperature remote sensing estimation method of the present invention performs time normalization on remotely sensed surface temperature data based on reanalysis data, effectively eliminating the uncertainty caused by observation time differences.

[0037] Second, the three-stage high-precision near-surface air temperature remote sensing estimation method of the present invention divides the study area into small areas based on the spatial clustering method, and constructs temperature estimation models in these smaller areas with more consistent environments. The adaptability of the model in different geographical environments is reduced through zoning modeling.

[0038] Third, the three-stage high-precision near-surface air temperature remote sensing estimation method of the present invention constructs an integrated model to integrate the estimation results of multiple single machine learning models. Through integration, the advantages of the estimation results of various basic models are complementary, which effectively reduces the deviation of a single model and generates more accurate estimation results. BRIEF DESCRIPTION OF THE DRAWINGS

[0039] Figure 1 This is a flow chart of the three-stage high-precision near-surface air temperature remote sensing estimation method of the present invention;

[0040] Figure 2 Comparison of MODIS LST before and after normalization on February 21, 2021: (a) is the original LST diagram, and (b) is the LST diagram after time normalization;

[0041] Figure 3 This is a schematic diagram of the partitioning results;

[0042] Figure 4a Scatter plot of the temperature estimated by the Cubist model and the measured temperature at the station;

[0043] Figure 4b Scatter plot of the temperature estimated by the Catboost model and the measured temperature at the station;

[0044] Figure 4c Scatter plot of the temperature estimated by the XGboost model and the measured temperature at the station;

[0045] Figure 4d Scatter plot of the temperature estimated by the RF model and the measured temperature at the station;

[0046] Figure 5 This is a scatter plot of the final estimated temperature and the measured temperature at the station;

[0047] Figure 6 This is the spatial distribution map of the average temperature in 2021 obtained by remote sensing estimation. DETAILED DESCRIPTION

[0048] The following examples may enable those skilled in the art to more fully understand the present invention, but are not intended to limit the present invention in any way.

[0049] The purpose of the present invention is to propose a three-stage high-precision near-surface air temperature remote sensing estimation method, which includes three stages: in the first stage, time normalization processing is performed on the remotely sensed surface temperature data to eliminate the uncertainty caused by the observation time difference; in the second stage, the spatial heterogeneity of the geographical environment within a large region is considered to divide it into multiple sub-regions with independent characteristics, and four machine learning models are used to estimate the temperature in each sub-region; in the third stage, the temperature estimation results of four single machine learning algorithms are integrated based on a generalized additive model to obtain high-precision near-surface air temperature.

[0050] 1. Time normalization of remotely sensed surface temperature

[0051] Because the MODIS remote sensing sensor swath exceeds 2,000 km and its orbital revisit period is 16 days, observed solar time varies significantly within the same swath and between adjacent swaths. Furthermore, within a 16-day revisit period, observed solar time for the same pixel on different days can vary significantly, with the difference reaching as much as two hours. Given that surface temperature is significantly affected by solar time, directly using raw remotely sensed surface temperatures to estimate near-surface air temperature is subject to significant uncertainty.

[0052] In order to eliminate the influence of different solar observation times in MODIS surface temperature data, it is necessary to perform time normalization processing. The time normalization equation for surface temperature is:

[0053] (1);

[0054] Where, T norm is the surface temperature normalized to the reference time, T ac is the actual surface temperature during the solar period, a and b are coefficients, t ac is solar time, t norm The Terra / MODIS daytime imaging time range is 9.9 ~ 12.8 (solar time), so t norm Set to 11.0.

[0055] The coefficients a and b in the time normalization equation are related to geographic location and environment and need to be solved pixel by pixel. Hourly surface temperatures are extracted based on the ERA5_land reanalysis meteorological data. Masking is performed using ERA5's total cloud cover information. Hourly surface temperatures during the clear sky phase are then composited monthly to obtain spatiotemporally continuous, monthly composite hourly ERA5_Land surface temperature data. ERA5_Land surface temperature data corresponds to UTC Coordinated Universal Time and needs to be converted to solar time using the following conversion formula:

[0056] (2);

[0057] Where, t solar_time is solar time, t UTC For Coordinated Universal Time, is the time difference between the standard time of the time zone where the study area is located and UTC. The Tzone_offset in China is 8, and Lon is the longitude value of the pixel.

[0058] Based on the monthly composite hourly ERA5_Land surface temperature data and the corresponding solar time, a least squares fit was performed pixel by pixel using formula (1) to obtain the coefficients a and b. The coefficients were then resampled to 1 km to match the resolution of the MODIS surface temperature data. Then, formula (1) was applied pixel by pixel to perform time normalization on the MODIS temperature product data, resulting in surface temperature data normalized to 11 o'clock.

[0059] 2. Temperature estimation based on regional modeling

[0060] The spatial auxiliary variables used in this paper are the vegetation index (NDVI), the water body index (MNDWI), albedo, elevation, downward solar radiation (DSR), the topographic index (CTI), cumulative day-per-year (DOY), longitude, and latitude. The MODIS reflectance product MOD09A1, MODIS albedo product MCD43A3, ERA5-Land, and SRTM / DEM data were resampled to a 1 km resolution to match the resolution of MODIS land surface temperature data. Extract the reflectance of the red band, near infrared band, green band, and shortwave infrared 1 band from the MODIS reflectance product MOD09A1, and use formula (3) and formula (4) to calculate NDVI and MNDWI respectively; extract albedo data and longitude and latitude data from the MODIS albedo product MCD43A3; extract downlink shortwave radiation DSR from ERA5-Land; obtain elevation data and CTI from SRTM / DEM. Among them, the calculation formulas of NDVI and MNDWI are:

[0061] (3);

[0062] (4);

[0063] Where, ρ NIR and ρ RED Represent the reflectance of the near-infrared band and the red band respectively, ρ Green is the green band reflectivity, ρ SWIR1 is the reflectivity of the shortwave infrared band 1.

[0064] Albedo, elevation, DSR, NDVI, and MNDWI were annually averaged and normalized to eliminate the influence of range differences during the classification process. Longitude and latitude were combined as classification features, and clustering and partitioning were performed using the iterative self-organizing data analysis model (ISODATA). The minimum and maximum number of clusters in the clustering process were set to 5 and 12, respectively. The study area was divided into several sub-areas through iteration.

[0065] By performing spatiotemporal matching of latitude, longitude, and time, we extracted spatial variables corresponding to the meteorological stations, including normalized surface temperature, NDVI, MNDWI, Albedo, elevation, DSR, CTI, DOY, Longitude, and Latitude. Using these spatial variables as independent variables and the measured temperature at the meteorological stations as the dependent variable, we constructed temperature estimation models for each subregion using four machine learning methods: Random Forest (RF), Extreme Gradient Boosting (XGBoost), Classification Boosting (CatBoost), and Cubist.

[0066] (5);

[0067] Where, represents the estimated temperature at station i; represents the spatial coordinates at site i; is the intercept of the equation; is the regression coefficient of the corresponding independent variable; are the normalized surface temperature, NDVI, MNDWI, Albedo, elevation, DSR, CTI, DOY, Longitude and Latitude at site i, respectively; is the random error term.

[0068] The four constructed machine learning models were applied to spatial independent variables such as surface temperature, NDVI, MNDWI, Albedo, elevation, DSR, CTI, DOY, Longitude and Latitude after time normalization to obtain spatial temperature estimation results.

[0069] (3) Temperature estimation based on integrated models

[0070] The estimated temperature obtained by each machine learning model in the second stage is used as the independent variable, and the measured temperature at the meteorological station is used as the dependent variable. The generalized additive model (GAM) is used to construct an integrated model:

[0071] (6);

[0072] Where, Ta To observe the near-surface air temperature; Ta model 1 、Ta model 2 、Ta model 3 、Ta model 4 are the estimated values of temperature by Cubist, XGBoost, Catboost and RF models respectively; S() represents the smooth function of the corresponding parameters; e i,j represents the residual error term for each observation.

[0073] The constructed integrated model was applied to the estimated temperature obtained by the four machine learning methods to obtain the final integrated temperature estimate.

[0074] Examples

[0075] The present invention is to estimate the near-surface temperature with high precision based on remote sensing data. The following are the specific implementation steps of the example. The technical process is shown in Figure 1 .

[0076] 1) The study area is located between 95° and 110° east longitude and 29° and 37° north latitude. Surface temperature and cloud cover information are extracted from the ERA5_Land and ERA5 data, respectively. Cloud mask processing and monthly synthesis are performed to generate spatiotemporally continuous, monthly-scale, hourly ERA5_Land surface temperature data. Equation (2) is used to convert the UTC time corresponding to the ERA5_Land surface temperature to solar time. The coefficients a and b in Equation (1) are then obtained by least squares fitting pixel by pixel. These values are then resampled to a 1 km resolution to match the MODIS surface temperature.

[0077] 2) Extract the clear sky surface temperature data and the corresponding solar time from the MODIS surface temperature product MOD11A1, and use formula (1) to time-normalize the MODIS surface temperature values pixel by pixel. Figure 2 A comparison of MODIS surface temperature before and after time normalization on February 21, 2021 is given. It can be seen that before time normalization (see Figure 2 In (a), the surface temperature at the edge of the scanning strip is significantly inconsistent due to the difference in satellite passing time. After time normalization correction (see Figure 2 (b) in the figure significantly reduces the inconsistency, and the spatial distribution of surface temperature is more continuous and uniform, effectively eliminating the temperature deviation caused by different imaging of the sun.

[0078] 3) The MODIS reflectance product MOD09A1, MODIS albedo product MCD43A3, MODIS albedo product MCD43A3, ERA5-Land, and SRTM / DEM data were resampled to 1 km resolution to match the MODIS surface temperature data resolution. The red, near-red, green, and shortwave infrared 1 band reflectances were extracted from the MODIS surface reflectance product MOD09A1, and the vegetation index NDVI and water body index MNDWI were calculated using formulas (3) and (4), respectively. Albedo and the corresponding longitude and latitude were extracted from the MODIS albedo product MCD43A3. Downward shortwave radiation DSR was extracted from ERA5-Land. Elevation and CTI data were extracted from SRTM / DEM.

[0079] 4) Annual average processing of Albedo, DEM, DSR, NDVI, and MNDWI was performed, and then used together with Longitude and Latitude as classification features, and the ISODATA iterative self-organizing method was used for classification to divide the study area into several sub-areas (see Figure 3 ). In the ISODATA clustering process, the minimum and maximum number of categories are set to 5 and 12 respectively.

[0080] 5) By performing spatiotemporal matching based on location, time, and region, we extracted the normalized surface temperature, NDVI, MNDWI, Albedo, DSR, CTI, elevation, longitude, latitude, and DOY values corresponding to the observed station temperature to construct a temperature estimation sample set. Using the measured station temperature as the dependent variable and the normalized surface temperature, NDVI, MNDWI, Albedo, DSR, CTI, elevation, longitude, latitude, and DOY values of the corresponding pixels as independent variables, we constructed a near-surface temperature estimation model for each subregion using the Cubist, XGBoost, RF, and Catboost methods (see Equation (5)). To optimize the model estimates, cross-validation was used to adjust the model parameter values during the fitting process. Figures 4a to 4d A scatter plot of the temperature estimation accuracy of four machine learning models is given.

[0081] 6) The four machine learning models constructed for each subregion were applied to the normalized surface temperature, NDVI, MNDWI, Albedo, DSR, CTI, elevation, longitude, latitude, and DOY spatial variables to calculate the estimated temperature for each subregion. The estimated temperatures for each subregion obtained by the four machine learning models were then mosaicked to obtain the estimated temperature for the entire study area using the four machine learning models.

[0082] 7) The estimated temperatures of the four machine learning models for the corresponding pixels of each meteorological station were extracted respectively. The measured temperature of the station was used as the dependent variable, and the estimated temperatures of the four machine learning models for the corresponding pixels were used as the independent variables. A generalized additive model was used to construct an integrated model for estimating temperature using the four machine learning methods (see formula (6)).

[0083] 8) Apply the constructed integrated model to the four machine learning models obtained in step 7 to estimate the temperature, and calculate the final remote sensing estimation results of the near-surface temperature in the study area. Figure 5 The verification scatter plot of the integrated results is given. The sample points are mainly distributed near the 1:1 line, and the determination coefficient R 2The average absolute error (MAE) is 0.99, the bias is approximately 0, and there is no significant overestimation or underestimation. Compared with the four single machine learning methods, the ensemble method has significantly improved the accuracy of temperature estimation. Figure 6 A spatial distribution map of the average temperature in the study area for 2021, obtained using the method presented in this paper, is provided. The annual average temperature ranged primarily from -8°C to 27°C, exhibiting significant spatial variability. Temperatures were higher in the southeast, typically between 15°C and 25°C, while temperatures were lower in the northwest, generally below 0°C.

[0084] The above are merely preferred embodiments of the present invention. The scope of protection of the present invention is not limited to the above embodiments. All technical solutions based on the principles of the present invention are within the scope of protection of the present invention. It should be noted that for those skilled in the art, various improvements and modifications that do not depart from the principles of the present invention should be considered within the scope of protection of the present invention.

Claims

1. A three-stage high-precision near-surface air temperature remote sensing estimation method, characterized in that: The method comprises the following steps: S1, time normalization of remotely sensed surface temperature data to eliminate imaging time differences; S2: The study area is divided into multiple sub-areas based on differences in natural conditions. Based on the time-normalized surface temperature data and the spatial auxiliary variables corresponding to the meteorological stations in the sub-areas, four machine learning models are used to estimate the temperature in each sub-area to obtain spatial temperature estimation results. S3, integrates the temperature estimation results of four machine learning models through a generalized additive model to obtain a more accurate near-surface temperature; Step S1 further includes the following steps: The surface temperature normalization formula is used to perform time normalization on the surface temperature data collected by the MODIS remote sensing sensor: T norm =T ac -a(t ac -t norm ) 2 +b(t ac -t norm ); Where, T norm is the surface temperature normalized to the reference time, T ac is the actual surface temperature during the solar period, a and b are coefficients, t ac is solar time, t norm is the base time; Hourly surface temperatures were extracted based on the ERA5_land reanalysis meteorological data. The total cloud cover information of ERA5 was used for masking. The clear-sky phase surface temperature was taken for hourly monthly synthesis to obtain the spatiotemporally continuous monthly composite hourly ERA5_Land surface temperature data. The following formula was used to convert the UTC coordinated universal time corresponding to the ERA5_Land surface temperature data to solar time: Where, t solar_time is solar time, t UTC is Coordinated Universal Time, Tzone_offset is the time difference between the standard time of the time zone where the study area is located and UTC, and Lon is the longitude value of the pixel; Based on the monthly composite hourly ERA5_Land surface temperature data and the corresponding solar time, the coefficients a and b are obtained by least squares fitting based on the surface temperature normalization formula for each pixel, and then resampled to a 1 km resolution to be consistent with the resolution of the MODIS surface temperature data; the MODIS surface temperature data is time-normalized by applying the surface temperature time normalization formula for each pixel to obtain surface temperature data normalized to 11 o'clock; in step S2, the process of clustering the study area into multiple sub-areas based on differences in natural conditions includes the following steps: Annual average and normalization processing was performed on albedo, elevation, downward solar radiation (DSR), vegetation index (NDVI), and water body index (MNDWI) to eliminate the influence of value range differences during the classification process. Longitude and latitude were combined as classification features, and an iterative self-organizing data analysis model was used for clustering and partitioning. The study area was divided into several sub-areas through iteration. In step S2, the process of estimating the temperature in each sub-area using the four machine learning models includes the following steps: By performing spatiotemporal matching of longitude and latitude and time, all spatial variables corresponding to the meteorological stations are extracted. Spatial variables include normalized surface temperature and spatial auxiliary variables including vegetation index NDVI, water body index MNDWI, albedo, elevation, downward solar radiation DSR, topographic index CTI, annual cumulative day DOY, longitude and latitude. With all spatial variables as independent variables and the measured temperature values at meteorological stations as dependent variables, four machine learning methods, random forest, extreme gradient boosting, classification boosting, and Cubist, were used to construct temperature estimation models in each sub-area: y i ~β0(μ i ,v i )+β1(μ i ,v i )x i1 +......+b 10 (m i ,v i )x i10 +e i ; Where y i represents the estimated temperature at station i; (μ i , v i ) represents the spatial coordinates of site i; β0(μ i , v i ) is the intercept of the equation; β1(μ i , v i )~β 10 (μ i , v i ) is the regression coefficient of the corresponding independent variable; x i1 ~x i10 are the normalized surface temperature, vegetation index NDVI, water body index MNDWI, albedo, elevation, downward solar radiation DSR, topographic index CTI, accumulated days per year DOY, longitude and latitude at site i; ε i is the random error term; The four machine learning models constructed were applied to the time-normalized surface temperature and spatial auxiliary variables including vegetation index NDVI, water body index MNDWI, albedo, elevation, downward solar radiation DSR, topographic index CTI, annual cumulative days DOY, longitude and latitude to obtain spatial temperature estimation results; In step S3, the temperature estimation results obtained by each machine learning model in step S2 are used as independent variables, and the measured temperature at the meteorological station is used as the dependent variable. The generalized additive model (GAM) is used to construct an integrated model: Yes(Yes model1 )+s(Ta model2 )+s(Ta model3 )+s(Ta model4 )+e i，j 4 Where Ta is the observed near-surface air temperature; Ta model 1 、Ta model 2 、Ta model 3 、Ta model 4 are the temperature estimates of the four machine learning models: random forest, extreme gradient boosting, classification enhancement, and Cubist; S() represents the smooth function of the corresponding parameters; e i,j represents the residual error term for each observation; The constructed ensemble model is applied to the estimated temperatures obtained by the four machine learning models to obtain the final integrated temperature estimate.

2. The three-stage high-precision near-surface air temperature remote sensing estimation method according to claim 1 is characterized in that: Resample the MODIS reflectance product MOD09A1, MODIS albedo product MCD43A3, ERA5-Land, and SRTM / DEM data to 1 km resolution to match the resolution of MODIS surface temperature data; extract the reflectance of the red band, near-infrared band, green band, and shortwave infrared 1 band from the preprocessed MODIS reflectance product MOD09A1; extract albedo data and longitude and latitude data from the preprocessed MODIS albedo product MCD43A3; extract solar downwelling radiation (DSR) from the preprocessed ERA5-Land; and extract elevation and terrain index (CTI) from the preprocessed SRTM / DEM data. The vegetation index NDVI and water body index MNDWI at a resolution of 500 m are calculated using the following formulas: Where, ρ NIR and ρ RED Represent the reflectance of the near-infrared band and the red band respectively; ρ Green is the green band reflectivity, ρ SWIR1 is the reflectivity of the shortwave infrared band 1.

Images