Gravity satellite groundwater storage calculation method based on hydrological model calibration coupling

Abstract

The application discloses a kind of gravity satellite groundwater storage calculation method based on hydrological model verification coupling, comprising: obtaining the equivalent water column height of gravity satellite, land soil moisture, meteorological precipitation temperature, remote sensing vegetation index, monitoring well water level and multi-source data of water intake opening mining quantity, satellite intermittent period is filled using meteorological driving LSTM;Extract the land water storage, soil water and crown layer water change sequence, obtain the initial groundwater storage change sequence by empirical mode decomposition and spatial regression correction;Construct a numerical model of groundwater flow, use the initial sequence as the observation field, and combine the weighted upscaled water supply to carry out ensemble prediction and dynamic calibration to obtain the parameter posterior distribution set;Update the model to generate a monthly full-grid water head field, which is converted into the final groundwater storage change result set;Based on the result set, construct a multilayer perception machine to simulate the combined scenario of climate and exploitation, and generate the dynamic response path of the storage. High precision of groundwater storage change inversion, provide help for management decision.

Description

Technical Field

[0001] This invention relates to the field of hydrogeology, and in particular to a method for calculating groundwater reserves using gravity satellites based on hydrological model verification coupling. Background Technology

[0002] Groundwater is a crucial component of water resources, and its over-exploitation can lead to geological and environmental problems such as land subsidence and seawater intrusion. Accurately understanding the spatiotemporal variations of regional groundwater reserves is a fundamental prerequisite for the scientific management and sustainable utilization of groundwater. Traditional groundwater monitoring relies on a network of monitoring wells, measuring water levels point-by-point using manual or automated equipment. This approach has inherent limitations, including poor spatial continuity, high maintenance costs, and difficulty in covering large-scale areas, especially in providing long-term, high-density isotropic observation data.

[0003] The Gravity Recovery and Climate Experiment Satellite (GRESE) measures inter-satellite distance changes using its onboard precision ranging system to invert the time-varying characteristics of Earth's gravity field, thereby obtaining changes in terrestrial water storage. Since changes in the gravity field are mainly caused by the distribution and movement of water, changes in groundwater storage can be indirectly obtained by subtracting surface water components from changes in terrestrial water storage. This method has the advantages of wide spatial coverage and good timeliness, and has been applied in groundwater storage assessment in some areas. However, relying solely on gravity satellite inversion suffers from problems such as coarse spatial resolution, signal leakage, and north-south striping errors, making it difficult to meet the needs of refined management in small- and medium-scale areas. Furthermore, existing fusion methods often involve simple comparisons or linear fitting between hydrological models and satellite inversion data, failing to achieve deep mutual constraints and calibration between remote sensing macroscopic observations and the physical mechanisms of hydrological models. This results in insufficient physical self-consistency and spatiotemporal continuity of the inversion results, limiting the practical application effectiveness of gravity satellite technology in the refined management of groundwater in key provincial over-extraction areas. Summary of the Invention

[0004] This invention overcomes the shortcomings of the prior art and provides a gravity satellite groundwater storage calculation method based on hydrological model verification coupling. Its important purpose is to improve the accuracy and spatiotemporal continuity of groundwater storage change inversion, and enhance the physical self-consistency and management decision support capabilities of the inversion results.

[0005] To achieve the above objectives, the first aspect of this invention provides a method for calculating groundwater storage using gravity satellites based on hydrological model verification coupling, comprising: We acquire high-grid data of gravity satellite equivalent water column, land surface soil moisture content and canopy water equivalent data, meteorological station precipitation and temperature data, remote sensing vegetation index data, monitoring well water level data and water intake metering data of the target area, and use a meteorological-driven long short-term memory network to fill in the satellite gaps to obtain a multi-source basic dataset. Based on the multi-source basic dataset, the land water storage change sequence and the soil water and canopy water change sequence are obtained. The signal is corrected by empirical mode decomposition into intrinsic mode function components and a spatial regression model is established to obtain the initial groundwater storage change sequence. A groundwater flow numerical model is constructed, with the initial groundwater storage change as the observation field. An ensemble forecast instance is generated based on the disturbance permeability coefficient, specific yield, and rainfall infiltration. The grid hydraulic head change is upscaled to the regional average simulated value with specific yield as the weight. The key parameters and states of the hydrological model are dynamically calibrated to obtain the parameter posterior distribution set. The groundwater flow numerical model is updated and run based on the posterior distribution set of the parameters. After the update, the model generates a monthly full-grid hydraulic head field sequence, and the hydraulic head field is converted into the final calculated value of groundwater storage change, forming a result dataset that has both macro trend and spatiotemporal fine features. Based on the aforementioned dataset, with the final change in groundwater reserves as the dependent variable, and precipitation, temperature, vegetation index, and mining disturbance index constructed from actual extraction volume as independent variables, a multilayer perceptron network is constructed to simulate the combined climate and mining scenarios and generate dynamic response paths for reserves.

[0006] In this scheme, the acquisition of gravity satellite equivalent water column high-grid data, land surface soil moisture content and canopy water equivalent data, meteorological station precipitation and temperature data, remote sensing vegetation index data, monitoring well water level data, and water intake metering data for the target area is achieved by using a meteorological-driven long short-term memory network to fill in satellite gaps, resulting in a multi-source basic dataset, specifically including: Acquire gravity satellite equivalent water column high grid data, land surface soil moisture and canopy water equivalent data, meteorological station precipitation and temperature observation records, vegetation index grid data extracted from remote sensing images, monthly measured water level data of monitoring wells, and monthly actual groundwater extraction sequence recorded by water intake metering facilities in the target area. All raster data are resampled to a unified grid through bilinear interpolation, and station data are interpolated to the same grid nodes through inverse distance weighted interpolation. The time dimension is uniformly aggregated into monthly averages to form an initial multi-source data stack. The multi-year monthly average and standard deviation are calculated for each variable sequence in each grid cell of the initial multi-source data stack. If the current month's value deviates from the corresponding monthly average by more than three times the standard deviation, it is marked as an outlier and set as a missing value. The same criteria for anomaly detection and removal are applied to the monitoring well water level and water intake extraction volume data on their original station sequences to obtain a multi-source sequence containing only true missing values ​​and satellite discontinuities. To address the long-term continuous discontinuities in the equivalent water column height sequence of gravity satellites caused by alternating satellite missions, grid cells are used as the filling objects. At least three consecutive periods before the discontinuity are selected as training samples. Monthly precipitation and monthly temperature are used as input features, and the corresponding equivalent water column height is used as the learning target. A long short-term memory regression prediction network with two hidden layers was constructed and iteratively trained to convergence using the mean squared error loss function and the adaptive moment estimation optimizer. The existing meteorological data during the intermission period was input into the trained network and predicted output month by month to complete the entire intermission period sequence. Simultaneously, short-term sporadic missing data in the vegetation index sequence, soil moisture content sequence, canopy water equivalent sequence, and monitoring well water level sequence were interpolated and completed grid-by-grid using a piecewise cubic spline interpolation method. After all data were processed, they were recombined according to a unified spatial grid and time index to form a multi-source basic dataset.

