Global climate long sequence data reconstruction method for heat stress risk assessment

Abstract

The invention discloses a global climate long-sequence data reconstruction method for heat stress risk assessment, and the method comprises the steps: carrying out the spatial interpolation of a meteorological data missing value of a target station through a ridge regression model based on the data of a spatially adjacent reference station; and time interpolation is carried out through a ridge regression model based on the physical relationship among different meteorological variables of the target station. According to the method, a ridge regression technology is adopted, coefficient estimation is stabilized by introducing a regularization term, multiple collinearity is effectively overcome, space-time collaborative interpolation is achieved, and the method is high in precision and self-adaptive.

Description

Technical Field

[0001] This invention relates to the field of meteorological data reconstruction technology, and in particular to a method for reconstructing long-series global climate data for heat stress risk assessment. Background Technology

[0002] In the context of global climate change, thermal stress risk assessment is crucial for protecting public health, guiding outdoor work, and developing adaptation policies. Wet-bulb black-bulb temperature (WBGT), as a composite indicator that comprehensively considers temperature and humidity, has been adopted by the International Organization for Standardization (ISO) as a core parameter for evaluating environmental heat load, and its long-term evolution trend forms the cornerstone of risk assessment.

[0003] Currently, the commonly used WBGT trend assessment and risk analysis methods in the industry are mainly based on two types of data sources: meteorological reanalysis data and raw meteorological station observation data. The typical implementation path is as follows: directly using the temperature and humidity fields provided by the above data sources, the historical daily WBGT series for each grid point or station is calculated; then, statistical methods (such as linear fitting) are used to analyze the long-term trend of WBGT changes, and the frequency and changes of days exceeding a preset risk threshold (e.g., WBGT ≥ 30°C) are calculated. Based on this, a qualitative or semi-quantitative judgment is made on the level and increasing / decreasing trend of regional heat risk.

[0004] However, the aforementioned existing technical solutions have the following inherent drawbacks: The spatiotemporal reliability of the assessment results is limited by the quality of the underlying data: While reanalysis data has the advantages of continuous spatiotemporal coverage and complete sequences, it is essentially a fusion product of numerical weather prediction models and assimilation systems on sparse observations. Systematic biases in the models themselves, changes in assimilation schemes, and the non-uniformity of the observation systems themselves all introduce "unnatural signals" into the reanalysis data, leading to uncertainty in its reproduction of the true climate state and thus affecting the reliability of the derived WBGT assessment conclusions. On the other hand, although ground-based station observations best reflect the local true climate state, the global observation network is extremely unevenly distributed (dense over land, sparse over ocean; dense in the Northern Hemisphere, sparse in the Southern Hemisphere), and historical data generally suffers from numerous missing, discontinuous, and various quality problems. Directly conducting trend statistics and risk assessments based on such incomplete and spatially sparse raw sequences is crucial for data-scarce regions (such as Africa). Against the backdrop of escalating global climate change, heat stress risk assessment has become a key support for protecting public health, regulating outdoor activities, and formulating scientific climate adaptation policies. Wet-bulb black-sphere temperature (WBGT), as a composite indicator that comprehensively considers multiple factors such as air temperature and humidity, has been established by the International Organization for Standardization (ISO) as a core parameter for heat stress assessment due to its accurate characterization of environmental heat load. Its long-term evolution patterns form the core foundation for conducting regional and even global heat risk assessments.

[0005] Therefore, the spatiotemporal reliability of the evaluation results is severely constrained, specifically in the following two aspects: On the one hand, while meteorological reanalysis data possesses advantages such as continuous spatiotemporal coverage and complete data sequences, it is essentially a product of the fusion of sparse global observational data by numerical weather prediction models and data assimilation systems. During the data generation process, systematic biases inherent in the numerical weather prediction models themselves, iterative updates to the assimilation scheme, and non-uniform variations in global observational systems at different times are all introduced into the reanalysis data in the form of "non-natural climate signals." This leads to inherent uncertainties in the accuracy of its reconstruction of the true climate state, and consequently directly affects the reliability of WBGT assessment conclusions derived from this type of data.

