A land water storage downscaling method and system based on monthly self-adaptation

Abstract

The land water storage anomaly downscaling method disclosed by the application can effectively generate high-precision monthly land water storage anomaly prediction results with significantly improved spatial resolution. The method intelligently integrates the advantages of different machine learning models, significantly enhances the prediction stability and adaptability in complex seasonal changes and variable geographical environments, and overcomes the limitations of insufficient fitting ability of traditional single models. The refined prediction data obtained finally can more accurately depict the local details and spatio-temporal patterns of water storage changes, providing more reliable data support for regional water resource accurate assessment, hydrological cycle process in-depth understanding, and timely warning of drought and flood disaster risks, and has important scientific value and wide application prospect.

Description

Technical Field

[0001] This invention relates to the interdisciplinary fields of hydrology, spatiotemporal data analysis, machine learning and hydrological remote sensing, and specifically to a method and system for downscaling terrestrial water storage based on monthly adaptive data. Background Technology

[0002] This patent relates to interdisciplinary technologies such as satellite hydrology, spatial data analysis, and machine learning, and in particular to a high-resolution downscaling method for terrestrial water storage based on a monthly adaptive hybrid model. Currently, the terrestrial water storage change observations provided by GRACE and its follow-up missions have become an important data source for global water resource monitoring, but their original spatial resolution is relatively coarse (approximately 400 kilometers), making it difficult to directly support refined water resource management at regional or local scales.

[0003] Although downscaling methods based on statistical relationships or multi-source data fusion have been extensively studied, most existing methods employ a uniform model architecture and fixed feature sets, failing to adequately consider the strong seasonal variations and monthly specificities in terrestrial water storage. Because the dominant hydrological and meteorological processes (such as summer evapotranspiration, winter freeze-thaw cycles, and rainy season runoff) differ significantly across months, using a single model often struggles to achieve consistently high accuracy throughout the year, especially during transitional seasons or extreme months where performance tends to degrade. Furthermore, traditional methods lack targeted strategies for handling the spatiotemporal heterogeneity, outliers, and modeling challenges related to small sample months, limiting the model's robustness and generalization ability.

[0004] Therefore, there is an urgent need to develop a downscaling framework that can adapt to monthly characteristics, dynamically select the optimal model architecture, and integrate targeted feature engineering and data optimization strategies to generate high-resolution, high-reliability monthly terrestrial water storage sequences, supporting refined water cycle research and water resource management. Summary of the Invention

[0005] To address the technical problems mentioned above, this invention provides a monthly adaptive method for downscaling terrestrial water storage, comprising the following steps: S1. Obtain monthly GRACE land water storage anomaly data and several high-resolution environmental prediction factor data for the study area within a preset time period. S2. Based on the monthly GRACE land water storage anomaly data and environmental prediction factor data, select the monthly feature factor set; before selection, the environmental prediction factor data needs to be resampled to the same resolution as the GRACE data. S3. Based on the monthly feature factor set and the corresponding monthly GRACE land water storage anomaly data, a hybrid downscaling model is constructed using a monthly independent modeling strategy. S4. Input the high spatial resolution environmental predictor data into the trained monthly hybrid downscaling model and output the high spatial resolution monthly land water storage anomaly prediction results.

[0006] Preferably, S1 includes: GRACE / GRACE-FO Mascon data from different institutions were acquired, and the average value was calculated to represent the terrestrial water storage anomaly in the study area. The acquired data were interpolated using singular spectrum interpolation to impute missing months.

[0007] Preferably, in step S2, the contribution of environmental predictive factor data is calculated monthly using partial least squares regression, and a monthly feature factor set is selected. The steps include: The environmental predictor data and the corresponding monthly GRACE terrestrial water storage anomaly data are input into the partial least squares regression model, and the variable projection importance index value of each environmental predictor data is calculated month by month; predictors with variable projection importance index values ​​greater than 0.8 are selected to form the feature factor set for that month.

[0008] Preferably, in step S3, for months in which the prediction performance of the CNN-LSTM model meets the preset requirements, the CNN-LSTM model is used for modeling to extract spatiotemporal features; for months in which the prediction performance of the CNN-LSTM model does not meet the preset requirements, the gradient boosting decision tree model is switched to be used for modeling to enhance the nonlinear fitting ability.

[0009] Preferably, S3 further includes: First, a CNN-LSTM model is used to initially model all months. Based on the modeling results, the months are divided into different performance levels. For months with good performance, the CNN-LSTM model is continued to be used. For months with poor performance, the gradient boosting decision tree model is used instead. For months with extremely poor performance, an enhanced gradient boosting decision tree model is used.

[0010] Preferably, S3 further includes: during model training, adopting differentiated dataset partitioning strategies, feature engineering processes, and data augmentation strategies for different months; and using dynamically changing learning rates and early stopping mechanisms to prevent overfitting.

[0011] Preferably, step S3 further includes: using test set data to evaluate the performance of the trained monthly downscaling model; and using the coefficient of determination R², root mean square error RMSE, and correlation coefficient CC to quantify the accuracy of the downscaling results.

[0012] The present invention also provides a monthly adaptive land water storage downscaling system, the system being used to implement the above method, comprising: a data acquisition module, a filtering module, a construction module, and an output module; The acquisition module is used to acquire monthly GRACE land water storage anomaly data and several high-resolution environmental prediction factor data of the study area within a preset time period. The filtering module is used to filter out monthly feature factor sets based on monthly GRACE land water storage anomaly data and environmental prediction factor data; before filtering, the environmental factor data needs to be resampled to the same resolution as the GRACE data. The construction module is used to build a hybrid downscaling model based on the monthly feature factor set and the corresponding monthly GRACE land water storage anomaly data, using a monthly independent modeling strategy. The output module is used to input high spatial resolution environmental prediction factor data into the trained monthly hybrid downscaling model and output high spatial resolution monthly land water storage anomaly prediction results.