[0007] In this scheme, the step of obtaining the land water storage change sequence and soil water and canopy water change sequence based on the multi-source basic dataset, and then using empirical mode decomposition to obtain intrinsic mode function components and establishing a spatial regression model for signal correction to obtain the initial groundwater storage change sequence specifically includes: Gravity satellite equivalent water column height sequences were extracted from multi-source basic datasets as a signal source for changes in terrestrial water storage. Simultaneously, soil moisture content sequences and canopy water equivalent sequences that had completed spatiotemporal registration were extracted as soil water change sequences and canopy water change sequences. Gaussian filtering is applied to the equivalent water column high sequence of gravity satellites. The time series of each grid is smoothed to suppress north-south striping errors and high-frequency noise, and the filtered land water storage change sequence is output. An integrated empirical mode decomposition algorithm was used to adaptively decompose the filtered land water storage change sequence and soil water change sequence into intrinsic mode function components arranged from high frequency to low frequency and residual trend terms reflecting long-term trends, resulting in two sets of intrinsic mode function components. Pairwise matching is performed on the components at the same frequency position in the two sets of intrinsic mode function components, and the Pearson correlation coefficient is calculated. Component pairs that exceed the preset significance threshold are judged as significantly correlated component pairs at the same frequency. A multivariate linear spatial regression correction model was established with monthly precipitation and monthly temperature as independent variables. The signal leakage and spatial attenuation deviation of the corresponding frequency bands were corrected on each grid. The components after frequency band correction of the regression model were linearly superimposed with the residual trend term to obtain the complete sequence of land water storage changes after time-frequency correction. The initial groundwater storage change sequence is output by subtracting the corresponding soil water change sequence and canopy water change sequence from the corrected terrestrial water storage change sequence month by month, with the equivalent water column height as the unit.

[0008] In this scheme, the construction of a groundwater flow numerical model and the dynamic calibration of key parameters and states of the hydrological model to obtain the posterior distribution set of parameters specifically include: Based on the hydrogeological map of the target area, the lithological distribution data of aquifers and the results of historical pumping tests, the permeability coefficient zone and the water yield zone were divided. The initial values ​​and a priori value ranges of parameters for each zone were determined. The MODFLOW software was used to discretize the study area into an orthogonal grid and set the boundary conditions. The rainfall infiltration recharge and the actual monthly extraction volume were input to complete the initial construction of the groundwater flow numerical model. Using the initial groundwater storage change sequence as the observation field, the Latin hypercube sampling method is used to randomly sample the permeability coefficient and specific yield values ​​of each parameter partition. At the same time, a normal distribution random perturbation is applied to the rainfall infiltration recharge. Several parallel model instances are generated through parameter combination to form an ensemble forecast field. At each assimilation time, all model instances in the set are run to obtain grid-by-grid head forecast values ​​and calculate head change values. The head change value of each grid is multiplied by the corresponding water supply value of the partition with water supply as the weight, and then the average is weighted by grid area. The result is then upscaled to convert the simulated value of regional average groundwater storage change. After all instances have completed the conversion process, the set of simulated values ​​is obtained. The forecast error covariance matrix is ​​calculated based on the set of simulated values. The observation error covariance matrix is ​​constructed by combining the initial groundwater storage change observation value with the preset observation error standard deviation. The Kalman gain matrix is ​​calculated and the state of the grid-by-grid water head forecast value for each instance is updated. At the same time, the permeability coefficient value and water yield value of each zone are updated. The complete process of ensemble operation, upscaling transformation, covariance calculation, gain calculation and state parameter update is carried out month by month along the time series. After all time series assimilation cycles are completed, the final converged permeability coefficient value and water yield value of each instance are extracted to generate the parameter posterior distribution set.

[0009] In this scheme, the groundwater flow numerical model is updated based on the posterior distribution set of the parameters. After the update, the model generates a monthly full-grid hydraulic head field sequence, which is then converted into the final calculated value of groundwater storage changes, forming a result dataset that combines macroscopic trends with fine spatiotemporal characteristics. Specifically, this includes: The posterior mean of the permeability coefficient and the posterior mean of the specific yield are extracted from the posterior distribution set of the parameters. The posterior mean is then used to replace the parameter values ​​of the corresponding partitions in the initially constructed groundwater flow numerical model. The grid division and boundary conditions remain unchanged, resulting in a high-fidelity groundwater flow numerical model calibrated with gravity satellite observation data. Driven by monthly rainfall infiltration recharge and monthly actual groundwater extraction, which are centrally covered by multi-source basic datasets throughout the entire study period, the updated groundwater flow numerical model is run month by month. The full-grid head distribution is calculated once at each time step, and the output is a monthly full-grid head field sequence covering the entire study period. The monthly full-grid hydraulic head field sequence is converted into the final calculated value of groundwater storage change. The hydraulic head difference between adjacent months of each grid is calculated and multiplied by the posterior mean of the water supply of the corresponding partition. After weighted averaging by grid area, it is converted into the final calculated value of groundwater storage change on a monthly scale in units of equivalent water column height. The final calculated value of groundwater storage change in the grid where the monitoring well is located is divided by the posterior mean of the water yield of the corresponding grid to obtain the predicted value of water level change. The predicted value of water level change is compared with the measured value of water level change over time, and the root mean square error and Pearson correlation coefficient are calculated to verify the inversion accuracy. The verified monthly groundwater storage change final calculation sequence, the monthly full-grid head field sequence, and the posterior mean of permeability coefficient and water yield in the parameter posterior set are output together to form a result dataset that combines the macroscopic trend accuracy of gravity satellites with the spatiotemporal fine features of distributed hydrological models.

[0010] In this scheme, the step of constructing a multilayer perceptron network based on the obtained dataset to simulate combined climate and mining scenarios and generate dynamic response paths for reserves specifically includes: The final calculated value of monthly groundwater storage change was extracted from the results dataset as the output target variable. The monthly precipitation, monthly temperature and vegetation index sequences of the same period were retrieved from the multi-source basic dataset as natural driving factors. The monthly mining disturbance intensity index was constructed by the historical monthly actual mining volume of the target area as a proxy variable for human activities. Time alignment and normalization were performed on all independent and dependent variables. A multilayer perceptron neural network containing an input layer, two hidden layers and an output layer was constructed. The input layer was set with four nodes corresponding to monthly precipitation, monthly average temperature, vegetation index and mining disturbance intensity index, respectively. The hidden layer adopted the modified linear unit activation function, and the output layer was the predicted value of groundwater storage change. The network uses the mean squared error loss function and the adaptive moment estimation optimizer to update the weights iteratively. The historical samples are divided into training set, validation set and test set according to time order. The training process monitors the value of the validation set loss function and implements an early stopping strategy. The training is terminated when the validation loss no longer decreases for several consecutive rounds. After the network training is completed, extract all connection weight matrices from the input layer to the hidden layer and from the hidden layer to the output layer and take their absolute values. For each input factor, calculate the sum of the products of the connection weights that lead to the output node through each hidden node. Divide the sum of the products of each factor by the sum of the products of all factors to obtain the normalized contribution score of each factor to the change in groundwater storage. Several sets of cross-combinations of climate scenarios and mining intensity scenarios are set. The climate scenarios are set according to the frequency level of precipitation, and the mining intensity scenarios are set according to the multiple of the mining disturbance index. The input vector of each scenario is fed into the trained network one by one. The network predicts and outputs the groundwater storage change value of the corresponding scenario on a monthly basis, generating dynamic response paths of storage under different combinations of climate conditions and mining intensity.

