A satellite-gravity collaborative forward modeling method for monitoring changes in confined water reserves

Abstract

The present invention discloses a satellite gravity collaborative forward simulation method for monitoring changes in confined water reserves, which comprises the following steps: S1: obtaining boundary and longitude and latitude information of a preset study area, and determining the boundary and longitude and latitude information of a preset background area according to a range that may be affected by an error; S2: obtaining monthly-scale monitoring well water level data, empirical values ​​of water supply degree in the background area, GRACE / GRACE-FO monthly-scale gravity satellite time-varying gravity field data, and non-groundwater component data based on the boundary and longitude and latitude information of the preset background area; S3: obtaining a monthly time series of groundwater reserve changes containing leakage errors; S4: obtaining a monthly time series of phreatic water reserve changes; S5: obtaining a monthly time series of confined water reserve changes on a 0.5° grid containing leakage errors; S6: running a preset iterative forward simulation algorithm; and S7: obtaining a corrected monthly time series of confined water reserve changes in the preset study area.

Description

Technical Field

[0001] The present invention relates to the intersection of gravity satellites and groundwater science, and in particular to a satellite gravity collaborative forward simulation method for monitoring changes in confined water reserves. Background Art

[0002] Groundwater, one of Earth's most important freshwater reserves, consists of both shallow and deep layers. It plays a crucial role in maintaining ecological balance, supporting the sustainable development of human society, and safeguarding human health. Particularly in arid and semi-arid regions and areas with water scarcity, groundwater provides a relatively stable water source, filling the gap left by insufficient surface water and meeting water needs across agriculture, industry, cities, and other sectors. Against the backdrop of strained water supply and demand, groundwater utilization is considered a key strategy for alleviating water stress.

[0003] Deep groundwater (confined water) is more stable than shallow groundwater. However, in areas with high water consumption, with the continuous expansion of human activities and the further intensification of water resource utilization, the long-term over-exploitation of deep groundwater has put groundwater under urgent pressure that cannot be ignored, and some areas have shown varying degrees of ground subsidence. Therefore, it is crucial to quantitatively analyze the total recoverable groundwater reserves in confined aquifer systems. In-depth research on the spatial distribution, recharge mechanism, and mutual relationship of deep groundwater with shallow groundwater will help to more effectively implement strategies for protection and rational utilization, and provide a reliable scientific basis for future water resource management and sustainable development.

[0004] Currently, the main methods for studying deep groundwater include monitoring using traditional monitoring well water level data, geophysical exploration methods, numerical modeling by simulating groundwater flow and recharge, and indirect monitoring using thermal infrared and microwave remote sensing technologies. The first two methods are more effective for small-scale monitoring, but the number and spatial distribution of monitoring points directly affect the research conclusions. While the latter two methods can achieve large-scale monitoring, model construction requires a large amount of prior knowledge and parameters, and the monitoring results are subject to high uncertainty. Combining the advantages and limitations of these methods, this paper proposes a deep groundwater storage monitoring method that combines monitoring well water level data with remote sensing technology.

[0005] In recent years, with the successful launch and continued development of the GRACE (Gravity Recovery and Climate Experiment) and GRACE-FO gravity satellites, jointly developed by the United States and Germany, scholars at home and abroad have conducted extensive and long-term research on monitoring changes in terrestrial and groundwater reserves. As remote sensing satellites currently capable of effectively monitoring changes in subsurface water reserves, they overcome the limitations of traditional ground-based observations and provide valuable data resources for groundwater science. These satellites have become effective tools for monitoring changes in water reserves globally and across large regions.

[0006] During the processing of GRACE / GRACE-FO data, operations such as truncation of spherical harmonic coefficients and filtering can cause varying degrees of inward or outward signal leakage, resulting in leakage errors. Therefore, the data must be corrected before conclusions are drawn to ensure the restoration of normal signals. Among the commonly used correction methods, the effectiveness of additive correction and scale factor correction is heavily dependent on the quantity and quality of hydrological model data. While iterative forward modeling overcomes the former's overreliance on prior knowledge, its ability to restore the spatial resolution of the signal is highly dependent on the uniformity of the spatial distribution of the true signal. Therefore, for groundwater, especially confined water, traditional iterative forward modeling methods struggle to achieve high-quality, high-precision signal recovery. Summary of the Invention

[0007] The present invention provides a satellite gravity collaborative forward simulation method for monitoring changes in pressurized water reserves, which is used to solve the problems existing in the above-mentioned prior art.

[0008] To achieve the above-mentioned object, the present invention provides a satellite gravity collaborative forward simulation method for monitoring changes in confined water reserves, which comprises:

[0009] S1: Obtain the boundary and longitude and latitude information of the preset study area, and determine the boundary and longitude and latitude information of the preset background area according to the range that the error may affect;

[0010] S2: Based on the boundaries and longitude and latitude information of the preset background area, monthly monitoring well water level data, empirical values ​​of background area water supply, GRACE / GRACE-FO monthly gravity satellite time-varying gravity field data, and non-groundwater component data are obtained. Non-groundwater component data include soil moisture data, snow water equivalent data, canopy water data, and surface water storage data.