[0006] On the other hand, while ground-based meteorological station observation data can most directly reflect the actual local climate conditions, the global meteorological observation network is significantly unevenly distributed, with denser observation stations in land areas and extremely sparser ones in ocean areas, and a much higher station density in the Northern Hemisphere than in the Southern Hemisphere. Furthermore, historical meteorological station data worldwide generally suffers from numerous quality issues such as missing data, discontinuous recordings, inconsistent formats, and errors. If trend statistics and thermal risk assessments are directly conducted based on this sparsely distributed, incomplete time-series raw data, the assessment results will exhibit extremely high uncertainty in data-scarce regions such as inland Africa, parts of South America, and vast oceans. Ultimately, this will result in a lack of sufficient robustness in global or large-regional thermal stress risk assessments, making it difficult to meet the needs of precise decision-making.

[0007] Given the inconsistent quality and significant gaps in historical meteorological station data worldwide, existing technological approaches have attempted to compensate for data deficiencies and improve sequence integrity through data interpolation. However, mainstream interpolation schemes still face key technical bottlenecks and have failed to fundamentally solve the problem. Specifically: The core principle of ordinary least squares regression imputation is to construct a linear regression model (Y=Xβ+ε) between the independent variable (such as the observation data of neighboring stations during the same period, and other relevant meteorological element data of this station) and the dependent variable (missing values ​​of the target variable of this station), estimate the regression coefficient β by minimizing the sum of squared residuals, and then use this model to predict missing values.

[0008] However, in meteorological data interpolation scenarios, the independent variables used to build models are often highly correlated. For example, temperature data from multiple meteorological stations within the same climate zone are affected by the same weather system, exhibiting strong synchronicity in their temporal changes; furthermore, meteorological variables such as temperature, dew point temperature, and air pressure at the same station inherently possess physical relationships. This high correlation among independent variables is known as multicollinearity, which directly leads to the following serious drawbacks: When severe multicollinearity exists, the (XᵀX) matrix in the model approaches a singular state, and the numerical calculation of its inverse matrix is ​​highly susceptible to fluctuations. This makes the estimated regression coefficient β extremely sensitive to even minor perturbations in the original data. In practical applications, this can lead to large fluctuations in the regression coefficients or even sign reversals, completely violating the physical laws governing meteorological elements. Regression models trained on data with severe multicollinearity exhibit extremely high variance in their extrapolation predictions (i.e., data interpolation), making it difficult to control the deviation between predicted and actual values, and resulting in extremely low reliability of the interpolation results.

[0009] In summary, directly applying ordinary least squares regression to interpolation of multi-site or multivariate meteorological data has been proven unreliable both in theoretical logic and in practical applications.