[0011] Compared with existing technologies, this solution has the following advantages: Regarding inversion accuracy, adaptive time-frequency decomposition and frequency-band correction of gravity satellite signals and soil water signals are performed by integrating empirical mode decomposition. This effectively solves the inversion bias caused by the forced mixing of signals at different scales in the traditional direct subtraction method, significantly improving the extraction accuracy of initial groundwater storage changes. Regarding physical consistency, an iterative ensemble Kalman filter assimilation mechanism is introduced. Using water yield as a weight, the distributed grid head changes are upscaled to regional average simulated values. Within a unified physical framework, the model's head state and key parameters are updated synchronously, allowing the model parameters to converge from prior values ​​dependent on empirical assignments to posterior optimal estimates constrained by remote sensing observations. This solves the problem of disconnect between empirical water yield assignments and scale transformation in traditional methods. Regarding spatiotemporal continuity, a long short-term memory network is used to fill in the gaps in satellite missions. Combined with the monthly full-grid simulation capability of the calibrated hydrological model, an uninterrupted, spatially continuous groundwater storage change dataset is generated, possessing the dual advantages of accurate macroscopic trends from gravity satellites and the spatiotemporal precision of distributed models. In terms of attribution analysis, the mining disturbance intensity index, constructed based on actual mining volume, was introduced as a technical proxy variable for human activities. An attribution algorithm system based on physically measurable parameters was established, providing objective and calculable decision support for the formulation of total groundwater extraction control indicators and mining intensity control thresholds. Attached Figure Description

[0012] To more clearly illustrate the technical solutions in the embodiments or examples of the present invention, the drawings used in the embodiments or examples will be 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 according to these drawings without creative effort.

[0013] Figure 1 A flowchart of a gravity satellite groundwater storage calculation method based on hydrological model verification coupling is provided as an embodiment of the present invention. Figure 2 A flowchart of a method for reconstructing and verifying groundwater storage change data according to an embodiment of the present invention; The realization of the objective, functional features and advantages of the present invention will be further explained in conjunction with the embodiments and with reference to the accompanying drawings. Detailed Implementation

[0014] To better understand the above-mentioned objectives, features, and advantages of the present invention, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments. It should be noted that, unless otherwise specified, the embodiments and features described in these embodiments can be combined with each other.

[0015] 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 therefore the scope of protection of the invention is not limited to the specific embodiments disclosed below.

[0016] Figure 1 A flowchart of a gravity satellite groundwater storage calculation method based on hydrological model verification coupling is provided as an embodiment of the present invention. like Figure 1 As shown, this invention provides a flowchart of a gravity satellite groundwater storage calculation method based on hydrological model verification coupling, including: We acquire high-grid data of gravity satellite equivalent water column, land surface soil moisture content and canopy water equivalent data, meteorological station precipitation and temperature data, remote sensing vegetation index data, monitoring well water level data and water intake metering data of the target area, and use a meteorological-driven long short-term memory network to fill in the satellite gaps to obtain a multi-source basic dataset. Based on the multi-source basic dataset, the land water storage change sequence and the soil water and canopy water change sequence are obtained. The signal is corrected by empirical mode decomposition into intrinsic mode function components and a spatial regression model is established to obtain the initial groundwater storage change sequence. A groundwater flow numerical model is constructed, with the initial groundwater storage change as the observation field. An ensemble forecast instance is generated based on the disturbance permeability coefficient, specific yield, and rainfall infiltration. The grid hydraulic head change is upscaled to the regional average simulated value with specific yield as the weight. The key parameters and states of the hydrological model are dynamically calibrated to obtain the parameter posterior distribution set. The groundwater flow numerical model is updated and run based on the posterior distribution set of the parameters. After the update, the model generates a monthly full-grid hydraulic head field sequence, and the hydraulic head field is converted into the final calculated value of groundwater storage change, forming a result dataset that has both macro trend and spatiotemporal fine features. Based on the aforementioned dataset, with the final change in groundwater reserves as the dependent variable, and precipitation, temperature, vegetation index, and mining disturbance index constructed from actual extraction volume as independent variables, a multilayer perceptron network is constructed to simulate the combined climate and mining scenarios and generate dynamic response paths for reserves.

[0017] Furthermore, in a preferred embodiment of the present invention, the acquisition of gravity satellite equivalent water column high-grid data, land surface soil moisture content and canopy water equivalent data, meteorological station precipitation and temperature data, remote sensing vegetation index data, monitoring well water level data, and water intake metering data of the target area, and the use of a meteorological-driven long short-term memory network to fill in satellite gaps to obtain a multi-source basic dataset, specifically includes: Acquire gravity satellite equivalent water column high grid data, land surface soil moisture and canopy water equivalent data, meteorological station precipitation and temperature observation records, vegetation index grid data extracted from remote sensing images, monthly measured water level data of monitoring wells, and monthly actual groundwater extraction sequence recorded by water intake metering facilities in the target area. All raster data are resampled to a unified grid through bilinear interpolation, and station data are interpolated to the same grid nodes through inverse distance weighted interpolation. The time dimension is uniformly aggregated into monthly averages to form an initial multi-source data stack. The multi-year monthly average and standard deviation are calculated for each variable sequence in each grid cell of the initial multi-source data stack. If the current month's value deviates from the corresponding monthly average by more than three times the standard deviation, it is marked as an outlier and set as a missing value. The same criteria for anomaly detection and removal are applied to the monitoring well water level and water intake extraction volume data on their original station sequences to obtain a multi-source sequence containing only true missing values ​​and satellite discontinuities. To address the long-term continuous discontinuities in the equivalent water column height sequence of gravity satellites caused by alternating satellite missions, grid cells are used as the filling objects. At least three consecutive periods before the discontinuity are selected as training samples. Monthly precipitation and monthly temperature are used as input features, and the corresponding equivalent water column height is used as the learning target. A long short-term memory regression prediction network with two hidden layers was constructed and iteratively trained to convergence using the mean squared error loss function and the adaptive moment estimation optimizer. The existing meteorological data during the intermission period was input into the trained network and predicted output month by month to complete the entire intermission period sequence. Simultaneously, short-term sporadic missing data in the vegetation index sequence, soil moisture content sequence, canopy water equivalent sequence, and monitoring well water level sequence were interpolated and completed grid-by-grid using a piecewise cubic spline interpolation method. After all data were processed, they were recombined according to a unified spatial grid and time index to form a multi-source basic dataset.

[0018] It is important to note that in monitoring regional groundwater storage changes, a unified spatiotemporal reference for multi-source data is a prerequisite for all subsequent calculations. Gravity satellite equivalent water column height data, soil moisture and canopy water equivalent data output by the land surface assimilation system, and remote sensing vegetation index data are typically stored with their own independent grid resolutions and projection methods. Meteorological station and monitoring well data, on the other hand, are discrete point records. Directly incorporating these into the inversion process would introduce systematic biases due to spatial misalignment and inconsistent temporal indexes. Therefore, this scheme employs bilinear interpolation to uniformly resample raster data to a preset grid, uses inverse distance weighted interpolation to spatialize station data to the same grid nodes, and aggregates the temporal resolution of all data into a monthly average, forming a spatially matched and temporally synchronized initial multi-source data stack. This operation eliminates the spatiotemporal mismatch problem between heterogeneous multi-source data at the data source, providing a unified reference framework for subsequent grid-by-grid, month-by-month signal decomposition, assimilation calibration, and attribution analysis.

