Gravity satellite groundwater storage variable vertical signal separation method based on physical model

Abstract

The invention discloses a gravity satellite groundwater storage variable vertical signal separation method based on a physical model, and the method comprises the steps: S1, building a physical response model of ground surface deformation and gravity change and groundwater storage variables, and enabling the physical response model to comprise a ground surface deformation response model and a gravity change response model; the method specifically comprises the following steps: step S2, constructing a gravity satellite and InSAR data joint inversion framework, and establishing an observation equation set; s3, regularization constraint conditions are introduced, and a weighted least square method is adopted for solving; and S4, carrying out groundwater reserve vertical separation and verification based on a mass conservation principle. The method overcomes the defect that an existing method depends on dense well observation data, fine monitoring of groundwater storage variables is achieved through physical model constraints, and the method can be widely applied to regional groundwater resource monitoring and management.

Description

Technical Field

[0001] The present invention relates to the technical field of groundwater resource monitoring, and in particular to a method for separating vertical signals of groundwater storage variables from gravity satellites based on physical model constraints. The method achieves fine vertical separation of groundwater storage variables by fusing gravity satellite and deformation monitoring data. Background Art

[0002] Groundwater resources are a vital foundation for maintaining national economy, people's livelihoods, and sustainable development. In recent years, with the impact of socioeconomic development and climate change, groundwater overexploitation has emerged in many regions of my country. In some areas, overexploitation of both shallow groundwater and deep confined water coexists. Refined monitoring and management of groundwater resources has become an urgent task.

[0003] The development of satellite gravity measurement technology has provided new technical means for large-scale groundwater monitoring. The GRACE satellite system, launched in 2002, and the GRACE Follow-On satellite system, launched in 2018, provide continuous observations of Earth's gravity field. These gravity satellites detect gravity field variations caused by changes in groundwater reserves by measuring intersatellite distances, thus enabling large-scale groundwater storage monitoring. Furthermore, InSAR (Interferometric Synthetic Aperture Radar) technology can indirectly provide information on groundwater reserves by observing surface deformation.

[0004] Currently, groundwater stratification monitoring relies primarily on well observation data. According to a study by Cao Jie et al. (2024), the integration of gravity satellite data and groundwater well data allows for the separate monitoring of shallow and deep groundwater in the North China Plain. Dai et al. (2024) used the generalized triangular hat method and Bayesian model averaging to fuse GRACE products, combined with high-density groundwater level observation data, to separate changes in shallow and deep groundwater reserves in the North China Plain. The results showed that between 2005 and 2016, the shallow groundwater loss rate (-13.7±4.6 mm / yr) was faster than that of the deep groundwater (-6.05±3.9 mm / yr). However, these well-based stratification methods struggle to obtain reliable vertical distribution information in areas with sparse or blank well sites, leading to significant uncertainty in groundwater resource assessment.

[0005] Gravity satellite signals reflect the overall changes in the groundwater system, including shallow and deep layers, and cannot directly distinguish the contributions of groundwater at different depths. This signal aliasing seriously restricts the refined management of groundwater resources. At the same time, the hydrogeological conditions and mining patterns in different regions vary significantly, and existing monitoring methods are difficult to adapt to complex and changing application scenarios. Although existing research has provided new ideas for gravity satellite groundwater stratification monitoring, these methods are still highly dependent on water well observation data and are difficult to apply to areas lacking a dense water well observation network. In addition, existing methods lack sufficient consideration of the physical mechanisms of the groundwater system.

[0006] Therefore, this field urgently needs to develop a vertical signal separation method for groundwater storage variables from gravity satellites based on physical models. The physical mechanism of the groundwater system should be fully considered, and multi-source observation data should be integrated to achieve refined monitoring of groundwater resources. The method should also have good applicability and promotion value. Summary of the Invention

[0007] The purpose of the present invention is to provide a method for separating vertical signals of groundwater storage variables from gravity satellites based on a physical model to solve the following technical problems existing in the prior art: 1) In view of the problem that existing gravity satellite technology cannot directly distinguish between shallow and deep groundwater storage variables, the present invention establishes a physical response model based on the effective stress principle and the crustal load elastic deformation theory to realize the vertical separation of groundwater signals from gravity satellites; 2) In view of the limitation of traditional methods that over-rely on dense water well observation data, the present invention proposes a joint inversion framework that integrates gravity satellite and InSAR deformation data, which can be extended to areas with sparse or blank water well sites; 3) In view of the problem that existing methods lack constraints on spatial continuity and temporal consistency, the present invention constructs constraints based on physical models to improve the reliability and stability of the separation results; 4) In view of the actual situation that hydrogeological conditions and mining modes in different regions vary significantly, the present invention establishes a systematic separation effect evaluation system to adapt to complex and changeable application scenarios.

[0008] By solving the above-mentioned technical problems, the method provided by the present invention can realize vertical and refined monitoring of groundwater storage variables over a large area, providing important technical support for groundwater resource management. The present invention is mainly used in the field of groundwater resource monitoring and management at the national and regional scales, and can be used for dynamic monitoring of groundwater over-exploitation areas, groundwater resource evaluation, and groundwater development and utilization planning. It is particularly suitable for large-scale areas with complex hydrogeological conditions, high groundwater exploitation intensity, and insufficient coverage by traditional monitoring methods. The present invention can also be promoted and applied to related fields such as earth system science research, water resource management, and ecological environmental protection, providing important technical support for global water cycle research, climate change research, and geological disaster early warning.

[0009] To achieve the above objectives, the present invention provides a method for separating vertical signals of groundwater storage variables from gravity satellites based on a physical model, which specifically includes the following steps:

[0010] Step S1: Constructing a physical response model of surface deformation, gravity change, and groundwater storage variables, wherein the physical response model includes a surface deformation response model and a gravity change response model; specifically, the following steps are involved:

[0011] Step S101: Based on the dual response characteristics of deep groundwater, i.e., pore elastic deformation and load elastic deformation, a surface deformation response model of deep groundwater is constructed;

[0012] Step S102: Based on the fact that non-deep groundwater mainly has load elastic deformation response characteristics, a surface deformation response model of non-deep groundwater is constructed;

[0013] Step S103: Based on the elastic load theory, a gravity change response model caused by deep and non-deep groundwater is constructed;

[0014] Step S2: Construct a joint inversion framework for gravity satellite and InSAR data and establish a set of observation equations; specifically, it includes:

[0015] Step S201: Decomposing the unknown parameter vector into two parts: deep groundwater reserve change and non-deep groundwater reserve change;

[0016] Step S202: constructing a joint observation vector, including InSAR surface deformation observations and gravity satellite gravity change observations;

[0017] Step S203: establishing a joint Green's function matrix, including the surface deformation response coefficients caused by deep and non-deep groundwater and the gravity change response coefficients;

[0018] Step S3: Introduce regularization constraints and use weighted least squares method to solve the problem; specifically, it includes:

[0019] Step S301: Setting a weight matrix to balance the contribution of different observation data;

[0020] Step S302: introducing a Laplace regularization term to impose a spatial smoothness constraint;

[0021] Step S303: solving the change of deep and non-deep groundwater reserves by minimizing the objective function;

[0022] Step S4: vertical separation and verification of groundwater reserves based on the principle of mass conservation, specifically including:

[0023] Step S401: Obtaining total water storage change information of the study area from gravity satellite data;

[0024] Step S402: extracting changes in surface water and soil water storage in combination with a land surface hydrological model;

[0025] Step S403: Calculating the change of deep groundwater reserves through a joint inversion framework;

[0026] Step S404: Calculating shallow groundwater storage changes according to the water balance equation;

[0027] Step S405: Verify the separation result using the monitoring well water level data.

[0028] In one embodiment of the present invention, the surface deformation response model in step S1 is expressed as:

[0029] L InSAR =(G poro +G elastic )·dH

[0030] Where, L InSAR is the InSAR surface deformation observation vector, G poro is the poroelastic deformation Green’s function, G elastic is the load elastic deformation Green's function, dH is the groundwater storage change vector, and is expressed in the form of equivalent thickness.

[0031] In one embodiment of the present invention, the poroelastic deformation Green's function G is defined as poro , then the deformation coefficient G caused by the change of unit water thickness of the j-th model grid cell on the surface above it is poro (j) is:

[0032]

[0033] Where B j is the Skempton coefficient, α j is the Biot-Willis coefficient, K j is the bulk modulus of the aquifer, ρ w is the water density, g is the acceleration due to gravity, b 0,j is the deformation depth;

[0034] Define load elastic deformation Green's function G elastic (i,j) is:

[0035] G elastic (i,j)=G e (i, j)·ρ w ·A cell

[0036] Where G e (i, j) is the load Lofgreen function describing the surface deformation caused by unit mass load, A cellis the area of ​​the model grid cell.

[0037] In one embodiment of the present invention, the calculation formula of the gravity change response model in step S1 is:

[0038] L T =G M ·dH

[0039] Where, L T is the gravity change observation vector of the gravity satellite, dH is the groundwater storage change vector, G M is the gravity change Green's function matrix, which is used to characterize the gravity change caused by the change of unit water thickness;

[0040] Among them, the Green function moment G of the gravity change per unit equivalent water height is M The calculation formula is:

[0041] G M =A G ·ρ w ·A cell

[0042] Where A G is the Green function matrix of unit mass gravity change, whose element A G (i, j) represents the contribution of the j-th grid cell unit mass change to the gravity change at the i-th gravity satellite observation point, ρ w is the water density, A cell is the area of ​​the model grid cell.

[0043] In one embodiment of the present invention, the observation equations of step S2 are:

[0044]

[0045] Where, L InSAR For InSAR surface deformation observation, L T is the gravity change observation of the gravity satellite, H 深层 and H 非深层 are the changes in the storage of deep groundwater and non-deep groundwater, respectively; G is the joint Green function matrix, which includes the deformation and gravity response coefficients of deep and non-deep groundwater, (G poro +G elastic ) 深层 G represents the complete deformation response coefficient of deep groundwater, including two deformation mechanisms: poroelasticity and load elasticity. elastic,非深层 It represents the deformation response coefficient of non-deep groundwater, mainly load elastic deformation, G M,深层 and G M,非深层 denote the response coefficients of gravity changes caused by deep and non-deep groundwater, respectively;

[0046] Among them, deep groundwater refers to groundwater in confined aquifers, and non-deep groundwater refers to shallow groundwater, surface water, soil water, etc.

[0047] In one embodiment of the present invention, the objective function introduced by the regularization constraint in step S3 is:

[0048]

[0049] Where, is the weighted residual term, W is the weight matrix, ||DH|| 2 is the regularization term, λ is the regularization parameter, and D is the Laplace operator;

[0050] Among them, the calculation formula of the weight matrix W is:

[0051]

[0052] Among them, Σ I and Σ T are the inverse matrices of the error covariance matrices of InSAR observations and gravity satellite observations, respectively; a and b are relative weight factors, and satisfy a+b=1.

[0053] In one embodiment of the present invention, the calculation formula for the solution obtained by the weighted least squares method in step S3 is:

[0054]

[0055] Where G T is the transpose of the joint Green function matrix G, W is the weight matrix, D T is the transpose of the Laplace operator matrix D, (·) -1 represents the matrix inversion, and the regularization parameter λ is determined by L-curve analysis or cross-validation method to balance the data fitting residual and the roughness of the solution.

[0056] In one embodiment of the present invention, the water balance relationship based on which the groundwater reserves are vertically separated in step S4 is:

[0057] ΔTWS=ΔSWS+ΔGWS 浅层 +ΔGWS 深层

[0058] Where ΔTWS is the change in total water storage in the study area, ΔSWS is the change in surface water and soil water storage, and ΔGWS is the change in total water storage in the study area. 浅层 is the change in shallow groundwater storage, ΔGWS 深层 Changes in deep groundwater reserves.

[0059] In one embodiment of the present invention, the verification of step S405 is specifically performed by calculating the correlation coefficient, root mean square error, Nash-Sutcliffe efficiency coefficient and the degree of agreement between the inversion result and the measured monitoring well water level data, and analyzing the spatial distribution characteristics.

[0060] In one embodiment of the present invention, the gravity satellite data uses GRACE or GRACE-FO satellite gravity field data, the InSAR data uses SAR image pairs of ascending and descending satellite orbits, and the method further includes a step of modifying parameters in the physical response model to make the method applicable to different hydrogeological conditions including karst areas, plain areas and mountainous areas.

[0061] The proposed method for separating vertical signals from groundwater storage variables using gravity satellites, based on a physical model, offers significant advantages in vertically isolating groundwater storage variations. Compared to traditional methods, this method avoids the difficulty of vertical separation when directly using GRACE data. By leveraging a physical model-driven joint inversion framework, it achieves the complementary advantages of multi-source remote sensing data, providing a scientific basis and technical support for the refined management and sustainable utilization of regional groundwater resources. BRIEF DESCRIPTION OF THE DRAWINGS