[0010] Current mainstream imputation methods generally adopt a single-dimensional approach, either performing only spatial imputation (such as using only data from neighboring stations during the same period to fill missing values) or only temporal imputation (such as building a time series model based on the station's historical data series to predict missing values). This separates the spatial and temporal dimensions of the data, leading to the following inherent drawbacks: Pure spatial interpolation completely ignores the rich information contained in the historical data of this station, such as the diurnal variation patterns of meteorological elements, seasonal cycle characteristics, and physical constraints between variables; while pure temporal interpolation abandons the strong external constraints provided by the spatial field and fails to make full use of the reference value of meteorological data in neighboring areas, resulting in the waste of potential data information.

[0011] When all neighboring stations in space experience data gaps during the target time period, pure spatial interpolation schemes will completely fail. Furthermore, when the target variable at the local station experiences long-term continuous data gaps, time-series models relying solely on historical data will rapidly decline in predictive power, and interpolation accuracy cannot be guaranteed. Interpolation methods that separate the spatial and temporal dimensions may produce physically unreasonable results. For example, the relative humidity calculated from the temperature value obtained through spatial interpolation and the dew point temperature observed at the local station may exceed 100%, violating fundamental principles of atmospheric thermodynamics. Such interpolated data cannot be used for subsequent climate trend analysis and risk assessment. Therefore, a method for reconstructing long-term global climate data for heat stress risk assessment is needed. Summary of the Invention

[0012] In view of the aforementioned existing problems, the present invention is proposed.

[0013] Therefore, this invention provides a global climate long-series data reconstruction method for heat stress risk assessment, which solves the problems of multicollinearity interference, fragmented use of spatiotemporal information, and poor physical consistency of interpolation results faced in multivariate meteorological data interpolation.

[0014] To solve the above-mentioned technical problems, the present invention provides the following technical solution: In a first aspect, the present invention provides a method for reconstructing long-series global climate data for heat stress risk assessment, which includes, for the missing meteorological data of the target station, firstly, spatially interpolating based on spatially adjacent reference station data through a ridge regression model, and then temporally interpolating based on the physical relationship between different meteorological variables of the target station through a ridge regression model.

[0015] Further, the spatial interpolation includes the following steps: grouping the daily-scale data of the target station and reference stations by variable and by month; using the effective observation sequence of the target variable for that month at the target station as the dependent variable Y, selecting the K highest-ranked reference stations' data of the same period and the same variable to form the independent variable matrix X, and constructing a ridge regression model. The regression coefficient calculation formula for the ridge regression model is as follows:

[0016] Where λ is the regularization parameter (ridge parameter) and I is the identity matrix; the constructed ridge regression model is used to imputate and predict the missing values ​​of the target variable for the target station in that month.

[0017] Furthermore, the time interpolation includes the following steps: grouping the spatially interpolated multivariate daily-scale data of the target station by month, wherein the multivariate data includes at least the target variable and key covariates; for each monthly group, constructing a bivariate ridge regression model using the physical relationships between the variables, wherein the bivariate ridge regression model uses the target variable as the dependent variable and at least one other meteorological variable as the independent variable, and its regression coefficients are... The calculation formula is:

[0018] in For the sequence of independent variables, For the dependent variable sequence, For regularization parameters, The identity matrix is ​​used; the constructed bivariate ridge regression model is used to perform secondary imputation on the target variable values ​​that are still missing after spatial imputation.

[0019] Furthermore, it also includes a model selection step: calculating the determination coefficients R for each ridge regression model used for spatial and temporal interpolation, respectively. 2 ; Set R 2 The minimum performance threshold, using only R 2 The output value of the model is higher than the preset threshold; when multiple qualified model output values ​​exist for the same missing data point, R is selected. 2 The highest output value is used as the final interpolation value for the missing data point.

[0020] Furthermore, the ridge parameter λ in the ridge regression model is determined automatically and adaptively using a random forest algorithm.

[0021] Furthermore, before spatial interpolation, a data preprocessing step is included: mapping all global meteorological stations to a regular latitude and longitude grid with a preset resolution; within each grid cell, a representative station is selected based on a preset comprehensive data integrity index, which is the number of days on which the target interpolation variable and the key covariate have valid observations on the same calendar date, and is weighted and scored in combination with the total number of valid record years of the station, and the station with the highest score is selected as the representative station of the grid.

[0022] Furthermore, the regular latitude and longitude grid is a global regular grid of 5°×5°.

[0023] Furthermore, the other meteorological variables include dew point temperature and local air pressure, and the key covariates include daily average dew point temperature and daily average local air pressure.

[0024] Furthermore, the steps for constructing the space reference station set are as follows: taking the grid where the target station is located as the center, all other stations in that grid and several geographically adjacent grids are initially selected as the space reference station candidate set; the comprehensive geographical distance index between the target station and each reference station candidate is calculated, the comprehensive geographical distance index is calculated based on spherical distance, latitude and longitude difference and altitude difference; all reference station candidates are sorted from smallest to largest according to the comprehensive geographical distance index, and the top K reference stations are selected to form the space reference station set, where K is set to 5-10 according to the data density.

[0025] Furthermore, when constructing the ridge regression model, the valid data for the corresponding month are divided into a training set and a validation set in an 8:2 ratio. The model is trained using the training set and its performance is validated using the validation set.

[0026] The beneficial effects of this invention are: The core of this invention employs ridge regression technology, which stabilizes coefficient estimation by introducing a regularization term, fundamentally solving the multicollinearity problem commonly found in meteorological data, making the interpolation model more robust and the results more reliable.

[0027] This invention employs a two-stage interpolation framework to maintain physical consistency of interpolation results across both spatiotemporal scales. This method does not rely on external data sources, utilizing only the inherent characteristics of station data, making it particularly suitable for processing global historical meteorological data. It ensures that only models meeting performance standards are adopted and automatically selects the optimal interpolation source, thus automating quality control in the interpolation process and significantly improving the accuracy of the final data product. Its standardized process facilitates automated batch processing, providing core technical support for constructing high-quality, long-sequence global meteorological reanalysis datasets. Attached Figure Description

[0028] To more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings used in the following description of the embodiments will be briefly introduced. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0029] Figure 1 A global distribution and data characteristic map of meteorological stations used in a global climate long-series data reconstruction method for heat stress risk assessment.

