A remote sensing air temperature reconstruction method based on kernel function time dimension filling and bias correction

Abstract

The application discloses a remote sensing air temperature reconstruction method based on kernel function time dimension filling and deviation correction, and comprises the following steps: taking the clear sky air temperature retrieved by thermal infrared remote sensing as the basis, adopting a time smoothing method based on kernel function weighting, comprehensively considering the annual periodic change and short-term dynamic fluctuation of air temperature, reconstructing the theoretical clear sky air temperature of a target date, and filling the missing time dimension of the remote sensing air temperature caused by cloud cover; quantifying the systematic deviation between the theoretical clear sky air temperature under the cloud and the actually observed cloud air temperature, combining a plurality of source environment variables to construct a deviation correction model to calculate the air temperature deviation term caused by the cloud disturbance; applying the deviation correction model to the pixel scale, correcting the theoretical clear sky air temperature under the cloud, obtaining the air temperature estimation value closer to the actual cloud condition, and generating continuous, seamless all-weather remote sensing air temperature data. The application can solve the technical problem of data missing of remote sensing air temperature estimation under cloud cover.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of thermal infrared remote sensing, and particularly relates to a remote sensing air temperature reconstruction method based on kernel function time dimension filling and bias correction. BACKGROUND

[0002] Near-surface air temperature is a core meteorological element in the process of surface-atmosphere energy exchange, and is widely used in climate change assessment, ecological environment monitoring, hydrological process simulation, agricultural production management and public health warning. Traditional air temperature observation mainly relies on ground meteorological stations, which has high observation accuracy and good time resolution, but due to the sparse distribution of stations, it is difficult to effectively reflect the spatial distribution characteristics of air temperature at regional scale and larger scale. Reanalysis products (such as ERA5, MERRA2) based on numerical model and multi-source observation assimilation technology can provide continuous, long-term gridded air temperature data, and have the advantages of global coverage and long-term stability. However, the resolution of reanalysis data is usually coarse (usually >10km), which makes it difficult to accurately capture local air temperature changes in areas with complex terrain or strong heterogeneity of underlying surface.

[0003] Thermal infrared remote sensing technology has become an important tool for studying regional scale surface characteristics due to its wide range, periodic and multi-temporal observation capability. Remote sensing inversion of surface or near-surface air temperature products has significant advantages in spatial resolution and coverage, providing an important supplement for high-precision spatio-temporal analysis of air temperature. However, thermal infrared remote sensing is limited by cloud cover and cannot obtain effective observations under cloudy conditions, resulting in a large number of missing values in air temperature products, which seriously restricts its application in continuous air temperature monitoring. To solve this problem, existing research has developed multi-source data fusion and spatio-temporal interpolation methods to try to fill the gaps in remote sensing air temperature data. Multi-source data fusion methods integrate reanalysis products or other satellite observations to realize air temperature gap filling by taking advantage of multiple sources. However, this method has problems such as mismatched spatial resolution and large differences in physical observation mechanisms, resulting in scale effects and insufficient physical consistency in the fusion results; spatio-temporal interpolation methods mainly estimate missing values based on the correlation of time series or spatial neighborhood information. Although it can alleviate the missing value problem to some extent, the results usually reflect the theoretical clear-sky air temperature and cannot effectively reflect the adjustment of cloud cover to the underlying surface radiation balance, resulting in systematic bias between the reconstruction results and the actual near-surface air temperature under cloudy or overcast conditions. SUMMARY

[0004] In view of the technical problem of data missing of remote sensing air temperature estimation under cloud coverage, the present application provides a remote sensing air temperature reconstruction method based on kernel function time dimension filling and bias correction, which breaks through the dependence of thermal infrared remote sensing on clear sky conditions, can effectively reconstruct the near-surface air temperature under the condition that the air temperature data is missing due to cloud cover, and significantly improves the spatiotemporal continuity of remote sensing air temperature products, providing reliable data support for climate monitoring, environmental management and agricultural production.

[0005] To achieve the above technical purpose, the technical scheme adopted by the present application is:

[0006] A remote sensing air temperature reconstruction method based on kernel function time dimension filling and bias correction, the method comprising the following steps:

[0007] S1, based on the clear sky air temperature retrieved by thermal infrared remote sensing, a time smoothing method based on kernel function weighting is used to consider the annual periodicity of air temperature and short-term dynamic fluctuations, and the seasonal air temperature of the same period in history and the effective air temperature of the adjacent date are used to reconstruct the theoretical clear sky air temperature of the target date by weighted average, to fill the time dimension of the missing remote sensing air temperature caused by cloud cover, and obtain the theoretical clear sky air temperature under cloud cover;

[0008] S2, based on the measured air temperature of the ground meteorological station, the systematic bias between the theoretical clear sky air temperature under cloud cover and the actual observed cloud cover air temperature is quantified, and a bias correction model is constructed combined with multi-source environmental variables to calculate the air temperature bias term caused by cloud disturbance;

