Method for inverting groundwater reserves by combining GNSS, GRACE and InSAR technologies

Abstract

The invention relates to the field of satellite gravity detection and application, in particular to a method for inverting groundwater reserves by combining GNSS, GRACE and InSAR technologies, which comprises the following steps: acquiring GNSS vertical displacement data, GRACE equivalent water height data and InSAR vertical deformation rate data; the GNSS vertical displacement data, the GRACE equivalent water height data and the InSAR vertical deformation rate data are associated to the total water reserve change, and a joint inversion equation set is constructed; the uncertainty of GNSS, InSAR and GRACE observation trends is quantified, weight coefficients corresponding to GNSS vertical displacement data, GRACE equivalent water height data and InSAR vertical deformation rate data are determined, and observation data vectors are constructed; laplacian roughness constraint is introduced, a joint inversion equation set is solved through step-by-step iterative calculation according to a weighted least square method by observing data vectors, and high-resolution groundwater reserve changes are determined. By combining GNSS, GRACE and InSAR technologies, the limitation of a traditional single technology is overcome, advantage complementation of multi-source data is realized, and the accuracy of underground water reserve inversion is remarkably improved.

Description

Technical Field

[0001] The present application relates to the field of satellite gravity detection and application, and more specifically, to a method for inverting groundwater reserves by combining GNSS, GRACE and InSAR technologies. Background Art

[0002] Against the backdrop of climate change and intensified human activities, accurate monitoring of river basin water resources is crucial for hydrological cycle research and water resources management. Currently, the Gravity Recovery and Climate Experiment (GRACE) / GRACE-FO satellites infer water storage from changes in the gravity field. However, their spatial resolution (approximately 330 km) and temporal resolution (monthly scale) limit their ability to monitor high-frequency or localized changes. Furthermore, filtering operations during data processing can lead to signal loss, and data gaps occur during satellite missions (e.g., the 11-month data gap between GRACE and GRACE-FO), further compromising continuity.

[0003] The Global Navigation Satellite System (GNSS), with its near-real-time, high-precision, and sensitivity to shortwave signals, can effectively complement GRACE. However, its inversion accuracy is significantly affected by the spatial distribution of stations. Meanwhile, Synthetic Aperture Radar (InSAR) technology offers the advantages of large-scale coverage and high spatial resolution, enabling the capture of detailed surface deformation. In recent years, integrating GNSS, GRACE, and InSAR has become a research hotspot in groundwater monitoring. However, existing methods often neglect the coupling effects of poroelastic signals, leading to uncertainty in inversion results. Summary of the Invention

[0004] In view of this, the present application provides a method for inverting groundwater reserves by combining GNSS, GRACE and InSAR technologies to solve the problems of low resolution, data discontinuity and missing poroelastic signals existing in traditional single technologies, and to achieve complementary advantages of multi-source data.

[0005] The technical solutions provided in this application are as follows:

[0006] A method for inverting groundwater reserves using GNSS, GRACE, and InSAR technologies, comprising:

[0007] Acquire GNSS vertical displacement data, GRACE equivalent water height data, and InSAR vertical deformation rate data;

[0008] The GNSS vertical displacement data, GRACE equivalent water height data, and InSAR vertical deformation rate data are correlated to the total water storage change using Green's function to construct a joint inversion equation set;

[0009] Quantify the uncertainties of GNSS, InSAR, and GRACE observation trends respectively, determine the weight coefficients corresponding to the GNSS vertical displacement data, GRACE equivalent water height data, and InSAR vertical deformation rate data, and construct an observation data vector;

[0010] The Laplace roughness constraint is introduced, and the joint inversion equations are solved step by step iterative calculation using the weighted least squares method through the observation data vector to determine the high-resolution groundwater storage changes.

[0011] In one possible implementation, after obtaining GNSS vertical displacement data, GRACE equivalent water height data, and InSAR vertical deformation rate data, the method further includes:

[0012] The GNSS vertical displacement data, equivalent water height data and vertical deformation rate data are preprocessed, including removing some stations of the GNSS vertical displacement data and projecting the InSAR vertical deformation rate data using an affine transformation; including converting the InSAR vertical deformation rate data into the WGS84 ellipsoidal height coordinate system using an affine transformation, aligning it with the GNSS vertical displacement, and controlling the InSAR vertical deformation rate data to remain consistent with the reference frame of the GNSS deformation.