[0030] Figure 1 (a) Spatial distribution of all stations; circle color indicates the number of data years (0-52 years); (b) The temporal evolution of the number of all sites and the number of sites with data years ≥ 35 years; (c) Probability density distribution of all sites and data years of sites with data years ≥ 35 years; dashed lines and dotted lines represent the mean and median, respectively.

[0031] Figure 2 This is a schematic diagram showing the spatial and latitudinal distribution of the wet-bulb black-sphere temperature trend from 1973 to 2024 in this implementation case.

[0032] Figure 3 This is a graph showing the temporal evolution of the global area-weighted average wet-bulb black-sphere temperature anomaly in this implementation case. Figure 3 (a) Global distribution pattern of WBGT linear trend (°C / decade); (b) WBGT trend chart distributed by latitudinal zone; Figure 4 This is a long-term trend chart (1973-2024) of the number of days of heat stress per year for the three heat stress categories defined by WBGT in this implementation case: (a) Schematic diagram of the number of days below 25°C (low heat stress); (b) Schematic diagram of the number of days between 25-30°C (moderate to high heat stress); (c) Schematic diagram of the number of days above 30°C (high to extremely high heat stress); (d) Schematic diagram of the distribution trend by latitude zone; (e) Schematic diagram of the probability density distribution of the trend. Detailed Implementation

[0033] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings.

[0034] Many specific details are set forth in the following description in order to provide a full understanding of the invention. However, the invention may also be practiced in other ways different from those described herein, and those skilled in the art can make similar extensions without departing from the spirit of the invention. Therefore, the invention is not limited to the specific embodiments disclosed below.

[0035] Secondly, the term "one embodiment" or "example" as used herein refers to a specific feature, structure, or characteristic that may be included in at least one implementation of the invention. The appearance of an embodiment in different places in this specification does not necessarily refer to the same embodiment, nor is it a single or selective embodiment that mutually excludes other embodiments.

[0036] This embodiment provides a method for reconstructing long-series global climate data for heat stress risk assessment, including the following steps: S1: Global Station Grid System and Representative Station Selection All global meteorological stations are mapped onto a latitude and longitude grid (e.g., 5°×5°) of a preset resolution. Within each grid cell, a representative station is selected based on a preset comprehensive data integrity index. The preferred comprehensive data integrity index is the number of days on which the target imputation variable (e.g., daily average temperature T) and key covariates (e.g., daily average dew point temperature Td, daily average local air pressure SP) have valid observations on the same calendar date. This score is then weighted and combined with the station's total number of valid record years, and the station with the highest score is selected as the representative station for that grid.

[0037] S2: Ridge Regression Interpolation in Spatial Dimensions For each representative station selected in step S1, the following spatial interpolation procedure is performed: S21: Construct a set of space reference stations: Using the representative station as the center, select all other stations in its grid and several geographically adjacent grids as the initial set of space reference stations.

[0038] S22: Calculate and rank geographic feature distances: Calculate the spherical distance between the representative station and each candidate reference station, and comprehensively consider the differences in latitude and longitude and altitude to calculate a comprehensive geographic distance index. Rank all candidate reference stations according to this index from smallest to largest.

[0039] S23: Grouping data by month: The daily-scale data of the representative station (dependent variable) and each sorted reference station (independent variable) are grouped by variable and by month. This step aims to eliminate the influence of seasonal cycles on the relationship between variables and to establish an independent regression model for each month.

[0040] S24: Constructing a Ridge Regression Model by Month: For each month of the target variable (e.g., temperature T), group the data by month. Using the effective observation sequence of the representative station for that month as the dependent variable Y, select the K highest-ranked reference stations (K can be set according to data density, e.g., the top 5-10) and their corresponding data to form the independent variable matrix X. Construct the following ridge regression model to estimate the regression coefficient β:

[0041] in, For regularization parameters (ridge parameters). The λI term is introduced to improve the condition number of the (XᵀX) matrix, thereby stabilizing the coefficient estimation and overcoming multicollinearity.

[0042] S25: Model Training and Imputation: Divide the valid data for the current month into a training set and a validation set (e.g., 8:2). Using the training set, determine the optimal ridge parameter λ through cross-validation or an automated algorithm (such as random forest optimization described later). Use the trained model to imputate and predict the missing values ​​of the target variable representing the station for that month.