[0009] S3, applying the bias correction model to the pixel scale to correct the theoretical clear sky air temperature under cloud cover, and obtaining the cloud air temperature estimate value closer to the actual cloud condition;

[0010] S4, fusing the clear sky air temperature retrieved by thermal infrared remote sensing and the corrected cloud air temperature estimate value to generate continuous and seamless all-weather remote sensing air temperature data.

[0011] Further, the method further comprises:

[0012] Obtaining remote sensing air temperature related data, including remote sensing data, reanalysis data and ground meteorological station observation data; wherein the remote sensing data includes remote sensing clear sky air temperature data set, MODIS land surface temperature data, MODIS reflectivity data and MODIS albedo data and SRTM digital elevation model data; the reanalysis data includes ERA5 data and derived ERA5-Land data set;

[0013] Extracting the clear sky air temperature retrieved by thermal infrared remote sensing in the research period range from the remote sensing clear sky air temperature data set;

[0014] The daily scale ground temperature data and its quality control band are extracted from the MODIS ground temperature data, cloud coverage pixels are identified by analyzing the QC information, and then the daily cloud continuous coverage days CCD of each pixel are calculated;

[0015] The red band, green band, near-infrared band and short-wave infrared band are extracted from the MODIS reflectivity data, which are used to calculate the normalized vegetation index NDVI and the normalized water body index MNDWI; the spatially complete normalized vegetation index NDVI and the normalized water body index MNDWI are obtained by using the monthly scale synthesis method;

[0016] The ground albedo data Albedo is extracted from the MODIS albedo data, and the cloud coverage influence is reduced by monthly synthesis;

[0017] The elevation information Elevation is extracted from the SRTM digital elevation model data, and the comprehensive terrain index CTI is calculated;

[0018] All remote sensing data are projected by Albers equal-area projection, and resampled to 1km spatial resolution by using bilinear interpolation method, so as to keep spatial consistency among data sources;

[0019] The hourly cloud base height data is extracted from the ERA5 data set, and the daily cloud base height CBH is obtained after daily averaging; the daily surface downward radiation data DSR is extracted from the ERA5-Land data set; the daily cloud base height and the daily surface downward radiation data are projected by Albers equal-area projection to resample to 1km spatial resolution, so as to keep consistency with remote sensing data;

[0020] The daily average air temperature of the station is extracted from the ground meteorological station observation data, and the extracted station average air temperature data is cleaned, de-redundant and quality detected.

[0021] Further, in step S1, the theoretical clear sky air temperature of the target date is reconstructed by using the following formula:

[0022] ;

[0023] In the formula, is the theoretical clear sky air temperature of the target day i, is the effective remote sensing air temperature observation value of the reference day j, and S is a set containing available reference days, is the weight assigned to the reference day j by using the kernel function weighting strategy;

[0024] ;

[0025] In the formula, is a seasonal kernel function for reflecting the annual periodic change characteristics of air temperature, the circular day order difference between the target day and the reference day in a year; is a local kernel function that captures short-term temporal similarity by evaluating the absolute temporal distance between the target date and the reference date, denotes the absolute temporal distance between the target day and the reference day; a is a dynamic adjustment coefficient that adaptively adjusts the relative contribution of the seasonal kernel and the local kernel in the weight based on a dynamic weighting mechanism of the density of valid observations and the intensity of short-term variations.

[0026] Further, the seasonal kernel function is:

[0027] ;

[0028] wherein, denotes the circular day order difference between the target day and the reference day in a year, is a bandwidth parameter of the seasonal kernel; denotes the annual accumulated day of the target day i, denotes the annual accumulated day of the reference day j.

[0029] The local kernel function is:

[0030] ;

[0031] wherein, denotes the absolute temporal distance between the target day and the reference day, is a bandwidth parameter of the local kernel; and denote the absolute time of the target day i and the reference day j, respectively.

[0032] The dynamic adjustment coefficient a is:

[0033] ;

[0034] ;

[0035] ;

[0036] wherein, is the density of valid observations, is the intensity of short-term variations, denotes the number of valid observations within a symmetric time window of size w, is the standard deviation of the near-surface air temperature within the symmetric time window.

[0037] The step S2 further comprises:

[0038] The multi-source environmental variables corresponding to the research area are screened from four aspects of cloud, solar radiation, surface characteristics and time characteristics, including: daily cloud continuous coverage days CCD and cloud base height CBH for reflecting cloud-related characteristics;

[0039] normalized vegetation index NDVI, normalized water body index MNDWI, surface albedo Albedo, elevation information Elevation, longitude Lon, latitude Lat and comprehensive terrain index CTI for reflecting surface characteristics;

[0040] annual accumulated day DOY for reflecting time characteristics;

[0041] XGBoost algorithm is used for regression modeling to build a bias correction model:

[0042] ;