[0013] In one possible implementation, the use of Green's function to correlate the GNSS vertical displacement data, GRACE equivalent water height data, and InSAR vertical deformation rate data to the total water storage change includes:

[0014] Based on the layered earth model, the vertical displacement response coefficient is determined to characterize the deformation and gravity changes caused by surface loads;

[0015] Calculate the poroelastic coefficient α and construct the poroelastic Green function matrix, which includes two sub-matrices G poro and G elastic ; Among them, the sub-matrix G poro It is used to describe the contribution of groundwater storage change ΔGWS to InSAR deformation and is expressed as:

[0016]

[0017] Submatrix G elastic It is used to describe the contribution of the total water storage change ΔTWS to the GNSS vertical displacement and is expressed as:

[0018]

[0019] Where r is the distance from the center of the grid cell to the observation point, R is the average radius of the earth, and ρ avgis the average density of the earth, α is the Biot-Willis coefficient, ν is the Poisson's ratio of the aquifer, E is the Young's modulus, ρ w is the water density, g is the acceleration of gravity, r0 is the impact radius, h n ' is the vertical displacement response coefficient;

[0020] The poroelastic Green's function matrix is ​​expressed as:

[0021]

[0022] Where I is the identity matrix, which is used to relate the total water storage TWS of GRACE to groundwater and other components.

[0023] In one possible implementation, constructing the joint observation equation includes:

[0024] According to the poroelastic Green's function matrix, the joint observation equation is constructed and expressed as:

[0025]

[0026] Where ΔTWS is the change in total water storage, which is decomposed into the change in groundwater storage ΔGWS and other components ΔOther.

[0027] In one possible implementation, quantifying the uncertainty of GNSS, InSAR, and GRACE observation trends includes:

[0028] Resample the GNSS vertical displacement data and add random noise with a mean of 0 and a standard deviation of 1 mm / year to each resampled data point. Fit a linear trend and record the slope. Repeat this process multiple times to obtain the distribution of trend values. Determine the uncertainty of the GNSS observation trend based on the standard deviation.

[0029] The InSAR vertical deformation rate data were spatially resampled. Random noise drawn from the error distribution was added to each pixel, and the regional average trend was calculated. This was repeated multiple times to determine the uncertainty of the InSAR observation trend.

[0030] The leakage effect is simulated using the GRACE observation operator. The simulated leakage signal is compared with the original signal, and the standard deviation of the leakage error is calculated. Random values ​​are extracted from the leakage error distribution and added to the GRACE trend. Multiple trend estimates affected by leakage are repeatedly generated to calculate the uncertainty of the GRACE observation trend.

[0031] In one possible implementation, determining weight coefficients corresponding to the GNSS vertical displacement data, the GRACE equivalent water height data, and the InSAR vertical deformation rate data includes:

[0032] Constructing corresponding covariance matrices based on the uncertainties of the GNSS, InSAR, and GRACE observation trends;

[0033] Based on the smoothness of the residuals and solutions of the L-curve equilibrium data fitting, the weight coefficient distribution is determined, which is expressed as:

[0034]

[0035] Where C GNSS , C InSAR , C GRACE are the covariance matrices, representing the uncertainties of GNSS, InSAR, and GRACE observations, respectively. The relative weights (a, b, c) are used to balance the strengths of GNSS, InSAR, and GRACE observations, respectively.

[0036] In one possible implementation, a Laplace roughness constraint is introduced, and the joint inversion equations are solved by step-by-step iterative calculation using the weighted least squares method using the observation data vector to determine the high-resolution groundwater storage change, including:

[0037] Construct the Laplace smoothing matrix, expressed as:

[0038]

[0039] According to the Laplace smoothing matrix, the Laplace operator D is introduced as a spatial smoothing constraint to construct the Laplace roughness constraint, which is expressed as:

[0040] min||G·xd|| 2 +λ||D·x|| 2

[0041] Where D is the Laplace operator matrix, G is the Green's function matrix, λ is the smoothing factor determined by cross-validation, d is the observed data vector, and x is the vector to be solved, which is used to determine high-resolution groundwater storage changes.