[0013] Compared with the prior art, the beneficial effects of the present invention are as follows: This invention significantly improves the overall accuracy and stability of spatial downscaling of terrestrial water storage anomalies, effectively overcoming the shortcomings of traditional single models in fitting complex seasonal changes. The resulting high spatial resolution predictions are more refined and reliable, clearly revealing local details of water storage changes, providing a more solid data foundation for precise regional water resource assessment and management, hydrological process mechanism research, and early warning of extreme droughts and floods. Furthermore, by intelligently integrating the advantages of different models, this method enhances its adaptability and robustness under various geographical and climatic conditions, demonstrating significant practical application value. Attached Figure Description

[0014] To more clearly illustrate the technical solution of the present invention, the drawings used in the embodiments are briefly introduced below. 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.

[0015] Figure 1 This is a schematic diagram of the method flow according to an embodiment of the present invention; Figure 2 This is a schematic diagram of the VIP values ​​of the study area in an embodiment of the present invention; Figure 2 In the table, (a) represents the VIP value in February; (b) represents the VIP value in May; (c) represents the VIP value in August; and (d) represents the VIP value in November. Figure 3 This is a time series diagram of monthly land water storage anomalies simulated by a monthly adaptive downscaling model in the study area of ​​this invention embodiment; Figure 4This is a spatial comparison map of terrestrial water storage anomalies before and after downscaling during typical periods in the study area of ​​this invention; wherein, Figure 4 (a) shows the comparison before and after downscaling in November 2003; (b) shows the comparison before and after downscaling in July 2010; (c) shows the comparison before and after downscaling in January 2016; and (d) shows the comparison before and after downscaling in August 2022. Detailed Implementation

[0016] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

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

[0018] Example 1 S1. Obtain monthly GRACE land water storage anomaly data and several high-resolution environmental prediction factor data for the study area within a preset time period.

[0019] The study acquired monthly GRACE land water storage anomaly data and various high-resolution environmental prediction factor data for the study area within a preset time period, and spatially resampled all data to a uniform resolution of 0.5°.

[0020] In this embodiment, terrestrial water storage anomaly information from GRACE data of multiple institutions within a preset historical period is extracted. Specifically, GRACE / GRACE-FO RL06.3 Mascon data published by CSR, GRACE / GRACE-FO RL06 v4 Mascon data published by JPL, and GRACE / GRACE-FO RL06 v2.0 Mascon data published by GSFC are used. To minimize uncertainty, this embodiment uses the average value of these three products to represent the terrestrial water storage anomaly (TWSA) in the study area, referred to as GRACE Average. The temporal resolution is monthly, and the spatial resolutions are 0.25°, 0.5°, and 0.5°, respectively. An area-weighted average method is used to resample the three data sets to 0.5°. Data are expressed in equivalent water column height (mm), and all grid values ​​represent surface quality deviation anomalies relative to the average baseline from January 2004 to December 2009. Due to problems with the GRACE sensor and the replacement of the GRACE-FO satellite, a total of 33 data periods have been lost since April 2002 (200206, 200207, 200306, 201101, 201106, 201205, 201210, 201303, 201308, 201309, 201402, 201407, 201412, 2). (01506, 201510, 201511, 201604, 201609, 201610, 201702, 201707, 201708, 201709, 201710, 201711, 201712, 201801, 201802, 201803, 201804, 201805, 201808, 201809). This embodiment uses Singular Spectral Interpolation (SSA) to impute missing data.

[0021] Singular Spectrum Analysis (SSA), as a non-parametric time series analysis method, exhibits unique advantages in multiple disciplines through its data-adaptive decomposition mechanism. The specific principles and methods of SSA are as follows: Given a discrete-time sequence X=[x1,x2,...,x] sampled at equal intervals. n ], where n represents the sequence length. Based on this sequence, its trajectory matrix Y can be constructed using the time-delay embedding method. The specific form is shown in the following formula. Each column of matrix Y consists of M consecutive sampling points from the original sequence, thus completely preserving the temporal characteristics of the sequence.

[0022] Based on the construction principle of the trajectory matrix, its column dimension L can be represented by the sequence length N and the embedding dimension M as L = N - M + 1. This matrix has a unique Hankel structure property, where the values ​​of elements on all ascending diagonal lines remain constant. Within the classic SSA methodology framework, the standard processing flow proposed by Broomhead and King involves constructing the lag covariance matrix Yᵀ and performing principal component analysis (PCA). However, from a computational efficiency perspective, a better solution is to directly perform Singular Value Decomposition (SVD) on the trajectory matrix Y, the mathematical expression of which is shown in the following equation: Where M and L are matrix indices representing the matrix dimensions, A is a diagonal matrix, UV is an orthogonal matrix, and V is the eigenvector of the lag covariance matrix. Principal component analysis generates two key outputs: Empirical Orthogonal Functions (EOFs) characterizing spatial modes, and Principal Components (PCs) describing time coefficients. Their mathematical relationships are shown below: .