[0043] Wherein is the difference between the theoretical clear sky temperature and the actual cloud temperature at each site, and f is the nonlinear relationship between the systematic bias of the site and the multi-source environmental variables corresponding to the research area screened out; CCD, CBH, NDVI, Albedo, DSR, MNDWI, Elevation, CTI, Lon, Lat and DOY are the values of cloud continuous coverage days CCD, cloud base height CBH, normalized vegetation index NDVI, surface albedo Albedo, surface downward radiation DSR, normalized water body index MNDWI, elevation Elevation, comprehensive terrain index CTI, longitude Lon, latitude Lat and annual accumulated day DOY, respectively;

[0044] The model parameters of the bias correction model are optimized by random search and five-fold cross-validation with the goal of minimizing the mean absolute error; During cross-validation, the ground meteorological stations are randomly divided into several subsets, and for each iteration, part of the subsets are used for model training, while the remaining subsets are used as the test set for model validation; Repeat the iteration to ensure that each subset is used as a validation set;

[0045] The model performance of the bias correction model is comprehensively considered by using various evaluation indexes including coefficient of determination R², mean absolute error MAE, root mean square error RMSE and bias Bias, and the model with the best performance is used as the final bias correction model.

[0046] Step S3 further comprises:

[0047] The established bias correction model is applied to the spatialized variable characteristics, and the air temperature bias term of all cloud-free pixels is estimated at the pixel level; The estimated air temperature bias term is superimposed on the theoretical clear sky temperature under the cloud, and the corrected cloud temperature estimation value is obtained:

[0048] ;

[0049] wherein, is the cloud-cleared air temperature estimated by cloud pixel correction, is the theoretical clear-sky air temperature obtained by time dimension reconstruction, is the air temperature bias term estimated at pixel scale.

[0050] Compared with the prior art, the present application has the following beneficial effects:

[0051] Firstly, the remote sensing air temperature reconstruction method based on kernel function time dimension filling and bias correction of the present application develops a time dimension filling method based on kernel function, innovatively fuses two complementary kernel functions of seasonal kernel and local kernel, effectively represents the periodic change and short-term dynamic characteristics of air temperature, and through time weighted smoothing, the effective filling of multi-temporal remote sensing air temperature data can be realized under the condition that the remote sensing air temperature is affected by cloud, and the spatio-temporal continuity of air temperature data is significantly improved.

[0052] Secondly, the remote sensing air temperature reconstruction method based on kernel function time dimension filling and bias correction of the present application aims at the limitation that the traditional spatio-temporal interpolation method can only estimate the theoretical clear-sky air temperature and is difficult to reflect the actual cloud-cleared air temperature, and based on the meteorological station measured data, a bias correction model is established. The model is applied to the spatial scale to realize the bias correction of the pixel-by-pixel air temperature, corrects the theoretical clear-sky air temperature to the real cloud-cleared air temperature, and thus the accuracy and physical interpretability of the reconstruction result are improved.

[0053] Thirdly, the remote sensing air temperature reconstruction method based on kernel function time dimension filling and bias correction of the present application has the characteristics of simple method and convenient implementation, and the reconstruction process does not need to rely on multi-source reanalysis data or complex physical model, but only needs to be based on the clear-sky air temperature inversed by remote sensing and the observation of meteorological station, combined with statistical modeling and kernel function weighting mechanism, so that high-precision air temperature reconstruction can be realized. The present application avoids the problems of resolution mismatch, inconsistent physical mechanism and difficult data acquisition in the traditional multi-source fusion method, and has strong independence and operability. BRIEF DESCRIPTION OF DRAWINGS

[0054] Figure 1 is the flow chart of the remote sensing air temperature reconstruction method based on kernel function time dimension filling and bias correction of the present application;

[0055] Figure 2 is the theoretical clear-sky air temperature and station observed air temperature scatter plot;

[0056] Figure 3 is the corrected air temperature and station observed air temperature scatter plot;

[0057] Figure 4 is the remote sensing clear-sky air temperature and all-weather air temperature spatial distribution map. DETAILED DESCRIPTION

[0058] The embodiments of the present application will be further described in detail with reference to the accompanying drawings.

[0059] The application discloses a remote sensing air temperature reconstruction method based on kernel function time dimension filling and bias correction.

[0060] S1, based on the clear-sky air temperature retrieved by thermal infrared remote sensing, a time smoothing method based on kernel function weighting is adopted, the annual periodic change and short-term dynamic fluctuation of air temperature are comprehensively considered, the historical seasonal air temperature of the same period and the effective air temperature of the adjacent date are utilized, the theoretical clear-sky air temperature of the target date is reconstructed by means of weighted average, the missing time dimension of the remote sensing air temperature caused by cloud cover is filled, and the theoretical clear-sky air temperature under the cloud is obtained.

[0061] S2, based on the measured air temperature of the ground meteorological station, the systematic deviation between the theoretical clear-sky air temperature under the cloud and the actually observed cloud-free air temperature is quantified, a bias correction model is constructed in combination with multi-source environmental variables to calculate the air temperature deviation caused by cloud disturbance.

[0062] S3, the bias correction model is applied to the pixel scale, the theoretical clear-sky air temperature under the cloud is corrected, and the cloud air temperature estimation value closer to the actual cloud condition is obtained.