[0042] Compared with the existing technology, the technical solution provided by this application has the following beneficial effects:

[0043] By integrating data from three different technologies, this method effectively overcomes key issues inherent in traditional single-source groundwater storage monitoring, such as low resolution, discontinuous data, and missing poroelastic signals. By leveraging Green's functions to correlate multi-source data with changes in total water reserves and constructing a joint inversion equation system, this method achieves deep data fusion and complementary advantages. Furthermore, by quantifying the uncertainty of each observation and introducing Laplace roughness constraints, the accuracy and reliability of the inversion results are further improved. This method enables accurate and high-resolution inversion of groundwater storage changes, providing more scientific and precise data support for water resources management, geological disaster warning, and regional hydrological research. BRIEF DESCRIPTION OF THE DRAWINGS

[0044] Figure 1 This is a flowchart of a method for inverting groundwater reserves by combining GNSS, GRACE and InSAR technologies, provided in Example 1 of the present application.

[0045] Figure 2 This is a flowchart of a method for inverting groundwater reserves by combining GNSS, GRACE and InSAR technologies, provided in Example 2 of the present application.

[0046] Figure 3 This is a flowchart of a method for removing low-correlation and irrelevant stations in GNSS vertical displacement data provided in Example 2 of the present application.

[0047] Figure 4 A flowchart of a method for constructing a joint inversion equation set through Green's function using a weighted least squares inversion framework provided in Example 2 of the present application. DETAILED DESCRIPTION

[0048] The following will combine the embodiments of the present application to clearly and completely describe the technical solutions in the embodiments of the present application. Obviously, the described embodiments are only part of the embodiments of the present application, not all of the embodiments. Based on the embodiments in the present application, all other embodiments obtained by ordinary technicians in this field without making creative work are within the scope of protection of this application.

[0049] Example 1

[0050] See also Figure 1 , is a flow chart of a method for inverting groundwater reserves using GNSS, GRACE, and InSAR technologies, as provided in Example 1 of this application. Figure 1 As shown in , the specific implementation steps of the above method include:

[0051] Step 101: Acquire GNSS vertical displacement data, GRACE equivalent water height data, and InSAR vertical deformation rate data.

[0052] Step 102: Use Green's function to correlate GNSS vertical displacement data, GRACE equivalent water height data, and InSAR vertical deformation rate data to the total water storage change, and construct a joint inversion equation group.

[0053] Step 103: quantify the uncertainties of the GNSS, InSAR, and GRACE observation trends respectively, determine the weight coefficients corresponding to the GNSS vertical displacement data, the GRACE equivalent water height data, and the InSAR vertical deformation rate data, and construct the observation data vector.

[0054] Step 104: introduce the Laplace roughness constraint, and solve the joint inversion equations by step-by-step iterative calculation using the weighted least squares method based on the observed data vector to determine the high-resolution groundwater storage change.

[0055] Compared with the prior art, the technical solution provided in Example 1 of the present application has the following beneficial effects:

[0056] This application uses a combination of GNSS, GRACE and InSAR technologies to invert groundwater reserves, effectively overcoming the limitations of traditional single technologies in groundwater reserve monitoring. By taking advantage of the high precision of GNSS vertical displacement data, the large-scale coverage of GRACE equivalent water height data and the high spatial resolution of InSAR vertical deformation rate data, the organic fusion and complementary advantages of multi-source data are achieved. By quantifying the uncertainty of each observation data and assigning corresponding weight coefficients, the observation data vector is constructed, which further improves the reliability of the inversion results. In addition, the introduction of Laplace roughness constraints and the use of weighted least squares method for step-by-step iterative solution can effectively reduce the influence of noise in the inversion process, significantly improve the inversion accuracy and spatial resolution of groundwater reserve changes, and provide more accurate, continuous and high-resolution results for the dynamic monitoring of groundwater reserves.

[0057] Example 2

[0058] See also Figure 2 , is a flow chart of a method for inverting groundwater reserves using GNSS, GRACE, and InSAR technologies, as provided in Example 2 of this application. Figure 2 As shown in , the specific implementation steps of the above method include:

[0059] Step 201: Acquire GNSS vertical displacement data, GRACE / GRACE-FO equivalent water height data, and InSAR vertical deformation rate data.