[0011] S3: The monthly-scale GRACE / GRACE-FO gravity satellite time-varying gravity field data are truncated to a certain order of spherical harmonic coefficients and processed. The average gravity field over a certain period of time is deducted to obtain a monthly time series of terrestrial water storage changes containing leakage errors. At the same time, the remaining non-groundwater component data are expanded into spherical harmonic coefficients and processed in the same way to obtain a monthly time series of non-groundwater component changes. The monthly time series of non-groundwater component changes is deducted from the monthly time series of terrestrial water storage changes containing leakage errors to obtain a monthly time series of groundwater storage changes containing leakage errors.

[0012] S4: Separate the monthly time series of phreatic well water level changes from the monthly time series of confined water well water level changes based on the monitoring well type. Then, multiply the phreatic well water level data with the empirical value of water supply obtained in S2 according to the latitude and longitude of the grid, expand it into spherical harmonic coefficients, and adopt the same data processing method as in S3 to obtain the monthly time series of phreatic reserve changes.

[0013] S5: The monthly time series of groundwater storage changes is deducted from the monthly time series of groundwater storage changes to obtain the monthly time series of confined water storage changes on the 0.5° grid with leakage errors;

[0014] S6: Based on the monthly time scale data of the water level change of the confined water well obtained in S4 and the monthly time series of the confined water storage change with leakage error obtained in S5, the preset iterative forward simulation algorithm is run by giving an arbitrary value within the range of hydrogeological conditions as the simulation initial value of the elastic water release coefficient;

[0015] S7: Based on the preset iterative forward simulation algorithm, the elastic water release coefficient value corresponding to the inversion value of the monthly time series of groundwater storage changes in the preset study area is obtained. The elastic water release coefficient value is multiplied with the monthly time scale data of the water level changes of the pressurized water wells obtained in S4 according to the spatial position to obtain the corrected monthly time series of the pressurized water storage changes in the preset study area.

[0016] In one embodiment of the present invention, in step S2, the background area water supply experience value is based on the spatial distribution map of the preset study area water supply experience value, and the water supply value at the 0.5° grid scale in the background area is calculated by an area-weighted average algorithm.

[0017] In one embodiment of the present invention, in step S3, the time-varying gravity field data of the GRACE / GRACE-FO gravity satellites are processed as follows:

[0018] Use spherical harmonic data truncated to order 60;

[0019] The value of C20 is replaced by the value measured by satellite laser ranging;

[0020] The spherical harmonic data were subjected to Gaussian filtering with a radius of 300 km and combined with P3M10 decorrelation processing, that is, the spherical harmonic coefficients of order 10 and above were fitted once with a third-order polynomial according to odd and even orders respectively and deducted.

[0021] In one embodiment of the present invention, in step S4, the monthly time series of the phreatic well water level change is interpolated within the background area using the inverse distance weighted interpolation method. The unknown point is estimated based on the weighted mean of the neighboring points. The distance between the interpolation point and the known sample point is used as the weight for weighted averaging. The sample point closer to the interpolation point is assigned a larger weight, thereby obtaining a monthly time series of the phreatic well water level change on a spatially continuous 0.5° grid scale within the background area.

[0022] The monthly time series of phreatic well water level changes and the water supply value at the 0.5° grid scale obtained in step S2 are multiplied according to the latitude and longitude of the grid to obtain the monthly time series of phreatic well reserve changes at the 0.5° grid scale. The calculation formula is as follows:

[0023] SGWSA=SGWLA×S y

[0024] Where: SGWSA is the monthly time series of phreatic water reserves change, SGWLA is the monthly time series of phreatic water well water level change, S y is the experience value of water supply degree;

[0025] The monthly time series of phreatic reserves changes is expanded into spherical harmonic coefficients, and forward simulation is performed using the same data processing method as in step S3 to obtain the monthly time series of phreatic reserves changes after forward simulation.

[0026] In one embodiment of the present invention, the calculation formula of step S5 is as follows:

[0027] DGWSA SH =GWSA SH -SGWSA SH

[0028] Where: DGWSA SH is the monthly time series of the confined water storage changes to be corrected, GWSA SH is the monthly time series of groundwater storage changes including leakage errors, SGWSA SH This is the monthly time series of the changes in phreatic reserves after forward simulation.

[0029] In one embodiment of the present invention, when the preset iterative forward simulation algorithm is executed in step S6, whether to stop the iteration is determined by comparing the difference between the monthly time series of the confined water storage change after the iterative simulation and the monthly time series of the confined water storage change containing the leakage error with a given threshold. The calculation formula is as follows:

[0030] DGWSA i =DGWLA×μ i

[0031] Where: DGWSA i is the monthly time series of changes in confined water reserves obtained in the i-th iteration, DGWLA is the monthly time series of changes in the confined water well water level in the study area, μ i is the elastic water release coefficient of the i-th iteration simulation,

[0032] μ i+1 =μ i ±λ