[0063] S4, the clear-sky air temperature retrieved by thermal infrared remote sensing and the corrected cloud air temperature estimation value are fused to generate continuous and seamless all-weather remote sensing air temperature data.

[0064] The application includes two core steps: (1) developing a time dimension filling method based on kernel function to preliminarily estimate the theoretical clear-sky air temperature under the cloud; (2) based on the measured air temperature of the meteorological station, a bias correction model is constructed on the basis of the estimated theoretical clear-sky air temperature under the cloud, and the bias correction model is applied to the pixel scale to correct the estimated theoretical clear-sky air temperature to the real air temperature under the cloud.

[0065] (I) Data preparation and pretreatment

[0066] The data used in the application includes three types: remote sensing data, reanalysis data and ground meteorological station observation data.

[0067] The remote sensing data included a clear-air temperature dataset from 2019 to 2023, MODIS land surface temperature product MOD11A1, MODIS reflectance product MOD09A1, MODIS albedo product MCD43A3, and SRTM DEM data. Daily-scale land surface temperature data and its quality control (QC) bands were extracted from MODIS land surface temperature product MOD11A1. Cloud cover pixels were identified by analyzing QC information, and the daily cloud cover duration (CCD) for each pixel was calculated. Red, green, near-infrared, and shortwave infrared bands were extracted from MODIS reflectance product MOD09A1 to calculate the Normalized Difference Vegetation Index (NDVI) and Normalized Difference Water Index (MNDWI). To mitigate the impact of cloud cover, a monthly-scale composite method was used to obtain spatially complete NDVI and MNDWI. Surface albedo data were extracted from the MODIS albedo product MCD43A3, and the impact of cloud cover was reduced through monthly composite processing. Elevation information was extracted from the SRTM digital elevation model (DEM) data, and the Comprehensive Topographic Index (CTI) was calculated. All the remote sensing data were subjected to Albers equal-area projection and resampled to a spatial resolution of 1 km using bilinear interpolation to ensure spatial consistency between data sources.

[0068] The reanalysis data primarily consisted of ERA5 data and the derived ERA5-Land dataset. Hourly cloud base height data were extracted from the ERA5 dataset and averaged daily to obtain the daily cloud base height (CBH). Daily surface downscattered radiation (DSR) data for the corresponding time periods were extracted from the ERA5-Land dataset. All these variables were also projected using Albers projection and uniformly resampled to a 1km spatial resolution to ensure consistency with the remote sensing data.

[0069] The meteorological station data comes from the China Meteorological Administration and provides daily average temperatures. The data undergoes cleaning, redundancy removal, and quality testing.

[0070] (ii) Kernel-based temporal dimension filling

[0071] To address the issue of missing remote sensing temperature data due to cloud cover, this invention proposes a kernel-weighted time-dimensional filling method. This method is based on the assumption that the temperature on a target date is influenced by both its annual seasonal cycle and the short-term dynamic changes of neighboring dates. Therefore, the theoretical clear-air temperature for the target date can be reconstructed by using historical seasonal temperatures from the same period and effective temperatures from neighboring dates through a weighted average. The basic expression of the reconstruction model is as follows:

[0072] (1);

[0073] In the formula, T is the theoretical clear sky air temperature of the target day i, S is the set of available reference days, and is the weight assigned to the reference day j.

[0074] Weight is assigned using a kernel function weighting strategy. Kernel function is a mathematical tool to measure similarity, which usually assigns higher weights to observations that are closer or more similar. Considering the annual periodicity and short-term variability of air temperature, two complementary kernel functions are designed in this invention, and the optimal weights are determined by weighted combination:

[0075] (2).

[0076] where, is the seasonal kernel, is the circular day order difference between the target day and the reference day within a year, is the local kernel, denotes the absolute time distance between the target day and the reference day, is the dynamic adjustment coefficient (0-1), which controls the relative contribution of the seasonal kernel and the local kernel in the weight.

[0077] ① Kernel function design

[0078] The seasonal kernel is used to reflect the annual periodicity of air temperature. The Gaussian kernel function is adopted:

[0079] (3).

[0080] where, denotes the circular day order difference between the target day and the reference day within a year, is the bandwidth parameter of the seasonal kernel, which is set to 30. Since the Gaussian kernel is non-zero over the entire domain, when is greater than , the seasonal kernel is set to 0 to ensure that only seasonally relevant observations contribute to the reconstruction process.

[0081] The local kernel captures short-term temporal similarity by evaluating the absolute time distance between the target date and the reference date. It assigns higher weights to observations close to the target date, also using the Gaussian kernel function expression:

[0082] (4).

[0083] where, is the absolute time distance between the target day and the reference day, is the bandwidth parameter of the local kernel, which is set to 5 in this embodiment to reflect the sensitivity of short-term changes.

[0084] ②Adaptive adjustment of dynamic adjustment coefficient α