[0019] Outliers inevitably exist in observational data due to instrument malfunctions, environmental interference, or human error in recording. Directly using these outliers for model calculations would severely reduce inversion accuracy and model stability. This step uses three standard deviations as a threshold to statistically discriminate each variable sequence in each grid cell. Data points exceeding the threshold are marked as missing, ensuring that outlier removal is statistically objective and repeatable, avoiding uncertainties introduced by subjective experience. To address the data gaps of over a year caused by the alternation of GRACE and GRACE-FO satellite missions, this step introduces a Long Short-Term Memory (LSTM) network for completion. This network excels at capturing the nonlinear temporal dependencies between meteorological drivers such as precipitation and temperature and changes in terrestrial water storage. By learning from at least three years of continuous historical data prior to the gap, a predictive model is built. After training with meteorological data from the gap period, the network can predict and output equivalent water column heights monthly, thus completing the gravity satellite sequence into a continuous, uninterrupted sequence. Compared to traditional linear interpolation or moving average imputation methods, long short-term memory (LSTM) network imputation can preserve the dynamic fluctuation characteristics of water storage changes, reduce trend deviation between the imputed segment and the true sequence, and provide a complete and reliable long-term input foundation for subsequent integrated empirical mode decomposition and ensemble Kalman filter assimilation. For short-term sporadic missing data in vegetation index, soil moisture, canopy water equivalent, and monitoring well water level sequences, piecewise cubic spline interpolation is used for grid-by-grid interpolation imputation. This method can maintain the second-order continuous smoothness of the interpolation curve at the segment nodes and avoid the interference of the broken line inflection points generated by linear interpolation on subsequent trend analysis. Finally, all data with outlier removal and missing data imputation are recombined according to a unified spatial grid and time index to form a multi-source basic dataset, providing standardized input for subsequent steps of terrestrial water storage change extraction and initial groundwater storage inversion.

[0020] Furthermore, in a preferred embodiment of the present invention, the step of obtaining the land water storage change sequence and soil water and canopy water change sequence based on the multi-source basic dataset, and then performing signal correction by empirical mode decomposition into intrinsic mode function components and establishing a spatial regression model to obtain the initial groundwater storage change sequence specifically includes: Gravity satellite equivalent water column height sequences were extracted from multi-source basic datasets as a signal source for changes in terrestrial water storage. Simultaneously, soil moisture content sequences and canopy water equivalent sequences that had completed spatiotemporal registration were extracted as soil water change sequences and canopy water change sequences. Gaussian filtering is applied to the equivalent water column high sequence of gravity satellites. The time series of each grid is smoothed to suppress north-south striping errors and high-frequency noise, and the filtered land water storage change sequence is output. An integrated empirical mode decomposition algorithm was used to adaptively decompose the filtered land water storage change sequence and soil water change sequence into intrinsic mode function components arranged from high frequency to low frequency and residual trend terms reflecting long-term trends, resulting in two sets of intrinsic mode function components. Pairwise matching is performed on the components at the same frequency position in the two sets of intrinsic mode function components, and the Pearson correlation coefficient is calculated. Component pairs that exceed the preset significance threshold are judged as significantly correlated component pairs at the same frequency. A multivariate linear spatial regression correction model was established with monthly precipitation and monthly temperature as independent variables. The signal leakage and spatial attenuation deviation of the corresponding frequency bands were corrected on each grid. The components after frequency band correction of the regression model were linearly superimposed with the residual trend term to obtain the complete sequence of land water storage changes after time-frequency correction. The initial groundwater storage change sequence is output by subtracting the corresponding soil water change sequence and canopy water change sequence from the corrected terrestrial water storage change sequence month by month, with the equivalent water column height as the unit.

[0021] It should be noted that the equivalent water column height sequence from gravity satellites reflects the overall change in regional terrestrial water storage. This signal simultaneously contains the superposition of multiple components, including groundwater, soil water, canopy water, and observational noise. Directly separating the groundwater signal from the original sequence would result in insufficient inversion accuracy due to noise interference and signal aliasing. Therefore, Gaussian filtering is required before processing the gravity satellite data. Gaussian filtering uses a convolution kernel to weight and smooth the time series of each grid, effectively suppressing north-south striping errors and high-frequency random noise. This results in a more accurate reflection of the actual fluctuation characteristics of hydrological processes in the filtered terrestrial water storage change sequence, providing a cleaner base signal for subsequent time-frequency decomposition.

[0022] Traditional direct subtraction of water balance data simply subtracts changes in terrestrial water storage from changes in soil water in the original time domain, ignoring the differences in the changing patterns of different water components over time. For example, soil water responses to precipitation events are concentrated in high-frequency monthly fluctuations, while groundwater changes are more reflected in low-frequency interannual trends. Directly subtracting the two in the original domain will mix high-frequency noise and low-frequency trends, generating spurious fluctuation signals. To address this issue, this step employs an integrated empirical mode decomposition algorithm to adaptively decompose the filtered terrestrial water storage change sequence and the soil water change sequence into time-frequency components. This algorithm calculates the envelope mean and iteratively sieves the data by adding white noise to the sequences multiple times. Without pre-setting basis functions, it can peel the signal layer by layer into intrinsic mode function components and residual trend terms arranged from high frequency to low frequency based on the time-scale characteristics of the data. After decomposition, each sequence obtains a set of components arranged in descending order of frequency. High-frequency components correspond to rapid fluctuations at the monthly scale, while low-frequency components correspond to slow trends at the interannual scale, achieving automatic separation of water storage change signals at different time scales.

[0023] Subsequently, components at the same frequency band position in the two component sets were matched pairwise, and the Pearson correlation coefficient between each pair of components was calculated. Component pairs with correlation coefficients exceeding a preset significance threshold were identified as significantly correlated components at the same frequency. The physical significance of this operation is that if there is a significant correlation between changes in terrestrial water storage and changes in soil water at the same time scale, it indicates that the fluctuations of both within that frequency band are dominated by common climate driving factors, and the signals are comparable. For the selected significantly correlated components at the same frequency, a multiple linear spatial regression correction model was established using monthly precipitation and monthly temperature as independent variables. This model corrected the deviation of gravity satellite signals within that frequency band caused by spatial leakage and attenuation on a grid-by-grid basis. Components with correlation coefficients below the threshold were identified as local high-frequency noise and directly removed. The components corrected by the regression model on a frequency band-by-frequency basis were linearly superimposed with the residual trend term to obtain the complete time-frequency corrected sequence of terrestrial water storage changes. By subtracting the corresponding soil water change series and canopy water change series from the corrected terrestrial water storage change series month by month, an initial groundwater storage change series with significantly improved accuracy can be output, providing more reliable observational constraints for subsequent hydrological model assimilation and calibration.