[0033] Where: μ i+1 is the elastic water release coefficient required for the i+1th iteration, and λ is the adjustment step size set by the preset iterative forward simulation algorithm. The judgment basis for adding or subtracting the step size here is the positive or negative difference between the monthly time series of confined water storage changes after the i-th iterative simulation and the monthly time series of confined water storage changes containing leakage errors. The entire iterative forward simulation algorithm ends when the difference is less than the given threshold.

[0034] In one embodiment of the present invention, in step S7, a harmonic analysis algorithm based on least squares fitting is used to decompose the signal to obtain the annual and semi-annual amplitudes, phases, and change trends. The calculation formula is as follows:

[0035]

[0036] Where: ΔH(t) is the monthly time series; t is time; a and b are the constant term and trend term respectively; A i and are the amplitude and phase respectively (i=1, for the anniversary; i=2, for the half-anniversary); ε(t) is the fitting residual.

[0037] The satellite gravity collaborative forward simulation method for monitoring changes in confined water reserves provided by the present invention performs iterative correction by combining gravity satellite and monitoring well water level data, combines macro-monitoring technology with micro-observation technology, and integrates the advantages of remote sensing technology and traditional hydrological monitoring methods to obtain higher-quality signals of confined water layers, providing a new idea for research on confined aquifers. BRIEF DESCRIPTION OF THE DRAWINGS

[0038] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments or the description of the prior art. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0039] Figure 1 This is a flow chart of a satellite gravity collaborative forward simulation method for monitoring changes in confined water reserves according to one embodiment of the present invention;

[0040] Figure 2 A flowchart for obtaining a monthly time series of changes in confined water reserves including leakage errors for this embodiment;

[0041] Figure 3 Flowchart of the iterative forward simulation algorithm of the embodiment based on the gravity satellite inversion of the confined water storage change monthly time series. DETAILED DESCRIPTION

[0042] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. All other embodiments obtained by ordinary technicians in this field based on the embodiments of the present invention without creative work are within the scope of protection of the present invention.

[0043] Figure 1 FIG. 1 is a flow chart of a satellite gravity collaborative forward simulation method for monitoring changes in confined water reserves according to an embodiment of the present invention. Figure 1 As shown, the present invention provides a satellite gravity collaborative forward simulation method for monitoring changes in confined water reserves, which includes:

[0044] S1: Obtain the boundary and longitude and latitude information of the preset study area, and determine the boundary and longitude and latitude information of the preset background area according to the range that the error may affect;

[0045] S2: Based on the boundaries and longitude and latitude information of the preset background area, monthly monitoring well water level data, empirical values ​​of background area water supply, GRACE / GRACE-FO monthly gravity satellite time-varying gravity field data, and non-groundwater component data are obtained. Non-groundwater component data include soil moisture data, snow water equivalent data, canopy water data, and surface water storage data.

[0046] S3: The monthly-scale GRACE / GRACE-FO gravity satellite time-varying gravity field data are truncated to a certain order of spherical harmonic coefficients and processed. The average gravity field over a certain period of time is deducted to obtain a monthly time series of terrestrial water storage changes containing leakage errors. At the same time, the remaining non-groundwater component data are expanded into spherical harmonic coefficients and processed in the same way to obtain a monthly time series of non-groundwater component changes. The monthly time series of non-groundwater component changes is deducted from the monthly time series of terrestrial water storage changes containing leakage errors to obtain a monthly time series of groundwater storage changes containing leakage errors.

[0047] S4: Separate the monthly time series of phreatic well water level changes from the monthly time series of confined water well water level changes based on the monitoring well type. Then, multiply the phreatic well water level data with the empirical value of water supply obtained in S2 according to the latitude and longitude of the grid, expand it into spherical harmonic coefficients, and adopt the same data processing method as in S3 to obtain the monthly time series of phreatic reserve changes.

[0048] S5: The monthly time series of groundwater storage changes is deducted from the monthly time series of groundwater storage changes to obtain the monthly time series of confined water storage changes on the 0.5° grid with leakage errors;

[0049] S6: Based on the monthly time scale data of the water level change of the confined water well obtained in S4 and the monthly time series of the confined water storage change with leakage error obtained in S5, the preset iterative forward simulation algorithm is run by giving an arbitrary value within the range of hydrogeological conditions as the simulation initial value of the elastic water release coefficient;

[0050] S7: Based on the preset iterative forward simulation algorithm, the elastic water release coefficient value corresponding to the inversion value of the monthly time series of groundwater storage changes in the preset study area is obtained. The elastic water release coefficient value is multiplied with the monthly time scale data of the water level changes of the pressurized water wells obtained in S4 according to the spatial position to obtain the corrected monthly time series of the pressurized water storage changes in the preset study area.

[0051] In one embodiment of the present invention, in step S2, the background area water supply experience value is based on the spatial distribution map of the preset study area water supply experience value, and the water supply value at the 0.5° grid scale in the background area is calculated by an area-weighted average algorithm.

[0052] In one embodiment of the present invention, in step S3, the time-varying gravity field data of the GRACE / GRACE-FO gravity satellites are processed as follows:

[0053] Use spherical harmonic data truncated to order 60;

[0054] The value of C20 is replaced by the value measured by satellite laser ranging;