[0085] In order to improve the adaptability of the model under different degrees of missing data, a dynamic weighting coefficient α is introduced to adaptively adjust the contribution of seasonal kernel and local kernel in the maximum weight. The application designs a dynamic weighting mechanism based on the density of effective observations and the intensity of short-term changes, and the expression of α is defined as:

[0086] (5);

[0087] (6);

[0088] (7);

[0089] In the formula, is the density of effective observations, is the intensity of short-term changes, is the number of effective observations in the symmetric time window of w size, is the standard deviation of Ta in the window.

[0090] (Three) Bias correction based on site

[0091] The multi-cloud pixel air temperature obtained in the foregoing time dimension filling step, although the missing values are reconstructed, still reflects the air temperature under the theoretical clear sky condition, and cannot reflect the actual adjustment effect of the cloud layer on the surface radiation balance and air temperature. Therefore, further bias correction is needed. The kernel function filling method has the advantages of simultaneously capturing seasonal and local fluctuations, and the reconstructed theoretical clear sky air temperature has strong smoothness and physical consistency in the time dimension. On this basis, the measured air temperature under the cloud at the site is used to capture the air temperature bias term, which can not only effectively reduce the interference of random noise, but also better depict the systematic influence of the cloud layer on the air temperature. Therefore, the application further quantifies the systematic bias between the theoretical clear sky air temperature under the cloud and the actual observed cloud air temperature based on the measured air temperature of the ground meteorological station, and constructs a bias correction model combined with multiple source environmental variables. The bias correction model is applied to the pixel scale, so as to realize the correction of the theoretical clear sky air temperature to the real air temperature under the cloud.

[0092] 1. Constructing a bias correction model

[0093] The meteorological station can provide all-weather (clear sky and under the cloud) near-surface air temperature observation. By using this characteristic, the systematic difference (denoted as ΔTa) between the theoretical clear sky air temperature estimated at the site scale and the measured air temperature under the cloud at the site can be quantified. The temperature difference reflects the influence of the cloud on the surface thermal conditions, and the bias correction model is developed by relating it to the environmental variables.

[0094] The weather station temperature observation under clear sky conditions is screened by using MODIS cloud mask, and then the theoretical clear sky temperature value of the weather station during cloudy period is estimated by using the time dimension filling method based on kernel function. The difference (ΔTa) between the theoretical clear sky temperature and the actual cloud temperature is calculated. This difference reflects the influence of cloud-induced disturbance on land-atmosphere radiation exchange and surface heat regulation. Therefore, it can be simulated by cloud, solar radiation, time characteristics and surface characteristics. The invention develops 11 bias correction models that can be freely selected according to the actual situation of the research area:

[0095] · Cloud-related variables: daily cloud continuous coverage days (CCD), daily cloud base height (CBH);

[0096] · Surface characteristic variables: normalized vegetation index (NDVI), normalized water body index (MNDWI), surface albedo (Albedo);

[0097] · Radiation variables: surface downward radiation (DSR);

[0098] · Topographic characteristic variables: elevation (Elevation), comprehensive terrain index (CTI);

[0099] · Spatial position variables: longitude, latitude;

[0100] · Time variable: annual accumulated day (DOY);

[0101] The bias correction model uses XGBoost algorithm for regression modeling. XGBoost is a high-efficiency gradient boosting tree model with strong non-linear fitting ability and can handle complex interactions between multiple variables. The model expression is as follows:

[0102] (8);

[0103] Wherein is the difference between the theoretical clear sky temperature and the actual cloud temperature of each site, and f is the nonlinear relationship between the site bias term and CCD (VCCD), CBH (VCBH), NDVI (VNDVI), MNDWI (VMNDWI), albedo (Valbedo), DSR (VDSR), elevation (VEle), CTI (VCTI), longitude (VLon), latitude (VLat), DOY (VDOY) and other 11 independent variables.

[0104] Model parameters are tuned by random search with five-fold cross-validation (CV) to minimize the mean absolute error (MAE) as the objective. During the cross-validation process, the weather stations are randomly divided into 5 subsets. For each iteration, four subsets (80% of the stations) are used for model training, while the remaining subset (20% of the stations) is used as the test set for model validation. This process is repeated five times to ensure that each subset is used as the validation set. Model performance is evaluated using a combination of metrics, including the coefficient of determination (R²), mean absolute error (MAE), root mean square error (RMSE), and bias (Bias). The model with the best performance is selected as the final bias correction model.

[0105] 2. Correcting the true air temperature under clouds

[0106] Based on the established bias correction model, it is applied to the spatialized variable characteristics to estimate the air temperature bias term of all cloud-free pixels at the pixel level. The estimated air temperature bias term is superimposed on the theoretical clear-sky air temperature to obtain the corrected cloud-free air temperature estimate, and the formula is:

[0107] (9);

[0108] In the formula, is the corrected cloud-free air temperature estimate of the cloud-free pixel, is the reconstructed theoretical clear-sky air temperature in the time dimension, is the estimated air temperature bias term at the pixel scale.