[0024] Furthermore, in a preferred embodiment of the present invention, the step of constructing a groundwater flow numerical model and dynamically calibrating the key parameters and states of the hydrological model to obtain a posterior distribution set of parameters specifically includes: Based on the hydrogeological map of the target area, the lithological distribution data of aquifers and the results of historical pumping tests, the permeability coefficient zone and the water yield zone were divided. The initial values ​​and a priori value ranges of parameters for each zone were determined. The MODFLOW software was used to discretize the study area into an orthogonal grid and set the boundary conditions. The rainfall infiltration recharge and the actual monthly extraction volume were input to complete the initial construction of the groundwater flow numerical model. Using the initial groundwater storage change sequence as the observation field, the Latin hypercube sampling method is used to randomly sample the permeability coefficient and specific yield values ​​of each parameter partition. At the same time, a normal distribution random perturbation is applied to the rainfall infiltration recharge. Several parallel model instances are generated through parameter combination to form an ensemble forecast field. At each assimilation time, all model instances in the set are run to obtain grid-by-grid head forecast values ​​and calculate head change values. The head change value of each grid is multiplied by the corresponding water supply value of the partition with water supply as the weight, and then the average is weighted by grid area. The result is then upscaled to convert the simulated value of regional average groundwater storage change. After all instances have completed the conversion process, the set of simulated values ​​is obtained. The forecast error covariance matrix is ​​calculated based on the set of simulated values. The observation error covariance matrix is ​​constructed by combining the initial groundwater storage change observation value with the preset observation error standard deviation. The Kalman gain matrix is ​​calculated and the state of the grid-by-grid water head forecast value for each instance is updated. At the same time, the permeability coefficient value and water yield value of each zone are updated. The complete process of ensemble operation, upscaling transformation, covariance calculation, gain calculation and state parameter update is carried out month by month along the time series. After all time series assimilation cycles are completed, the final converged permeability coefficient value and water yield value of each instance are extracted to generate the parameter posterior distribution set.

[0025] It should be noted that the construction of the groundwater flow numerical model is the physical carrier for achieving deep integration of gravity satellite observations and hydrological processes. During model building, the study area was divided into different permeability coefficient and specific yield zones based on hydrogeological maps, aquifer lithology distribution, and historical pumping test results. This is because these parameters exhibit significant spatial heterogeneity, and using a single uniform value would fail to characterize the differences in hydraulic response between different geological units. The MODFLOW software was selected to discretize the study area into orthogonal grids and set constant head and impermeable boundaries. Simultaneously, rainfall infiltration recharge and monthly actual extraction volume were input as source and sink terms, forming an initial numerical model that reflects the basic characteristics of groundwater flow in the study area. However, the permeability coefficient and specific yield in the model at this stage are only initial estimates based on prior experience or limited pumping test data, and contain significant uncertainties. Directly using them for prediction will lead to significant deviations between the simulation results and actual groundwater dynamics.

[0026] To correct for uncertainties in the model parameters, the initial groundwater storage change sequence output from the previous stage is used as an observation constraint, and an iterative ensemble Kalman filter assimilation method is introduced. First, Latin hypercube sampling is used to perform stratified random sampling of the permeability coefficient and specific yield for each parameter partition. Simultaneously, a normal random perturbation is applied to the rainfall infiltration recharge. Dozens of parallel model instances are generated through parameter combinations to form an ensemble forecast field. Latin hypercube sampling can achieve more uniform sample coverage within the parameter value range, avoiding the local clustering problem that may occur with traditional Monte Carlo sampling, and obtaining an effective representation of the parameter space with a smaller number of samples. A key problem to be solved during the assimilation process is the mismatch between the equivalent water column height retrieved from gravity satellites and the grid head output by the model in terms of physical dimensions and spatial scale. This scheme defines a physical transformation operator with specific yield as the weight, multiplies the head change of each grid by the specific yield value of the corresponding partition, and then averages it by grid area, upscaling it to convert it into a simulated equivalent water column height value of the regional average groundwater storage change. The significance of this transformation lies in the fact that water yield is itself a physical parameter linking changes in hydraulic head and water volume. Upscaling with water yield as the weight has a clear physical basis, rather than being a simple empirical regression fit. After obtaining the set of simulated values ​​at each assimilation time, the set dispersion is calculated to obtain the forecast error covariance. Combined with the observation error covariance, the Kalman gain matrix is ​​calculated, and the gain matrix is ​​used to synchronously update the hydraulic head state and parameter values ​​for each instance. As the assimilation process progresses month by month along the time series, the model parameters gradually converge from prior values ​​to a posterior distribution constrained by remote sensing observations, ultimately outputting a posterior distribution set of parameters. This provides a calibrated physical model foundation for subsequently generating a high-precision groundwater storage change dataset.

[0027] It is worth mentioning that the mathematical expression for the physical transformation operator is as follows: , In the formula, This represents the simulated value of the regional average groundwater storage change, expressed in equivalent water column height. This represents the total number of grid cells within the study area. For the first The area of ​​each grid; For the first The water yield value of the parameter partition where each grid is located; For the first Each grid at the current moment Head value; For the first The grid at the previous moment The head value.

[0028] Figure 2 A flowchart of a method for reconstructing and verifying groundwater storage change data according to an embodiment of the present invention; like Figure 2 As shown, this invention provides a flowchart of a method for reconstructing and verifying groundwater storage change data, including: The posterior mean of the permeability coefficient and the posterior mean of the specific yield are extracted from the posterior distribution set of the parameters. The posterior mean is then used to replace the parameter values ​​of the corresponding partitions in the initially constructed groundwater flow numerical model. The grid division and boundary conditions remain unchanged, resulting in a high-fidelity groundwater flow numerical model calibrated with gravity satellite observation data. Driven by monthly rainfall infiltration recharge and monthly actual groundwater extraction, which are centrally covered by multi-source basic datasets throughout the entire study period, the updated groundwater flow numerical model is run month by month. The full-grid head distribution is calculated once at each time step, and the output is a monthly full-grid head field sequence covering the entire study period. The monthly full-grid hydraulic head field sequence is converted into the final calculated value of groundwater storage change. The hydraulic head difference between adjacent months of each grid is calculated and multiplied by the posterior mean of the water supply of the corresponding partition. After weighted averaging by grid area, it is converted into the final calculated value of groundwater storage change on a monthly scale in units of equivalent water column height. The final calculated value of groundwater storage change in the grid where the monitoring well is located is divided by the posterior mean of the water yield of the corresponding grid to obtain the predicted value of water level change. The predicted value of water level change is compared with the measured value of water level change over time, and the root mean square error and Pearson correlation coefficient are calculated to verify the inversion accuracy. The verified monthly groundwater storage change final calculation sequence, monthly full-grid head field sequence, and the posterior mean of permeability coefficient and posterior mean of water yield from the parameter posterior set are output together to form a result dataset that combines the macroscopic trend accuracy of gravity satellites with the spatiotemporal fine features of distributed hydrological models.

[0029] It should be noted that after ensemble Kalman filtering assimilation, the permeability coefficients and water yield of each parameter partition in the model have converged from a priori range dependent on empirical assignments to a posterior distribution constrained by gravity satellite observation data. This step extracts the posterior mean of the permeability coefficient and the posterior mean of the water yield for each partition from the parameter posterior distribution set and directly replaces the corresponding parameter values ​​in the initial model. Because the parameter values ​​are dynamically adjusted at each time step during assimilation as the Kalman gain is updated, directly using the parameter state at a certain moment in the assimilation stage as the final result will deviate from the optimal estimate due to the instantaneous influence of the observation error at that moment. Taking the posterior mean after the entire time series converges can integrate all historical observation information, obtain the most statistically stable parameter estimate, and enable the updated model to truly possess "high-fidelity" characteristics.