[0055] The spherical harmonic data were subjected to Gaussian filtering with a radius of 300 km and combined with P3M10 decorrelation processing, that is, the spherical harmonic coefficients of order 10 and above were fitted once with a third-order polynomial according to odd and even orders respectively and deducted.

[0056] In one embodiment of the present invention, in step S4, the monthly time series of the phreatic well water level change is interpolated within the background area using the inverse distance weighted interpolation method. The unknown point is estimated based on the weighted mean of the neighboring points. The distance between the interpolation point and the known sample point is used as the weight for weighted averaging. The sample point closer to the interpolation point is assigned a larger weight, thereby obtaining a monthly time series of the phreatic well water level change on a spatially continuous 0.5° grid scale within the background area.

[0057] The monthly time series of phreatic well water level changes and the water supply value at the 0.5° grid scale obtained in step S2 are multiplied according to the latitude and longitude of the grid to obtain the monthly time series of phreatic well reserve changes at the 0.5° grid scale. The calculation formula is as follows:

[0058] SGWSA=SGWLA×S y

[0059] Where: SGWSA is the monthly time series of phreatic water reserves change, SGWLA is the monthly time series of phreatic water well water level change, S y is the experience value of water supply degree;

[0060] The monthly time series of phreatic reserves changes is expanded into spherical harmonic coefficients, and forward simulation is performed using the same data processing method as in step S3 to obtain the monthly time series of phreatic reserves changes after forward simulation.

[0061] In one embodiment of the present invention, the calculation formula of step S5 is as follows:

[0062] DGWSA SH =GWSA SH -SGWSA SH

[0063] Where: DGWSA SH is the monthly time series of the confined water storage changes to be corrected, GWSA SH is the monthly time series of groundwater storage changes including leakage errors, SGWSA SH This is the monthly time series of the changes in phreatic reserves after forward simulation.

[0064] In one embodiment of the present invention, when the preset iterative forward simulation algorithm is executed in step S6, whether to stop the iteration is determined by comparing the difference between the monthly time series of the confined water storage change after the iterative simulation and the monthly time series of the confined water storage change containing the leakage error with a given threshold. The calculation formula is as follows:

[0065] DGWSA i =DGWLA×μ i

[0066] Where: DGWSA iis the monthly time series of changes in confined water reserves obtained in the i-th iteration, DGWLA is the monthly time series of changes in the confined water well water level in the study area, μ i is the elastic water release coefficient of the i-th iteration simulation,

[0067] μ i+1 =μ i ±λ

[0068] Where: μ i+1 is the elastic water release coefficient required for the i+1th iteration, and λ is the adjustment step size set by the preset iterative forward simulation algorithm. The judgment basis for adding or subtracting the step size here is the positive or negative difference between the monthly time series of confined water storage changes after the i-th iterative simulation and the monthly time series of confined water storage changes containing leakage errors. The entire iterative forward simulation algorithm ends when the difference is less than the given threshold.

[0069] In one embodiment of the present invention, in step S7, a harmonic analysis algorithm based on least squares fitting is used to decompose the signal to obtain the annual and semi-annual amplitudes, phases, and change trends. The calculation formula is as follows:

[0070]

[0071] Where: ΔH(t) is the monthly time series; t is time; a and b are the constant term and trend term respectively; A i and are the amplitude and phase respectively (i=1, for the anniversary; i=2, for the half-anniversary); ε(t) is the fitting residual.

[0072] The present invention is described below with a specific embodiment:

[0073] Step S1: Obtain the boundary and longitude and latitude information of the study area, and determine the boundary and longitude and latitude information of the background area according to the range that the error may affect;

[0074] In this embodiment, based on the distribution of monitoring wells, the boundaries and longitude and latitude information of the study area are preset to be the longitude and latitude grid data of 0.5° within the North China Plain, totaling 28 grids. The boundaries and longitude and latitude information of the background area are set to the longitude and latitude range of 29.75°N-42.25°N and 108.25°E-124.75°E, representing the entire background area. This setting range is for illustration only and does not limit the scope of implementation of the present invention.

[0075] Step S2: Based on the boundaries and longitude and latitude information of the preset background area, monthly monitoring well water level data, empirical values ​​of background area water supply, monthly GRACE / GRACE-FO gravity satellite time-varying gravity field data, and non-groundwater component data are obtained. The non-groundwater component data includes soil moisture data, snow water equivalent data, canopy water data, and surface water storage data.

[0076] This step primarily involves preprocessing the data required for the implementation of this method. The GRACE / GRACE-FO monthly-scale time-varying gravity field data from the gravity satellites can be directly downloaded from the Internet. Based on previous studies of uncertainty in models released by various data production agencies worldwide and in China, this example prefers the GRACE Science Data System (SDS) version 6 Earth Gravity Model product (CSR RL06 Level-2) released by the Center for Space Research (CSR) at the University of Texas at Austin. Other data agencies available in this system include the Jet Propulsion Laboratory (JPL) of NASA and the Deutschen GeoForschungsZentrum (GFZ) in Potsdam, Germany.