[0060] In the embodiments of the present application, GNSS data is obtained from public networks such as the China Crustal Observation Network CMONOC and IGS, solid tides, atmospheric pressure loading, and instrument drift are removed, and the vertical displacement time series is extracted. At the same time, the GNSS data is eliminated from noise interference, and stations with high signal-to-noise ratios are retained to obtain the above-mentioned GNSS vertical displacement data. The equivalent water height data of GRACE / GRACE-FO is derived from JPL Mascon RL06 and is obtained by filling in the data gaps through spatiotemporal interpolation. The InSAR vertical deformation rate data is derived from the Sentinel-1 ascending / descending orbit data. The above-mentioned InSAR vertical deformation rate data is obtained by inverting the deformation time series based on SBAS-InSAR and estimating the linear deformation rate.

[0061] Step 202: Preprocess the GNSS vertical displacement data, equivalent water height data, and vertical deformation rate data. This includes removing some stations from the GNSS vertical displacement data and projecting the vertical rate data obtained in step 201 using an affine transformation to keep it consistent with the GNSS deformation reference frame, thereby obtaining the final valid data.

[0062] The projection transformation of the vertical rate data specifically includes converting the InSAR deformation rate into the WGS84 ellipsoid high coordinate system using affine transformation and aligning it with the GNSS vertical displacement.

[0063] For removing some stations from GNSS vertical displacement data, see Figure 3 , is a flow chart of a method for removing less relevant and irrelevant stations in GNSS vertical displacement data provided in the second embodiment of the present application. Figure 3 As shown in , the above method specifically includes:

[0064] Step 2021: Use GRACE groundwater storage change data to forward model deformation, remove the impact of crustal movement caused by external far-field water mass changes from the time series of GNSS vertical displacement data, and only retain GNSS time series with a high completion rate.

[0065] Step 2022: By forward modeling the vertical land motion of the GRACE groundwater storage changes and calculating the correlation coefficients between the time series, remove the stations with correlation less than 0.3.

[0066] Specifically, the above correlation coefficient is the correlation coefficient between the GNSS vertical displacement and the GRACE mascon water storage change. In the embodiment of the present application, the retained GNSS stations are used for inversion, and the excluded GNSS stations (such as BJFS and TJN1) are used for InSAR vertical rate solution.

[0067] Step 2023: retain the stations where the aquifer poroelastic deformation is the main feature, and use them for InSAR vertical velocity determination in subsequent data processing.

[0068] Step 203: Use Green's function to correlate GNSS vertical displacement data, GRACE equivalent water height data, and InSAR vertical deformation rate data to the total water storage change, and construct a joint inversion equation group.

[0069] The present embodiment establishes a weighted least squares inversion framework to correlate GNSS vertical displacement, GRACE equivalent water height, and InSAR deformation rate with total water storage changes. A system of equations is constructed based on poroelasticity and elastic Green's functions.

[0070] Specifically, see Figure 4 , is a flow chart of a method for constructing a joint inversion equation set by using a weighted least squares inversion framework through Green's function provided in Example 2 of this application. Figure 4 As shown in , the above method specifically includes:

[0071] Step 2031: Using a layered earth model, assuming that the earth is an elastic, isotropic, spherically symmetric layered structure, calculate the Load-Love number to characterize the deformation and gravity changes caused by surface loads (such as changes in water storage). For order n≤60, the formula is as follows:

[0072] h n ′=0.603-0.011n+0.0002n 2 (Applicable to land areas)

[0073] Where h n ′ is the vertical displacement response coefficient.

[0074] Step 2032: Use Biot theory to calculate the poroelastic coefficient α and construct a poroelastic Green's function matrix. The poroelastic Green's function matrix includes two sub-matrices G poro and G elastic .

[0075] In the embodiment of the present application, the dimension of the Green function matrix G is (p+q+s)×2m, where m is the number of spatial grid cells, p is the number of GNSS stations, q is the number of InSAR effective pixels, and s is the number of GRACEmascon data points. The Green function matrix contains two sub-matrices, where the sub-matrix G poro Describes the contribution of groundwater change (ΔGWS) to InSAR deformation, with dimensions of p×m. Submatrix G elastic Describes the contribution of total water storage change (ΔTWS) to GNSS vertical displacement, with a dimension of q×m. The above submatrix G poro and submatrix Gelastic They can be expressed as:

[0076]

[0077] Where r is the distance from the center of the grid cell to the observation point (m), R is the average radius of the earth (6.371×10 6 m), ρ avg The average density of the earth (5517 kg / m 3 ). α is the Biot-Willis coefficient (dimensionless, generally 0.7 to 0.95), calculated using the Biot theoretical equation. ν is the Poisson's ratio of the aquifer, E is the Young's modulus (Pa), ρ w is the water density, ρ w =1000kg / m 3 , gravitational acceleration g = 9.81 m / s 2 , r0 is the influence radius, which is usually taken as the thickness of the aquifer.

[0078] The above poroelastic Green's function matrix can be expressed as:

[0079]

[0080] Where I is the identity matrix, which is used to relate GRACE's TWS to groundwater and other components.

[0081] Step 2033: Based on the above-mentioned poroelastic Green's function matrix, a joint observation equation is constructed, which is expressed as:

[0082]

[0083] Where ΔTWS is the total water storage change (Total Water Storage Change), which is decomposed into groundwater storage change (Groundwater Storage Change; hereinafter referred to as: ΔGWS) and other components (ΔOther), and G is the coefficient matrix constructed by the Green's function.

[0084] Step 204: Use the bootstrapping method to quantify the uncertainty of GNSS, InSAR and GRACE observation trends and generate a covariance matrix.

[0085] In the embodiment of the present application, assuming that the errors of each observation data are spatially independent (non-correlated), the above covariance matrix is ​​a diagonal matrix, and the diagonal elements are the variances of each site.

[0086] In the embodiment of the present application, a large number of "pseudo-data sets" are generated by repeatedly sampling with replacement from the original data through the bootstrap method, and then the distribution of statistics (such as mean, trend) and their uncertainty are estimated. Specifically, for GNSS observations (standard deviation = 1mm / year), it is assumed that the GNSS time series error obeys the normal distribution (Gaussian noise) with a standard deviation of 1mm / year. The original GNSS time series is resampled, and random noise with a mean of 0 and a standard deviation of 1mm / year is added to each resampled data set. Fit a linear trend, record the slope value, and repeat multiple times to obtain the distribution of the trend value, whose standard deviation is the uncertainty of the GNSS observation trend.

[0087] For InSAR observations (standard deviation = 7 mm / year), the InSAR vertical deformation rate was spatially block resampled (preserving spatial correlation). Random noise drawn from the error distribution was added to each pixel, and the regional average trend was calculated. This was repeated multiple times to obtain the uncertainty.

[0088] For GRACE observations (with a leakage error standard deviation of 2 cm), we simulated the true mass change through forward modeling and applied the model output to GRACE observation operators (such as spherical harmonic truncation and filtering) to simulate leakage effects. Comparing the simulated leakage signal with the original signal revealed a statistically significant leakage error standard deviation of 2 cm. We then extracted random values ​​from the leakage error distribution (with a standard deviation of 2 cm) and added them to the GRACE trend. This process repeatedly generated multiple leakage-affected trend estimates to calculate the uncertainty in the GRACE observation trend.

[0089] Step 205: Determine the weight coefficients of GNSS, InSAR, and GRACE using the L-curve method, and construct the observation data vector d, which is expressed as:

[0090]

[0091] The equivalent water height change ΔTWS provided by GRACE represents the total water storage change derived from the inversion; the vertical displacement ΔGNSS provided by GNSS represents processed station data; and the deformation rate ΔInSAR provided by InSAR represents processed deformation rate data. These data are organized into a column vector when constructing the d vector, which serves as the input for the inversion problem.

[0092] In the embodiment of the present application, the above-mentioned weight coefficients are estimated by the corresponding GNSS, InSAR and GRACE observation uncertainty covariance matrices, and are expressed as:

[0093]

[0094] Among them, C GNSS , C InSAR , C GRACEare covariance matrices representing the uncertainties of GNSS, InSAR, and GRACE observations, respectively, and are evaluated by the bootstrapping method described in step 204. The relative weights (a, b, c) are used to balance the strengths of GNSS, InSAR, and GRACE observations, respectively, satisfying: a + b + c = 1.