[0030] After the model is updated, the updated groundwater flow numerical model is run month by month, driven by monthly rainfall infiltration recharge and monthly actual groundwater extraction, which are centrally covered by a multi-source basic dataset throughout the entire study period. At each time step, the full-grid head distribution is calculated, outputting a monthly full-grid head field sequence covering the entire study period. Unlike the traditional method of directly using the updated head values ​​during the assimilation process, the significance of rerunning the model here lies in the fact that the head field during the assimilation stage is instantaneously corrected by the Kalman gain at each moment, and the corrected head sequence may exhibit discontinuous jumps in the time dimension. By rerunning the model with the parameters calibrated to the optimal posterior mean, a smooth and continuous head field sequence that conforms to both physical laws and water balance constraints can be generated. After obtaining the monthly full-grid head field, the physical transformation operator defined and used in the assimilation stage is used to convert the head field into the final calculated value of groundwater storage change. Specifically, the head difference between adjacent months in each grid is calculated and multiplied by the posterior mean of the corresponding zone's water yield. This is then weighted by grid area and upscaled to convert it into the final calculated value of monthly groundwater storage change in equivalent water column units. To verify the inversion accuracy, the final calculated value of the grid containing the monitoring well is divided by the posterior mean of the corresponding grid's water yield to obtain the predicted water level change value. This is then compared with the measured water level change value over time, and the root mean square error and Pearson correlation coefficient are calculated. The significance of this conversion verification lies in the fact that the monitoring well water level is the only independent data source that can be directly verified on the ground. The conversion using the water yield mean within the parameter posterior set achieves the closed-loop verification logic of "using the assimilated and calibrated model parameters to verify the assimilated and calibrated model results." Finally, all the results data that have been verified for accuracy will be output together to form a results dataset that combines the macroscopic trend accuracy of gravity satellites with the spatiotemporal fine features of distributed hydrological models, providing complete data support for subsequent attribution analysis and management decisions.

[0031] Furthermore, in a preferred embodiment of the present invention, the step of constructing a multilayer perceptron network based on the obtained dataset to simulate combined climate and mining scenarios and generate dynamic response paths for reserves specifically includes: The final calculated value of monthly groundwater storage change was extracted from the results dataset as the output target variable. The monthly precipitation, monthly temperature and vegetation index sequences of the same period were retrieved from the multi-source basic dataset as natural driving factors. The monthly mining disturbance intensity index was constructed by the historical monthly actual mining volume of the target area as a proxy variable for human activities. Time alignment and normalization were performed on all independent and dependent variables. A multilayer perceptron neural network containing an input layer, two hidden layers and an output layer was constructed. The input layer was set with four nodes corresponding to monthly precipitation, monthly average temperature, vegetation index and mining disturbance intensity index, respectively. The hidden layer adopted the modified linear unit activation function, and the output layer was the predicted value of groundwater storage change. The network uses the mean squared error loss function and the adaptive moment estimation optimizer to update the weights iteratively. The historical samples are divided into training set, validation set and test set according to time order. The training process monitors the value of the validation set loss function and implements an early stopping strategy. The training is terminated when the validation loss no longer decreases for several consecutive rounds. After the network training is completed, extract all connection weight matrices from the input layer to the hidden layer and from the hidden layer to the output layer and take their absolute values. For each input factor, calculate the sum of the products of the connection weights that lead to the output node through each hidden node. Divide the sum of the products of each factor by the sum of the products of all factors to obtain the normalized contribution score of each factor to the change in groundwater storage. Several sets of cross-combinations of climate scenarios and mining intensity scenarios are set. The climate scenarios are set according to the frequency level of precipitation, and the mining intensity scenarios are set according to the multiple of the mining disturbance index. The input vector of each scenario is fed into the trained network one by one. The network predicts and outputs the groundwater storage change value of the corresponding scenario on a monthly basis, generating dynamic response paths of storage under different combinations of climate conditions and mining intensity.

[0032] It should be noted that, based on the high-precision groundwater storage change sequence obtained in the aforementioned steps, further quantifying the impact of different driving factors on storage changes is a crucial leap from "understanding the patterns of change" to "supporting management decisions." This step extracts the final calculated monthly groundwater storage change value from the results dataset as the output target variable, and simultaneously retrieves the monthly precipitation, monthly temperature, and vegetation index from the multi-source basic dataset as natural driving factors. Regarding the impact of human activities, this scheme introduces a monthly mining disturbance intensity index constructed based on actual extraction volume as a technical proxy variable for human activities. Its core concept is to transform the physical consequences of administrative control policies, i.e., changes in actual groundwater extraction, into continuously measurable technical parameters, so that the attribution analysis is based on purely physically measurable data.

[0033] After preparing the factor data, a multilayer perceptron neural network consisting of an input layer, two hidden layers, and an output layer was constructed to learn the nonlinear mapping relationship between multiple factors and groundwater storage changes. The significance of choosing an artificial neural network instead of a traditional linear regression model lies in the fact that the response of groundwater storage changes to precipitation and extraction typically exhibits significant threshold effects and lag characteristics. The impact of extraction may be masked during periods of abundant precipitation, while consecutive dry years coupled with high-intensity extraction can trigger a nonlinear accelerated decline. A simple linear superposition relationship cannot characterize this complex hydrological response mechanism. The network input layer has four nodes corresponding to monthly precipitation, monthly average temperature, vegetation index, and extraction disturbance intensity index, respectively. The hidden layer uses a modified linear unit activation function to introduce nonlinear expressive power, and the output layer is the predicted value of groundwater storage changes. During training, a mean squared error loss function and an adaptive moment estimation optimizer are used. Historical samples are divided chronologically, and an early stopping strategy is implemented to prevent overfitting and ensure the network's generalization prediction ability for future samples.

[0034] After network training, a connection weight-based algorithm is used to quantify the contribution of each input factor. The absolute values ​​of all connection weight matrices are extracted, and the sum of the weight products of each input factor's path through hidden nodes to the output node is calculated. After normalization, the percentage contribution of each factor is obtained. This contribution quantification process reveals the primary and secondary relationships of precipitation, temperature, vegetation cover, and human extraction on groundwater storage changes in the study area from a data-driven perspective. Finally, multiple combinations of climate scenarios and extraction intensity scenarios are set. The input vectors of each scenario are fed into the trained network one by one, and the corresponding groundwater storage changes are predicted and output monthly, generating dynamic response paths for reserves under different combinations of climate conditions and extraction intensities. This scenario simulation capability allows managers to intuitively compare the expected groundwater recovery effects under various combinations of schemes, such as strict enforcement of extraction control in normal water years and routine management in dry years, providing quantitative support for formulating total groundwater extraction control indicators and extraction intensity management thresholds.