[0077] It should be noted that the acquisition of the Earth's gravity field data in the present invention can also select other gravity satellite products according to actual scientific research needs, such as Tianqin-1, including but not limited to GRACE / GRACE-FO gravity satellites. This setting condition is only for illustration and does not limit the implementation conditions of the present invention.

[0078] The processing of other required data mainly includes the following steps:

[0079] Step S201: Monthly monitoring well water level data are acquired based on the pre-set background area boundaries and longitude and latitude information, and data for missing months are interpolated. In this step, considering the uncertainty introduced by conventional interpolation methods, a least-squares fit with the GWS component of the CLSM2.2 model is used to minimize the sum of squared errors. Data are then filled in, and monitoring wells with missing values ​​greater than half of the total number of months are removed to achieve quality control of the monitoring well water level data. After deducting the average water level from January 2003 to December 2016, a monthly time series of water level changes for all monitoring wells is obtained.

[0080] Step S202: Based on the spatial distribution map of the empirical value of water supply degree in the North China Plain, the water supply degree value at the 0.5° grid scale in the background area is calculated by using an area-weighted average algorithm;

[0081] Step S203: The non-groundwater components referred to herein include, but are not limited to, soil water, snow water equivalent, and canopy water. The soil moisture data can be sourced from the Common Land Model (CLM) soil moisture dataset of the Global Land Data Assimilation System (GLDAS), or from measured data or other hydrological model data. Surface water storage data can be derived from reservoir storage data published in the China Water Resources Bulletin and relevant river basin water resources bulletins. The water supply value within the study area can be derived from the comprehensive water supply value in the Survey and Evaluation of Sustainable Groundwater Utilization in the North China Plain. The present invention does not limit the above data, and the selected data can be adjusted based on the specific conditions of the study area.

[0082] Preferably, in this embodiment, the arithmetic mean of the multi-model simulation values ​​is selected for this part of the data to obtain a monthly time series of non-groundwater component reserve changes at a 0.5° grid scale. The models used include: Noah of GLDAS-2.1, Variable Infiltration Capacity (VIC), and monthly solution data of the Comprehensive Land Surface Model (CLSM).

[0083] Of course, the definition of non-groundwater components needs to be adjusted according to actual data calculation requirements, and should not be limited to the contents listed in this embodiment. Those skilled in the art can also make appropriate expansions, selections, and combinations based on the needs of specific model calculations.

[0084] Step S3: The monthly-scale GRACE / GRACE-FO gravity satellite time-varying gravity field data are truncated to spherical harmonic coefficients of a certain order and processed. The average gravity field over a certain period of time is deducted to obtain a monthly time series of terrestrial water storage changes containing leakage errors. At the same time, the remaining non-groundwater component data are expanded into spherical harmonic coefficients and processed in the same way to obtain a monthly time series of non-groundwater component changes. The monthly time series of non-groundwater component changes is deducted from the monthly time series of terrestrial water storage changes containing leakage errors to obtain a monthly time series of groundwater storage changes containing leakage errors.

[0085] The gravity potential coefficient, C20, reflects the gravity potential perturbation caused by the Earth's oblateness. However, the GRACE / GRACE-FO gravity satellites are insensitive to this perturbation, resulting in poor solution accuracy. Therefore, it is necessary to replace it. Furthermore, since gravity satellites primarily monitor medium- and long-wavelength signals, their accuracy in identifying short-wavelength signals is low. Therefore, in addition to the true signal, the signals acquired by the GRACE / GRACE-FO gravity satellites also contain a significant amount of noise. As the spatial scale of the study area decreases, the impact of high-order signal noise becomes increasingly pronounced. Therefore, filtering and smoothing the data are necessary.

[0086] Preferably, in order to better weaken the systematic error of the GRACE / GRACE-FO gravity satellite and remove the "stripe" feature in the north-south radial direction, in addition to the filtering operation, a decorrelation operation can also be adopted, that is, a combination of odd and even order coefficients in different ways is used, and the related errors are fitted with a sliding polynomial and deducted from the original coefficients.

[0087] Generally, the conventional process for processing the time-varying gravity field of the GRACE / GRACE-FO gravity satellite mainly includes the following steps: expansion and truncation of spherical harmonic coefficients, removal of low-order terms, replacement of C20 terms, filtering and decorrelation processing, where filtering methods include but are not limited to Gaussian filtering, DDK filtering, FAN filtering, etc.

[0088] As a preferred approach, GRACE / GRACE-FO gravity satellite data can be processed using spherical harmonics truncated to 60th order. To minimize noise and preserve the true signal, subsequent data processing is performed: the C20 term is replaced with values ​​obtained from satellite laser ranging. The spherical harmonics are then Gaussian filtered with a radius of 300 km and combined with P3M10 decorrelation to reduce the impact of high-order noise. Specifically, spherical harmonic coefficients of order 10 and above are fitted once with a third-order polynomial for both odd and even orders, and then subtracted. The Gaussian filter is essentially a spatially averaged weighted function, a low-pass filter that effectively suppresses high-order noise. The filtered spherical harmonic coefficients are converted into equivalent water height grid data, and finally, the average gravity field from January 2003 to December 2016 is subtracted to obtain a monthly time series of terrestrial water storage changes at a 0.5° grid scale, including leakage errors. Of course, other existing processing algorithms can also be used to reduce noise on the truncated data, tailored to the specific characteristics of the study area, but this will not be elaborated here.