[0023] Each mode Zi has the same magnitude as Y. Therefore, a new time series can be represented by the average of the diagonal elements, called the reconstructed component (RCs, denoted as θ), as shown in the formula: In this equation, all elements in Zi satisfy j + k = p + 1, where the superscripts p, j, and k are position indices. The sum of all RCs equals the original input time series, as shown in the formula: .

[0024] RCs are sorted in descending order based on the magnitude of their corresponding singular values, which directly reflects the energy contribution of each component in the original signal. The first-order RCs typically contain long-term trend terms and significant periodic oscillation components, while the subsequent-order RCs mainly characterize noise information. Based on this characteristic, in practical applications, the total number of original components K is often not used directly. Instead, an appropriate truncation threshold K′ is set to selectively retain RCs with significant physical meaning, thereby achieving effective signal denoising and feature extraction.

[0025] This embodiment employs the iterative SSA imputation algorithm developed by Kondrashov and Ghil. Its core innovation lies in introducing a dual-loop mechanism into the traditional SSA framework: the inner loop optimizes data imputation by dynamically updating the gap estimate Xgap, terminating the iteration when Xgap converges to a stable state; the outer loop increases the complexity of sequence reconstruction by gradually increasing the number of modes K involved in the reconstruction. The algorithm initializes with Xgap=0 and retains the optimization parameters after each iteration for use in subsequent loops. Addressing the missing data characteristics of GRACE satellite data, the imputation problem is divided into two difficulty levels: SSA-filling-a handles discontinuous missing data within the GRACE mission period, while SSA-filling-b specifically addresses the 11-month data connection issue between the GRACE and GRACE-FO missions. This categorized approach effectively reduces the uncertainty caused by different missing data types.

[0026] S2. Based on monthly GRACE land water storage anomaly data and environmental prediction factor data, select monthly characteristic factor sets.

[0027] In this embodiment, before screening, the environmental factor data needs to be resampled to the same resolution as the GRACE data. Then, the partial least squares regression method is used to calculate the contribution of each environmental prediction factor to the anomaly of terrestrial water storage on a monthly basis, and the monthly feature factor set is screened according to the preset contribution threshold (VIP>0.8).

[0028] Monthly hydrological data from ERA5-Land with a spatial resolution of 0.1°, collected concurrently with the GRACE gravity satellite data, were selected in this embodiment. Eleven environmental predictive factors with high correlation to terrestrial water storage were chosen: snow depth, surface thermal radiation, soil temperature, total precipitation, temperature, total evaporation, runoff, soil volumetric water content, bare soil evaporation, vegetation transpiration, and normalized difference vegetation index (NDVI). These data were resampled to a 0.5° grid, consistent with the spatial resolution of the GRACE Average data. Specifically, snow depth, surface thermal radiation, soil temperature, total precipitation, temperature, total evaporation, runoff, soil volumetric water content, bare soil evaporation, and vegetation transpiration were from the ERA5-Land dataset; while the selected NDVI data came from MOD13C2 version 6, with a monthly scale, a spatial resolution of 0.05°, and a time span from January 2003 to December 2024. MODIS NDVI data is global monthly vegetation index data provided by the Moderate Resolution Imaging Spectroradiometer (MODIS). It is an important indicator characterizing vegetation growth and biomass, with values ​​typically ranging from -1 to 1. Values ​​close to 0 or negative generally represent non-vegetated areas, while values ​​close to 1 indicate areas with dense vegetation cover. This data is generated based on a 16-day composite algorithm, which effectively reduces cloud interference and provides clear monthly vegetation images, thus having wide applications in climate change monitoring, agricultural production assessment, and ecological research.

[0029] The processed environmental predictor factors and the corresponding monthly Grace Average land water storage anomaly data are input into the partial least squares regression model. The variable projection importance index (VIP value) of each factor is calculated monthly. The principle of partial least squares regression is as follows: The mathematical principle of partial least squares regression (PLSR) can be stated as follows: Given the original feature matrix X0 (n×m dimensional) and response matrix Y0 (n×n dimensional), X and Y are obtained after standardization. The algorithm first solves for the first pair of principal component axis vectors. a 1 (m×1 dimension) and b 1 (n×1 dimension), thus obtaining the first pair of principal components of X and Y. u 1 and v 1. The extraction of this pair of principal components follows the principle of maximizing covariance, reflecting the characteristics of PLSR in integrating principal component analysis and canonical correlation analysis.

[0030] It should be noted that the projection vector defined here... a 1 and b 1 represents parameters to be determined under the unit norm constraint, and their specific values ​​need to be obtained through subsequent optimization algorithms.

[0031] The core algorithm of Canonical Correlation Analysis (CCA) lies in finding the optimal projection direction between two sets of variables through linear transformation, thereby maximizing the correlation coefficient between the extracted canonical variables. .

[0032] The core algorithm principle of Principal Component Analysis (PCA) is to find a new set of coordinate axes through orthogonal transformation, so that the projection variance along each principal component direction is maximized sequentially, that is: ; Combining the ideas of CCA and PCA yields the solution approach for PLSR, namely: .

[0033] In mathematical terms, this can be expressed as follows: Maximize<Xa1,Yb1> ; .

[0034] for a 1 and b To solve for 1, this embodiment introduces Lagrange multipliers: ; To each a 1 and b 1. Taking the partial derivative, we get: ; ; After sorting, we can obtain: ; ; objective function The requirements are the highest.

[0035] Eigenvalue decomposition can be used to determine that the eigenvector corresponding to the largest eigenvalue of the symmetric matrix XTYYTX is the desired eigenvector. a 1. Similarly, the leading eigenvector of matrix YTXXTY is b 1. Based on this solution, the first pair of canonical correlation components of X and Y can be calculated. u 1, v 1. This pair of components has a maximized cross-covariance property, forming the most significant correlation pattern between the two sets of variables.