[0035] Furthermore, the method for calculating groundwater storage based on hydrological model verification coupling provided by this invention also includes the following steps: The monthly full-grid head field sequence and the calibrated permeability coefficient field and water yield field are extracted from the results dataset. The unit extraction pulse is applied to each grid using the calibrated groundwater flow numerical model and the head response value of the whole area is recorded. The unit extraction pulse response matrix characterizing the head influence relationship between arbitrary grids is constructed. The final calculated value of monthly groundwater storage change is deducted from the climate-driven natural change component predicted by the multilayer perceptron network. The residual component of storage change caused by human mining activities in each grid is separated. The residual component of storage change in each grid is combined with the response matrix to establish the mining intensity inversion equation. The non-negative least squares algorithm is used to solve the mining intensity inversion equation month by month, and the mining intensity estimate of each grid is constrained to be no less than zero. Finally, the monthly independent mining intensity estimate sequence of each grid is output, and the mining intensity estimate of the grid corresponding to the grid whose reserve change residual is less than the preset residual threshold is reduced to zero. Calculate the multi-year average and standard deviation of the mining intensity estimate for each grid. Grids whose multi-year average exceeds one standard deviation of the regional average mining intensity are marked as grids with high mining intensity. Adjacent grids with high mining intensity are merged into over-extraction hot spot clusters according to spatial connectivity. Extract the boundary range, center location and average mining intensity of each hot spot cluster. Spatial overlay analysis of the spatial extent of each hot spot cluster with land use classification map is performed. The type of mining entity in the hot spot cluster is determined based on the land use type with the highest proportion within the coverage area of ​​the hot spot cluster. Kuness analysis and seasonal index calculation are performed on the monthly variation curve of mining intensity of each hot spot cluster to distinguish between perennial over-extraction areas, seasonal over-extraction areas and intermittent over-extraction areas. The spatial location, mining entity type, temporal pattern characteristics, and average mining intensity of each hot spot cluster are integrated into a spatiotemporal source map of mining intensity. Differentiated control strategies are generated for different types of over-extraction hot spots and pushed out together with the spatiotemporal source map of mining intensity.

[0036] It is important to note that understanding the spatial distribution and temporal patterns of extraction activities is a prerequisite for precise control in the management of regional groundwater over-extraction. Traditional extraction statistics rely on intake metering data, which cannot effectively cover scattered extraction activities that are not metered or have incomplete metering. This step first applies unit extraction pulses to each grid using a calibrated groundwater flow numerical model and records the hydraulic head response values ​​for the entire region, constructing a response matrix characterizing the hydraulic head influence relationship between any grids. The physical significance of this operation is that extraction activities at any location will be transmitted to the surrounding area through the aquifer. The response matrix quantifies this spatial hydraulic transmission relationship into a computable matrix form, providing a physical mapping basis for back-analyzing the independent extraction intensity of each grid from the regional total signal. Subsequently, the final calculated value of monthly groundwater storage changes is subtracted from the climate-driven natural change component predicted by the multilayer perceptron network, separating the residual component of storage changes caused by human extraction. The inversion equation is solved monthly using a non-negative least squares algorithm. The algorithm automatically shrinks the extraction intensity of grids with near-zero residuals to zero, achieving adaptive differentiation between extraction and non-extraction areas. After obtaining the grid-by-grid mining intensity estimate, the grids with consistently high mining intensity are marked and clustered into over-extraction hotspots by calculating the multi-year average and standard deviation. This cluster is then overlaid with the land use classification map to determine the type of mining entity. At the same time, the kurtosis and seasonality index of the monthly variation curves are analyzed to distinguish between perennial, seasonal, and intermittent over-extraction characteristics. The resulting spatiotemporal source map of mining intensity provides management departments with a data-driven technical tool for accurately locating the source of over-extraction, distinguishing the causes of over-extraction, and grasping the timing of over-extraction.

[0037] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in the present invention should be included within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. A method for calculating groundwater storage using gravity satellites based on hydrological model verification coupling, characterized in that, include: We acquire high-grid data of gravity satellite equivalent water column, land surface soil moisture content and canopy water equivalent data, meteorological station precipitation and temperature data, remote sensing vegetation index data, monitoring well water level data and water intake metering data of the target area, and use a meteorological-driven long short-term memory network to fill in the satellite gaps to obtain a multi-source basic dataset. Based on the multi-source basic dataset, the land water storage change sequence and the soil water and canopy water change sequence are obtained. The signal is corrected by empirical mode decomposition into intrinsic mode function components and a spatial regression model is established to obtain the initial groundwater storage change sequence. A groundwater flow numerical model is constructed, with the initial groundwater storage change as the observation field. An ensemble forecast instance is generated based on the disturbance permeability coefficient, specific yield, and rainfall infiltration. The grid hydraulic head change is upscaled to the regional average simulated value with specific yield as the weight. The key parameters and states of the hydrological model are dynamically calibrated to obtain the parameter posterior distribution set. The groundwater flow numerical model is updated and run based on the posterior distribution set of the parameters. After the update, the model generates a monthly full-grid hydraulic head field sequence, and the hydraulic head field is converted into the final calculated value of groundwater storage change, forming a result dataset that has both macro trend and spatiotemporal fine features. Based on the aforementioned dataset, with the final change in groundwater reserves as the dependent variable, and precipitation, temperature, vegetation index, and mining disturbance index constructed from actual extraction volume as independent variables, a multilayer perceptron network is constructed to simulate the combined climate and mining scenarios and generate dynamic response paths for reserves.

2. The method for calculating groundwater storage based on hydrological model verification coupling according to claim 1, characterized in that, The acquisition of gravity satellite equivalent water column high-grid data, land surface soil moisture content and canopy water equivalent data, meteorological station precipitation and temperature data, remote sensing vegetation index data, monitoring well water level data, and water intake metering data for the target area, along with the use of a meteorological-driven long short-term memory network to fill in satellite gaps, yields a multi-source basic dataset, specifically including: Acquire gravity satellite equivalent water column high grid data, land surface soil moisture and canopy water equivalent data, meteorological station precipitation and temperature observation records, vegetation index grid data extracted from remote sensing images, monthly measured water level data of monitoring wells, and monthly actual groundwater extraction sequence recorded by water intake metering facilities in the target area. All raster data are resampled to a unified grid through bilinear interpolation, and station data are interpolated to the same grid nodes through inverse distance weighted interpolation. The time dimension is uniformly aggregated into monthly averages to form an initial multi-source data stack. The multi-year monthly average and standard deviation are calculated for each variable sequence in each grid cell of the initial multi-source data stack. If the current month's value deviates from the corresponding monthly average by more than three times the standard deviation, it is marked as an outlier and set as a missing value. The same criteria for anomaly detection and removal are applied to the monitoring well water level and water intake extraction volume data on their original station sequences to obtain a multi-source sequence containing only true missing values ​​and satellite discontinuities. To address the long-term continuous discontinuities in the equivalent water column height sequence of gravity satellites caused by alternating satellite missions, grid cells are used as the filling objects. At least three consecutive periods before the discontinuity are selected as training samples. Monthly precipitation and monthly temperature are used as input features, and the corresponding equivalent water column height is used as the learning target. A long short-term memory regression prediction network with two hidden layers was constructed and iteratively trained to convergence using the mean squared error loss function and the adaptive moment estimation optimizer. The existing meteorological data during the intermission period was input into the trained network and predicted output month by month to complete the entire intermission period sequence. Simultaneously, short-term sporadic missing data in the vegetation index sequence, soil moisture content sequence, canopy water equivalent sequence, and monitoring well water level sequence were interpolated and completed grid-by-grid using a piecewise cubic spline interpolation method. After all data were processed, they were recombined according to a unified spatial grid and time index to form a multi-source basic dataset.