[0043] S3: Ridge Regression Interpolation in Time Dimension After completing the spatial interpolation in step S2, further temporal interpolation is performed on the partially interpolated data sequence for each representative station: S31: Single-station data grouped by month: The multivariate daily-scale data of this station after spatial interpolation (including at least the target variable T and covariates Td and SP) are grouped by month again.

[0044] S32: Construct a bivariate time-ridge regression model: Group by month and construct a bivariate ridge regression model using the known physical relationships between variables. For example, a model could be constructed as follows: Model M1: Temperature T is the dependent variable, and dew point temperature Td is the independent variable.

[0045] Model M2: with temperature T as the dependent variable and local air pressure SP as the independent variable.

[0046] Model M3: Dew point temperature Td is the dependent variable, and local air pressure SP is the independent variable (can be adjusted as needed).

[0047] Each model is constructed in the same way as S24. Right now: : Where X is the sequence of independent variables and Y is the sequence of dependent variables.

[0048] It should be noted that the β of spatially interpolated temperature T represents the contribution weight of the T data from the five neighboring stations to the T data of the target station; the β of temporally interpolated temperature T represents... This refers to the dew point temperature (Td) data of this station, and its correlation weight with the T data of this station.

[0049] S33: Model Training and Secondary Imputation: Similarly, divide the training / validation sets and determine the optimal λ. Use the trained temporal ridge regression model to perform secondary imputation on the target variable values ​​that are still missing after spatial imputation.

[0050] S4: Fusion of Model Performance Control and Optimal Interpolation Values Strict quality control is implemented during the model building process in steps S24 and S32: The determination coefficient R0 for each trained model is calculated using the reserved validation set. 2 .

[0051] Set an acceptable minimum performance threshold (e.g., R). 2 >0.5). Only R 2 Only models exceeding this threshold are allowed for actual interpolation.

[0052] For the same missing data point, there may be both qualified spatial interpolation models and multiple temporal interpolation models that can provide interpolated values. In this case, select R from all qualified models. 2 The value output by the highest-performing model is used as the final imputation value for the missing point. This selection mechanism ensures that the final data product consists of the best prediction results available at that time and place.

[0053] The dataset used is the Global Surface Daily Summary Dataset released by the National Oceanic and Atmospheric Administration of a certain country. This dataset is derived from the Integrated Surface Hourly Dataset, which integrates global data archived at the Federal Climate Integration Center by the Air Force Climatology Center and the National Environmental Information Center of a certain country. This dataset provides daily meteorological elements available at each station, such as mean temperature, dew point temperature, mean sea level pressure, and mean local pressure. All data are reported and summarized in Greenwich Mean Time, consistent with the original weather / hourly observation standards.

[0054] Global multi-factor daily GSOD data is continuously updated. Between 1947 and 2024, data was recorded from 27,890 stations, but some of these sites have inherent data quality issues. Referring to the quality control module implemented in the RClimDex software, a comprehensive QC procedure was applied to station metadata and variables including temperature, dew point temperature, and local atmospheric pressure. This process included the following steps: spatial location verification, temporal repeatability checks, threshold-based verification, and manual inspection and correction.

[0055] Spatial location verification: Stations with missing longitude or latitude, or whose coordinates are both reported as zero, have been removed. Valid longitude and latitude ranges are required to be -180° to 180° and -90° to 90°, respectively.

[0056] Time repeatability check: Merge records from the same site identifier, identify entries with repeating dates, and reduce them to a single observation for each calendar day.

[0057] Threshold-based validation: Values ​​outside the physically reasonable range were marked as missing. The acceptable ranges for each variable were defined as follows: air temperature was limited to -89.4°C to 57.8°C, consistent with officially recorded global extreme values; dew point temperature was limited to -98.2°C to 50°C to reflect thermodynamic limits and observational records; and both station pressure and sea level pressure were limited to the range of 300 hPa to 1100 hPa.

[0058] The analysis utilized quality control station data, encompassing a total of 24,651 stations. These stations were divided into 5°×5° latitude and longitude grid cells. Within each grid cell, a representative station was selected based on optimal data integrity, specifically the availability of synchronous temperature, dew point temperature, and local air pressure on the same day, as well as the highest number of valid records. This process resulted in 1,164 representative stations globally. Although these representative stations exhibited better data integrity compared to other stations within their respective grid cells, significant data gaps and periods of low-quality data still existed. Therefore, a key objective of this invention is to apply an optimal method to impute missing values ​​as completely as possible.