[0062] 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.

[0063] Figure 1 It is a schematic diagram of the overall flow of the technical solution of the present invention.

[0064] Figure 2 Schematic diagram of the physical response model of surface deformation and groundwater storage volume in an embodiment of the present invention.

[0065] Figure 3 This is a structural diagram of the joint inversion framework of gravity satellite and InSAR data in an embodiment of the present invention. DETAILED DESCRIPTION

[0066] 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.

[0067] Figure 1This is a schematic diagram of the overall process of the technical solution of the present invention, showing the complete process of the vertical signal separation method of groundwater storage variables from gravity satellites based on physical models, including a data input module, a physical model construction module, a joint inversion module and a result evaluation module, as well as the data flow relationship between each module. Figure 2 This is a schematic diagram of the physical response model of surface deformation and groundwater storage in an embodiment of the present invention. It describes the physical relationship between surface deformation, effective stress and groundwater pressure changes, and demonstrates the impact mechanism of groundwater changes at different depths on surface deformation. Figure 3 This is a structural diagram of the joint inversion framework of gravity satellite and InSAR data in an embodiment of the present invention, which illustrates how gravity change observations and deformation observations jointly constrain the vertical separation process of groundwater storage variables, and demonstrates the mechanism of spatial continuity and temporal consistency constraints. Figures 1 to 3 As shown, the present invention provides a method for separating vertical signals of groundwater storage variables from gravity satellites based on a physical model, which is a joint inversion method for gravity satellite and InSAR data based on a physical model. In view of the limitations of existing groundwater monitoring technologies in vertical signal separation, the present invention aims to achieve high-precision monitoring of regional groundwater storage variables, especially deep and shallow groundwater storage variables. The core innovations of the present invention are: 1) a differentiated physical model of deep and shallow groundwater responses is constructed, revealing the differences in the deformation response mechanisms of groundwater at different depths; 2) a joint inversion framework for gravity satellite and InSAR data based on differences in physical response mechanisms is proposed, which realizes the effective separation of vertical groundwater signals by relying on only two remote sensing data sources; 3) an adaptive regularization solution algorithm is developed to overcome the ill-conditioned solution problem in traditional methods and improve the stability and reliability of the solution.

[0068] In establishing the physical model, the present invention clearly distinguishes between two main types of groundwater reserves: deep groundwater refers specifically to groundwater located in the main confined aquifer system in the region, which is surrounded by relatively impermeable layers above and below, has strong pressure-bearing capacity, and changes in water reserves cause surface deformation through both pore elastic deformation and load elastic deformation mechanisms; shallow groundwater refers to groundwater that has no stable impermeable layer on the upper part or has weak impermeable performance, weak pressure-bearing capacity, and changes in water reserves mainly cause surface deformation through the load elastic deformation mechanism.

[0069] The specific implementation process of the present invention includes the following steps:

[0070] Step S1: constructing a physical response model of surface deformation, gravity change, and groundwater storage variables, wherein the physical response model includes a surface deformation response model and a gravity change response model;

[0071] This step aims to establish a mathematical relationship between surface InSAR deformation observations, gravity satellite gravity change observations, and groundwater storage variables, providing a physical model foundation for subsequent joint inversion. The physical response model consists of two parts: the surface deformation response model and the gravity change response model, specifically including:

[0072] Step S101: Based on the dual response characteristics of deep groundwater, i.e., pore elastic deformation and load elastic deformation, a surface deformation response model of deep groundwater is constructed;

[0073] Step S102: Based on the fact that non-deep groundwater mainly has load elastic deformation response characteristics, a surface deformation response model of non-deep groundwater is constructed;

[0074] Surface deformation response model: This model comprehensively considers two groundwater-related deformation mechanisms: (1) poroelastic deformation and (2) elastic loading deformation.

[0075] In this embodiment, the surface deformation response model can be expressed as:

[0076] L InSAR =(G poro +G elastic )·dH (1)

[0077] Where, L InSAR is the InSAR surface deformation observation vector (unit of length), G poro is the pore elastic deformation Green's function, which is used to describe the pore compression / expansion deformation caused by the change of groundwater pressure. elastic is the load-elastic deformation Green’s function, which is used to describe the surface elastic deformation caused by groundwater mass loading / unloading, and dH is the groundwater storage change vector (mm).

[0078] For ease of understanding, the following will explain the detailed principles of the poroelastic deformation mechanism. First, according to the Terzaghi effective stress principle, the effective stress σ' is the difference between the total stress σ and the pore pressure P, which can be calculated as:

[0079] σ'=σ-P (2)

[0080] When the pore pressure decreases (e.g., groundwater extraction), the effective stress increases, causing the aquifer skeleton to compress and deform, manifesting as surface subsidence, where the volumetric strain dε v The following relationships exist with the effective stress change dσ', the total stress change dσ and the pore pressure change dP:

[0081]

[0082] Where dε v is the volume strain (dimensionless), db is the change in aquifer thickness (unit of length), b is the initial aquifer thickness (unit of length), dσ' is the effective stress change (unit of pressure), K is the bulk modulus of the aquifer skeleton (Bulk Modulus, pressure unit), dσ is the total stress change (unit of pressure), α is the Biot-Willis coefficient (dimensionless), and dP is the pore pressure change (unit of pressure). The Biot-Willis coefficient α characterizes the amplitude of deformation caused by pore pressure changes, and its value ranges from 0 to 1. Its value is related to the pressure degree of the aquifer. The higher the pressure degree, the closer the value is to 1. In this field, this parameter needs to be determined through field tests or laboratory tests based on the hydrogeological conditions of the actual study area, rather than simply using a fixed value. Therefore, the present invention does not limit its value.

[0083] Among them, the total stress change dσ and the pore pressure P have the following relationship:

[0084]

[0085] Where σ is the total stress (pressure unit), P is the pore pressure (pressure unit), and B is the Skempton's coefficient (dimensionless), ranging from 0 to 1. The Skempton's coefficient B needs to be determined based on the characteristics of the aquifer in the specific study area.

[0086] Substituting formula (4) into formula (3), we can get the volume strain dε v The relationship with the pore pressure change dP is:

[0087]