[0089] Forward modeling was performed on the water storage data of the non-groundwater components in S2. Here, the spherical harmonic coefficients were expanded to the 60th order cutoff. The same processing method as the above data was adopted and summed to obtain the monthly time series of non-groundwater component storage changes at the 0.5° grid scale after forward modeling.

[0090] The monthly time series of non-groundwater component reserves changes after forward simulation is deducted from the monthly time series of terrestrial water storage changes containing leakage errors. Preferably, a difference process is directly adopted here to obtain the monthly time series of groundwater storage changes on the 0.5° grid scale containing leakage errors. The calculation formula is as follows:

[0091] GWSA=TWSA-SMSA-SWEA-CWSA

[0092] Where TWSA is the monthly time series of terrestrial water storage changes obtained by inverting GRACE / GRACE-FO gravity satellite data; SMSA, SWEA, and CWSA are the monthly time series of soil water storage changes, snow water equivalent changes, and canopy water storage changes after forward simulation, respectively; and GWSA is the monthly time series of groundwater storage changes separated from the GRACE / GRACE-FO gravity satellite time-varying signals.

[0093] It should be noted that, depending on the different preset study areas, the monthly time series of groundwater storage changes in the background area outside the preset study area can use hydrological model data, measured data or data released in the bulletin.

[0094] Step S4: Separate the monthly time series of phreatic well water level changes from the monthly time series of confined water well water level changes based on the monitoring well type. Then, multiply the phreatic well water level data by the empirical value of water supply degree obtained in S2 according to the latitude and longitude of the grid, expand it into spherical harmonic coefficients, and adopt the same data processing method as in S3 to obtain the monthly time series of phreatic reserve changes.

[0095] In this embodiment, due to the limited distribution range of monitoring wells, the inverse distance weighted (IDW) interpolation method is adopted to interpolate the monthly time series of submersible well water level changes within the background area. That is, the unknown point is estimated based on the weighted mean of the neighboring points, and the distance between the interpolation point and the known sample point is used as the weight for weighted averaging. The closer the sample point is to the interpolation point, the greater the weight is assigned to it, so as to obtain the monthly time series of submersible well water level changes on a spatially continuous 0.5° grid scale within the background area.

[0096] The above interpolation processing may adopt, for example, the Kriging interpolation method, the sub-region averaging method, etc., or other interpolation methods in the prior art.

[0097] Then, the monthly time series of the phreatic well water level change is multiplied by the water supply value at the 0.5° grid scale obtained in S2 according to the latitude and longitude of the grid, and the calculation formula for the monthly time series of the phreatic well reserve change at the 0.5° grid scale is obtained as follows:

[0098] SGWSA=SGWLA×S y

[0099] Where: SGWSA is the monthly time series of phreatic water reserves change, SGWLA is the monthly time series of phreatic water well water level change, S y It is the experience value of water supply degree.

[0100] Then, the monthly time series of phreatic reserves change is expanded into spherical harmonic coefficients, and the forward simulation is performed using the same data processing method as in S3 to obtain the monthly time series of phreatic reserves change after forward simulation.

[0101] S5: The monthly time series of groundwater storage changes is deducted from the monthly time series of groundwater storage changes to obtain the monthly time series of confined water storage changes on a 0.5° grid scale that includes leakage errors;

[0102] like Figure 2 The flowchart of this embodiment is to obtain the monthly time series of changes in confined water reserves containing leakage errors. In this embodiment, the difference between the two is directly adopted to obtain the monthly time series of changes in confined water reserves on a 0.5° grid scale containing leakage errors in the background area. The calculation formula is as follows:

[0103] DGWSA SH =GWSA SH -SGWSA SH

[0104] Where: DGWSA SH is the monthly time series of the confined water storage changes to be corrected, GWSA SH is the monthly time series of groundwater storage changes including leakage errors, SGWSA SH is the monthly time series of phreatic reserves changes after forward simulation;

[0105] S6: Based on the monthly time scale data of water level changes in the confined water wells in the study area obtained in S4 and the monthly time series of confined water storage changes with leakage errors obtained in S5, the preset iterative forward simulation algorithm is run by giving any value within the range of hydrogeological conditions as the simulation initial value of the elastic water release coefficient;

[0106] like Figure 3 The flowchart of the iterative forward simulation algorithm of this embodiment based on the gravity satellite inversion of the confined water storage change monthly time series is shown. The difference between the confined water storage change monthly time series after iterative simulation and the confined water storage change monthly time series containing leakage error is compared with the size of a given threshold to determine whether to stop the iteration. The main calculation formula involved is as follows:

[0107] DGWSA i =DGWLA×μ i

[0108] Where: DGWSAi is the monthly time series of changes in confined water reserves obtained in the i-th iteration, DGWLA is the monthly time series of changes in the confined water well water level in the study area, μ i is the elastic water release coefficient of the i-th iteration simulation.

[0109] μ i+1 =μ i ±λ