3. The method for calculating groundwater storage based on hydrological model verification coupling according to claim 1, characterized in that, The process involves obtaining land water storage change sequences and soil and canopy water change sequences based on the multi-source basic dataset, performing empirical mode decomposition into intrinsic mode function components, and establishing a spatial regression model for signal correction to obtain the initial groundwater storage change sequence. Specifically, this includes: Gravity satellite equivalent water column height sequences were extracted from multi-source basic datasets as a signal source for changes in terrestrial water storage. Simultaneously, soil moisture content sequences and canopy water equivalent sequences that had completed spatiotemporal registration were extracted as soil water change sequences and canopy water change sequences. Gaussian filtering is applied to the equivalent water column high sequence of gravity satellites. The time series of each grid is smoothed to suppress north-south striping errors and high-frequency noise, and the filtered land water storage change sequence is output. An integrated empirical mode decomposition algorithm was used to adaptively decompose the filtered land water storage change sequence and soil water change sequence into intrinsic mode function components arranged from high frequency to low frequency and residual trend terms reflecting long-term trends, resulting in two sets of intrinsic mode function components. Pairwise matching is performed on the components at the same frequency position in the two sets of intrinsic mode function components, and the Pearson correlation coefficient is calculated. Component pairs that exceed the preset significance threshold are judged as significantly correlated component pairs at the same frequency. A multivariate linear spatial regression correction model was established with monthly precipitation and monthly temperature as independent variables. The signal leakage and spatial attenuation deviation of the corresponding frequency bands were corrected on each grid. The components after frequency band correction of the regression model were linearly superimposed with the residual trend term to obtain the complete sequence of land water storage changes after time-frequency correction. The initial groundwater storage change sequence is output by subtracting the corresponding soil water change sequence and canopy water change sequence from the corrected terrestrial water storage change sequence month by month, with the equivalent water column height as the unit.

4. The method for calculating groundwater storage based on hydrological model verification coupling according to claim 1, characterized in that, The construction of the groundwater flow numerical model, and the dynamic calibration of the key parameters and states of the hydrological model to obtain the posterior distribution set of parameters, specifically includes: Based on the hydrogeological map of the target area, the lithological distribution data of aquifers and the results of historical pumping tests, the permeability coefficient zone and the water yield zone were divided. The initial values ​​and a priori value ranges of parameters for each zone were determined. The MODFLOW software was used to discretize the study area into an orthogonal grid and set the boundary conditions. The rainfall infiltration recharge and the actual monthly extraction volume were input to complete the initial construction of the groundwater flow numerical model. Using the initial groundwater storage change sequence as the observation field, the Latin hypercube sampling method is used to randomly sample the permeability coefficient and specific yield values ​​of each parameter partition. At the same time, a normal distribution random perturbation is applied to the rainfall infiltration recharge. Several parallel model instances are generated through parameter combination to form an ensemble forecast field. At each assimilation time, all model instances in the set are run to obtain grid-by-grid head forecast values ​​and calculate head change values. The head change value of each grid is multiplied by the corresponding water supply value of the partition with water supply as the weight, and then the average is weighted by grid area. The result is then upscaled to convert the simulated value of regional average groundwater storage change. After all instances have completed the conversion process, the set of simulated values ​​is obtained. The forecast error covariance matrix is ​​calculated based on the set of simulated values. The observation error covariance matrix is ​​constructed by combining the initial groundwater storage change observation value with the preset observation error standard deviation. The Kalman gain matrix is ​​calculated and the state of the grid-by-grid water head forecast value for each instance is updated. At the same time, the permeability coefficient value and water yield value of each zone are updated. The complete process of ensemble operation, upscaling transformation, covariance calculation, gain calculation and state parameter update is carried out month by month along the time series. After all time series assimilation cycles are completed, the final converged permeability coefficient value and water yield value of each instance are extracted to generate the parameter posterior distribution set.

5. The method for calculating groundwater storage based on hydrological model verification coupling according to claim 1, characterized in that, The process involves updating the groundwater flow numerical model based on the posterior distribution set of the parameters. After the update, the model generates a monthly full-grid hydraulic head field sequence. This hydraulic head field is then converted into the final calculated value of groundwater storage changes, forming a dataset that combines macroscopic trends with fine spatiotemporal characteristics. Specifically, this dataset includes: The posterior mean of the permeability coefficient and the posterior mean of the specific yield are extracted from the posterior distribution set of the parameters. The posterior mean is then used to replace the parameter values ​​of the corresponding partitions in the initially constructed groundwater flow numerical model. The grid division and boundary conditions remain unchanged, resulting in a high-fidelity groundwater flow numerical model calibrated with gravity satellite observation data. Driven by monthly rainfall infiltration recharge and monthly actual groundwater extraction, which are centrally covered by multi-source basic datasets throughout the entire study period, the updated groundwater flow numerical model is run month by month. The full-grid head distribution is calculated once at each time step, and the output is a monthly full-grid head field sequence covering the entire study period. The monthly full-grid hydraulic head field sequence is converted into the final calculated value of groundwater storage change. The hydraulic head difference between adjacent months of each grid is calculated and multiplied by the posterior mean of the water supply of the corresponding partition. After weighted averaging by grid area, it is converted into the final calculated value of groundwater storage change on a monthly scale in units of equivalent water column height. The final calculated value of groundwater storage change in the grid where the monitoring well is located is divided by the posterior mean of the water yield of the corresponding grid to obtain the predicted value of water level change. The predicted value of water level change is compared with the measured value of water level change over time, and the root mean square error and Pearson correlation coefficient are calculated to verify the inversion accuracy. The verified monthly groundwater storage change final calculation sequence, the monthly full-grid head field sequence, and the posterior mean of permeability coefficient and water yield in the parameter posterior set are output together to form a result dataset that combines the macroscopic trend accuracy of gravity satellites with the spatiotemporal fine features of distributed hydrological models.

6. The method for calculating groundwater storage based on hydrological model verification coupling according to claim 1, characterized in that, Based on the aforementioned dataset, a multilayer perceptron network is constructed to simulate combined climate and mining scenarios, generating dynamic response paths for reserves. Specifically, this includes: The final calculated value of monthly groundwater storage change was extracted from the results dataset as the output target variable. The monthly precipitation, monthly temperature and vegetation index sequences of the same period were retrieved from the multi-source basic dataset as natural driving factors. The monthly mining disturbance intensity index was constructed by the historical monthly actual mining volume of the target area as a proxy variable for human activities. Time alignment and normalization were performed on all independent and dependent variables. A multilayer perceptron neural network containing an input layer, two hidden layers and an output layer was constructed. The input layer was set with four nodes corresponding to monthly precipitation, monthly average temperature, vegetation index and mining disturbance intensity index, respectively. The hidden layer adopted the modified linear unit activation function, and the output layer was the predicted value of groundwater storage change. The network uses the mean squared error loss function and the adaptive moment estimation optimizer to update the weights iteratively. The historical samples are divided into training set, validation set and test set according to time order. The training process monitors the value of the validation set loss function and implements an early stopping strategy. The training is terminated when the validation loss no longer decreases for several consecutive rounds. After the network training is completed, extract all connection weight matrices from the input layer to the hidden layer and from the hidden layer to the output layer and take their absolute values. For each input factor, calculate the sum of the products of the connection weights that lead to the output node through each hidden node. Divide the sum of the products of each factor by the sum of the products of all factors to obtain the normalized contribution score of each factor to the change in groundwater storage. Several sets of cross-combinations of climate scenarios and mining intensity scenarios are set. The climate scenarios are set according to the frequency level of precipitation, and the mining intensity scenarios are set according to the multiple of the mining disturbance index. The input vector of each scenario is fed into the trained network one by one. The network predicts and outputs the groundwater storage change value of the corresponding scenario on a monthly basis, generating dynamic response paths of storage under different combinations of climate conditions and mining intensity.

Images