[0088] Formula (5) shows that the volume strain dε v It is linearly positively correlated with the pore pressure change dP, and the proportional coefficient is Among them, this coefficient is related to the poroelastic parameters B, α, and K of the aquifer, which must be determined according to the specific hydrogeological conditions of the study area. When the pore pressure decreases (dP<0), the volume strain dε v <0, indicating volume compression.

[0089] It is further assumed that the change in pore pressure dP and the change in groundwater storage dH approximately satisfy the hydrostatic pressure relationship:

[0090] dP=ρ w gdH (6)

[0091] Where, ρ wis the water density, g is the acceleration of gravity, and dH is the change in groundwater storage. Substituting the hydrostatic pressure relationship into formula (5), we can obtain the volume strain dε v The relationship with groundwater storage change dH is:

[0092]

[0093] Where, ρ w is the water density (mass / volume unit), g is the acceleration due to gravity (length / time squared unit), and dH is the change in groundwater storage (length unit). Formula (7) establishes a quantitative relationship between volume strain and groundwater storage change. When the water storage decreases (dH<0), the volume strain dε v <0, indicating volume compression.

[0094] The vertical deformation db of the ground surface can be approximated as the volume strain dε v The integral over the deformation depth b0 is:

[0095]

[0096] Where db is the vertical deformation of the ground surface (where subsidence is negative and uplift is positive, and the units are in length units), and b0 is the deformation depth, which is approximately equal to the aquifer thickness (in length units). Equation (8) establishes a linear relationship between vertical deformation of the ground surface and changes in groundwater storage.

[0097] To simplify the expression, in this embodiment, the poroelastic deformation Green's function G is defined as poro , then the deformation coefficient G caused by the unit water storage change of the j-th model grid cell on the surface above it is poro (j) is:

[0098]

[0099] Where B j is the Skempton coefficient, α j is the Biot-Willis coefficient, K j is the bulk modulus of the aquifer, ρ w is the water density, g is the acceleration due to gravity, b 0,j is the deformation depth. In this embodiment, G poro The size of the aquifer is related to the poroelastic parameters (B, α, K), water density (ρ w ), gravitational acceleration (g) and deformation depth (b0). Formula (9) shows that the poroelastic deformation Green's function G poro Calculations can be performed based on aquifer medium parameters. It should be emphasized that these parameter values ​​should be determined based on the actual hydrogeological conditions of the study area rather than directly using empirical fixed values.

[0100] Load elastic deformation response model: Changes in water mass will add or reduce loads on the Earth's surface, causing elastic deformation of the Earth's crust. Based on elastic load theory, the relationship between the vertical deformation of the surface δh and the change in mass load per unit area δm can be expressed as:

[0101] δh(r)=∫G e (rr′)·δm(r')dr′ (10)

[0102] Where r and r' are the plane coordinates of the observation point and the load position, respectively, G e is the load elastic deformation Green's function. In the spherical harmonic function domain, the vertical deformation of the surface caused by the surface load can be expressed as:

[0103]

[0104] Where a is the radius of the earth, m E is the mass of the Earth, h' l is the load Love number, ΔC lm is the spherical harmonic coefficient of mass loading change, Y lm (θ, φ) are spherical harmonics. Load Love number h' l It describes the elastic response characteristics of the earth to loads, and its value is related to the internal elastic structure of the earth (such as the PREM model).

[0105] In the discretized model grid, the vertical deformation caused by the mass load change of the jth grid cell at the i-th observation point can be expressed as:

[0106]

[0107] The relationship between groundwater quality change δm and water storage change dH is:

[0108] δm=ρ w ·A cell ·dH (13)

[0109] Substituting formula (13) into formula (12), the load elastic deformation caused by the change of water storage is obtained as follows:

[0110]

[0111] In this embodiment, the load elastic deformation Green's function G is defined as elastic (i,j) is:

[0112] G elastic (i,j)=G e (i,j)·ρ w ·A Gell(15)

[0113] Where G e (i, j) is the load Lofgreen function describing the surface deformation caused by unit mass load, A cell is the area of ​​the model grid cell.

[0114] The load elastic deformation can be simplified as:

[0115] δh elastic =G elastic ·dH (16)

[0116] Due to the difference in deformation response mechanisms between deep and shallow groundwater, that is, the deformation response mechanisms of deep and shallow groundwater are significantly different. Deep groundwater (confined aquifer) has both poroelastic deformation and load-elastic deformation responses. Therefore, the total deformation is:

[0117] δh 深层 =(G poro +G elastic )·dH 深层 (17)

[0118] Non-deep groundwater (such as phreatic aquifers and soil water) mainly exhibits load-elastic deformation response and hardly produces poroelastic deformation. Therefore:

[0119] δh 非深层 =G elastic ·dH 非深层 (18)

[0120] This different deformation response mechanism is the important physical basis for the method of this invention to achieve vertical separation of groundwater. It is important to emphasize that the determination of whether the groundwater in the study area belongs to the deep or shallow type should be based on the actual hydrogeological conditions and quantitative assessment of aquifer characteristic parameters (such as pressure bearing capacity and Biot-Willis coefficient), rather than simply dividing it according to a fixed depth.

[0121] Step S103: Based on the elastic load theory, a gravity change response model caused by deep and non-deep groundwater is constructed;

[0122] Gravity change response model: This model is based on the elastic load theory and describes the changes in the earth's gravity field caused by changes in groundwater quality. Changes in the earth's gravity field are caused by a variety of factors, including changes in groundwater quality, changes in soil moisture, glacier melting, tectonic movements, etc. The method of the present invention mainly focuses on gravity changes related to groundwater. Under the assumption of the elastic earth model, changes in the mass load on the earth's surface will cause the earth to deform and the gravity field to change. This change in the earth's gravity field caused by the surface mass load can be described by the elastic load theory. In the spherical harmonic function domain, the gravity potential U caused by the change in the earth's surface mass load dmg It can be expressed as:

[0123]

[0124] Formula (19) is the spherical harmonic expansion of the gravitational potential caused by mass load, where G is the universal gravitational constant, M E is the mass of the Earth, a is the average radius of the Earth, r is the distance from the center of the Earth, h' l and k' l are the load Love number, which is related to the elastic structure of the earth (such as the PREM model); ΔC lm is the spherical harmonic coefficient of mass loading variation, which is related to the mass loading distribution; is a spherical harmonic function; l is the order, m is the degree; θ and are the co-latitude and longitude respectively. Gravity change dg and gravity position U g The relationship is Taking the radial derivative of Equation (19) yields the spherical harmonic expression for the gravity variation dg. In practical applications, the spherical harmonic series often need to be truncated, for example, to a maximum order of 60 or 90 to match the spatial resolution of the gravity satellite data.

[0125] In order to link the gravity change model with the discretized model grid cells and gravity satellite observation points, it is necessary to discretize the above continuous spherical harmonic function expression. Discretize the study area into N grid cells, assuming that the mass change of the jth grid cell is Δm j , then the gravity change L of the i-th gravity satellite observation point (such as the mascon center) T (i) It can be expressed as the linear superposition of the mass change contributions of all grid cells:

[0126]

[0127] Where L T (i) is the gravity change observation value of the i-th gravity satellite observation point (unit of gravity acceleration), Δm j is the mass change of the jth model grid cell (mass unit), A G (i,j) is the gravity variation Green's function matrix A G The element represents the contribution of the mass change of the jth grid cell to the gravity change of the i-th gravity satellite observation point (gravitational acceleration / mass unit), and N is the total number of model grid cells. Formula (20) establishes the gravity satellite observation L T and the model grid unit mass change Δm j The linear relationship between them. Green's function matrix A G Element A G(i, j) can be calculated based on elastic load theory and spherical harmonic expansion.

[0128] Formula (20) can be simplified in matrix form as:

[0129] L T =A G ·ΔM (21)

[0130] Where, L T is the gravity change observation vector of the gravity satellite (N dimension), ΔM=[Δm1,Δm2,...,Δm N ] T is the model grid unit mass change vector (N dimension), A G is the gravity change Green function matrix (M x N dimensions, M is the number of gravity satellite observation points). Formula (21) is the gravity change response model, which takes the gravity satellite gravity change observation L T It is linked to the model grid cell mass change ΔM and is also an important component of the joint inversion framework.

[0131] Groundwater quality change ΔM G The change in groundwater storage dH and the area of ​​the model grid unit A cell and water density ρ w The relationship is:

[0132] ΔM G =ρ w ·A cell ·dH (22)

[0133] Where, ΔM G is the groundwater quality change (mass unit), ρ w is the water density (mass / volume unit), A cell is the area of ​​the model grid cell (unit of area), and dH is the change in groundwater storage in the model grid cell (unit of length). Formula (22) establishes the relationship between groundwater quality change and water storage change. It should be noted that Formula (22) assumes that the change in groundwater quality is uniformly distributed over the entire area A of the model grid cell. cell , and vertically per unit depth. A more accurate calculation of groundwater quality changes may require consideration of parameters such as aquifer thickness and specific water storage rate. To simplify the model, the present invention uses formula (22) to express the above relationship.

[0134] If it is assumed that the gravity change is mainly caused by the change of groundwater quality, then substitute formula (22) into formula (21) to obtain the gravity change response model:

[0135] L T =A G ·ρ w·A cell dH=G M ·dH (23)

[0136] Where G M =A G ·ρ w ·A cell is the gravity change Green's function matrix (gravitational acceleration / length unit), which represents the gravity change caused by the change of unit water storage. Formula (23) converts the gravity change observation L of the gravity satellite into T It is directly linked to the groundwater storage change dH, providing another observation equation for joint inversion.

[0137] In this embodiment, the calculation formula of the gravity change response model is:

[0138] L T =G M ·dH

[0139] Where, L T is the gravity change observation vector of the gravity satellite, dH is the groundwater storage change vector, G M is the gravity change Green's function matrix, which is used to characterize the gravity change caused by unit water storage change;

[0140] Among them, the Green function moment G of the gravity change per unit equivalent water height is M The calculation formula is:

[0141] G M =A G ·ρ w ·A cell

[0142] Where A G is the Green function matrix of unit mass gravity change, whose element A G (i, j) represents the contribution of the j-th grid cell unit mass change to the gravity change at the i-th gravity satellite observation point, ρ w is the water density, A cell is the area of ​​the model grid cell.

[0143] Step S2: Construct a joint inversion framework of gravity satellite and InSAR data and establish a set of observation equations;

[0144] This step aims to establish a joint inversion framework for InSAR surface deformation observations and gravity satellite gravity change observations. Based on the physical response model established in step S1, the two observation data are linked to groundwater storage changes and the corresponding observation equations are established. Based on the complete surface deformation response model (Equation (1)) and gravity change response model (Equation (23)) constructed in step S1, the joint inversion observation equations are established.

[0145] Specifically include:

[0146] To achieve vertical separation of deep and shallow groundwater storage variables, step S201: decompose the unknown parameter vector into two parts: deep groundwater storage change and non-deep groundwater storage change;

[0147]

[0148] Where H 深层 Represents the change in deep groundwater reserves, usually referring to the change in groundwater reserves in confined aquifers, H 非深层 Represents changes in shallow groundwater reserves, including changes in shallow groundwater, soil water, surface water, etc.

[0149] Step S202: constructing a joint observation vector, including InSAR surface deformation observations and gravity satellite gravity change observations;

[0150] Step S203: establishing a joint Green's function matrix, including the surface deformation response coefficients caused by deep and non-deep groundwater and the gravity change response coefficients;

[0151] In this embodiment, the observation equations for the joint inversion are established as follows:

[0152]

[0153] Where, L InSAR For InSAR surface deformation observation, L T is the gravity change observation of the gravity satellite, H 深层 and H 非深层 are the changes in the storage of deep groundwater and non-deep groundwater, respectively; G is the joint Green function matrix, which includes the deformation and gravity response coefficients of deep and non-deep groundwater, (G poro +G elastic ) 深层 G represents the complete deformation response coefficient of deep groundwater, including two deformation mechanisms: poroelasticity and load elasticity. elastic,非深层 It represents the deformation response coefficient of non-deep groundwater, mainly load elastic deformation, G M,深层 and G M,非深层 denote the response coefficients of gravity changes caused by deep and non-deep groundwater, respectively;

[0154] Among them, deep groundwater refers to groundwater in confined aquifers, and non-deep groundwater refers to shallow groundwater, surface water, soil water, etc.

[0155] The above set of observation equations reflects a key difference in the response mechanisms of deep and shallow groundwater: deep groundwater affects surface deformation through both poroelastic and load-elastic mechanisms, while shallow groundwater primarily affects surface deformation through load-elastic deformation. It is precisely this difference in physical response mechanisms that enables the present invention to achieve vertical groundwater signal separation through joint inversion. To highlight the core role of InSAR and gravity satellite data, this invention primarily uses these two data types for joint inversion.

[0156] Step S3: Introduce regularization constraints and use weighted least squares method to solve the problem; specifically, it includes:

[0157] Step S301: Setting a weight matrix to balance the contribution of different observation data;

[0158] Step S302: introducing a Laplace regularization term to impose a spatial smoothness constraint;

[0159] Step S303: solving the change of deep and non-deep groundwater reserves by minimizing the objective function;

[0160] Since the inversion problem is usually an ill-posed problem, in order to obtain a stable and reliable solution, it is necessary to introduce regularization constraints. At the same time, considering that the error characteristics of InSAR and gravity satellite data may be different, the weighted least squares method is used to solve it. In order to improve the stability and physical rationality of the inversion results, regularization constraints, such as spatial smoothness constraints, are introduced. The present invention uses Laplace regularization to impose spatial smoothness constraints, and the regularized objective function Φ(H) is the weighted sum of the weighted residual term and the regularization term. In this embodiment, the objective function introduced by the regularization constraint in step S3 is:

[0161]

[0162] Where, is the weighted residual term, and W is the weight matrix, which reflects the weight of the observation data error. It is used to balance the contribution of the two observation data of InSAR and gravity satellite, and reduce the impact of data with large errors on the inversion results. ||DH|| 2is the regularization term, λ is the regularization parameter, which controls the regularization intensity, and D is the Laplace operator, which is used to impose spatial smoothness constraints. Laplace regularization assumes that the change in groundwater reserves is smooth in space. By minimizing the sum of the squares of the second-order derivatives (curvature) of the solution, the high-frequency oscillation of the solution can be suppressed and the stability of the solution can be improved. The selection of the regularization parameter λ is crucial. If λ is too small, it will lead to insufficient regularization and unstable solution; if λ is too large, it will lead to excessive smoothing and loss of solution resolution. The present invention adopts methods such as L-curve analysis or cross-validation to select the optimal regularization parameter λ to balance the data fitting residual and the roughness of the solution. The goal of formula (26) is to obtain a spatially smooth solution (minimize the regularization term) as much as possible while ensuring the accuracy of data fitting (minimize the weighted residual term).

[0163] Among them, the weight matrix W can be expressed in detail as:

[0164]

[0165] Among them, Σ I and Σ T are the inverse matrices of the error covariance matrices for InSAR and gravity satellite observations (results of the GRACE mascon solution), respectively. The diagonal elements of these matrices are the inverse of the variance of each observation, representing the uncertainty of each observation. The larger the variance (higher the uncertainty), the smaller its inverse, and the lower its weight in the weighting.

[0166] a and b are relative weighting factors used to balance the contributions of InSAR and GRACE-FO data in the inversion. They are scalars, satisfying a + b = 1. These weights are determined through an optimization method to achieve the best data fit and solution smoothness. The weighting factors a and b are primarily determined using the L-curve method, which is achieved by running a series of inversions with different values ​​of a, b, and the smoothing factor λ. The goal is to find the optimal combination of weights so that the inversion results can both fit the observed data well (small residuals) and maintain the spatial smoothness of the solution (avoiding overfitting noise). The L-curve method is used to evaluate the residual norm and the roughness norm of the solution. The optimal solution is usually located at the "corner" of the L-curve, indicating that the best balance has been achieved between data fit and solution smoothness.

[0167] By appropriately setting the weight matrix W, the joint inversion framework can effectively utilize the complementary advantages of the two remote sensing data sources, improve the reliability and resolution of the inversion results, and achieve effective separation of deep and shallow groundwater storage variables.

[0168] In this embodiment, in step S3, by minimizing the objective function Φ(H), the optimal estimation solution of the groundwater storage change H can be obtained. That is, the least squares solution, the optimal estimate solution It can be expressed as:

[0169]

[0170] Formula (28) gives the analytical solution to the regularized least squares problem, where G T is the transpose of the joint Green function matrix G, W is the weight matrix, D T is the transpose of the Laplace operator matrix D, (·) -1 represents the matrix inversion. The regularization parameter λ is determined through L-curve analysis or cross-validation to balance the data fitting residual and the roughness of the solution. In actual calculations, when the number of model grid cells is large, directly solving Equation (28) is computationally intensive. An iterative optimization algorithm, such as the Gauss-Newton method or the conjugate gradient method, can be used to iteratively solve Equation (23) for the minimum value and obtain the optimal estimate of the water storage change, H.

[0171] Step S4: vertical separation and verification of groundwater reserves based on the principle of mass conservation, specifically including:

[0172] Step S401: Obtaining total water storage change information of the study area from gravity satellite data;

[0173] Step S402: extracting surface water and soil storage changes in combination with the land surface hydrological model;

[0174] Step S403: Calculating the change of deep groundwater reserves through a joint inversion framework;

[0175] Step S404: Calculating shallow groundwater storage changes according to the water balance equation;

[0176] Step S405: Verify the separation result using the monitoring well water level data.

[0177] After the joint inversion obtains the changes in deep and non-deep groundwater reserves, the present invention further realizes the precise vertical separation of groundwater storage variables based on the principle of conservation of mass. The theoretical basis for the vertical separation of groundwater storage variables is the regional water storage balance relationship, that is, the regional total water storage change (ΔTWS) is equal to the sum of the surface water, soil water storage change (ΔSWS) and groundwater storage change, and the groundwater storage change can be further decomposed into shallow groundwater storage change and deep groundwater storage change. In this embodiment, the water balance relationship based on which the vertical separation of groundwater reserves in step S4 is:

[0178] ΔTWS=ΔSWS+ΔGWS 浅层 +ΔGWS深层 (29)

[0179] Where ΔTWS is the change in total water storage in the study area, ΔSWS is the change in surface water and soil water storage, and ΔGWS is the change in total water storage in the study area. 浅层 is the change in shallow groundwater storage, ΔGWS 深层 Changes in deep groundwater reserves.

[0180] Based on the above water balance relationship, the change in shallow groundwater storage can be derived by the following formula:

[0181] ΔGWS 浅层 =ΔTWS-ΔSWS-ΔGWS 深层 (30)

[0182] In practice, the authors first obtain information on changes in total water reserves in the study area from GRACE / GRACE-FO gravity satellite data. Then, combined with land surface hydrological models (such as GLDAS), they extract changes in surface and soil water reserves. Simultaneously, they use the proposed joint inversion framework to directly calculate changes in deep groundwater reserves. Finally, shallow groundwater reserves are calculated using the water balance equation. This approach leverages the sensitivity of gravity satellite data to total water reserves and the ability of InSAR data to monitor groundwater-induced deformation, effectively separating changes in deep and shallow groundwater reserves.

[0183] To ensure the scientificity and reliability of the separation results, the present invention designed a systematic verification scheme, which specifically includes:

[0184] (1) Quantitative indicator verification: Calculate the Pearson correlation coefficient (r), root mean square error (RMSE), and Nash-Sutcliffe efficiency coefficient (NSE) between the inversion results and the actual well data, requiring r>0.7 and NSE>0.5 as the basic verification standards;

[0185] (2) Spatial consistency verification: Calculate the consistency index (SC) between the inversion results and the spatial distribution of multiple monitoring wells to evaluate the degree of consistency of spatial distribution characteristics, including hot spot area identification and change gradient;

[0186] (3) Time series feature verification: Wavelet analysis is used to compare the seasonal and long-term trend characteristics of the inversion results with the observed data to verify whether the inversion results can accurately capture the temporal dynamic characteristics of groundwater changes.

[0187] Deep groundwater verification is primarily conducted using time series data from deep groundwater monitoring wells within the study area. Measured water level changes are converted into reserves and systematically compared with the inversion results. Similarly, shallow groundwater verification is conducted using data from monitoring wells in the corresponding layers. The verification process focuses not only on numerical consistency but also on analyzing the degree of agreement between spatial distribution characteristics. Furthermore, time series analysis assesses the temporal evolution of the inversion results based on seasonal variation patterns and long-term trends to ensure consistency with actual observations.

[0188] In this embodiment, the gravity satellite data uses GRACE or GRACE-FO satellite gravity field data, the InSAR data uses SAR image pairs of ascending and descending satellite orbits, and the method may further include a step of modifying parameters in the physical response model to make the method applicable to different hydrogeological conditions including karst areas, plain areas and mountainous areas.

[0189] The summary of parameters, variables, etc. involved in the present invention is shown in the following table:

[0190]

[0191]

[0192]

[0193] The parameters in the table are categorized as constants, variables, observables, calculated parameters, and calculated results, depending on their role in the inversion process. Physical parameters are typically determined based on regional geological conditions, while weights and regularization parameters are determined through optimization methods to balance data fit and solution smoothness.

[0194] The following embodiment is a scheme for vertical separation of groundwater storage volume using a certain area in the North China Plain as an example. This embodiment selects a typical over-exploited area in the North China Plain as the research area. The aquifer system in this area is mainly composed of Quaternary loose sediments and has a multi-layer structure. The lithology of the shallow aquifer is mainly sand and gravel, and the deep aquifer is mainly medium and fine sand. Over-exploitation of shallow and deep groundwater coexists in the region, and the groundwater level continues to decline. It is one of the most typical groundwater over-exploitation areas in my country.

[0195] During the data preparation phase, GRACE monthly gravity field data were selected as the primary data source. Using the RL06 series of data products, the influence of non-hydrological signals such as the atmosphere and ocean was first removed. Then, Gaussian filtering was applied to the data, and finally, gravity field changes were converted into water-equivalent height. SAR data from the Sentinel-1 satellite were also acquired, and image pairs from ascending and descending orbits were selected. For deformation field processing, SBAS was used to extract surface deformation information and perform orbit error and atmospheric delay corrections. Furthermore, hydrogeological data for the study area were collected, including basic information such as aquifer structure, stratigraphic lithology, and hydrological parameters.

[0196] In constructing the physical response model, model parameters were set based on the hydrogeological conditions of the study area. Formation compressibility was determined based on the aquifer's lithologic characteristics, and the elastic storage coefficient was determined based on pumping test data. The calculation of the Green's function also considered the study area and data resolution requirements.

[0197] During the model solution process, the characteristics of groundwater extraction over many years in the study area were considered when selecting the initial iteration parameters. During the joint inversion, the weights of gravity observations and deformation observations should be appropriately set based on their respective observation accuracies to reflect the role of gravity data in constraining the total volume.

[0198] In terms of constraint setting, spatial continuity constraints consider the integrity of the hydrogeological units in the study area. Temporal consistency constraints focus on seasonal variations, and constraint weights can be determined based on the dynamic characteristics of groundwater levels at seasonal scales. Reasonable convergence criteria are set during the iterative solution process.

[0199] A person skilled in the art can implement the above technical solution to achieve vertical separation of groundwater storage variables. The specific effect may vary depending on the conditions of the study area and the parameter selection. A person skilled in the art can make equivalent substitutions or appropriate adjustments to the above implementation scheme based on specific application scenarios. As long as they do not violate the core technical concept of the present invention, they shall fall within the scope of protection of the present invention.

[0200] The physical model-based vertical signal separation method for groundwater storage variables from gravity satellites provided by the present invention has the following significant advantages over existing technologies:

[0201] 1) This paper establishes a physical response model for surface deformation and groundwater storage based on the effective stress principle, clarifying the quantitative relationship between the two and providing a reliable theoretical basis for vertical separation of gravity satellite groundwater signals. By incorporating the theory of elastic deformation under crustal load, a complete physical constraint system is established, ensuring the physical rationality of the separation results.