[0095] The embodiment of the present application adopts the L curve to balance the data fitting residual and the smoothness of the solution, and adjusts the weight distribution in combination with the experience of the hydrological model. Optionally, the second embodiment of the present application sets the weight coefficient of the GNSS time series to a = 0.17, the weight coefficient of the InSAR vertical deformation rate to b = 0.30, and the weight coefficient of the equivalent water height data of GRACE / GRACE-FO to c = 0.53 to balance high resolution and long period signals. Of course, the above weight coefficients can also be set according to actual conditions, and the embodiment of the present application does not make specific limitations.

[0096] Step 206: Introduce the Laplace roughness constraint and calculate the high-resolution groundwater storage change (ΔGWS) by step-by-step iterative calculation using the weighted least squares method.

[0097] In the embodiment of the present application, a Laplace smoothing matrix is ​​constructed, which is expressed as:

[0098]

[0099] According to the above Laplace smoothing matrix, the Laplace operator D is introduced as the spatial smoothing constraint to construct the Laplace roughness constraint, which is expressed as:

[0100] min||G·xd|| 2 +λ||D ij ·x|| 2

[0101] Wherein, G is the Green function matrix, λ is the smoothing factor determined by cross-validation, and in the embodiment of the present application, λ is set to 10 3 , d is the observation data vector, including the GRACE equivalent water height, GNSS vertical displacement, and InSAR vertical deformation rate, and x is the vector to be solved.

[0102] In this embodiment, the number of iterations and the convergence tolerance are set, and the optimal groundwater reserve change is iteratively solved using the least squares method. Optionally, this application sets the least squares method to 50 iterations and the convergence tolerance to 1e-6. Of course, the number of iterations and the convergence tolerance can be set based on actual conditions and are not specifically limited in this embodiment.

[0103] Compared with the prior art, the technical solution provided in Example 2 of the present application has the following beneficial effects:

[0104] By integrating data from three different technologies, the team effectively overcomes many of the challenges inherent in traditional single-source groundwater storage monitoring. Furthermore, through the fusion of multi-source data, data discontinuity is addressed, making monitoring of groundwater storage changes more continuous and stable. Furthermore, by introducing Laplace roughness constraints and combining iterative calculations with weighted least squares, the team further improves the reliability and accuracy of the inversion results, enabling a more realistic reflection of groundwater storage changes.

[0105] Example 3

[0106] Experimental verification shows that the method proposed in this application can improve spatial resolution to a 10km grid (compared to approximately 250km for the traditional GRACE method). Compared with groundwater monitoring well data, the correlation coefficient reached 0.79 (compared to 0.6-0.7 for the traditional method), and the root mean square error (RMSE) was 1.9cm / year (compared to 2.5-4cm / year for the traditional method), outperforming the traditional GRACE method.

[0107] Specifically, in the data preparation stage, Example 3 of the present application selects a major groundwater over-exploitation area in a plain area, performs regional scope and grid division, and divides the area into 200 10km×10km grid units (dimension m=200).

[0108] During the data collection and preprocessing stage, 20 public GNSS continuous monitoring data were collected, non-hydrological signals such as solid tide and atmospheric pressure loading were removed, and vertical displacement time series were extracted (time resolution: daily).

[0109] The Mascon product (RL06 version) provided by NASA JPL fills the 11-month data gap between GRACE and GRACE-FO (using spatiotemporal interpolation) and is converted to equivalent water height (EWH) with a spatial resolution of 300 km downscaled to 10 km (based on hydrological model constraints).

[0110] Based on the ESA Sentinel-1 satellite orbit raising data (2015-2024), about 600 images were processed by SBAS-InSAR to generate vertical deformation rate maps with a spatial resolution of 20m and resampled to a 10km grid.

[0111] After obtaining the above data, the data preprocessing includes: (1) GNSS site screening; (2) InSAR data projection.

[0112] GNSS site screening involved calculating the correlation coefficient between GNSS vertical displacement and GRACEmascon water storage changes, and removing sites with correlation coefficients < 0.3 (ultimately retaining 16 sites, p = 16). Poroelastic modeling was performed on the InSAR vertical velocities at the removed sites (e.g., BJFS and TJN1) for inversion.