[0059] The ridge regression model is constructed as follows: The core algorithm of ordinary least squares (OLS) requires that the matrix X^TX be invertible. However, when variables are highly correlated, X^TX approaches singularity, leading to highly unstable estimations of regression coefficients and sensitivity to small changes in the data. This is a major limitation of OLS. To address this issue, ridge regression adds a positive constant λ to the diagonal of X^TX, i.e.: , where I is the identity matrix.

[0060] As the ridge parameter λ increases, the ridge regression coefficient β shrinks towards zero. The trajectory of β as a function of λ is called the ridge trace. The value of λ is determined when the ridge trace stabilizes. This method alleviates multicollinearity and generally produces more reliable estimates than ordinary least squares. Based on this method, this invention performs regression-based interpolation in both spatiotemporal dimensions to fill in missing values ​​in the daily records of representative stations within each 5°×5° grid.

[0061] Spatial interpolation uses neighboring stations within the same grid as reference sequences for each representative station. 1) Calculate geospatial features: Calculate spherical distance metrics (longitude, latitude, and altitude differences) between the representative station and each reference station, then sort the reference stations in order of proximity; 2) Group data by month: To account for seasonal variations in temperature, dew point temperature, and local air pressure, data are grouped by variable and month; 3) Construct spatial models and interpolation: Ridge regression models are constructed for the same variables within the same month using data from the nearest reference station. Missing values ​​are imputed as much as possible to generate daily series with spatial gap filling for the representative station.

[0062] Subsequently, time interpolation was applied to each representative station, and the spatially interpolated series was further processed by utilizing the temporal relationships between variables. The data were again grouped by variable and month. Ridge regression was used to construct bivariate time regression models, specifically the relationships between daily temperature and dew point temperature, temperature and local air pressure, and dew point temperature and local air pressure within the same month. These models were used to impute the remaining missing values.

[0063] In spatiotemporal interpolation, data are divided into training and test sets in an 8:2 ratio. Model performance is evaluated using the coefficient of determination; only R-squared values ​​are accepted. 2 For models with a value >0.5, when multiple models are available for the same missing value, choose R. 2 The model with the highest score. The optimal ridge parameters were determined using a random forest method, and the model fit requires at least five years of data.

[0064] A comprehensive spatiotemporal interpolation framework was used to obtain a globally representative dataset consisting of daily records from 1,032 stations, characterized by its wide temporal coverage and high data quality. Given the improved data quality after 1973, this invention focuses on the period from 1973 to 2024. Stations with at least 35 years of valid data were selected for analysis, where a valid month was defined as having fewer than 7 days of missing data, and a valid year was defined as having fewer than 1 month of missing data. The spatial distribution of the selected stations is as follows: Figure 1 As shown.

[0065] Wet-bulb temperature (WBGT) is used as a quantitative measure of thermal comfort in ISO standards. The formula for calculating WBGT is as follows: (2) (3) (4) In the formula, T is the temperature (in °C), e is the vapor pressure (in hPa), and SP is the local atmospheric pressure (in hPa). If SP is missing, it is estimated using the sea level pressure using the formula: (5) H represents the station's altitude (in meters).

[0066] Annual averages of temperature, vapor pressure, and WBGT were calculated for each station. The stations were then gridded onto a uniform global 5°×5° grid by calculating the annual averages for all stations within each grid cell. For each grid cell, a decade-long linear trend over the study period (1973–2024) was calculated using ordinary least squares regression. The statistical significance of each trend was assessed using the nonparametric Mann-Kendall test, and only significant trends at the p<0.05 level were retained for further analysis to ensure robustness.

[0067] This invention employs an area-weighted average method to calculate the regional average WBGT time series. The steps include: (1) calculating the annual average WBGT for each station; (2) multiplying the value of each station by the cosine of its latitude and summing the weighted values ​​for all stations; (3) summing the latitudinal cosine values ​​for all stations; the regional WBGT time series is obtained by dividing the result of step (2) by the result of step (3). This method effectively considers the decrease in representativeness per unit area with increasing latitude, minimizing the spatial bias in constructing the regional average.