[0202] 2) This invention innovatively constructs a joint inversion framework for gravity satellite and InSAR data, overcoming the existing methods' reliance on densely populated water well observations. By constraining spatial continuity and temporal consistency, the reliability and stability of the separation results are improved, making the method applicable to areas with sparse or no water well observations.

[0203] 3) The technical solution proposed in this invention is highly adaptable. By introducing different parameters and correction factors, it can be applied to different hydrogeological conditions, such as karst areas, plain areas, and mountainous areas. This invention can provide important technical support for groundwater resource management at the national and regional levels, especially in dynamic monitoring of groundwater overexploitation areas, water resource assessment, and development and utilization planning.

[0204] Compared to existing technologies, this method not only solves the difficulty of separating vertical signals in groundwater monitoring using gravity satellites, but also improves the reliability of the separation results through physical model constraints. This method enables refined monitoring of groundwater resources over a large area, providing a scientific basis for groundwater resource management and sustainable utilization.

[0205] 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.

[0206] 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.

[0207] 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 method for separating vertical signals of groundwater storage variables from gravity satellites based on a physical model, characterized in that: The following steps are involved: Step S1: Constructing a physical response model of surface deformation, gravity change, and groundwater storage variables, wherein the physical response model includes a surface deformation response model and a gravity change response model; specifically, the following steps are involved: Step S101: Based on the dual response characteristics of deep groundwater, i.e., pore elastic deformation and load elastic deformation, a surface deformation response model of deep groundwater is constructed; Step S102: Based on the fact that non-deep groundwater mainly has load elastic deformation response characteristics, a surface deformation response model of non-deep groundwater is constructed; Step S103: Based on the elastic load theory, a gravity change response model caused by deep and non-deep groundwater is constructed; Step S2: Construct a joint inversion framework for gravity satellite and InSAR data and establish a set of observation equations; specifically, it includes: Step S201: Decomposing the unknown parameter vector into two parts: deep groundwater and non-deep groundwater reserve changes; Step S202: constructing a joint observation vector, including InSAR surface deformation observations and gravity satellite gravity change observations; Step S203: establishing a joint Green's function matrix, including the surface deformation response coefficients caused by deep and non-deep groundwater and the gravity change response coefficients; Step S3: Introduce regularization constraints and use weighted least squares method to solve the problem; specifically, it includes: Step S301: Setting a weight matrix to balance the contribution of different observation data; Step S302: introducing a Laplace regularization term to impose a spatial smoothness constraint; Step S303: solving the change of deep and non-deep groundwater reserves by minimizing the objective function; Step S4: vertical separation and verification of groundwater reserves based on the principle of mass conservation, specifically including: Step S401: Obtaining total water storage change information of the study area from gravity satellite data; Step S402: extracting changes in surface water and soil water storage in combination with a land surface hydrological model; Step S403: Calculating the change of deep groundwater reserves through a joint inversion framework; Step S404: Calculating shallow groundwater storage changes according to the water balance equation; Step S405: Verify the separation result using the monitoring well water level data.