[0109] Finally, the reconstructed results under clear-sky conditions are fused with the corrected cloud-free air temperature results to generate continuous, seamless all-weather remote sensing air temperature data products, providing high spatial and temporal resolution data support for related application research.

[0110] The present application is based on the all-weather near-surface air temperature reconstruction of remote sensing data, and the following are the specific implementation steps of the example. The technical process is shown in Figure 1 .

[0111] (1) The study area is between 100-115 °E and 30-40 °N. The remote sensing clear-sky air temperature dataset from 2019 to 2023 is prepared. The daily surface temperature data and its quality control (QC) band are extracted from the MODIS land surface temperature product MOD11A1. By analyzing the QC information, cloud-covered pixels are identified, and the daily cloud cover duration (CCD) of each pixel is calculated. The red, green, near-infrared, and shortwave infrared bands are extracted from the MODIS reflectance product MOD09A1. The 500-meter resolution normalized vegetation index (NDVI) and normalized water index (MNDWI) are obtained by band operation. The monthly synthesis of NDVI and MNDWI is performed to obtain the monthly scale normalized vegetation index NDVI and normalized water index MNDWI. The 500-meter resolution albedo data (Albedo) is extracted from the MODIS albedo product MCD43A3, and the monthly synthesis is performed to obtain the monthly scale albedo Albedo. The elevation data is extracted from the SRTM DEM data, and the comprehensive terrain index CTI is calculated.

[0112] The 0.25 ° resolution hourly cloud base height data is extracted from the ERA5 data, and the daily average cloud base height CBH is obtained. The 0.01 ° resolution daily surface downward radiation data DSR is extracted from the ERA5-Land data. The above data is projected to Albers projection and resampled to 1 km spatial resolution using bilinear interpolation.

[0113] (2) The kernel function-based time dimension filling method (equations (1) to (7)) is used to estimate the theoretical clear-sky air temperature under cloud cover on a pixel-by-pixel basis. Figure 2 The scatter plot between the theoretical clear-sky air temperature and the station observed air temperature is drawn. The results show that the filling accuracy is: R² between 0.89 and 0.92, the mean absolute error (MAE) between 2.50 °C and 2.73 °C, the root mean square error (RMSE) between 3.16 °C and 3.48 °C, and the bias (Bias) between 1.32 °C and 1.82 °C. The positive bias of the consistency indicates that the theoretical clear-sky air temperature generally overestimates the actual cloud air temperature, because the filling value reflects the air temperature under clear-sky conditions, and the cooling effect of the cloud layer is not considered.

[0114] (3) Based on the quality control (QC) band in the MODIS land surface temperature product, the clear-sky observation records of the weather station are identified. For the dates with invalid observations due to cloud cover, the same kernel function-based time dimension filling method is used to interpolate the cloud air temperature at the station scale, and the theoretical clear-sky air temperature on that date is obtained. The difference between the estimated theoretical clear-sky air temperature and the measured cloudy or overcast air temperature at the same date is calculated to obtain the air temperature deviation term (ΔTa) at the station scale.

[0115] (4) Based on the longitude and latitude of the weather stations and the date, the spatio-temporal matching is performed to extract the environmental variable information corresponding to the station temperature deviation term from the spatial feature dataset. According to the geographical location of the station, the daily cloud continuous cover days (CCD) of the corresponding pixel are extracted, which represents the cloud cover frequency in a period of time; the daily cloud base height (CBH) is extracted, which reflects the vertical distribution characteristics of the cloud layer in the atmosphere; the monthly scale synthesized normalized vegetation index (NDVI) and normalized water body index (MNDWI) are extracted, which reflect the surface coverage characteristics; the surface albedo (Albedo) is obtained, which indicates the surface energy reflection ability; the daily surface downward radiation (DSR) is extracted, which represents the available solar radiation; the elevation information (Elevation) and comprehensive terrain index (CTI) are extracted, which represent the terrain characteristics; finally, the annual cumulative day (DOY) is calculated according to the observation date, which represents the seasonal variation characteristics of the temperature.

[0116] (5) The XGBoost model is inputted with the station scale temperature deviation term as the dependent variable and the eleven environmental characteristic variables (CCD, CBH, NDVI, MNDWI, Albedo, DSR, Elevation, CTI, Longitude, Latitude, DOY) extracted above as the independent variables, and the model training is performed (formula (8)). The samples are randomly divided into five subsets, which are alternately used as the training set and the validation set. The optimal hyperparameter combination is selected to determine the best deviation correction model by minimizing the mean absolute error (MAE).