[0068] The time series of global mean wet-bulb black-sphere temperatures from 1973 to 2024 shows a significant warming trend and considerable interannual variability, such as... Figure 2 As shown, the 0.36°C increase per decade (p<0.05) indicates a significant intensification of global heat stress over the past fifty years. The ten-year moving average smooths out interannual fluctuations, revealing more clearly the persistence and acceleration of warming. The annual global WBGT values ​​are relatively concentrated, with a global average of 0.098°C, ranging from -1.15°C to 1.36°C, and a standard deviation of 0.60°C. This narrow distribution suggests that despite a clear long-term upward trend, heat stress remains within a relatively consistent range. Significant interannual variations exist, with interannual differences ranging from -0.70°C to 0.73°C and a standard deviation of 0.27°C, reflecting significant variability between certain years, which may be related to internal climate modes.

[0069] Spatial and latitudinal analysis of the WBGT trend from 1973 to 2024 reveals significant spatial heterogeneity and strong latitudinal dependence in the intensification of global heat stress. Figure 3As shown, the vast majority of the global region (82.6% of the grids) experienced a statistically significant increase in WBGT, particularly in the tropics and high-latitude continental regions of the Northern Hemisphere. Insignificant increases were observed in 11.6% of the grids, primarily in the oceanic and coastal zones. In contrast, WBGT declines were less common, with insignificant and significant declines covering only 4.0% and 1.8% of the global total, respectively, and mainly confined to areas around 30°S. Latitudinally, the increase in global WBGT increased from south to north. The Northern Hemisphere exhibited considerable variability. The Arctic region (>70°N) showed the strongest warming, with a median trend exceeding 0.40°C / decade, consistent with the polar amplification effect. Variation was more limited in the Southern Hemisphere, highlighting the moderating role of maritime climate. These results emphasize that heat stress is most severe in the tropics and Arctic, indicating an urgent need for latitude-specific public health and adaptation strategies.

[0070] Systemic restructuring of global heat stress risk: The global upward trend in WBGT has led to significant regional variations in the risk of different levels of heat stress. Figure 4 This section analyzes the long-term trends in the annual frequency of days with temperatures below 25°C (low heat stress), days between 25 and 30°C (moderate to high heat stress), and days above 30°C (high to extremely high heat stress) as defined by the WBGT, from 1973 to 2024. The results reveal a coordinated shift in global heat conditions characterized by a systematic reduction in the number of low-stress days and a significant redistribution of heat exposure to more dangerous conditions.

[0071] A widespread and significant reduction in the number of days with WBGT below 25°C was observed, particularly across the mid-latitudes. This contraction of low heat stress conditions was most pronounced in Europe, East Asia, southeastern North America, southern Africa, and parts of South America and Australia. Latitudinal analysis revealed the most significant losses occurring between 10°N and 50°N and between 10°S and 30°S. This trend implies a shortening of the duration of the thermal comfort period, which traditionally buffers against seasonal heat.

[0072] The trends within the range of moderate to high heat stress show a clear latitudinal differentiation. A significant increase was detected between 30°N and 60°N and south of 20°S, particularly in the northern continental regions, indicating that the period of moderate to high heat stress is expanding and intensifying towards the poles. Conversely, a significant decrease was generally observed in the core tropical and subtropical regions (20°S to 30°N). This pattern suggests a thermodynamic shift in the tropics, with an increasing number of days within this moderate range transitioning to higher risk categories (≥30°C), resulting in a net decrease.

[0073] The number of days with WBGT exceeding 30°C has shown a widespread and statistically significant increase globally, with the strongest trend concentrated in tropical and subtropical regions. The largest increases have occurred in already vulnerable areas such as South Asia, the Arabian Peninsula, the Sahara, northern Australia, and the Caribbean. Latitudinal analysis confirms that the strongest increases are concentrated between 30°N and 20°S, highlighting the accelerating risk of extreme heat stress in areas with high background humidity and temperature.

[0074] Synergistic trends across these three WBGT thresholds depict a systemic, large-scale reorganization of the global heat stress pattern. Tropical regions are experiencing a shift directly converting days of moderate stress risk into days of extreme heat risk, while mid-latitudes bear a double burden: a significant decrease in low-stress-risk days and a substantial increase in moderate-stress days. Simultaneously, probability density functions show a significant shift from low-risk to high-risk categories. The changing characteristics of different risks are more pronounced in the Northern Hemisphere than in the Southern Hemisphere. Box plots show the median and mean trends for each latitudinal zone, revealing significant hemispheric asymmetry, consistent with established Arctic amplification effects and different climate response patterns across hemispheres. This latitudinal redistribution highlights the pervasive increase in heat risk and its polar expansion, underscoring the urgent need for tailored public health strategies and adaptation policies to specific trajectories of heat exposure changes across different geographic regions.