2. The method according to claim 1, characterized in that The surface deformation response model in step S1 is expressed as: L InSAR =(G poro +G elastic )·dH Where, L InSAR is the InSAR surface deformation observation vector, G poro is the poroelastic deformation Green’s function, G elastic is the load elastic deformation Green's function, dH is the groundwater storage, which is expressed in the form of equivalent thickness.

3. The method according to claim 2, characterized in that Define the poroelastic deformation Green's function G poro , then the deformation coefficient G caused by the unit water storage change of the j-th model grid cell on the surface above it is poro (j) is: Where B j is the Skempton coefficient, α j is the Biot-Willis coefficient, K j is the bulk modulus of the aquifer, ρ w is the water density, g is the acceleration due to gravity, b 0,j is the deformation depth; Define load elastic deformation Green's function G elastic (i,j) is: G elastic (i,j)=G e (i,j)·ρ w ·A cell Where G e (i, j) is the load Lofgreen function describing the surface deformation caused by unit mass load, A cell is the area of ​​the model grid cell.

4. The method according to claim 1, wherein The calculation formula of the gravity change response model in step S1 is: L T =G M ·dH Where, L T is the gravity change observation vector of the gravity satellite, dH is the groundwater storage change vector, G M is the gravity change Green's function matrix, which is used to characterize the gravity change caused by unit water storage change; Among them, the Green function moment G of the gravity change per unit equivalent water height is M The calculation formula is: G M =A G ·r w ·A cell Where A G is the Green function matrix of unit mass gravity change, whose element A G (i, j) represents the contribution of the j-th grid cell unit mass change to the gravity change at the i-th gravity satellite observation point, ρ w is the water density, A cell is the area of ​​the model grid cell.