[0036] ; .

[0037] Based on the principal component regression approach, X and Y can be correlated with their respective principal components. u 1, v 1. Establish a regression model: It should be noted that, although the load vector p 1, q 1 and projection vector a 1, b While they differ in their mathematical definitions, they are interconnected through specific transformation relationships. Based on the first pair of principal components... u 1 and v The linear mapping relationship between 1 and 2 can be expressed as the principal components of the explanation matrix X. u The function form of 1, and the residual terms generated during this transformation are recorded in residual matrices E1 and E2 respectively: Calculated using the least squares method p 1, q 1, r 1: From the results obtained above, we can deduce that... a 1, p The relationship between 1 and 2 is: Based on the above analysis, a 1 represents u 1 is the projection direction vector in X space, and p 1 represents the coefficient vector that minimizes the residual matrix E1 obtained through least squares regression; the two have different mathematical meanings. During the algorithm iteration, the current unexplained X residual E1 is used as the X input for the next round, while the Y residual E2 is used as the new Y input, and the above regression steps are repeated. This process continues until any of the following termination conditions are met: (1) The residual E2 reaches the preset accuracy threshold; (2) The number of principal components extracted reaches the upper limit determined by the rank of the X matrix.

[0038] Assuming the final algorithm extractsk If there are principal components, then a series of vectors can be represented as a 1, a 2, ..., ak ; u 1, u 2, ..., UK ; v 1, v 2, ..., vk ; p 1, p 2, ..., pk ; r 1, r 2, ..., r k, where UI , uj ( i ≠ j They are mutually orthogonal. ai , aj ( i ≠ j It is also orthogonal, but pi , pj They are generally not orthogonal. X and Y It can ultimately be expressed as: ; .

[0039] use (at this time i=j ), (At this point, i≠j) The relationship can be expressed in matrix form as follows: ; .

[0040] That is, the regression equation X→Y, where D=ART. During the algorithm process, the calculated values ​​of A and R are collected and then PLSR can be used for prediction.

[0041] Therefore, for a newly input data x First, we use A to calculate each principal component, that is: Then substitute: Or directly substitute: The predicted value of vector y can then be obtained.

[0042] Predictors with a projected importance index value greater than 0.8 are selected to form the feature factor set for the month.

[0043] S3. Based on the monthly feature factor set and the corresponding monthly GRACE land water storage anomaly data, a hybrid downscaling model is constructed using a monthly independent modeling strategy.

[0044] Based on the monthly feature factor set and the corresponding monthly GRACE land water storage anomaly data, a hybrid downscaling model is constructed using a monthly independent modeling strategy. Specifically, for months in which the prediction performance of the CNN-LSTM model (convolutional neural network-long short-term memory network model) meets the preset requirements, this model is used to extract spatiotemporal features. For months in which the prediction performance of the CNN-LSTM model does not meet the preset requirements, the gradient boosting decision tree model (GBDT) is switched to enhance the nonlinear fitting ability.

[0045] First, a CNN-LSTM model is used to initially model all months. The modeling results identify poorly performing months, and the 12 months are then divided into three performance levels: months with good performance continue using the CNN-LSTM model, months with poor performance are replaced with a Gradient Boosting Tree (GBDT) model, and months with extremely poor performance are treated with an Enhanced Gradient Boosting Tree (EBR) model. During training, the system employs a refined hierarchical training strategy, performing differentiated preprocessing procedures for different months. Regarding dataset partitioning, month-specific filtering is performed based on the month labels of the original data. Separate training and test sets are constructed for each month, with the test set proportion dynamically adjusted according to the month's performance level (10% for poorly performing months and 20% for normally performing months). A stratified sampling strategy based on the quintiles of the target variable is used to ensure consistent data distribution. For feature and model settings, the system configures specific feature engineering processes for different months. For low-performing months, data augmentation strategies including supplementing with data from neighboring months, cyclic sampling, and controlled noise injection are implemented, combined with progressively relaxed outlier handling thresholds to optimize data quality.

[0046] During model training, the GBDT model captures complex nonlinear relationships through a large number of trees (400-700 trees) and a relatively deep tree structure (7-10 layers), and employs weak regularization to avoid overfitting. The CNN-LSTM model extracts local features through convolutional layers, captures temporal dependencies through LSTM layers, and integrates multi-scale information by combining spatial feature fusion layers. Although the model training results may vary due to random noise in data augmentation, training set partitioning, network weight initialization, and random subsampling by GBDT, the system ensures the reproducibility of the results by fixing the random seed (e.g., using the month as the random_state), controlling the noise amplitude, and applying an early stopping mechanism. In the prediction phase, the saved model is used to perform deterministic forward propagation, resulting in consistent outputs under the same input. The entire system ultimately generates high-precision spatiotemporal prediction results with a resolution of 0.1°, significantly improving the prediction accuracy, stability, and reliability in complex seasonal variation scenarios.

[0047] During model training, a dynamically varying learning rate is used, ranging from 3×10⁻⁵ to 1×10⁻⁴. An early stopping mechanism is employed to prevent overfitting, with the number of early stopping rounds set between 12 and 40 rounds based on the characteristics of each month. Depending on the differences in data characteristics across months, either the Huber loss function or the SmoothL1 loss function is selected as the optimization objective for model training.

[0048] In this embodiment, the model performance evaluation includes: The performance of the monthly downscaling model was evaluated using test set data. Before evaluation, the model output was inversely standardized. The accuracy of the downscaling results was quantitatively analyzed using three indicators: coefficient of determination R2, root mean square error RMSE, and correlation coefficient CC.

[0049] (1) Coefficient of determination R2 The coefficient of determination represents the proportion of the total variation in observed values ​​that can be explained by the model's simulation. It ranges from [0, 1]. The closer the value is to 1, the better the model's performance and the more reliable the predictions. The mathematical expression is as follows: .

[0050] (2) Correlation coefficient CC The correlation coefficient measures the strength and direction of the linear relationship between two variables, reflecting the tightness of the linear association between simulated values ​​and observed values. Its value ranges from -1 to 1; the closer the absolute value is to 1, the stronger the linear correlation. A positive value indicates a positive correlation, and a negative value indicates a negative correlation. The mathematical expression is as follows: .

[0051] (3) Root Mean Square Error (RMSE) The root mean square error (RMSE) represents the "average error magnitude" between the simulated and observed values. A smaller RMSE is better, ideally approaching 0 (no error). Its mathematical expression is as follows: in, Xreal , i Indicates the first i Groundwater observation data from each station, This represents the average value of groundwater observation data from all stations. Xpre , i The model simulation represents the first i Groundwater data from each station, This represents the average groundwater data from all sites simulated by the model. n The number of observation stations.

[0052] S4. Input the high spatial resolution environmental predictor data into the trained monthly hybrid downscaling model and output the high spatial resolution monthly land water storage anomaly prediction results.

[0053] The environmental predictor data with a high spatial resolution of 0.1° is input into the trained monthly hybrid downscaling model. After inverse standardization, the monthly land water storage anomaly prediction results with a high spatial resolution of 0.1° are output.

[0054] Example 2 The following will describe in detail, with reference to this embodiment, how the present invention solves the technical problems in practical work, and the process is as follows: Figure 1 As shown.

[0055] Taking a region in my country as an example, GRACE data, ERA5-Land monthly hydrological data, and Normalized Difference Vegetation Index (NDVI) data from January 2003 to December 2024 were used. ERA5-Land is a global land surface reanalysis dataset developed by the European Centre for Medium-Range Weather Forecasts (ECMWF). This dataset provides high spatiotemporal resolution meteorological data covering the global land surface since 1981, including 35 variables such as surface temperature, bare soil evaporation, canopy top evaporation, vegetation evaporation, and runoff. NDVI data comes from MOD13C2 version 6, with a monthly scale, a spatial resolution of 0.05°, and a time span from January 2003 to December 2024. MODIS NDVI data is global monthly vegetation index data provided by the Moderate Resolution Imaging Spectroradiometer (MODIS), an important indicator characterizing vegetation growth and biomass, with values ​​typically ranging from -1 to 1. Values ​​close to 0 or negative generally represent non-vegetated areas, while values ​​close to 1 indicate areas with dense vegetation cover. This data is generated based on a 16-day composite algorithm, which effectively reduces cloud interference and provides clear monthly vegetation images, thus having wide applications in climate change monitoring, agricultural production assessment, and ecological research.

[0056] This embodiment uses the average value of the three products mentioned in Embodiment 1 to represent the terrestrial water storage anomaly (TWSA) in the study area, referred to as GRACE Average. The temporal resolution is monthly, and the spatial resolutions are 0.25°, 0.5°, and 0.5°, respectively. These three data points are resampled to 0.5°. The data are expressed in equivalent water column height (mm), and all grid values ​​represent surface quality deviation anomalies relative to the average baseline from January 2004 to December 2009. Due to GRACE sensor problems and the replacement of the GRACE-FO satellite, a total of 31 data periods were lost between January 2003 and December 2024 (200306, 201101, 201106, 201205, 201210, 201303, 201308, 201309, 201402, 201407, 201412, 2015). (06, 201510, 201511, 201604, 201609, 201610, 201702, 201707, 201708, 201709, 201710, 201711, 201712, 201801, 201802, 201803, 201804, 201805, 201808, 201809). This study uses Singular Spectral Interpolation (SSA) to interpolate missing data to obtain complete satellite data for the period from January 2003 to December 2024.

[0057] Furthermore, partial least squares regression was used to calculate the contribution of each environmental predictor to the anomaly in terrestrial water storage on a monthly basis. Based on a preset contribution threshold (VIP > 0.8), a set of monthly characteristic factors was selected, including: Monthly hydrological data from ERA5-Land with a spatial resolution of 0.1° from January 2003 to December 2024 were selected. Based on literature review, 11 environmental predictor data with high correlation to terrestrial water storage were chosen: snow depth, surface thermal radiation, soil temperature, total precipitation, temperature, total evaporation, runoff, soil volumetric water content, bare soil evaporation, vegetation transpiration, and NDVI. These data were resampled to a 0.5° grid with the same spatial resolution as the GRACE Average data. Snow depth, surface thermal radiation, soil temperature, total precipitation, temperature, total evaporation, runoff, soil volumetric water content, bare soil evaporation, and vegetation transpiration were from the ERA5-Land dataset; while the selected NDVI data came from MOD13C2 version 6, with a monthly scale, a spatial resolution of 0.05°, and a time span from January 2003 to December 2024. MODIS NDVI data is global monthly vegetation index data provided by the Moderate Resolution Imaging Spectroradiometer (MODIS). It is an important indicator characterizing vegetation growth and biomass, with values ​​typically ranging from -1 to 1. Values ​​close to 0 or negative generally represent non-vegetated areas, while values ​​close to 1 indicate areas with dense vegetation cover. This data is generated based on a 16-day composite algorithm, which effectively reduces cloud interference and provides clear monthly vegetation images, thus having wide applications in climate change monitoring, agricultural production assessment, and ecological research.

[0058] The processed environmental predictive factors and the corresponding monthly Grace Average land water storage anomaly data were input into a partial least squares regression model. The variable projection importance index (VIP value) for each factor was calculated monthly. The results show that the dominant influencing factors differ significantly across months (Table 1). Figure 2 As shown, taking February, May, August and November as examples, the normalized vegetation index has the most significant impact in February and August (VIP=1.7, VIP=2.1), while temperature makes a significant contribution in May and November (VIP=1.9, VIP=1.8).

[0059] Table 1 In Table 1, the feature factor set for each month is arranged from largest to smallest according to the VIP value.

[0060] Based on the screening criterion of VIP>0.8, the optimal combination of predictive factors for each month was finally determined, as shown in Table 2. This result indicates that the driving mechanism of changes in terrestrial water storage has obvious seasonal characteristics.

[0061] Table 2 .

[0062] Furthermore, a hybrid downscaling model is constructed using a monthly independent modeling strategy, including: First, a CNN-LSTM model is used to initially model all months. The criteria for judging months with good performance are as follows: (1) Good performance: R²≥0.55 indicates that the model can explain more than half of the data variation, CC≥0.75 indicates strong linear correlation, and RMSE≤75 indicates that the error is within an acceptable range. The data in Table 3 for January (R²=0.54, close to the threshold), April (R²=0.64), and June (R²=0.50, close to the threshold) meet this characteristic.

[0063] (2) Poor performance: R² is between 0.40 and 0.55, indicating that the model has some explanatory power but is insufficient; CC is between 0.65 and 0.75, indicating moderate linear correlation; RMSE is between 75 and 90, indicating significant error but not extreme. February, May, July, August, September, and October in Table 3 meet this characteristic.

[0064] (3) Extremely poor performance: R²≤0.40 indicates that the model has very weak explanatory power, CC≤0.65 indicates poor linear correlation, and RMSE≥90 indicates that the prediction error is very large. In Table 3, March (R²=0.37), November (R²=0.26), and December (R²=0.37) meet this characteristic, especially November, which has the largest improvement, confirming the judgment that its initial performance was extremely poor.

[0065] During training, the system employs a refined hierarchical training strategy, performing differentiated preprocessing procedures for different months. Regarding dataset partitioning, month-specific filtering is performed based on the month labels of the original data, constructing separate training and test sets for each month. The proportion of the test set is dynamically adjusted according to the month's performance level (10% for poor-performing months and 20% for normal-performing months), and a stratified sampling strategy based on the quintiles of the target variable is used to ensure consistent data distribution. In terms of feature and model settings, the system configures specific feature engineering processes for different months. For low-performing months, data augmentation strategies including supplementing with data from neighboring months, cyclic sampling, and controlled noise injection are implemented, combined with progressively relaxed outlier handling thresholds to optimize data quality. During model training, the GBDT model captures complex nonlinear relationships through a large number of trees (400-700 trees) and a relatively deep tree structure (7-10 layers), and uses weak regularization to avoid overfitting; the CNN-LSTM model extracts local features through convolutional layers, captures temporal dependencies through LSTM layers, and integrates multi-scale information by combining spatial feature fusion layers.

[0066] Although model training results may vary due to random noise in data augmentation, training set partitioning, network weight initialization, and GBDT random subsampling, the system ensures reproducibility by fixing the random seed, controlling noise amplitude, and applying an early stopping mechanism. During the prediction phase, deterministic forward propagation is performed using the saved model, resulting in consistent outputs under the same input. The entire system ultimately generates high-precision spatiotemporal predictions with a resolution of 0.1°, significantly improving prediction accuracy, stability, and reliability in complex seasonal scenarios.

[0067] During model training, a dynamically varying learning rate is used, ranging from 3×10⁻⁵ to 1×10⁻⁴, and an early stopping mechanism is employed to prevent overfitting. The number of patience rounds for the early stopping mechanism is set between 12 and 40 rounds based on the characteristics of the month. Depending on the differences in data characteristics of the month, either the Huber loss function or the SmoothL1 loss function is selected as the optimization objective for model training.

[0068] Table 3 Table 4 .

[0069] From the model performance evaluation results, the monthly adaptive downscaling model achieved a comprehensive improvement in downscaling performance compared to the CNN-LSTM model (as shown in Tables 3 and 4): the R² of the monthly adaptive downscaling model on the test set was higher than that of the CNN-LSTM model in all months, with the most significant improvement in November, from 0.26 to 0.86. The R² improvements in February, August, and September also exceeded 0.2, and the R² exceeded 0.6 in most months, reaching a good fit level. Regarding the Pearson correlation coefficient (CC), which reflects the consistency of spatial trends, the monthly adaptive downscaling model showed a significant improvement in performance compared to the CNN-LSTM model. The monthly adaptive downscaling model achieved a CC score exceeding 0.8 on all test sets, with November reaching a strong correlation level of 0.96, significantly outperforming the optimal value (0.85) of the CNN-LSTM model. This indicates that its predictions have a better linear correlation with the true values. In terms of the root mean square error (RMSE), which measures prediction accuracy, the monthly adaptive downscaling model achieved varying degrees of reduction in RMSE on all test sets. The RMSE in November decreased from 97.80 to 39.54, while the RMSE in February, August, and other months decreased to around 50, indicating a significant reduction in prediction error.

[0070] The monthly adaptive downscaling model achieves a comprehensive improvement in the prediction performance of data for each month. It constructs a hierarchical strategy of "monthly feature adaptation + model optimization + full-process fine-tuning": First, for months adapted to CNN-LSTM, its spatiotemporal feature modeling advantages are retained, while adaptively adjusting network depth (deeper convolutional / LSTM layers), hierarchical sampling, and robust standardization improve data quality. This is combined with month-specific hyperparameters (learning rate, batch size), learning rate scheduling, and early stopping mechanisms to amplify the modeling capabilities of CNN-LSTM. For months where CNN-LSTM performs poorly, it switches to gradient boosting trees, which are more suitable for small samples and more resistant to noise, and specifically optimizes the number and depth of trees. The system dynamically adjusts outlier detection thresholds and adds data augmentation strategies to address issues such as complex monthly data distribution. Simultaneously, seasonal feature engineering is integrated into all months, supplementing them with physically meaningful derived features (e.g., the difference between 2-meter air temperature and surface temperature represents the near-surface-surface heat exchange process, and the difference between precipitation and actual evapotranspiration represents the regional water balance). Combined with flexible outlier handling and dynamic threshold adjustment preprocessing techniques, the system aligns data input, model selection, and training strategies with the data patterns of each month. This not only further improves the adaptability of CNN-LSTM to each month but also compensates for the shortcomings of gradient boosting trees in adapting to each month, achieving overall optimization of prediction performance metrics for all months.

[0071] Figure 3 This paper presents a time series comparison between the 0.5°×0.5° spatial resolution GRACE monthly land water storage anomaly (Original data) from 2003 to 2024 and the 0.1°×0.1° high-resolution results obtained after downscaling using a monthly adaptive downscaling model. The two curves show a highly consistent upward and downward trend for most of the period, with good agreement on fluctuations between 2003 and 2005, changes around 2014, and the overall evolution pattern in the later period. This indicates that the monthly adaptive downscaling model can effectively reconstruct the temporal variation characteristics of land water storage. The coefficient of determination (R²) between the two is 0.80, the root mean square error (RMSE) is 43.07 mm, and the correlation coefficient (CC) is 0.94. These indicators collectively demonstrate that the downscaled results have high interpretability, small mean bias, and a very strong linear correlation with the original data. The results show that the monthly adaptive downscaling model has a significant advantage in improving the spatial resolution of land water storage data and can provide reliable high-resolution data support for subsequent research.

[0072] Furthermore, the high spatial resolution environmental predictor data is input into the trained monthly hybrid downscaling model, including: Monthly environmental predictors with a spatial resolution of 0.1° (i.e., snow depth, surface thermal radiation, soil temperature, total precipitation, temperature, total evaporation, runoff, soil volumetric water content, bare soil evaporation, vegetation transpiration, and NDVI data) are input into the corresponding monthly downscaling model. The model outputs monthly terrestrial water storage anomaly raster data with the same spatial resolution (0.1°). The original resolution of the NDVI data is 0.05°; therefore, the NDVI data is first resampled to 0.1° to maintain consistency with other 0.1° high-resolution environmental predictors before being input into the downscaling model to predict 0.1° terrestrial water storage anomaly data.

[0073] Given the large time span and significant seasonality of the study, to avoid repetitive analysis across the entire time period, this invention focuses on comparing typical months with significant model differences and key hydrological implications. Based on time series (… Figure 3 Based on the characteristics of GRACE inversion, the spatial distribution characteristics of the original terrestrial water storage anomaly (TWSA) data and the downscaling results based on the monthly adaptive downscaling model were compared and analyzed in November 2003, July 2010, January 2016 and August 2022 as representative periods.

[0074] In November 2003, the initial period of the study, the TWSA (Transient Water Scale) was in a relatively low range due to the influence of winter climate and data foundation. This helped to test the model's spatial simulation capability under low-value scenarios with weak data support and stable hydrological processes. By comparing the original data with the downscaling results, the model's adaptability at the start of the time series could be assessed. July 2010 was at the peak of the high-water season, with the TWSA curve fluctuating dramatically. Concentrated summer rainfall led to rapid accumulation and complex redistribution of terrestrial water storage. This month was a key node for testing the model's ability to capture dynamic changes in high-value areas and reconstruct the spatial heterogeneity of complex hydrological processes. In January 2016, the TWSA was in a significantly low period in the series, and the separation between the original data and the downscaled model curve was high, making it suitable for analyzing the model's performance under low-value anomaly scenarios. August 2022, as the final high-water month of the study, allowed for the assessment of the model's adaptability and stability to recent hydrological environmental characteristics by comparing data before and after downscaling.

[0075] like Figure 4As shown, the results demonstrate that the monthly adaptive downscaling model effectively improves the spatial resolution of the data. The TWSA ranges for the study area in the four time periods are: November 2003 (-160~2mm), July 2010 (-240~6mm), January 2016 (-300~-60mm), and August 2022 (-250~20mm). Significant differences in the spatial distribution of TWSA exist between different years. Specifically, in November 2003, the original GRACE data, due to resolution limitations, did not adequately depict regional details; the monthly adaptive downscaling model, however, exhibited richer spatial details and a higher degree of fit to reality. In July 2010, the downscaling model effectively reflected regional differences, and the spatial variations conformed to actual hydrological characteristics. In January 2016, the model more naturally and accurately depicted the low-to-high value transition characteristics in the western part of the study area and other regions. In August 2022, the model demonstrated good spatial distribution accuracy and matching with the actual scene in the north and southwest, and overall showed significant advantages in accuracy and matching. In summary, the monthly adaptive downscaling model exhibits good reliability in simulating high-precision monthly changes in land water storage in this region.

[0076] Example 3 This embodiment also provides a monthly adaptive land water storage downscaling system, including: an acquisition module, a filtering module, a construction module, and an output module; the acquisition module is used to acquire monthly GRACE land water storage anomaly data and several high-resolution environmental prediction factor data of the study area within a preset time period; the filtering module filters out monthly feature factor sets based on the monthly GRACE land water storage anomaly data and environmental prediction factor data; before filtering, the environmental factor data needs to be resampled to the same resolution as the GRACE data; the construction module is used to construct a hybrid downscaling model based on the monthly feature factor set and the corresponding monthly GRACE land water storage anomaly data using a monthly independent modeling strategy; the output module is used to input the high spatial resolution environmental prediction factor data into the trained monthly hybrid downscaling model and output the high spatial resolution monthly land water storage anomaly prediction results.

[0077] The embodiments described above are merely preferred embodiments of the present invention and are not intended to limit the scope of the present invention. Various modifications and improvements made to the technical solutions of the present invention by those skilled in the art without departing from the spirit of the present invention should fall within the protection scope defined by the claims of the present invention.

Claims

1. A monthly adaptive method for downscaling terrestrial water storage, characterized in that, Includes the following steps: S1. Obtain monthly GRACE land water storage anomaly data and several high-resolution environmental prediction factor data for the study area within a preset time period. S2. Based on the monthly GRACE land water storage anomaly data and environmental prediction factor data, select the monthly feature factor set; before selection, the environmental prediction factor data needs to be resampled to the same resolution as the GRACE data. S3. Based on the monthly feature factor set and the corresponding monthly GRACE land water storage anomaly data, a hybrid downscaling model is constructed using a monthly independent modeling strategy. S4. Input the high spatial resolution environmental predictor data into the trained monthly hybrid downscaling model and output the high spatial resolution monthly land water storage anomaly prediction results.

2. The method for downscaling terrestrial water storage based on monthly adaptive scaling according to claim 1, characterized in that, S1 includes: GRACE / GRACE-FO Mascon data from different institutions were acquired, and the average value was calculated to represent the terrestrial water storage anomaly in the study area. The acquired data were interpolated using singular spectrum interpolation to impute missing months.

3. The method for downscaling terrestrial water storage based on monthly adaptive scaling according to claim 1, characterized in that, In step S2, the contribution of environmental prediction factor data is calculated monthly using partial least squares regression, and a monthly feature factor set is selected. The steps include: The environmental predictor data and the corresponding monthly GRACE terrestrial water storage anomaly data are input into the partial least squares regression model, and the variable projection importance index value of each environmental predictor data is calculated month by month; predictors with variable projection importance index values ​​greater than 0.8 are selected to form the feature factor set for that month.

4. The method for downscaling terrestrial water storage based on monthly adaptive scaling according to claim 1, characterized in that, In S3, for months in which the prediction performance of the CNN-LSTM model meets the preset requirements, the CNN-LSTM model is used for modeling to extract spatiotemporal features. For months in which the CNN-LSTM model's prediction performance fails to meet the preset requirements, the gradient boosting decision tree model is switched to enhance nonlinear fitting capabilities.

5. The method for downscaling terrestrial water storage based on monthly adaptive scaling according to claim 4, characterized in that, S3 further includes: First, a CNN-LSTM model is used to initially model all months. Based on the modeling results, the months are divided into different performance levels. For months with good performance, the CNN-LSTM model is continued to be used. For months with poor performance, the gradient boosting decision tree model is used instead. For months with extremely poor performance, an enhanced gradient boosting decision tree model is used.

6. The method for downscaling terrestrial water storage based on monthly adaptive scaling according to claim 4, characterized in that, S3 also includes: during model training, adopting differentiated dataset partitioning strategies, feature engineering processes, and data augmentation strategies for different months; and using dynamically changing learning rates and early stopping mechanisms to prevent overfitting.

7. The method for downscaling terrestrial water storage based on monthly adaptive scaling according to claim 4, characterized in that, S3 further includes: evaluating the performance of the trained monthly downscaling model using test set data; and using the coefficient of determination R0. 2 The root mean square error (RMSE) and correlation coefficient (CC) are used to quantify the accuracy of the downscaling results.

8. A monthly adaptive terrestrial water storage downscaling system, said system being used to implement the method according to any one of claims 1-7, characterized in that, include: The module includes a data acquisition module, a filtering module, a data construction module, and an output module. The acquisition module is used to acquire monthly GRACE land water storage anomaly data and several high-resolution environmental prediction factor data of the study area within a preset time period. The filtering module is used to filter out monthly feature factor sets based on monthly GRACE land water storage anomaly data and environmental prediction factor data; before filtering, the environmental factor data needs to be resampled to the same resolution as the GRACE data. The construction module is used to build a hybrid downscaling model based on the monthly feature factor set and the corresponding monthly GRACE land water storage anomaly data, using a monthly independent modeling strategy. The output module is used to input high spatial resolution environmental predictor data into the trained monthly hybrid downscaling model and output high spatial resolution monthly land water storage anomaly prediction results.

Images