[0075] It should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of the present invention without departing from the spirit and scope of the technical solutions of the present invention, and all such modifications or substitutions should be covered within the scope of the claims of the present invention.

Claims

1. A method for reconstructing long-series global climate data for heat stress risk assessment, characterized in that: For missing meteorological data at the target station, spatial interpolation is first performed using a ridge regression model based on spatially adjacent reference station data, and then temporal interpolation is performed using a ridge regression model based on the physical relationships between different meteorological variables at the target station.

2. The method for reconstructing long-series global climate data for heat stress risk assessment according to claim 1, characterized in that, The spatial interpolation includes the following steps: grouping the daily-scale data of the target station and reference stations by variable and by month; using the effective observation sequence of the target variable for that month at the target station as the dependent variable Y, selecting the K highest-ranked reference stations' data of the same period and the same variable to form the independent variable matrix X, and constructing a ridge regression model. The regression coefficients of the ridge regression model are calculated using the following formula: ; Where λ is the regularization parameter (ridge parameter) and I is the identity matrix; the constructed ridge regression model is used to imputate and predict the missing values ​​of the target variable for the target station in that month.

3. The method for reconstructing long-series global climate data for heat stress risk assessment according to claim 1 or 2, characterized in that, The time interpolation includes the following steps: grouping the spatially interpolated multivariate daily-scale data of the target station by month, wherein the multivariate variables include at least a target variable and key covariates; for each month group, constructing a bivariate ridge regression model using the physical relationships between variables, wherein the bivariate ridge regression model uses the target variable as the dependent variable and at least one other meteorological variable as the independent variable, and its regression coefficients are... The calculation formula is: ; in For the sequence of independent variables, For the dependent variable sequence, For regularization parameters, The identity matrix is ​​used; the constructed bivariate ridge regression model is used to perform secondary imputation on the target variable values ​​that are still missing after spatial imputation.

4. The method for reconstructing long-series global climate data for heat stress risk assessment according to claim 1, characterized in that, It also includes a model selection step: calculating the determination coefficients R for each ridge regression model used for spatial and temporal interpolation, respectively. 2 ; Set R 2 The minimum performance threshold, using only R 2 The output value of the model is higher than the preset threshold; when multiple qualified model output values ​​exist for the same missing data point, R is selected. 2 The highest output value is used as the final interpolation value for the missing data point.

5. The method for reconstructing long-series global climate data for heat stress risk assessment according to claim 1, characterized in that, The ridge parameter λ in the ridge regression model is determined automatically and adaptively using the random forest algorithm.

6. The method for reconstructing long-series global climate data for heat stress risk assessment according to claim 1, characterized in that, Before spatial interpolation, a data preprocessing step is also included: mapping all global meteorological stations to a regular latitude and longitude grid with a preset resolution; within each grid cell, a representative station is selected based on a preset comprehensive data integrity index, which is the number of days on which the target interpolation variable and the key covariate have valid observations on the same calendar date, and is weighted and scored in combination with the total number of valid record years of the station, and the station with the highest score is selected as the representative station of the grid.

7. The method for reconstructing long-series global climate data for heat stress risk assessment according to claim 6, characterized in that, The regular latitude and longitude grid is a global regular grid of 5°×5°.

8. The method for reconstructing long-series global climate data for heat stress risk assessment according to claim 3, characterized in that, The other meteorological variables include dew point temperature and local air pressure, and the key covariates include daily average dew point temperature and daily average local air pressure.

9. The method for reconstructing long-series global climate data for heat stress risk assessment according to claim 2, characterized in that, The steps for constructing a set of space reference stations are as follows: taking the grid where the target station is located as the center, all other stations in that grid and several geographically adjacent grids are initially selected as the space reference station candidate set; the comprehensive geographical distance index between the target station and each candidate reference station is calculated, and the comprehensive geographical distance index is calculated based on spherical distance, latitude and longitude difference and altitude difference; all candidate reference stations are sorted from smallest to largest according to the comprehensive geographical distance index, and the top K reference stations are selected to form the space reference station set, where K is set to 5-10 according to the data density.

10. The method for reconstructing long-series global climate data for heat stress risk assessment according to claim 2 or 3, characterized in that, When constructing a ridge regression model, the valid data for the corresponding month are divided into a training set and a validation set in an 8:2 ratio. The model is trained using the training set and its performance is validated using the validation set.

Images