5. The method according to claim 1, wherein The observation equations of step S2 are: Where, L InSAR For InSAR surface deformation observation, L T is the gravity change observation of the gravity satellite, H 深层 and H 非深层 are the changes in the storage of deep groundwater and non-deep groundwater, respectively; G is the joint Green function matrix, which includes the deformation and gravity response coefficients of deep and non-deep groundwater, (G poro +G elastic ) 深层 G represents the complete deformation response coefficient of deep groundwater, including two deformation mechanisms: poroelasticity and load elasticity. elastic,非深层 It represents the deformation response coefficient of non-deep groundwater, mainly load elastic deformation, G M,深层 and G M,非深层 denote the response coefficients of gravity changes caused by deep and non-deep groundwater, respectively; Among them, deep groundwater refers to groundwater in confined aquifers, and non-deep groundwater refers to shallow groundwater, surface water and soil water.

6. The method according to claim 1, characterized in that The objective function introduced by the regularization constraint in step S3 is: Where, is the weighted residual term, W is the weight matrix, ||DH|| 2 is the regularization term, λ is the regularization parameter, and D is the Laplace operator; Among them, the calculation formula of the weight matrix W is: Among them, Σ I and Σ T are the inverse matrices of the error covariance matrices of InSAR observations and gravity satellite observations, respectively; a and b are relative weight factors, and satisfy a+b=1.

7. The method according to claim 1, characterized in that The calculation formula for the solution obtained by the weighted least squares method in step S3 is: Where G T is the transpose of the joint Green function matrix G, W is the weight matrix, D T is the transpose of the Laplace operator matrix D, (·) -1 represents the matrix inversion, and the regularization parameter λ is determined by L-curve analysis or cross-validation method to balance the data fitting residual and the roughness of the solution.

8. The method according to claim 1, characterized in that The water balance relationship based on which the vertical separation of groundwater reserves in step S4 is performed is: ΔTWS=ΔSWS+ΔGWS 浅层 +ΔGWS 深层 Where ΔTWS is the change in total water storage in the study area, ΔSWS is the change in surface water and soil water storage, and ΔGWS is the change in total water storage in the study area. 浅层 is the change in shallow groundwater storage, ΔGWS 深层 Changes in deep groundwater reserves.

9. The method according to claim 1, characterized in that The verification of step S405 is specifically performed by calculating the correlation coefficient, root mean square error, Nash-Sutcliffe efficiency coefficient and the degree of agreement between the inversion result and the measured monitoring well water level data, and analyzing the spatial distribution characteristics.

10. The method according to claim 1, characterized in that The gravity satellite data uses GRACE or GRACE-FO satellite gravity field data, and the InSAR data uses SAR image pairs of ascending and descending satellite orbits. The method also includes a step of modifying parameters in a physical response model to make the method applicable to different hydrogeological conditions including karst areas, plain areas, and mountainous areas.

Images