[0113] InSAR data projection involves converting the InSAR deformation rates into the WGS84 ellipsoidal high-level coordinate system using an affine transformation, which is aligned with the GNSS vertical displacement.

[0114] Then, the joint inversion equations are constructed, which includes the following steps:

[0115] (1) Green's function calculation:

[0116] Elastic response: Calculate the Load-Love number (h n ′=0.609), and generate the elastic deformation Green's function matrix of each grid cell (dimension: 16×200).

[0117] Poroelastic response: Biot theory is used to calculate the poroelastic coefficient (α = 0.85), r0 is taken as 60m, An approximate value of -5.3 is taken at the 10 km scale to construct the poroelastic Green's function matrix (dimension: 200 × 200).

[0118] (2) Construct the joint observation equation, which is expressed as:

[0119]

[0120] According to the records in the above embodiments 1 and 2, uncertainty assessment, weight allocation, roughness constraint and inversion solution are performed, specifically including:

[0121] (1) Bootstrapping method to assess uncertainty:

[0122] GNSS: normal random noise (standard deviation = 1 mm / year) was added;

[0123] InSAR: Generate error distribution (standard deviation = 7 mm / year) based on coherence threshold (>0.3);

[0124] GRACE: Leakage error was estimated by forward modeling (standard deviation = 2 cm).

[0125] (2) L-curve method to determine weights:

[0126] Weight coefficients: a = 0.17 (GNSS), b = 0.30 (InSAR), c = 0.53 (GRACE), balancing high resolution and long period signals.

[0127] (3) Construct the Laplace smoothing matrix:

[0128]

[0129] (4) Construct the objective function:

[0130] min||G·xd|| 2 +λ||D·x|| 2

[0131] Among them, λ is 10 3 .

[0132] (5) Solution algorithm: Least squares iteration 50 times, convergence tolerance set to 1e-6.

[0133] Through comparative tests, this application obtained the following results:

[0134] (1) Groundwater storage changes (ΔGWS) show that the average annual decline rate in the North China Plain is -2.5±0.6 cm / ;

[0135] (2) The spatial resolution is increased to 10 km grid (the traditional GRACE method is about 250 km);

[0136] (3) Compared with the groundwater monitoring well data, the correlation coefficient reached 0.79 (0.6-0.7 for the traditional method), and the root mean square error (RMSE) was 1.9 cm / year (RMSE for the traditional method was 2.5-4 cm / year), which is better than the traditional GRACE method.

[0137] Although the embodiments of the present application have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and variations may be made to these embodiments without departing from the principles and spirit of the present application, and the scope of the present application is defined by the appended claims and their equivalents.

Claims

1. A method for inverting groundwater reserves by combining GNSS, GRACE and InSAR technologies, characterized in that: include: Acquire GNSS vertical displacement data, GRACE equivalent water height data, and InSAR vertical deformation rate data; The GNSS vertical displacement data, GRACE equivalent water height data, and InSAR vertical deformation rate data are correlated to the total water storage change using Green's function to construct a joint inversion equation set; Quantify the uncertainties of GNSS, InSAR, and GRACE observation trends respectively, determine the weight coefficients corresponding to the GNSS vertical displacement data, GRACE equivalent water height data, and InSAR vertical deformation rate data, and construct an observation data vector; The Laplace roughness constraint is introduced, and the joint inversion equations are solved step by step iterative calculation using the weighted least squares method through the observation data vector to determine the high-resolution groundwater storage changes.

2. The method for inverting groundwater reserves by combining GNSS, GRACE and InSAR technologies according to claim 1 is characterized in that: After obtaining GNSS vertical displacement data, GRACE equivalent water height data, and InSAR vertical deformation rate data, the method further includes: The GNSS vertical displacement data, equivalent water height data and vertical deformation rate data are preprocessed, including removing some stations of the GNSS vertical displacement data and projecting the InSAR vertical deformation rate data using an affine transformation; including converting the InSAR vertical deformation rate data into the WGS84 ellipsoidal height coordinate system using an affine transformation, aligning it with the GNSS vertical displacement, and controlling the InSAR vertical deformation rate data to remain consistent with the reference frame of the GNSS deformation.