[0117] (6) The optimal XGBoost deviation correction model obtained by training is applied to the spatial feature variables in the study area to estimate the temperature deviation term pixel by pixel. Then, the deviation term is superimposed with the theoretical clear-sky air temperature obtained by time dimension filling to obtain the corrected air temperature under cloudy conditions (formula (9)). Figure 3 The scatter plot of the corrected air temperature and the station observed air temperature is given, and the sample points are mainly distributed near the 1:1 line, indicating good consistency between the two. The filling accuracy of the corrected air temperature is R² stable to 0.98, the mean absolute error (MAE) is between 0.94°C and 0.98°C, and the root mean square error (RMSE) is between 1.25°C and 1.29°C. It is worth noting that the bias (Bias) is always close to 0°C, indicating that there is no obvious systematic overestimation or underestimation in all years. Compared with the theoretical clear-sky air temperature (as shown in Figure 2 , the corrected air temperature has higher accuracy and smaller deviation. This improvement benefits from the fact that the deviation correction model effectively compensates for the cloud cooling effect that the theoretical clear-sky air temperature cannot reflect, making the reconstruction result have higher accuracy.

[0118] (7) Fuse the corrected air temperature with the clear-sky air temperature to obtain the spatially continuous all-weather remote sensing air temperature. Figure 4 The spatial distribution results of the clear-sky air temperature and the all-weather remote sensing air temperature on the 15th of each month in 2020 are shown. The method of the application can effectively fill all air temperature gaps under different cloud cover conditions, and successfully generate an air temperature field with high spatiotemporal continuity. Compared with the clear-sky air temperature, the all-weather air temperature not only accurately restores the large-scale air temperature distribution pattern in China, but also reflects the fine-scale air temperature changes caused by local topography and surface type in detail.

[0119] Although the preferred embodiments of the application have been described, those skilled in the art can make further changes and modifications to the embodiments once they know the basic inventive concept. Therefore, the appended claims are intended to be interpreted as including all the preferred embodiments and all changes and modifications falling within the scope of the application.

[0120] Obviously, those skilled in the art can make various modifications and variations to the application without departing from the spirit and scope of the application. Thus, if these modifications and variations of the application fall within the scope of the claims of the application and their equivalent technologies, the application also intends to include these modifications and variations.

Claims

1. A remote sensing air temperature reconstruction method based on kernel function time dimension filling and bias correction, characterized in that, The method comprises the following steps: S1, based on the clear-sky air temperature retrieved by thermal infrared remote sensing, a time smoothing method based on kernel function weighting is adopted to comprehensively consider the annual periodicity and short-term dynamic fluctuation of air temperature, and the seasonal air temperature of the same period in history and the effective air temperature of the adjacent date are used to reconstruct the theoretical clear-sky air temperature of the target date by weighted average, so as to fill the time dimension of the missing remote sensing air temperature caused by cloud cover, and obtain the theoretical clear-sky air temperature under the cloud; S2, based on the measured air temperature of the ground meteorological station, the systematic deviation between the theoretical clear-sky air temperature under the cloud and the actually observed cloud air temperature is quantified, and a deviation correction model is constructed combined with multi-source environmental variables to calculate the air temperature deviation caused by cloud disturbance; S3, the deviation correction model is applied to the pixel scale to correct the theoretical clear-sky air temperature under the cloud, and the cloud air temperature estimation value closer to the actual cloud condition is obtained; S4, the clear-sky air temperature retrieved by thermal infrared remote sensing and the corrected cloud air temperature estimation value are fused to generate continuous and seamless all-weather remote sensing air temperature data; In step S1, the following formula is used to reconstruct the theoretical clear-sky air temperature of the target date: where T theo (i) is the theoretical clear sky air temperature for the target day i, T j is the effective remote sensing air temperature observation for the reference day j, S is the set of available reference days, w ij is the weight assigned to the reference day j using a kernel function weighting strategy; w ij = a * K seasonal (DOY ij ) + (1 - a) * K local (time ij ); where K seasonal (DOY ij ) is a seasonal kernel function for reflecting the annual periodic variation characteristics of air temperature, DOY ij is the circular day order difference between the target day and the reference day in a year; K local (time ij ) is a local kernel function, which captures short-term temporal similarity by evaluating the absolute time distance between the target date and the reference date, time ij denotes the absolute time distance between the target day and the reference day; a is a dynamic adjustment coefficient, which is based on a dynamic weighting mechanism of the density of valid observations and the intensity of short-term changes to adaptively adjust the relative contribution of the seasonal kernel and the local kernel in the weight. The seasonal kernel function is: where DOY ij = min(|DOY i - DOY j |, 365 - |DOY i - DOY j |) represents the circular day order difference between the target day and the reference day, σ s is the bandwidth parameter of the seasonal kernel; DOY i represents the year-accumulated day of the target day i, DOY j represents the year-accumulated day of the reference day j.

2. The remote sensing air temperature reconstruction method based on kernel function time dimension filling and bias correction according to claim 1, characterized in that, The method further comprises: Obtaining remote sensing air temperature related data, including remote sensing data, reanalysis data and ground meteorological station observation data; wherein the remote sensing data includes remote sensing clear-sky air temperature data set, MODIS land surface temperature data, MODIS reflectivity data and MODIS albedo data, and SRTM digital elevation model data; the reanalysis data includes ERA5 data and derived ERA5-Land data set; Extracting the clear-sky air temperature retrieved by thermal infrared remote sensing in the research period range from the remote sensing clear-sky air temperature data set; Extracting the daily scale land surface temperature data and its quality control band from the MODIS land surface temperature data, identifying the cloud coverage pixels by analyzing the QC information, and then calculating the daily cloud continuous coverage days CCD of each pixel; Extracting the red band, green band, near-infrared band and short-wave infrared band from the MODIS reflectivity data for calculating the normalized vegetation index NDVI and the normalized water body index MNDWI; obtaining the spatially complete normalized vegetation index NDVI and the normalized water body index MNDWI by using the monthly scale synthesis method; Extracting the surface albedo data Albedo from the MODIS albedo data to reduce the influence of cloud cover by monthly synthesis; Extracting the elevation information Elevation from the SRTM digital elevation model data, and calculating the comprehensive terrain index CTI; all the remote sensing data are projected by Albers equal-area projection, and resampled to 1km spatial resolution by using the bilinear interpolation method to maintain the spatial consistency between the data sources; Extracting the hourly cloud base height data from the ERA5 data set, and obtaining the daily cloud base height CBH after daily averaging; Extracting the daily surface downward radiation data DSR from the ERA5-Land data set; the daily cloud base height and the daily surface downward radiation data are projected by Albers equal-area projection to resample to 1km spatial resolution to maintain consistency with the remote sensing data; Daily station average temperature is extracted from ground meteorological station observation data, and the extracted station average temperature data is cleaned, de-redundant and quality detected. 3.The remote sensing air temperature reconstruction method based on kernel function time dimension filling and bias correction according to claim 1, characterized in that, The local kernel function is: where time ij is the absolute time distance between the target day and the reference day, σ t is the bandwidth parameter of the local kernel; time i and time j denote the absolute time of the target day i and the reference day j, respectively.

4. The remote sensing air temperature reconstruction method based on kernel function time dimension filling and bias correction according to claim 1, characterized in that, The dynamic adjustment coefficient alpha is: where p is the density of valid observations, is the intensity of short-term variations, N valid denotes the number of valid observations within a symmetric time window of size w, σ local is the standard deviation of near-surface air temperature within the symmetric time window.

5. The remote sensing air temperature reconstruction method based on kernel function time dimension filling and bias correction according to claim 1, characterized in that, Step S2 further includes: From the four angles of cloud, solar radiation, surface characteristics and time characteristics, the multi-source environmental variables corresponding to the research area are screened, including: day cloud continuous coverage days CCD and day cloud bottom height CBH for reflecting cloud related characteristics; ground surface downward radiation DSR for representing solar radiation characteristics; Normalized vegetation index NDVI, normalized water body index MNDWI, surface albedo Albedo, elevation information Elevation, longitude Lon, latitude Lat and comprehensive terrain index CTI for reflecting surface characteristics; Annual accumulated day DOY for reflecting time characteristics; An XGBoost algorithm is used for regression modeling to build a bias correction model: ΔT a_s = f(V CCD , V CBH , V NDVI , V albedo , V DSR , V MNDWI , V Ele , V CTI , V Lon , V Lat , V DOY ) where ΔT a_s is the difference between the theoretical clear-sky temperature and the actual under-cloud temperature at each site, and f is the nonlinear relationship between the systematic bias of the site and the multi-source environmental variables corresponding to the selected study area. V CCD ,V CBH ,V NDVI , V albedo ,V DSR ,V MNDWI ,V Ele ,V CTI , V Lon , V Lat ,V DOY are the values of the daily cloud cover (CCD), the daily cloud base height (CBH), the normalized difference vegetation index (NDVI), the albedo (Albedo), the downward surface radiation (DSR), the modified normalized difference water index (MNDWI), the elevation (Elevation), the compound topographic index (CTI), the longitude (Lon), the latitude (Lat), and the day of year (DOY), respectively. With the minimum mean absolute error as the target, the model parameters of the bias correction model are optimized through random search and five-fold cross-validation; in the cross-validation process, the ground meteorological stations are randomly divided into several subsets, and for each iteration, part of the subsets are used for model training, and the remaining subsets are used as the test set for model validation; repeated iterations are performed to ensure that each subset is used as a validation set; Various evaluation indexes including determination coefficient R 2 , mean absolute error MAE, root mean square error RMSE and bias Bias are adopted to comprehensively consider the model performance of the bias correction model, and the model with the best performance is taken as the final bias correction model.

6. The remote sensing air temperature reconstruction method based on kernel function time dimension filling and bias correction according to claim 1, characterized in that, Step S3 further includes: The established bias correction model is applied to the spatialized variable characteristics, and the temperature bias term of all cloud-free pixels is estimated at the pixel level; the estimated temperature bias term is superimposed on the theoretical clear sky temperature under the cloud to obtain the corrected temperature estimation value under the cloud: T a = T theo + ΔT a_p ; In the formula, T a is the estimated value of the cloud-free air temperature after cloud pixel correction, T theo is the theoretical clear-sky air temperature under the cloud after reconstruction in the time dimension, ΔT a_p is the estimated air temperature deviation term at the pixel scale.

Images