A high-resolution spatialization method for near-surface air temperature in glacierized regions

Abstract

This invention discloses a high-resolution spatialization method for near-surface air temperature in glacial regions. This method integrates measured data from high-altitude meteorological stations, satellite remote sensing data, and machine learning-optimized spatial interpolation techniques to construct a continuous temperature field under complex terrain conditions. First, by fusing station-measured and satellite remote sensing data, a data complementarity method suitable for sparsely observed areas is established. Based on a machine learning optimization algorithm, spatialized glacial temperature records from meteorological station observations are used as training data, while satellite remote sensing data with broader coverage is employed. Surface temperature retrieved through a radiative transfer model is used as input features to correct the surface temperature of the study area, generating a spatial distribution of air temperature that better matches ground observations. Finally, by fusing altitude and vertical temperature lapse rate as core constraint variables and combining benchmark verification based on station observation data, the accuracy of near-surface air temperature spatialization calculations in high-altitude complex terrain areas is improved.

Description

Technical Field

[0001] This invention relates to the field of near-surface air temperature multi-source data fusion modeling technology, and in particular to a high-resolution spatialization method for near-surface air temperature in glacial regions. Background Technology

[0002] Against the backdrop of global warming, the accelerated melting of glaciers makes the accurate acquisition of near-surface air temperature and its spatial distribution data particularly important. This data is crucial not only for analyzing glacier melting dynamics and predicting future cryosphere changes, but also for understanding the Asian monsoon system and studying global climate change. As a key input parameter for hydrological and glacier models, obtaining near-surface air temperature data faces severe challenges in high-altitude regions such as the Tibetan Plateau. The sparse distribution of meteorological observation stations in these areas, coupled with the harsh natural environment, results in a scarcity of in-situ observational data on glacier temperatures, rendering traditional air temperature monitoring methods inadequate for practical needs.

[0003] Current temperature monitoring technologies face several challenges. First, while ground-based observations provide highly accurate data, the extreme environmental conditions at high altitudes make station deployment extremely difficult, resulting in insufficient spatial representativeness of the data. For high-altitude areas with complex terrain, spatializing near-surface temperature data in the absence of observational data makes accuracy difficult to guarantee. Second, while satellite remote sensing technology can achieve large-scale temperature monitoring across geographic areas, it has numerous limitations. For example, the temporal resolution of Landsat series satellite surface temperature products is constrained by satellite orbital parameters, resulting in long revisit periods. Furthermore, data loss due to cloud cover and other factors exists. In addition, the unique radiation characteristics of highly reflective surfaces such as ice and snow, as well as the shading effect caused by mountainous terrain, introduce significant uncertainties into temperature retrieval.

[0004] Currently, although some studies have attempted to combine ground observations, satellite remote sensing, and numerical simulations to obtain temperature data, these multi-source data fusion methods still have significant shortcomings. On the one hand, satellite-based surface temperature retrieval is limited to clear-sky conditions, and effective data cannot be obtained under cloudy weather. On the other hand, existing methods often ignore the dynamic changes in the vertical temperature lapse rate over time and space at high altitudes. Particularly noteworthy is that, within the framework of machine learning algorithms, few studies have been able to systematically utilize local temperature data with clearly defined geographic coordinates as constraints to optimize temperature retrieval results based on remote sensing data. This problem is particularly prominent in glacial regions, as most existing frameworks fail to adequately consider the unique topographic and climatic characteristics of glacial areas.

[0005] Therefore, existing technologies still need to be improved and enhanced. Summary of the Invention

[0006] In view of the shortcomings of the prior art, the purpose of this invention is to provide a high-resolution spatialization method for near-surface air temperature in glacial regions, which aims to solve the problems of traditional techniques relying on sparse station interpolation or single satellite data, resulting in significant errors in complex terrain, and lacking dynamic constraints on the vertical temperature lapse rate, leading to low spatialization accuracy of air temperature.

[0007] To achieve the above objectives, the present invention adopts the following technical solution: A high-resolution spatialization method for near-surface air temperature in glacial regions includes: Step 1: Collect near-surface temperature observation data and corresponding station coordinates provided by meteorological observation stations in the glacier area. Establish a near-surface temperature data sample set in the glacier area according to the time series. Mark the values ​​of 0 and invalid data in the sample set as missing values ​​and retain high-quality data. Step 2: Based on the near-surface temperature data sample set of the glacier region, establish a preliminary temperature distribution field. Use the nearest neighbor interpolation method to assign values ​​to the grid points of the preliminary temperature distribution field using high-quality data, and output the preliminary spatial distribution of temperature. Step 3: Use triangulation-based interpolation to reconstruct the preliminary temperature distribution field into a high-resolution temperature distribution field for the missing value data; Step 4: Based on the near-surface air temperature observation data provided by the meteorological observation station and the satellite remote sensing surface temperature data, perform time synchronization and spatial alignment to refine the high-resolution temperature distribution field and obtain a spatialized model of near-surface air temperature in the glacier area; Step 5: Integrate the near-surface temperature data measured by meteorological stations with the satellite remote sensing surface temperature data that has been calibrated with latitude and longitude coordinates and matched with band conversion through the ground-satellite collaborative framework, and build a benchmark training set in combination with the digital elevation model to train the spatialization model of near-surface temperature in the glacier area. Step 6: Correct the error of the trained spatial model of near-surface air temperature in the glacier area by calculating the vertical temperature lapse rate and using the inverse distance weighted interpolation method for elevation variables, and generate the spatial distribution field of near-surface air temperature by combining the air temperature data of the glacier area.