3. The method for inverting groundwater reserves by combining GNSS, GRACE and InSAR technologies according to claim 1, characterized in that: The method of using Green's function to relate the GNSS vertical displacement data, GRACE equivalent water height data, and InSAR vertical deformation rate data to the total water storage change includes: Based on the layered earth model, the vertical displacement response coefficient is determined to characterize the deformation and gravity changes caused by surface loads; Calculate the poroelastic coefficient α and construct the poroelastic Green function matrix, which includes two sub-matrices G poro and G elastic ; Among them, the sub-matrix G poro It is used to describe the contribution of groundwater storage change ΔGWS to InSAR deformation and is expressed as: Submatrix G elastic It is used to describe the contribution of the total water storage change ΔTWS to the GNSS vertical displacement and is expressed as: Where r is the distance from the center of the grid cell to the observation point, R is the average radius of the earth, and ρ avg is the average density of the earth, α is the Biot-Willis coefficient, ν is the Poisson's ratio of the aquifer, E is the Young's modulus, ρ w is the water density, g is the acceleration of gravity, r0 is the impact radius, h n ' is the vertical displacement response coefficient; The poroelastic Green's function matrix is ​​expressed as: Where I is the identity matrix, which is used to relate the total water storage TWS of GRACE to groundwater and other components.

4. The method for inverting groundwater reserves by combining GNSS, GRACE and InSAR technologies according to claim 3, characterized in that: Constructing the joint observation equation includes: According to the poroelastic Green's function matrix, the joint observation equation is constructed and expressed as: Where ΔTWS is the change in total water storage, which is decomposed into the change in groundwater storage ΔGWS and other components ΔOther.

5. The method for inverting groundwater reserves by combining GNSS, GRACE and InSAR technologies according to claim 1, characterized in that: The uncertainties in quantified trends in GNSS, InSAR, and GRACE observations include: Resample the GNSS vertical displacement data and add random noise with a mean of 0 and a standard deviation of 1 mm / year to each resampled data point. Fit a linear trend and record the slope. Repeat this process multiple times to obtain the distribution of trend values. Determine the uncertainty of the GNSS observation trend based on the standard deviation. The InSAR vertical deformation rate data were spatially resampled. Random noise drawn from the error distribution was added to each pixel, and the regional average trend was calculated. This was repeated multiple times to determine the uncertainty of the InSAR observation trend. The leakage effect is simulated using the GRACE observation operator. The simulated leakage signal is compared with the original signal, and the standard deviation of the leakage error is calculated. Random values ​​are extracted from the leakage error distribution and added to the GRACE trend. Multiple trend estimates affected by leakage are repeatedly generated to calculate the uncertainty of the GRACE observation trend.

6. The method for inverting groundwater reserves by combining GNSS, GRACE and InSAR technologies according to claim 1, characterized in that: Determining the weight coefficients corresponding to the GNSS vertical displacement data, GRACE equivalent water height data, and InSAR vertical deformation rate data includes: Constructing corresponding covariance matrices based on the uncertainties of the GNSS, InSAR, and GRACE observation trends; Based on the smoothness of the residuals and solutions of the L-curve equilibrium data fitting, the weight coefficient distribution is determined, which is expressed as: Where C GNSS , C InSAR , C GRACE are the covariance matrices, representing the uncertainties of GNSS, InSAR, and GRACE observations, respectively. The relative weights (a, b, c) are used to balance the strengths of GNSS, InSAR, and GRACE observations, respectively.

7. The method for inverting groundwater reserves using GNSS, GRACE, and InSAR technologies according to claim 1, characterized in that: The Laplace roughness constraint is introduced, and the joint inversion equations are solved step by step iterative calculation using the weighted least squares method using the observation data vector to determine the high-resolution groundwater storage change, including: Construct the Laplace smoothing matrix, expressed as: According to the Laplace smoothing matrix, the Laplace operator D is introduced as a spatial smoothing constraint to construct the Laplace roughness constraint, which is expressed as: min||G·xd|| 2 +λ||D ij ·x|| 2 Where G is the Green's function matrix, λ is the smoothing factor determined by cross-validation, d is the observed data vector, and x is the vector to be solved, which is used to determine the high-resolution groundwater storage changes.

Images