[0110] Where: μ i+1 is the elastic water release coefficient required for the i+1th iteration, and λ is the adjustment step size set by the preset iterative forward simulation algorithm. The judgment basis for adding or subtracting the step size here is the positive or negative difference between the monthly time series of confined water storage changes after the i-th iterative simulation and the monthly time series of confined water storage changes containing leakage errors. The entire iterative forward simulation algorithm ends when the difference is less than the given threshold.

[0111] In this example, with reference to the parameters of the confined aquifer in the North China Plain, the elastic water release coefficient is initially set to a range of 10⁻³ to 10⁻², with 0.005 used as the initial value for the elastic water release coefficient simulation, and a step size of 0.0001. It should be noted that the above-given range, initial value, and step size may be slightly adjusted depending on the geological conditions of the study area. These settings are for illustrative purposes only and do not limit the implementation conditions of the present invention.

[0112] S7: Based on the preset iterative forward simulation algorithm, the elastic water release coefficient value corresponding to the inversion value of the monthly time series of groundwater storage changes in the preset study area is obtained. The elastic water release coefficient value is multiplied by the monthly time scale data of the water level changes of the confined water wells obtained in S4 according to the spatial position to obtain the corrected monthly time series of the confined water storage changes in the preset study area;

[0113] Based on the preset iterative forward simulation algorithm, when the difference in S6 is less than the given threshold, the iteration ends, and the simulated value of the monthly time series of the change in the preset confined water reserves in the preset study area is used as the inversion value of the monthly time series of the change in the preset confined water reserves in the preset study area, and the simulated value of the elastic water release coefficient in the preset study area is used as the water release coefficient value corresponding to the inversion value of the groundwater reserves change trend in the preset study area.

[0114] Finally, the elastic water release coefficient value corresponding to the inversion value of the monthly time series of the change in the confined water storage in the preset study area is multiplied by the monthly time series of the change in the confined water level obtained by S4 to obtain the inverted monthly time series of the change in the confined water storage in the preset study area.

[0115] In a preferred embodiment, in order to better evaluate the temporal changes of each component, including long-term trends, interannual fluctuations, cyclical changes and other signals, a harmonic analysis algorithm based on least squares fitting can be used to decompose the signal to obtain the annual and semi-annual amplitudes, phases and change trends. The calculation formula is as follows:

[0116]

[0117] Where: ΔH(t) is the monthly time series; t is time; a and b are the constant term and trend term respectively; A i and are the amplitude and phase respectively (i=1, for the anniversary; i=2, for the half-anniversary); ε(t) is the fitting residual.

[0118] The satellite gravity collaborative forward simulation method for monitoring changes in confined water reserves provided by the present invention performs iterative correction by combining gravity satellite and monitoring well water level data, combines macro-monitoring technology with micro-observation technology, and integrates the advantages of remote sensing technology and traditional hydrological monitoring methods to obtain higher-quality signals of confined water layers, providing a new idea for research on confined aquifers.

[0119] Those skilled in the art will appreciate that the accompanying drawings are merely schematic diagrams of an embodiment, and the modules or processes in the accompanying drawings are not necessarily required to implement the present invention.

[0120] Those skilled in the art will appreciate that the modules in the apparatuses of the embodiments may be distributed in the apparatuses of the embodiments as described in the embodiments, or may be located in one or more apparatuses different from the embodiments with corresponding changes. The modules in the above embodiments may be combined into one module or further divided into multiple sub-modules.

[0121] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit it. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the aforementioned embodiments, or make equivalent replacements for some of the technical features therein. However, these modifications or replacements do not deviate the essence of the corresponding technical solutions from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A satellite gravity collaborative forward modeling method for monitoring changes in confined water reserves, characterized in that: include: S1: Obtain the boundary and longitude and latitude information of the preset study area, and determine the boundary and longitude and latitude information of the preset background area according to the range that the error may affect; S2: Based on the boundaries and longitude and latitude information of the preset background area, monthly monitoring well water level data, empirical values ​​of background area water supply, GRACE / GRACE-FO monthly gravity satellite time-varying gravity field data, and non-groundwater component data are obtained. Non-groundwater component data include soil moisture data, snow water equivalent data, canopy water data, and surface water storage data. S3: The monthly-scale GRACE / GRACE-FO gravity satellite time-varying gravity field data are truncated to a certain order of spherical harmonic coefficients and processed. The average gravity field over a certain period of time is deducted to obtain a monthly time series of terrestrial water storage changes containing leakage errors. At the same time, the remaining non-groundwater component data are expanded into spherical harmonic coefficients and processed in the same way to obtain a monthly time series of non-groundwater component changes. The monthly time series of non-groundwater component changes is deducted from the monthly time series of terrestrial water storage changes containing leakage errors to obtain a monthly time series of groundwater storage changes containing leakage errors. S4: Separate the monthly time series of phreatic well water level changes from the monthly time series of confined water well water level changes based on the monitoring well type. Then, multiply the phreatic well water level data with the empirical value of water supply obtained in S2 according to the latitude and longitude of the grid, expand it into spherical harmonic coefficients, and adopt the same data processing method as in S3 to obtain the monthly time series of phreatic reserve changes. S5: The monthly time series of groundwater storage changes is deducted from the monthly time series of groundwater storage changes to obtain the monthly time series of confined water storage changes on the 0.5° grid with leakage errors; S6: Based on the monthly time scale data of the water level change of the confined water well obtained in S4 and the monthly time series of the confined water storage change with leakage error obtained in S5, the preset iterative forward simulation algorithm is run by giving an arbitrary value within the range of hydrogeological conditions as the simulation initial value of the elastic water release coefficient; S7: Based on the preset iterative forward simulation algorithm, the elastic water release coefficient value corresponding to the inversion value of the monthly time series of groundwater storage changes in the preset study area is obtained. The elastic water release coefficient value is multiplied with the monthly time scale data of the water level changes of the pressurized water wells obtained in S4 according to the spatial position to obtain the corrected monthly time series of the pressurized water storage changes in the preset study area.