[0008] Furthermore, in step 1, the Isolation Forest algorithm is used to mark values ​​of 0 and invalid data in the sample set as missing values. The Isolation Forest algorithm then extracts data from the data sample set without replacement. f A default-sized subsample is used, and an isolation tree is trained for each subsample. For each tree, features and splitting values ​​are randomly selected. Features include time, near-surface air temperature observation data, and observation point coordinates. The data space is recursively partitioned until a data point is isolated. The maximum depth of the isolation tree is set to [value missing]. The rounded-up value is used to avoid overfitting. A tree is constructed by randomly splitting the data. Outliers, because they deviate from the population, can be isolated to shallow nodes of the tree with only a few partitions. Anomaly scores are calculated. ; in h This represents the path length, which is the number of edges from the root node to the isolated data point. E ( h ( x )) is a point x Average path length across all trees f Number of subsamples ， c( f The normalization factor for path length is 0 to 1. The standardized score range is 0 to 1. The threshold is dynamically adjusted in combination with the glacial climate background. If the abnormal score is higher than the threshold, it is judged as abnormal.

[0009] Furthermore, in step 1, the threshold is dynamically adjusted in combination with the glacier climate background. If the abnormal score is higher than the threshold, it is judged as abnormal. Specifically, for spatial anomalies, when the temperature of a certain observation point in the glacier observation station area deviates significantly from that of the surrounding stations, it is marked as abnormal even if the score does not exceed the threshold, and missing values ​​are marked. For glacier areas with extreme climate disturbances and significant diurnal temperature variations, adaptive optimization is performed by combining multi-source data. For cold wave scenarios, spatial verification is performed by introducing remote sensing reanalysis data. If the low temperature is consistent over a large area, it is determined to be a real event and the data is retained. For stations with large diurnal temperature variations, sliding window standardization is adopted, and the mean and standard deviation are calculated with a 7-day window. The current value is marked as an anomaly when it exceeds the mean ± 3 times the standard deviation, so as to avoid misjudging normal diurnal fluctuations.

[0010] Furthermore, in step 2, the coverage area of ​​the weather stations is divided into grids. Based on the nearest neighbor algorithm, the observation data of each weather station is used as the initial near-surface temperature distribution estimate for the corresponding grid. The temperature data used comes from automatic weather stations and temperature recorders in different watersheds. The coordinates of the weather stations within the target coverage area are organized into a set in logical order: X (station) ={1,…, X j ,…, X n},in X (station) Represents the set of coordinates of weather stations. X j Represents the coordinates of a single weather station. n Indicates the total number of stations; The coverage area of ​​the weather station is divided into grids: X (grid) ={1,…, X i ,…, Xm},in X (grid) Represents the set of grid coordinates. X i Represents the center coordinates of a single grid cell. m Indicates the number of grid points; Assign a unique index to each site X (station) ={1,…, X j ,…, X n} is used for subsequent matching, which calculates the distance between grid points and stations based on the meteorological station index number and performs nearest neighbor matching; Based on nearest neighbor learning, a brute-force search is used to calculate the distance between each station and each grid center point: ; L This indicates that Euclidean distance is used as the metric, which is the straight-line distance between two points in the standard coordinate system. When the station and the grid have the same index, the near-surface temperature value of the station is directly assigned to the corresponding grid. Within the coverage area of ​​the meteorological station, the observation values ​​will be distributed to each neighboring grid to form an initial temperature distribution field. The meteorological observation data is used to assign values ​​to the surrounding station areas to ensure that the temperature data of each grid comes from the actual observation values.

[0011] Furthermore, in step 3, based on the missing value data from step 1, the mean and standard deviation of the data sample set are calculated, using a 3x3 matrix. s The criteria include reassigning extreme temperature observations exceeding ±3 standard deviations, cleaning and standardizing the data, and recording the spatial distribution characteristics of all outliers, including altitude and latitude / longitude information. Then, a Voronoi diagram is constructed based on station coordinates, using the stations as the center points of the map elements to achieve spatial partitioning. The actual coordinates of the stations are the centers of the Voronoi units. Furthermore, a Delaunay triangulation is constructed using the intersections of Voronoi units and station coordinates, ensuring full coverage of the study area. An irregular triangulation is constructed, and the initial near-surface air temperature distribution estimation results are iteratively corrected through triangulation. The residuals between the interpolated grid node temperatures and the measured data are calculated. Nodes with residuals greater than 2°C undergo secondary interpolation, iterative residual optimization is performed to improve local accuracy, and adaptive threshold optimization is used to ensure that the reconstructed temperature field conforms to the actual climate distribution patterns.

[0012] Furthermore, in the interpolation process of step 3, to avoid errors in single-station observation data leading to overestimation or underestimation of surrounding areas, a smoothing function needs to be applied to the interpolated near-surface air temperature distribution within the station network coverage area: A 5×5 sliding window is used to calculate the median of the temperature values ​​within the window and replace the center pixel value to effectively remove noise and outliers. The window size is dynamically adjusted, and when the proportion of valid data points within the window is detected to be lower than the threshold, the window is expanded to ensure sufficient samples. When performing median filtering, special attention is paid to key areas such as the glacier-bare rock transition zone and the moraine region to avoid excessive smoothing that could lead to the loss of temperature transition information. Filter intensity thresholds are set separately for glacier surfaces and non-glacier surfaces to ensure that small-scale temperature gradients are preserved. The local gradient is calculated for the filtered temperature field. If the temperature difference between adjacent grid cells exceeds a threshold, the standard deviation of the Gaussian kernel in that region is reduced. s Values ​​are optimized to preserve the characteristics of climate abrupt changes and to achieve spatial continuity. Cross-validation is used to calculate the root mean square error and standard deviation to ensure that the interpolated temperature field conforms to the observation trend. For areas with large residuals, corrections are made by combining satellite remote sensing surface temperature data to adjust the near-surface air temperature distribution and output a high-resolution continuous near-surface temperature field.

[0013] Furthermore, in step 4, time synchronization and spatial alignment are performed based on near-surface air temperature observation data provided by meteorological observation stations and satellite remote sensing surface temperature data, including: Step 41: Temperature data collected by meteorological observation stations are usually recorded at fixed time periods, with a certain time interval between adjacent observation times. When the satellite's transit time falls between two consecutive observation times at a meteorological station, the station observation times in the remote sensing data fusion are interpolated. Based on all available station observation data, a third-order uniform B-spline curve is fitted to interpolate the station temperatures to the precise remote sensing imaging timestamp, ensuring that the temperature data is interpolated in time and matched with the satellite center's transit time. This process connects two adjacent observation times within each station. Meteorological records for different time periods are divided into subdivisions. n Sub-intervals The step size of each subinterval can be expressed as: ; Where, Δ t This indicates the time step, which is the duration of the interval between two time points. This is the end time of the period. The start time of the time period. n This indicates the number of sub-intervals after the subdivision. to Any time between can be represented as , i = 0, 1,2, ..., n ; Step 42: Use a third-order uniform B-spline curve to... and Interpolation is performed on the temperature observation records corresponding to each time point. Each basis function is a cubic polynomial, and the spline curves are uniformly distributed in the parameter node space. Each control point corresponds to one basis function. For two time points... arrive Observation data points between Each of them represent The temperature observation value at time, the cubic spline interpolation function f( t The following conditions must be met: ; in, The start time of the time period. n This indicates the number of sub-intervals after further subdivision; to Any time between is represented as , i = 0, 1, 2, ..., n Each paragraph All are cubic polynomials calculated according to the following rules: ; For B-spline representation, the interpolation function is expressed as: ; Where, β j ( t ) represents the cubic B-spline basis functions. α j To be based on data constraints f i ( t i )= y i The calculated weighting coefficients are used as the f( of the cubic polynomial). t The cubic spline interpolation function is obtained by multiplying the basis functions and the control point weights and then summing the results to obtain the expression of the interpolation curve. It can obtain the temperature value at any satellite center time, achieve time synchronization and spatial alignment, improve the high-resolution temperature distribution field, and obtain a spatialized model of near-surface temperature in glacier areas.

[0014] Step 43: When the image time of the remote sensing satellite is obtained, but no meteorological observation data is recorded at that time, the meteorological observation information at the next time point is predicted. Then, the temperature is interpolated between the previous time point and the current time point to locate and align with the remote sensing image time to a specific time point. The temperature interpolation result at this precise time point is used to represent the near-surface temperature distribution of the meteorological station. The near-surface temperature at the next time point of each meteorological station is predicted using support vector regression. Finally, the optimized model is obtained by minimizing the total loss and maximizing the interval. Support Vector Regression (SVR) builds a regression model by finding the optimal hyperplane, minimizing the error between the predicted values ​​and the observed values ​​of the training samples. For nonlinear models, SVR uses a kernel function to map to the feature space before performing regression. The loss function of SVR is expressed as: ; in, y This represents the temperature observation value recorded by the weather station. or The predicted value of near-surface air temperature is within the tolerance range of regression error. If the deviation is within the range of ε, it will not be considered an error. y , or ) is the loss function, ignoring minor inaccuracies in its calculation process.

[0015] Furthermore, in step 4, in order to effectively analyze and utilize the residuals of support vector regression in predicting air temperature and their impact on temperature prediction, a diverse dataset is used for model training, and measures are set to avoid overfitting.

[0016] Furthermore, in step 5, the high-resolution near-surface air temperature spatialization data generated in steps 1-4, the surface temperature retrieved from remote sensing features after radiometric calibration and atmospheric correction, and the digital elevation model are used as variables to prepare input data. Resampling is used to ensure that the elevation data and the spatial resolution of the temperature spatialization grid are consistent. All data were resampled to a uniform grid with the same resolution as the temperature field, maintaining spatial alignment, and the feature matrix was... X The grid points contain features including the spatialized near-surface air temperature distribution of glaciers, elevation variables, surface temperature retrieved from remote sensing data, and target variables. y Interpolated temperature values ​​corresponding to grid points; The training and test sets were divided to ensure uniform spatial distribution and avoid regional bias. Five-fold cross-validation was used to optimize the model hyperparameters. Random forest regression was selected as the candidate model, and decision tree regression was used as the baseline comparison model. Multilayer perceptron was used to verify the ability to capture nonlinear relationships. During the training of the random forest model, each tree is trained independently based on randomly selected subsamples. Hyperparameter optimization adopts a grid search method, and the parameters adjusted include the number of trees, the maximum depth of the trees, and the minimum number of samples required for node splitting. The model performance is evaluated by the mean squared error to assess the magnitude of the prediction bias, and the coefficient of determination is used to measure the model's explanatory power for temperature changes. Through feature importance analysis, the influence of variables such as surface temperature and digital elevation model on temperature prediction is determined.

[0017] Furthermore, in step 6, the regional temperature vertical lapse rate is calculated using meteorological station data to correct for the glacier area. When a station is located in the glacier area, a local temperature vertical lapse rate is used to correct each target grid point based on its altitude and the temperature of the reference station, performing temperature field altitude correction. If multiple stations affect the target point, the adjusted temperature is combined by weighted average of multiple stations. The temperature of the target point is calculated by weighting the observation data of surrounding meteorological stations using the inverse distance weighted interpolation method. In areas with sparse data, the search range is expanded or virtual stations are added. Combined with the adjustment rate of temperature change with altitude, a spatial distribution field of near-surface temperature is generated. Inverse distance weighted interpolation is performed on all grid points, and then altitude correction is performed using the temperature vertical lapse rate. In the glacier area, a local lapse rate is used for further optimization. Finally, a cross-validation method is used, with some station data as the validation set, to calculate the root mean square error and coefficient of determination before and after correction. At the same time, the residuals are analyzed to verify the error improvement in high-altitude areas.

[0018] The technical solution adopted in this invention has the following beneficial effects: The near-surface air temperature multi-source data fusion method provided by this invention has significant innovative advantages and application effects. This method innovatively integrates multi-source data, including ground station observations, satellite remote sensing data, and digital elevation models, overcoming the limitations of single data sources, specifically addressing the unique environmental conditions of high-altitude glacier regions. By establishing an optimized vertical temperature lapse rate correction algorithm for glaciers, it corrects the spatial heterogeneity bias of near-surface air temperature under complex terrain conditions. Furthermore, the proposed elevation-constrained inverse distance weighted interpolation technique adjusts the weight allocation for different regions, reducing the errors generated by single interpolation methods in areas with complex terrain.

[0019] In terms of practical effectiveness, this method improves the accuracy of spatializing near-surface air temperature in high-altitude glacier regions. Test data shows that after processing with this method, the root mean square error of air temperature data is reduced by an average of 3.2 degrees Celsius, and the standard deviation is reduced by an average of 1.8 degrees Celsius, significantly outperforming single linear regression algorithms. During glacier surface temperature correction, by integrating near-surface observation data, elevation information, and a vertical lapse rate correction model, unreasonable temperature abrupt changes between glacier and non-glacier areas are effectively reduced, while maintaining the continuity of air temperature spatial distribution. Furthermore, the generated high-resolution air temperature products can take into account the influence of the unique topography of glacier regions on air temperature distribution, providing data support for glacier change research.

[0020] This method also exhibits good applicability and scalability. Adaptive optimization can be achieved by dynamically adjusting the vertical temperature lapse rate parameter to suit the data characteristics of different regions and time periods. Systematic bias correction based on machine learning techniques further enhances the model's stability under various climatic conditions. This method not only meets the needs of scientific research but also provides reliable data for practical applications such as site selection for cold-region engineering projects and cryosphere hydrological simulation. Because it significantly reduces reliance on ground station data, this method is particularly suitable for temperature monitoring in high-altitude glacier regions such as the Qinghai-Tibet Plateau, providing new model data support and technical means for glacier dynamics research in the context of global climate change. Attached Figure Description

[0021] Figure 1 A flowchart illustrating a high-resolution spatialization method for near-surface air temperature in glacial regions provided by this invention; Figure 2 This is a spatial distribution map of meteorological stations and Voronoi units in an embodiment of the present invention; Figure 3 This is a flowchart illustrating the data fusion process in an embodiment of the present invention. Detailed Implementation

[0022] The present invention will be further explained below with reference to specific implementation schemes, but is not limited to the present invention. The present invention provides a method for spatializing near-surface temperature of glaciers, such as... Figure 1 The spatialization method of near-surface temperature of glaciers shown and Figure 3 The data fusion process shown includes the following steps: S1: Data Preparation and Preprocessing S1 specifically includes the following steps: S11: Data Sources, Data Standardization, and Missing Value Labeling Collect near-surface air temperature data from meteorological stations in the glacier area ( T The data and corresponding station coordinates (longitude, latitude, altitude) are used to form a time series sample set, and the land surface temperature (LST) retrieved from satellite remote sensing data is integrated as an auxiliary feature. Temperature data is standardized, and values ​​of 0 or invalid data are marked as missing values ​​for subsequent imputation. Meteorological station coordinates (spatial features) and timestamps (temporal features) are merged to construct a multidimensional input matrix X=[time step, temperature value, longitude, latitude, altitude].

[0023] S12: Anomaly Detection Based on Isolation Forest The Isolation Forest algorithm is used to construct an iTree set of isolated trees through random sampling and recursive space partitioning for model training. The Isolation Forest algorithm is used to mark 0 and invalid data in the sample set as missing values. The Isolation Forest algorithm then extracts data from the data sample set without replacement. f A default-sized subsample is used, and an isolation tree is trained for each subsample. For each tree, features and splitting values ​​are randomly selected. Features include time, near-surface air temperature observation data, and observation point coordinates. The data space is recursively partitioned until a data point is isolated. The maximum depth of the isolation tree is set to [value missing]. The rounded-up value is used to avoid overfitting. A tree is constructed by randomly splitting the data. Outliers, because they deviate from the population, can be isolated to shallow nodes of the tree with only a few partitions. Anomaly scores are calculated. ; in h This represents the path length, which is the number of edges from the root node to the isolated data point. E ( h ( x )) is a point x Average path length across all trees f Number of subsamples ， c( f The normalization factor for path length is 0 to 1. The standardized score range is 0 to 1. The threshold is dynamically adjusted in combination with the glacial climate background. If the abnormal score is higher than the threshold, it is judged as abnormal.

[0024] Time step (observation time series), temperature value ( T The spatial coordinates (longitude, latitude, altitude) are set as features.

[0025] An anomaly score is generated by calculating the path length of each sample. The score ranges from [0,1]. The closer the score is to 1, the higher the probability of an anomaly (such as extreme high / low temperature or sudden drop in temperature events). Normal data points have scores close to 0. This method can effectively identify anomalous low temperature abrupt events during the summer glacier melting period.

[0026] By setting an outlier threshold using quantile methods (such as the 95th percentile) or manual verification (such as a score > 0.6 being considered outlier), outliers are automatically marked, thresholds are set, and data filtering is implemented to retain high-quality data for subsequent data entry.

[0027] In response to spatial anomalies, if the temperature at a certain observation point in the glacier observation station area deviates significantly from that of the surrounding stations, it is marked as an anomaly and missing values ​​are marked. For glacier areas with extreme climate disturbances and significant diurnal temperature variations, adaptive optimization is performed by combining multi-source data. For cold wave scenarios, spatial verification is performed by introducing remote sensing reanalysis data. If the low temperature is consistent over a large area, it is determined to be a real event and the data is retained. For stations with large diurnal temperature variations, sliding window standardization is adopted, and the mean and standard deviation are calculated with a 7-day window. The current value is marked as an anomaly when it exceeds the mean ± 3 times the standard deviation, so as to avoid misjudging normal diurnal fluctuations.

[0028] S2 is constructed based on the preliminary temperature distribution field interpolated from the nearest neighbor (NN): S21: Meteorological station data processing and target area grid division Input data includes station coordinates (latitude and longitude or projected coordinate system, such as WGS84 / UTM) and near-surface temperature observations; Data processing includes filling in missing data (such as assigning values ​​to sensor fault records) and ensuring that coordinate units are consistent (such as unifying them to the coordinate unit of meters or decimal degrees). The grid parameter definition includes generating the grid resolution (e.g., 30m) and determining the grid range (determined by the boundary of the target region, such as the glacier region of the Qinghai-Tibet Plateau). Generate grid point coordinate matrix according to resolution X (grid) ={1,…, X i ,…, X m The initial distribution of} is calculated, including the site index and minimum interval calculation (brute force search), the nearest neighbor distance between sites is calculated, the minimum distance between each site is determined to evaluate the interpolation reliability, the nearest neighbor distance Min_Distance of each site is output, and the site density insufficient area is identified (if Min_Distance>threshold, a warning should be issued and processed).

[0029] That is, the coverage area of ​​the weather stations is divided into grids. Based on the nearest neighbor algorithm, the observation data of each weather station is used as the initial near-surface temperature distribution estimate for the corresponding grid. The temperature data used comes from automatic weather stations and temperature recorders in different watersheds. The coordinates of the weather stations within the target coverage area are arranged into a set in logical order. X (station) ={1,…, X j ,…, X n},in X (station) Represents the set of coordinates of weather stations. X j Represents the coordinates of a single weather station. nIndicates the total number of stations; The coverage area of ​​the weather station is divided into grids: X (grid) ={1,…, X i ,…, X m},in X (grid) Represents the set of grid coordinates. X i Represents the center coordinates of a single grid cell. m This indicates the number of grid points.

[0030] Assign a unique index to each site X (station) ={1,…, X j ,…, X n} is used for subsequent matching, which calculates the distance between grid points and stations based on the meteorological station index number and performs nearest neighbor matching; Calculate the Euclidean distance from the grid point to all stations, assign the temperature observation value of the nearest station to the corresponding grid point, output the preliminary spatial distribution of temperature, the data structure is a regular grid temperature matrix, and draw a near-surface temperature distribution map to verify whether it reflects the spatial distribution characteristics of the stations, and mark the grids without station coverage (such as the Distance>threshold area). Based on nearest neighbor learning, a brute-force search is used to calculate the distance between each station and each grid center point: ; L This indicates that Euclidean distance is used as the metric, which is the straight-line distance between two points in the standard coordinate system. When the station and the grid have the same index, the near-surface temperature value of the station is directly assigned to the corresponding grid. Within the coverage area of ​​the meteorological station, the observation values ​​will be distributed to each neighboring grid to form an initial temperature distribution field. The meteorological observation data is used to assign values ​​to the surrounding station areas to ensure that the temperature data of each grid comes from the actual observation values.

[0031] During the spatialization of temperature data from weather stations, there may be a problem of no station coverage on the edge grid. This can be addressed by setting an upper limit for the search radius or explicitly marking empty grid points. If uneven station density leads to temperature abrupt changes in sparse areas, inverse distance weighting or Kriging interpolation can be introduced for smoothing optimization. When the representativeness of high-altitude stations is insufficient, a single station is difficult to reflect temperature changes under complex terrain. Digital elevation model (DEM) data can be combined with interpolation through methods such as temperature-elevation correction to improve the rationality of spatial distribution.

[0032] S22: Spatialization process of near-surface air temperature based on nearest neighbor (NN) interpolation Collect observational data from meteorological stations, including station coordinates (latitude, longitude, and projected coordinates) and temperature values. Check and remove missing or outlier values ​​to ensure data quality. Standardize the coordinate system (such as WGS84 latitude and longitude or UTM projection) to avoid subsequent calculation errors. Set the spatial range (minimum and maximum latitude, longitude, or projected coordinates) and determine the grid resolution, which can be adjusted according to research needs. Calculate the Euclidean distance between each weather station and all other stations, record the distance and index number of the nearest neighbor station, estimate the minimum distance between stations, statistically analyze the distribution of the minimum distance between stations across the entire region, identify sparse areas of stations (such as cases where the distance between stations is too large), mark areas with insufficient stations, and subsequently consider different interpolation strategies or supplementary data to evaluate the spatial distribution using the nearest neighbor distance calculation. According to the defined resolution, the target area is divided into regular grids, and the coordinates of the center point of each grid are calculated. For each grid point, its Euclidean distance to all weather stations is calculated, the nearest station is found, the nearest neighbor is assigned a temperature value, and the temperature value of the station is directly assigned to the corresponding grid. Regular grid division and nearest neighbor matching generate grid center points. If multiple grid points match the same station, the same observation value is used. If a grid point has no valid station coverage (distance exceeds the threshold), it can be marked as invalid data (NaN). For densely populated areas, the grid values ​​basically reflect the actual measurements at the stations, and no excessive additional smoothing is needed to generate a rasterized temperature field. Each grid cell contains the temperature values ​​of the nearest neighbor stations. Record metadata for information such as interpolation method, grid resolution, and data source; Draw a temperature distribution map for visual inspection, observe the matching between station and grid values, and ensure its rationality; Cross-validation (e.g., interpolating after removing a certain site and then comparing it with the original data) is used to assess error at some sites. Subsequently, inverse distance weighting or Rickin interpolation can be used to smooth abrupt change areas, and terrain correction can be performed by adjusting high-altitude station data in conjunction with elevation data to improve the near-surface temperature distribution in sparsely distributed observation stations.

[0033] S3 Reconstruction of Glacier Temperature Data Based on Staged Quality Control: S31: Interpolation based on triangulation Invalid records in the original meteorological observation data (such as missing values ​​due to sensor failure or communication interruption) are uniformly marked as 0 values ​​(or missing value symbols). The mean and standard deviation of the dataset are calculated, and a 3x3 matrix is ​​used. sThe criteria include reassigning values ​​to extreme temperature observations that exceed ±3 standard deviations, cleaning and standardizing the data, and recording the spatial distribution characteristics of all outliers (such as altitude and latitude / longitude information) to facilitate subsequent analysis of their potential environmental factors.

[0034] After handling anomalies, valid station data is used as grid nodes. A Voronoi diagram is constructed based on station coordinates, making the stations the center points of the primitives to achieve spatial partitioning. The actual coordinates of the stations are the centers of the Voronoi cells. Further, a Delaunay triangulation is constructed using the intersections of the Voronoi cells and the station coordinates, performing Delaunay triangulation to ensure full coverage of the study area. An irregular triangulation (Delaunay Triangulation) is then constructed, such as... Figure 2 As shown; The residuals between the interpolated grid node temperatures and the measured data are calculated. For nodes with residuals >2°C, a second interpolation is performed. Iterative residual optimization is carried out to improve local accuracy. Adaptive threshold optimization is adopted to ensure that the reconstructed temperature field conforms to the actual climate distribution pattern.

[0035] S32: Smoothing of near-surface air temperature distribution To avoid errors in single-station observation data leading to overestimation or underestimation of surrounding area values, a smoothing function needs to be applied to the interpolated near-surface air temperature distribution within the station network coverage area.

[0036] A 5×5 (or adaptive size) sliding window is used to calculate the median of the temperature values ​​within the window and replace the center pixel value to effectively remove noise and outliers. The window size is dynamically adjusted (e.g., 3×3→5×5→7×7). If the proportion of valid data points within the window is found to be lower than the threshold (e.g., <50%), the window is expanded to ensure sufficient samples. When performing median filtering, special attention is paid to key areas such as the glacier-bare rock transition zone and glacial till region to avoid excessive smoothing that would lead to the loss of temperature transition information. Filter intensity thresholds are set for glacier surface and non-glacier surface respectively to ensure that small-scale temperature gradients are preserved and high-frequency noise components are filtered out. Noise control is achieved through phased quality control. Calculate the local gradient for the filtered temperature field. If the temperature difference between adjacent grid cells exceeds a threshold, reduce the standard deviation of the Gaussian kernel in that region. s Values ​​were selected to preserve the characteristics of climate abrupt changes and to optimize spatial continuity, ensuring the accuracy of temperature field reconstruction and the continuity of temperature spatial distribution. Cross-validation was used to calculate the root mean square error and standard deviation to ensure that the interpolated temperature field conforms to the observation trend. For areas with large residuals, corrections were made by combining satellite remote sensing surface temperature data to adjust the near-surface air temperature distribution. Integrate all optimization steps to output a high-resolution continuous near-surface temperature field and record quality control (outlier distribution, error analysis, etc.).

[0037] S4: Strictly align the temperature observations from meteorological stations with the time of remote sensing images to ensure consistent spatiotemporal resolution.

[0038] S41: Station observation time interpolation Temperature data collected by meteorological observation stations are usually recorded at fixed time periods, with a certain time interval between adjacent observation times. When the satellite's transit time is between two consecutive observation times at a meteorological station, the station observation time in the remote sensing data fusion is interpolated. Based on all available station observation data, a third-order uniform B-spline curve is fitted to interpolate the station temperature to the precise remote sensing imaging timestamp, so that the temperature data is interpolated in time and matches the transit time of the satellite center.

[0039] Extracting meteorological stations at adjacent observation times The temperature sequence is evenly divided into n Sub-intervals are divided into time intervals and time steps are defined. , where Δ t This indicates the time step, which is the duration of the interval between two time points. This is the end time of the period. The start time of the time period. n This indicates the number of sub-intervals after the subdivision. to Any time between can be represented as , i = 0, 1, 2, ..., n Ensure that the remote sensing imaging time falls within a certain sub-interval; Based on the site Based on the observed temperature values, fit a third-order B-spline curve f( t The curve passes through all observation points. ,Right now Each paragraph Given a cubic polynomial, ensure smoothness and output the interpolation function f( t It can calculate the temperature value at any satellite center time.

[0040] Using a third-order uniform B-spline curve to and Interpolation is performed on the temperature observation records corresponding to each time point. Each basis function is a cubic polynomial, and the spline curves are uniformly distributed in the parameter node space. Each control point corresponds to one basis function. For two time points... arrive Observation data points between Each of them represent The temperature observation value at time, the cubic spline interpolation function f( t The following conditions must be met: ; in, The start time of the time period. n This indicates the number of sub-intervals after further subdivision; to Any time between is represented as , i = 0, 1, 2, ..., n Each paragraph All are cubic polynomials calculated according to the following rules: ; For B-spline representation, the interpolation function is expressed as: ; Where, β j ( t ) represents the cubic B-spline basis functions. α j To be based on data constraints f i ( t i )= y i The calculated weighting coefficients are used as the f( of the cubic polynomial). t The cubic spline interpolation function is obtained by multiplying the basis functions and the control point weights and then summing the results to obtain the expression of the interpolation curve. It can obtain the temperature value at any satellite center time, achieve time synchronization and spatial alignment, improve the high-resolution temperature distribution field, and obtain a spatialized model of near-surface temperature in glacier areas.

[0041] S42: Remote sensing image time matching If there are no station observation records for the remote sensing center time, the temperature at the time of remote sensing imaging is obtained directly by interpolation. If the remote sensing time is outside the station's observation range, proceed to step S43 to predict the temperature at future times, and then interpolate.

[0042] S43: Missing Time Temperature Prediction When the image time of a remote sensing satellite is acquired, but no meteorological observation data is recorded at that time, the meteorological observation information for the next time point is predicted. Then, the temperature is interpolated between the previous time point and the current time point to locate and align with the remote sensing image time to a specific time point. The temperature interpolation result at this precise time point is used to represent the near-surface temperature distribution of the meteorological station. The near-surface temperature of each meteorological station at the next time point is predicted using support vector regression. Finally, the optimized model is obtained by minimizing the total loss and maximizing the interval. Extracting historical temperature series , From preceding observations ,..., y i-1 Generated by random perturbation , , For the first i The real temperature at a given moment, predicting the next moment. y m+1 ; The hyperplane is trained using the support vector regression algorithm to minimize the prediction error, and the current time step is input. temperature ,predict t m+1 Time value y pred Using the B-spline curve of S42 in Interval interpolation is used to obtain the temperature at the satellite's center at that time.

[0043] Support Vector Regression (SVR) builds a regression model by finding the optimal hyperplane, minimizing the error between the predicted values ​​and the observed values ​​of the training samples. For nonlinear models, SVR uses a kernel function to map to the feature space before performing regression. The loss function of SVR is expressed as: ; in, y This represents the temperature observation value recorded by the weather station. or The predicted value of near-surface air temperature is within the tolerance range of regression error. If the deviation is within the range of ε, it will not be considered an error. y , or ) is the loss function, ignoring minor inaccuracies in its calculation process.

[0044] S44: Model Validation and Error Control The model is computed by dividing the test set into independent sets. The deviation between the predicted and the true values ​​is quantified by the root mean square error (RMSE). The explanatory power of the model is evaluated by the coefficient of determination (R²). The residual distribution is analyzed to identify systematic errors. Cross-validation is used to adjust the hyperparameters, limit the model complexity, and avoid overfitting to noisy data.

[0045] To effectively analyze and utilize the residuals of support vector regression in predicting air temperature and their impact on temperature prediction, diverse datasets were used for model training, and measures were implemented to avoid overfitting, as shown below: ; in, y This represents the temperature data observed by the weather station. m This represents the end point of the observation record. i Iterate from the initial time point 0 to the current time point. m At a certain point in time, m +1 represents the total number of observation points required for model training, while y i Indicates the first i Observation records at each time point; ; ; ; ; in, X It is a matrix of temperature observation data y The generated ( m +1)×( k +1) dimensional matrix, I for( m +1)×( m +1) dimensional identity matrix, matrix X any element in x ij Represented as: ; in, d ij It is a random number between -1 and 1, that is, -1 < d ij <1, m This is the end time point of the observation record. i From the initial time point 0 to the... m The sequence number of each time point. k Representation matrix X The largest column index, j The column index is 0, and its value ranges from 0 to 1. k , k +1 is a matrix X The total number of columns; Divide the training sample set D = {( x 1, y 1), …, ( xi , y i ), …, ( x m , y m The data generation rules are as follows: x 1. From initial observations y 0 generation, specifically satisfying: x 1j = ( d 0j + 1) y 0, where j = 0,1,…, k ; This dataset is based on meteorological records and uses a regression algorithm to generate the training set. X train yes m ×( k +1) dimension matrix is ​​derived from temperature observation matrix y Generate, its elements x ij Corresponding to the 1st to the 2nd m Time point matrix y train Also covering numbers 1 to 2 m The target for predicting the observations at the () time point is the () m The remote sensing image of the temperature value at time point +1) is located at the [number]th time point. m With the ( m +1) Between time points, it is represented as: ; ; ; No. m The training data at each time point was generated through calculation, where x m+1 =(1+ d mj ) y m , d mj For random disturbance terms, y m It is the first m One training data value, y m+1 As the true value of the test data, a near-surface air temperature was established based on a regression algorithm. T The prediction model and its calculation process are as follows: ; Near-surface air temperature at the moment of remote sensing imaging is processed by a cubic spline function. y m+1 Compared with the predicted value y prediction For interpolation calculations, in a model trained based on observation time steps and their corresponding near-surface air temperature records, the data representation used for training is assumed to be: ; ; In the m At a certain point in time, y m It is the first m One training data value, y prediction Yes m The predicted value at time +1; the air temperature at the time of remote sensing imaging was obtained using a cubic spline function. y i and y prediction Interpolation estimation is performed, and the remote sensing image center acquires time points. y m+1 The interpolation results can be used as test data. The model is evaluated on an independent test set, and the predictive performance of different residual subgroups is quantified by the root mean square error and coefficient of determination to verify its stability.

[0046] S45: Data Integration and Spatiotemporal Rasterization The interpolated temperature data from meteorological stations and remote sensing images are combined with a unified coordinate system and grid resolution to ensure that the temperature value of each grid point is aligned with the center of the remote sensing pixel, thus generating a spatiotemporally synchronized raster dataset. The near-surface temperature field has the same timestamp and spatial resolution as the remote sensing image.

[0047] To improve matching accuracy, time synchronization and spatial alignment are required, and resampling is used to ensure that the resolution of the surface temperature data is consistent with the interpolated data from the observation stations. A time-series temperature interpolation method that combines third-order uniform B-spline interpolation and support vector regression is adopted to improve the time matching accuracy between remote sensing data and ground meteorological observations. Based on discrete station temperature observation data, a continuous temperature change model is constructed using third-order uniform B-spline curves. The temperature data is accurately interpolated to the transit time of the remote sensing image. For cases where the remote sensing imaging time exceeds the most recent observation time, a support vector regression model is used to predict the temperature at the next moment. Combined with historical data, the temperature value at the satellite transit time is generated through spline interpolation, achieving precise spatiotemporal alignment. By enhancing the diversity of training samples through data perturbation and verifying model stability based on independent test sets and quantitative indicators, the predictive ability to resist overfitting is ensured. This step combines the advantages of physical interpolation and machine learning, effectively supporting the integrated application of remote sensing and meteorological observation.

[0048] S5: Construction of Spatialized Training Set for Glacier Temperature Based on Multi-Source Data and Model Training: S51: Data Integration and Feature Matrix Construction The high-resolution near-surface air temperature spatialization data generated by steps S1-S4, the surface temperature retrieved from remote sensing features after radiometric calibration and atmospheric correction, and the digital elevation model (DEM) are used as variables to prepare the input data. Resampling is used to ensure that the elevation data is consistent with the spatial resolution of the temperature spatialization grid. All data were resampled to a uniform grid with the same resolution as the temperature field, maintaining spatial alignment, and the feature matrix was... X The grid points contain features including the spatialized near-surface air temperature distribution of glaciers, elevation variables, surface temperature retrieved from remote sensing data, and target variables. y Interpolated temperature values ​​corresponding to grid points.

[0049] S52: Training Set Partitioning and Model Selection The training and test sets were divided to ensure uniform spatial distribution and avoid regional bias. Five-fold cross-validation was used to optimize the model hyperparameters. Random forest regression was selected as the candidate model because it has strong anti-overfitting ability. Decision tree regression was used as the baseline comparison model. Multilayer perceptron was used to verify the ability to capture nonlinear relationships.

[0050] S53: Model Training and Hyperparameter Optimization During the training of the random forest model, each tree is trained independently based on randomly selected subsamples. Hyperparameter optimization adopts the grid search method. The main parameters adjusted include the number of trees, the maximum depth of the trees, and the minimum number of samples required for node splitting. The model performance is evaluated by the mean squared error to assess the magnitude of the prediction bias, and the coefficient of determination is used to measure the model's explanatory power for temperature changes. Through feature importance analysis, the influence of variables such as surface temperature and digital elevation model on temperature prediction can be determined. To evaluate the performance of different models, we compared the performance of random forest regression, decision tree regression, and multilayer perceptron on the test set. Furthermore, we analyzed the spatial distribution of the prediction errors of the models through residual analysis. We found that there are biases caused by the cooling effect in high-altitude glacier regions, which need to be corrected in subsequent steps. We also analyzed the characteristics of the model output. The optimal model is applied to predict grid points across the entire region, generating a near-surface temperature distribution map of the glacier area. The output results are fused with the temperature lapse rate data from step S6. Uncertainty is assessed through quantile regression and outlier labeling, generating a prediction interval of 5%-95% quantiles, quantifying the confidence range, and prompting manual verification for areas with residuals greater than 2 times the standard deviation.

[0051] A multi-source data fusion model was constructed, which integrates near-surface temperature data measured by meteorological stations with satellite remote sensing surface temperature data that has been calibrated with latitude and longitude coordinates and matched with band conversion through a ground-satellite collaborative framework. By combining ground observation data with digital elevation models, a spatialized near-surface temperature field was established in the glacier region as a benchmark training set. Surface temperature data was retrieved using high-resolution remote sensing images to make up for the insufficient observation coverage caused by the sparse stations in the plateau region.

[0052] S6: Spatialization process of glacier temperature based on altitude dependence optimization, data fusion and storage as follows Figure 3 As shown.

[0053] S61: Calculation and Adjustment of Vertical Temperature Lag Rate The regional temperature vertical lapse rate (LR, unit: ℃ / 100m) is calculated using meteorological station data to correct for glacier areas. If the station is located in a glacier area, the local temperature vertical lapse rate is used. Each target grid point is corrected according to its altitude and the temperature of the reference station to perform temperature field altitude correction. If multiple stations affect the target point, the adjusted temperature is combined by weighted average of multiple stations.

[0054] S62: Inverse Distance Weighted Interpolation and Smooth Transition Based on the inverse distance weighted interpolation method, the temperature of the target point is calculated by weighting the observation data of surrounding meteorological stations. The search range is expanded or virtual stations are added in areas with sparse data. Combined with the adjustment rate of temperature change with altitude, the spatial distribution field of near-surface temperature is finally generated. First, inverse distance weighted interpolation is performed on all grid points to obtain the initial temperature field. Then, the vertical temperature lapse rate is used for elevation correction, and the local lapse rate is used for further optimization in the glacier region. To assess the accuracy of the results, a cross-validation method was used, with a subset of station data serving as the validation set. The root mean square error and coefficient of determination before and after correction were calculated, and the residuals were analyzed to examine the error improvement in high-altitude areas. During the calculation process, the exponent of the inverse distance weighted interpolation can be adjusted according to different glacier regions. The final output temperature field is consistent with the resolution of the digital elevation model and is optimized for glacier regions. Sensitivity analysis of the adjustment of interpolation weights and descent rates is performed to quantify the uncertainty of the results and facilitate subsequent applications.

[0055] By optimizing the temperature spatialization process based on altitude dependence and using the regional temperature vertical lapse rate to perform gradient adjustment on the original temperature records, interpolation errors in high-altitude areas are effectively reduced. Combined with the inverse distance weighting method, a smooth transition of the temperature field is achieved in sparse observation intervals. For the special environment of glacier areas (such as cooling effects), station-specific temperature lapse rates are used for differentiated fitting to enhance the local adaptability of high-altitude areas, thereby improving the accuracy of glacier temperature input. The reliability and applicability of the temperature spatialization model under complex terrain are comprehensively considered.

[0056] Following the steps outlined above, the grid division and spatialization of near-surface air temperature data for the glacier region have been completed. This method ensures that even with insufficient or completely missing observational data, near-surface air temperature at any location within the glacier region can still be assigned and corrected. This technical approach is not only applicable to temperature estimation in glacier regions but can also be extended to temperature field construction in complex terrain and areas with sparse meteorological stations, providing more meteorological data support for high-altitude cryosphere research.

[0057] To address the challenges of near-surface air temperature monitoring in high-altitude glacial regions, this invention proposes a high-resolution spatialization method for near-surface air temperature based on multi-source data fusion. Traditional techniques rely on sparse station interpolation or single satellite data, which results in significant errors under complex terrain and lacks dynamic constraints on the vertical temperature lapse rate.

[0058] Practical data experiments were conducted on multi-source data fusion. By integrating measured data from multiple glaciers on the Qinghai-Tibet Plateau, surface temperature products from Landsat 8 / 9 satellites, and high-precision digital elevation models, a "glacier near-surface air temperature optimization fusion algorithm" specifically for glacier areas was developed. This algorithm uses the Isolation Forest (IForest) algorithm to identify anomalous temperature records and mark them as missing or outliers. It then uses the Delaunay triangulation method to interpolate and reconstruct missing values ​​in the original data and improves data quality through noise filtering techniques.

[0059] Regarding the correction of vertical temperature lapse rate, based on the temperature-altitude relationship model established in previous studies, a zonal correction for spatial heterogeneity of air temperature was achieved by introducing a quantitative description of the glacier cooling effect. This cooling effect is defined as the difference between the measured temperature at the glacier surface and the temperature estimated based on the vertical temperature lapse rate in non-glacier areas at the same altitude. In practical applications, spatiotemporally continuous temperature field products have been constructed for four typical glaciers: Guliya, Aru, Namunani, and Dunde.

[0060] To improve the accuracy of spatialized models of near-surface air temperature, this invention innovatively integrates multiple machine learning algorithms. Random forests, multilayer perceptrons, and decision tree regression models are trained using spatially processed station observation data, effectively correcting systematic biases in large-scale surface temperature retrieval results. By combining elevation data and vertical temperature lapse rate information, high-resolution air temperature products with a spatial resolution of 30 meters can be generated. Cross-validation of experimental data results was performed, and the test results on 58 independent validation samples showed that, compared with linear regression, this method reduced the root mean square error by an average of 3.2℃, significantly improving the consistency and continuity of the air temperature fields between glacial and non-glacial regions. Test results show that after processing with this algorithm, compared with linear regression, the root mean square error of air temperature data can be reduced by 0.3 to 7.6℃ (an average improvement of 3.2℃), and the standard deviation can be reduced by 0.3 to 5.0℃ (an average reduction of 1.8℃).

[0061] To further optimize the spatial interpolation process, this study proposes an elevation-constrained inverse distance-weighted interpolation method. This method effectively addresses the problems of topographic bias and sparse station data in high-altitude glacier regions by fusing station observations calibrated with remote sensing data and elevation information. By dynamically adjusting the weighting coefficients of elevation and vertical lapse rate in the interpolation, temperature prediction bias in high-altitude areas can be reduced.

[0062] Practical application of spatialized near-surface air temperature data from multiple glaciers on the Qinghai-Tibet Plateau demonstrates that the method established in this invention effectively compensates for the low accuracy of air temperature spatialization caused by the scarcity of observational data in glacial areas. Model sensitivity analysis shows that the performance of this method is correlated with the density of station data, suggesting that future efforts should focus on strengthening in-situ observations in glacial areas to better quantify and study the complex feedback mechanisms between glaciers and climate. The high spatial resolution air temperature products generated by the experiment will provide reliable data support for practical applications such as cold-region engineering planning and cryosphere hydrological model construction.

[0063] The specific embodiments of this invention are written in a progressive manner, emphasizing the differences between the various implementation schemes, and similar parts can be referred to each other. This method has been successfully applied in several experiments on spatializing near-surface air temperature in glaciers, and the relevant datasets have been published by the World Data Center for Climate (WDCC).

[0064] Multiple public datasets have been released (included in WDCC / DKRZ): Wang, Tianyun: Spatial distribution of air temperature in high-elevation glacierized regions: fromobservations in four catchments on theTibetan Plateau, World Data Center for Climate (WDCC) at DKRZ, https: / / doi.org / 10.26050 / WDCC / HGR, 2025a. Wang, Tianyun: Spatial Distribution of Near-surface Air Temperaturein the Glacierized Areas of theTibetan Plateau, World Data Center for Climate(WDCC) at DKRZ, https: / / doi.org / 10.26050 / WDCC / GATP, 2025b. Wang, Tianyun: Near-surface air temperature dataset for the Qinghai-Tibet Plateau (2019)derived from thermal infrared remote sensing andelevation-constrained modeling, World Data Center for Climate (WDCC) at DKRZ,https: / / doi.org / 10.26050 / WDCC / QTPTIR, 2025c. Wang, Tianyun: Spatialization of near-surface air temperature andupdating based on thermalinfrared remote sensing information in the Qinghai-Tibet Plateau, World Data Center for Climate (WDCC) at DKRZ, https: / / doi.org / 10.26050 / WDCC / QTP, 2025d. The spatial interpolation method used in this invention for correcting model training is described in detail in the specification. The complete Python code required for generating spatialized near-surface air temperature products for the Tibetan Plateau glacial region is available at: Wang, Tianyun (2025). Untitled ItemSpatialization for Near-surfaceAir Temperature in Glacierized Regions of the Qinghai-Tibet Plateau.figshare.Dataset. https: / / figshare.com / articles / dataset / _ / 30888905.

Claims

1. A method for spatializing near-surface air temperature in high-resolution glacial regions, characterized in that, include: Step 1: Collect near-surface temperature observation data and corresponding station coordinates provided by meteorological observation stations in the glacier area. Establish a near-surface temperature data sample set in the glacier area according to the time series. Mark the values ​​of 0 and invalid data in the sample set as missing values ​​and retain high-quality data. Step 2: Based on the near-surface temperature data sample set of the glacier region, establish a preliminary temperature distribution field. Use the nearest neighbor interpolation method to assign values ​​to the grid points of the preliminary temperature distribution field using high-quality data, and output the preliminary spatial distribution of temperature. Step 3: Use triangulation-based interpolation to reconstruct the preliminary temperature distribution field into a high-resolution temperature distribution field for the missing value data; Step 4: Based on the near-surface air temperature observation data provided by the meteorological observation station and the satellite remote sensing surface temperature data, perform time synchronization and spatial alignment to refine the high-resolution temperature distribution field and obtain a spatialized model of near-surface air temperature in the glacier area; Step 5: Integrate the near-surface temperature data measured by meteorological stations with the satellite remote sensing surface temperature data that has been calibrated with latitude and longitude coordinates and matched with band conversion through the ground-satellite collaborative framework, and build a benchmark training set in combination with the digital elevation model to train the spatialization model of near-surface temperature in the glacier area. Step 6: Correct the error of the trained spatial model of near-surface air temperature in the glacier area by calculating the vertical temperature lapse rate and using the inverse distance weighted interpolation method for elevation variables, and generate the spatial distribution field of near-surface air temperature by combining the air temperature data of the glacier area.

2. The high-resolution spatialization method for near-surface air temperature in glacial regions according to claim 1, characterized in that, In step 1, the Isolation Forest algorithm is used to mark values ​​of 0 and invalid data in the sample set as missing values. The Isolation Forest algorithm then extracts data from the data sample set without replacement. φ A default-sized subsample is used, and an isolation tree is trained for each subsample. For each tree, features and splitting values ​​are randomly selected. Features include time, near-surface air temperature observation data, and observation point coordinates. The data space is recursively partitioned until a data point is isolated. The maximum depth of the isolation tree is set to [value missing]. The rounded-up value is used to avoid overfitting. A tree is constructed by randomly splitting the data. Outliers, because they deviate from the population, can be isolated to shallow nodes of the tree with only a few partitions. Anomaly scores are calculated. ； in h This represents the path length, which is the number of edges from the root node to the isolated data point. E ( h ( x )) is a point x Average path length across all trees φ Number of subsamples ， c( φ The normalization factor for path length is 0 to 1. The standardized score range is 0 to 1. The threshold is dynamically adjusted in combination with the glacial climate background. If the abnormal score is higher than the threshold, it is judged as abnormal.

3. The high-resolution spatialization method for near-surface air temperature in glacial regions according to claim 2, characterized in that, In step 1, the threshold is dynamically adjusted in combination with the glacier climate background. If the abnormal score is higher than the threshold, it is judged as abnormal. Specifically, for spatial anomalies, when the temperature of a certain observation point in the glacier observation station area deviates significantly from that of the surrounding stations, it is marked as abnormal even if the score does not exceed the threshold, and missing values ​​are marked. For glacier areas with extreme climate disturbances and significant diurnal temperature variations, adaptive optimization is performed by combining multi-source data. For cold wave scenarios, spatial verification is performed by introducing remote sensing reanalysis data. If the low temperature is consistent over a large area, it is determined to be a real event and the data is retained. For stations with large diurnal temperature variations, sliding window standardization is adopted, and the mean and standard deviation are calculated with a 7-day window. The current value is marked as an anomaly when it exceeds the mean ± 3 times the standard deviation, so as to avoid misjudging normal diurnal fluctuations.

4. The high-resolution spatialization method for near-surface air temperature in glacial regions according to claim 1, characterized in that, In step 2, the coverage area of ​​the weather stations is divided into grids. Based on the nearest neighbor algorithm, the observation data of each weather station is used as the initial near-surface temperature distribution estimate for the corresponding grid. The temperature data used comes from automatic weather stations and temperature recorders in different watersheds. The coordinates of the weather stations within the target coverage area are organized into a set in logical order: X (station) ={1,…, X j ,…, X n },in X (station) Represents the set of coordinates of weather stations. X j Represents the coordinates of a single weather station. n Indicates the total number of stations; The coverage area of ​​the weather station is divided into grids: X (grid) ={1,…, X i ,…, X m },in X (grid) Represents the set of grid coordinates. X i Represents the center coordinates of a single grid cell. m Indicates the number of grid points; Assign a unique index to each site X (station) ={1,…, X j ,…, X n } is used for subsequent matching, which calculates the distance between grid points and stations based on the meteorological station index number and performs nearest neighbor matching; Based on nearest neighbor learning, a brute-force search is used to calculate the distance between each station and each grid center point: ； L This indicates that Euclidean distance is used as the metric, which is the straight-line distance between two points in the standard coordinate system. When the station and the grid have the same index, the near-surface temperature value of the station is directly assigned to the corresponding grid. Within the coverage area of ​​the meteorological station, the observation values ​​will be distributed to each neighboring grid to form an initial temperature distribution field. The meteorological observation data is used to assign values ​​to the surrounding station areas to ensure that the temperature data of each grid comes from the actual observation values.

5. The high-resolution spatialization method for near-surface air temperature in glacial regions according to claim 1, characterized in that, In step 3, based on the missing value data from step 1, the mean and standard deviation of the data sample set are calculated, using a 3x3 matrix. σ The criteria include reassigning extreme temperature observations exceeding ±3 standard deviations, cleaning and standardizing the data, and recording the spatial distribution characteristics of all outliers, including altitude and latitude / longitude information. Then, a Voronoi diagram is constructed based on station coordinates, using the stations as the center points of the map elements to achieve spatial partitioning. The actual coordinates of the stations are the centers of the Voronoi units. Furthermore, a Delaunay triangulation is constructed using the intersections of Voronoi units and station coordinates, ensuring full coverage of the study area. An irregular triangulation is constructed, and the initial near-surface air temperature distribution estimation results are iteratively corrected through triangulation. The residuals between the interpolated grid node temperatures and the measured data are calculated. Nodes with residuals greater than 2°C undergo secondary interpolation, iterative residual optimization is performed to improve local accuracy, and adaptive threshold optimization is used to ensure that the reconstructed temperature field conforms to the actual climate distribution patterns.

6. The high-resolution spatialization method for near-surface air temperature in glacial regions according to claim 5, characterized in that, In the interpolation process of step 3, to avoid errors in single-station observation data leading to overestimation or underestimation of surrounding areas, a smoothing function needs to be applied to the interpolated near-surface air temperature distribution within the station network coverage area: A 5×5 sliding window is used to calculate the median of the temperature values ​​within the window and replace the center pixel value to effectively remove noise and outliers. The window size is dynamically adjusted, and when the proportion of valid data points within the window is detected to be lower than the threshold, the window is expanded to ensure sufficient samples. When performing median filtering, special attention is paid to key areas such as the glacier-bare rock transition zone and the moraine region to avoid excessive smoothing that could lead to the loss of temperature transition information. Filter intensity thresholds are set separately for glacier surfaces and non-glacier surfaces to ensure that small-scale temperature gradients are preserved. The local gradient is calculated for the filtered temperature field. If the temperature difference between adjacent grid cells exceeds a threshold, the standard deviation of the Gaussian kernel in that region is reduced. σ Values ​​are optimized to preserve the characteristics of climate abrupt changes and to achieve spatial continuity. Cross-validation is used to calculate the root mean square error and standard deviation to ensure that the interpolated temperature field conforms to the observation trend. For areas with large residuals, corrections are made by combining satellite remote sensing surface temperature data to adjust the near-surface air temperature distribution and output a high-resolution continuous near-surface temperature field.

7. The high-resolution spatialization method for near-surface air temperature in glacial regions according to claim 1, characterized in that, Step 4 involves time synchronization and spatial alignment based on near-surface air temperature observation data provided by meteorological observation stations and satellite remote sensing surface temperature data. Step 41: Temperature data collected by meteorological observation stations are usually recorded at fixed time periods, with a certain time interval between adjacent observation times. When the satellite's transit time falls between two consecutive observation times at a meteorological station, the station observation times in the remote sensing data fusion are interpolated. Based on all available station observation data, a third-order uniform B-spline curve is fitted to interpolate the station temperatures to the precise remote sensing imaging timestamp, ensuring that the temperature data is interpolated in time and matched with the satellite center's transit time. This process connects two adjacent observation times within each station. Meteorological records for different time periods are divided into subdivisions. n Sub-intervals The step size of each subinterval can be expressed as: ； Where, Δ t This indicates the time step, which is the duration of the interval between two time points. This is the end time of the period. The start time of the time period. n This indicates the number of sub-intervals after the subdivision. to Any time between can be represented as ; Step 42: Use a third-order uniform B-spline curve to... and Interpolation is performed on the temperature observation records corresponding to each time point. Each basis function is a cubic polynomial, and the spline curves are uniformly distributed in the parameter node space. Each control point corresponds to one basis function. For two time points... arrive Observation data points between Each of them represent The temperature observation value at time, the cubic spline interpolation function f( t The following conditions must be met: ； in, The start time of the time period. n This indicates the number of sub-intervals after further subdivision; to Any time between is represented as , i = 0, 1, 2, ..., n Each paragraph All are cubic polynomials calculated according to the following rules: ； For B-spline representation, the interpolation function is expressed as: ； Where, β j ( t ) represents the cubic B-spline basis functions. α j To be based on data constraints f i ( t i )= y i The calculated weighting coefficients are used as the f( of the cubic polynomial). t The cubic spline interpolation function is obtained by multiplying the basis functions and the control point weights and then summing the results to obtain the expression of the interpolation curve. It can obtain the temperature value at any satellite center time, achieve time synchronization and spatial alignment, improve the high-resolution temperature distribution field, and obtain a spatialized model of near-surface temperature in glacier areas. Step 43: When the image time of the remote sensing satellite is obtained, but no meteorological observation data is recorded at that time, the meteorological observation information at the next time point is predicted. Then, the temperature is interpolated between the previous time point and the current time point to locate and align with the remote sensing image time to a specific time point. The temperature interpolation result at this precise time point is used to represent the near-surface temperature distribution of the meteorological station. The near-surface temperature at the next time point of each meteorological station is predicted using support vector regression. Finally, the optimized model is obtained by minimizing the total loss and maximizing the interval. Support Vector Regression (SVR) builds a regression model by finding the optimal hyperplane, minimizing the error between the predicted values ​​and the observed values ​​of the training samples. For nonlinear models, SVR uses a kernel function to map to the feature space before performing regression. The loss function of SVR is expressed as: ； in, y This represents the temperature observation value recorded by the weather station. η The predicted value represents the near-surface air temperature, and is within the tolerance range of regression error. If the deviation is within... ε If the value is within the specified range, it will not be considered an error, L( y , η ) is the loss function, ignoring minor inaccuracies in its calculation process.

8. The high-resolution spatialization method for near-surface air temperature in glacial regions according to claim 7, characterized in that, In step 4, to effectively analyze and utilize the residuals of support vector regression in predicting air temperature and their impact on temperature prediction, a diverse dataset is used for model training, and measures are implemented to avoid overfitting, as shown below: ； in, y This represents the temperature data observed by the weather station. m This represents the end point of the observation record. i Iterate from the initial time point 0 to the current time point. m At a certain point in time, m +1 represents the total number of observation points required for model training, while y i Indicates the first i Observation records at each time point; ； ； ； ； in, X It is a matrix of temperature observation data y The generated ( m +1)×( k +1) dimensional matrix, I for( m +1)×( m +1) dimensional identity matrix, matrix X any element in x ij Represented as: ； in, δ ij It is a random number between -1 and 1, that is, -1 < δ ij < 1, m This is the end time point of the observation record. i From the initial time point 0 to the... m The sequence number of each time point. k Representation matrix X The largest column index, j The column index is 0, and its value ranges from 0 to 1. k , k +1 is a matrix X The total number of columns; Divide the training sample set D = {( x 1, y 1), …, ( x i , y i ), …, ( x m , y m The data generation rules are as follows: x 1. From initial observations y 0 generation, specifically satisfying: x 1j = ( δ 0j + 1) y 0, where j = 0,1,…, k ; This dataset is based on meteorological records and uses a regression algorithm to generate the training set. X train yes m ×( k +1) dimension matrix is ​​derived from temperature observation matrix y Generate, its elements x ij Corresponding to the 1st to the 2nd m Time point matrix y train Also covering numbers 1 to 2 m The target for predicting the observations at the () time point is the () m The remote sensing image of the temperature value at time point +1) is located at the [number]th time point. m With the ( m +1) Between time points, it is represented as: ； ； ； No. m The training data at each time point was generated through calculation, where x m+1 =(1+ δ mj ) y m , δ mj For random disturbance terms, y m It is the first m One training data value, y m+1 As the true value of the test data, a near-surface air temperature was established based on a regression algorithm. T The prediction model and its calculation process are as follows: ； Near-surface air temperature at the moment of remote sensing imaging is processed by a cubic spline function. y m+1 Compared with the predicted value y prediction For interpolation calculations, in a model trained based on observation time steps and their corresponding near-surface air temperature records, the data representation used for training is assumed to be: ； ； In the m At a certain point in time, y m It is the first m One training data value, y prediction Yes m The predicted value at time +1; the air temperature at the time of remote sensing imaging was obtained using a cubic spline function. y i and y prediction Interpolation estimation is performed, and the remote sensing image center acquires time points. y m+1 The interpolation results can be used as test data. The model is evaluated on an independent test set, and the predictive performance of different residual subgroups is quantified by the root mean square error and coefficient of determination to verify its stability.

9. The high-resolution spatialization method for near-surface air temperature in glacial regions according to claim 1, characterized in that, In step 5, the high-resolution near-surface air temperature spatialization data generated in steps 1-4, the surface temperature retrieved from remote sensing features after radiometric calibration and atmospheric correction, and the digital elevation model are used as variables to prepare the input data. Resampling is used to ensure that the elevation data and the spatial resolution of the temperature spatialization grid are consistent. All data were resampled to a uniform grid with the same resolution as the temperature field, maintaining spatial alignment, and the feature matrix was... X The grid points contain features including the spatialized near-surface air temperature distribution of glaciers, elevation variables, surface temperature retrieved from remote sensing data, and target variables. y Interpolated temperature values ​​corresponding to grid points; The training and test sets were divided to ensure uniform spatial distribution and avoid regional bias. Five-fold cross-validation was used to optimize the model hyperparameters. Random forest regression was selected as the candidate model, and decision tree regression was used as the baseline comparison model. Multilayer perceptron was used to verify the ability to capture nonlinear relationships. During the training of the random forest model, each tree is trained independently based on randomly selected subsamples. Hyperparameter optimization adopts a grid search method, and the parameters adjusted include the number of trees, the maximum depth of the trees, and the minimum number of samples required for node splitting. The model performance is evaluated by the mean squared error to assess the magnitude of the prediction bias, and the coefficient of determination is used to measure the model's explanatory power for temperature changes. Through feature importance analysis, the influence of variables such as surface temperature and digital elevation model on temperature prediction is determined.

10. The high-resolution spatialization method for near-surface air temperature in glacial regions according to claim 1, characterized in that, In step 6, the vertical temperature lapse rate of the region is calculated using meteorological station data to correct for the glacier area. When a station is located in the glacier area, the local vertical temperature lapse rate is used to correct each target grid point according to its altitude and the temperature of the reference station, performing temperature field altitude correction. If multiple stations affect the target point, the adjusted temperature is combined by weighted average of multiple stations. The temperature of the target point is calculated by weighting the observation data of surrounding meteorological stations using the inverse distance weighted interpolation method. In areas with sparse data, the search range is expanded or virtual stations are added. Combined with the adjustment rate of temperature change with altitude, a spatial distribution field of near-surface temperature is generated. Inverse distance weighted interpolation is performed on all grid points, and then altitude correction is performed using the vertical temperature lapse rate. In the glacier area, the local lapse rate is used for further optimization. Finally, the cross-validation method is used, with some station data as the validation set, to calculate the root mean square error and coefficient of determination before and after correction. At the same time, the residuals are analyzed to verify the error improvement in high-altitude areas.

Images