2. The satellite gravity collaborative forward simulation method for monitoring changes in confined water reserves according to claim 1 is characterized in that: In step S2, the empirical value of water supply degree in the background area is based on the spatial distribution map of the empirical value of water supply degree in the preset study area, and the water supply degree value at the 0.5° grid scale in the background area is calculated by the area weighted average algorithm.

3. The satellite gravity collaborative forward simulation method for monitoring changes in confined water reserves according to claim 1 is characterized in that: In step S3, the time-varying gravity field data of the GRACE / GRACE-FO gravity satellites are processed as follows: Use spherical harmonic data truncated to order 60; The value of C20 is replaced by the value measured by satellite laser ranging; The spherical harmonic data were subjected to Gaussian filtering with a radius of 300 km and combined with P3M10 decorrelation processing, that is, the spherical harmonic coefficients of order 10 and above were fitted once with a third-order polynomial according to odd and even orders respectively and deducted.

4. The satellite gravity collaborative forward simulation method for monitoring changes in confined water reserves according to claim 2 is characterized in that: In step S4, the monthly time series of the phreatic well water level change is interpolated within the background area using the inverse distance weighted interpolation method. The unknown point is estimated based on the weighted mean of the neighboring points. The distance between the interpolation point and the known sample point is used as the weight for weighted averaging. The sample point closer to the interpolation point is given a larger weight, so as to obtain the monthly time series of the phreatic well water level change on a spatially continuous 0.5° grid scale within the background area. The monthly time series of phreatic well water level changes and the water supply value at the 0.5° grid scale obtained in step S2 are multiplied according to the latitude and longitude of the grid to obtain the monthly time series of phreatic well reserve changes at the 0.5° grid scale. The calculation formula is as follows: SGWSA=SGWLA×S y Where: SGWSA is the monthly time series of phreatic water reserves change, SGWLA is the monthly time series of phreatic water well water level change, S y is the experience value of water supply degree; The monthly time series of phreatic reserves changes is expanded into spherical harmonic coefficients, and forward simulation is performed using the same data processing method as in step S3 to obtain the monthly time series of phreatic reserves changes after forward simulation.

5. The satellite gravity collaborative forward simulation method for monitoring changes in confined water reserves according to claim 1 is characterized in that: The calculation formula of step S5 is as follows: DGWSA SH =GWSA SH -SGWSA SH Where: DGWSA SH is the monthly time series of the confined water storage changes to be corrected, GWSA SH is the monthly time series of groundwater storage changes including leakage errors, SGWSA SH This is the monthly time series of the changes in phreatic reserves after forward simulation.

6. The satellite gravity collaborative forward simulation method for monitoring changes in confined water reserves according to claim 1 is characterized in that: When running the preset iterative forward simulation algorithm in step S6, whether to stop the iteration is determined by comparing the difference between the monthly time series of the confined water storage changes after iterative simulation and the monthly time series of the confined water storage changes containing leakage errors with a given threshold. The calculation formula is as follows: DGWSA i =DGWLA×μ i Where: DGWSA i is the monthly time series of changes in confined water reserves obtained in the i-th iteration, DGWLA is the monthly time series of changes in the confined water well water level in the study area, μ i is the elastic water release coefficient of the i-th iteration simulation, m i+1 =μ i ±λ Where: μ i+1 is the elastic water release coefficient required for the i+1th iteration, and λ is the adjustment step size set by the preset iterative forward simulation algorithm. The judgment basis for adding or subtracting the step size here is the positive or negative difference between the monthly time series of confined water storage changes after the i-th iterative simulation and the monthly time series of confined water storage changes containing leakage errors. The entire iterative forward simulation algorithm ends when the difference is less than the given threshold.

7. The satellite gravity collaborative forward simulation method for monitoring changes in confined water reserves according to claim 1 is characterized in that: Step S7: Decompose the signal based on the harmonic analysis algorithm of least squares fitting to obtain the annual and semi-annual amplitudes, phases and change trends. The calculation formula is as follows: Where: ΔH(t) is the monthly time series; t is time; a and b are the constant term and trend term respectively; A i and are amplitude and phase respectively, i=1 is the anniversary, i=2 is the half-anniversary; ε(t) is the fitting residual.

Images