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

Abstract

The application discloses a gravity satellite underground water storage variable vertical signal separation method based on a physical model, and comprises the following steps: S1, constructing a physical response model of surface deformation and gravity change and underground water storage variable, wherein the physical response model comprises a surface deformation response model and a gravity change response model; S2, constructing a joint inversion framework of a gravity satellite and InSAR data, and establishing an observation equation set; S3, introducing a regularization constraint condition, and solving by using a weighted least square method; and S4, performing underground water storage vertical separation and verification based on the principle of mass conservation. The application overcomes the dependence of the prior art method on dense water well observation data, realizes fine monitoring of the underground water storage variable through the physical model constraint, and can be widely applied to regional underground water resource monitoring and management.

Description

Technical Field

[0001] This invention relates to the field of groundwater resource monitoring technology, and more specifically, to a method for vertical signal separation of groundwater storage variables based on physical model constraints using gravity satellite data. This method achieves precise vertical separation of groundwater storage variables by fusing gravity satellite and deformation monitoring data. Background Technology

[0002] Groundwater resources are a vital foundation for national welfare and sustainable development. In recent years, with socio-economic development and the impact of climate change, many regions in my country have experienced groundwater over-extraction. In some areas, both shallow groundwater and deep confined aquifers are over-extracted. Therefore, refined monitoring and management of groundwater resources has become an urgent priority.

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

[0004] Currently, groundwater stratification monitoring mainly relies on well observation data. According to Cao et al. (2024), fusing gravity satellite data and groundwater well data can achieve separate monitoring of shallow and deep groundwater in the North China Plain. Dai et al. (2024) used the generalized triangular cap method and Bayesian model averaging method to fuse GRACE products, combined with high-density groundwater level observation data, to achieve separation of shallow and deep groundwater storage changes in the North China Plain. The results showed that during 2005-2016, the shallow groundwater depletion rate (-13.7±4.6 mm / yr) was faster than that of deep groundwater (-6.05±3.9 mm / yr). However, these well-data-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 groundwater systems, including shallow and deep layers, but cannot directly distinguish the contributions of groundwater at different depths. This signal aliasing severely restricts the refined management of groundwater resources. Furthermore, significant differences in hydrogeological conditions and extraction patterns across different regions make existing monitoring methods ill-suited to complex and ever-changing application scenarios. Although current research provides new insights into gravity satellite-based groundwater stratification monitoring, these methods still heavily rely on well observation data, making them difficult to apply to areas lacking dense well observation networks. In addition, existing methods lack sufficient consideration of the physical mechanisms of groundwater systems.

[0006] Therefore, there is an urgent need in this field to develop a method for separating vertical signals of groundwater storage variables based on physical models from gravity satellites. This method should fully consider the physical mechanisms of the groundwater system, integrate multi-source observation data, and achieve refined monitoring of groundwater resources. It should also have good applicability and promotional value. Summary of the Invention

[0007] The purpose of this invention is to provide a physical model-based method for vertical signal separation of groundwater storage variables from gravity satellites, addressing the following technical problems in existing technologies: 1) To address the inability of existing gravity satellite technology to directly distinguish between shallow and deep groundwater storage variables, this invention establishes a physical response model based on the effective stress principle and the theory of elastic deformation of crustal load, achieving vertical separation of groundwater signals from gravity satellites; 2) To address the limitation of traditional methods relying excessively on dense well observation data, this invention proposes a joint inversion framework integrating gravity satellite and InSAR deformation data, which can be extended to areas with sparse or blank well sites; 3) To address the lack of constraints on spatial continuity and temporal consistency in existing methods, this invention constructs constraint conditions based on a physical model, improving the reliability and stability of the separation results; 4) Considering the significant differences in hydrogeological conditions and mining modes in different regions, this invention establishes a systematic evaluation system for the separation effect, adapting to complex and varied application scenarios.

[0008] By addressing the aforementioned technical problems, the method provided by this invention enables vertically refined monitoring of groundwater storage variables over a large area, providing crucial technical support for groundwater resource management. This invention is primarily applied to groundwater resource monitoring and management at national and regional scales, and can be used for dynamic monitoring of groundwater over-extraction areas, groundwater resource assessment, and groundwater development and utilization planning. It is particularly suitable for areas with complex hydrogeological conditions, high groundwater extraction intensity, and insufficient coverage by traditional monitoring methods. This invention can also be extended 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, this invention provides a method for separating vertical signals of groundwater storage variables from gravity satellites based on a physical model, specifically including the following steps:

[0010] Step S1: Construct a physical response model of surface deformation, gravity change, and groundwater storage change. This physical response model includes a surface deformation response model and a gravity change response model; specifically:

[0011] Step S101: Based on the dual response characteristics of deep groundwater, which simultaneously possesses pore elastic deformation and load elastic deformation, a surface deformation response model for deep groundwater is constructed.

[0012] Step S102: Based on the load elastic deformation response characteristics of non-deep groundwater, construct a surface deformation response model for non-deep groundwater;

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

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

[0015] Step S201: Decompose the unknown parameter vector into two parts: changes in deep groundwater storage and changes in non-deep groundwater storage;

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

[0017] Step S203: Establish a joint Green's function matrix, which includes 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 solve using the weighted least squares method; specifically including:

[0019] Step S301: Set the weight matrix to balance the contributions of different observation data;

[0020] Step S302: Introduce a Laplace regularization term to apply a spatial smoothing constraint;

[0021] Step S303: Solve for the changes in deep and non-deep groundwater storage by minimizing the objective function;

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

[0023] Step S401: Obtain information on changes in total water storage in the study area from gravity satellite data;

[0024] Step S402: Extract changes in surface water and soil water storage by combining the land surface hydrological model;

[0025] Step S403: Calculate the changes in deep groundwater storage using the joint inversion framework;

[0026] Step S404: Calculate the change in shallow groundwater storage based on the water balance equation;

[0027] Step S405: Verify the separation results using 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] In the formula, L InSAR G is the InSAR surface deformation observation vector. poro G is the Green's function for pore elastic deformation. elastic dH is the Green's function for load elastic deformation, and dH is the vector of groundwater storage change, expressed in the form of equivalent thickness.

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

[0032]

[0033] In the formula, B j Let α be the Skempton coefficient. j K represents the Biot-Willis coefficient. j ρ is the bulk modulus of the aquifer. w Let g be the density of water, g be the acceleration due to gravity, and b be the acceleration due to gravity. 0,j Deformation depth;

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

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

[0036] In the formula, G e (i,j) is the Lovegreen function describing the surface deformation caused by a unit mass load, A cellThis represents the area of ​​the model's mesh cells.

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

[0038] L T =G M ·dH

[0039] In the formula, L T dH is the gravity change observation vector from the gravity satellite, dH is the groundwater storage change vector, and G is the groundwater storage change vector. M This is the Green's function matrix for gravity variation, used to characterize the change in gravity caused by a unit change in water thickness;

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

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

[0042] In the formula, A G Let A be the Green's function matrix representing the change in gravitational force of a single unit mass, and let its element A be the Green's function matrix. G (i,j) represents the contribution of the unit mass change of the j-th grid cell to the gravity change at the i-th gravity satellite observation point, ρ w Let A be the density of water. cell This represents the area of ​​the model's mesh cells.

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

[0044]

[0045] In the formula, L InSAR For InSAR surface deformation observation, L T For gravity satellite observation of gravity changes, H 深层 and H 非深层 Let G represent the changes in deep and shallow groundwater reserves, respectively, and G be the joint Green's function matrix, containing the deformation and gravity response coefficients of deep and shallow groundwater. poro +G elastic ) 深层 G represents the complete deformation response coefficient of deep groundwater, encompassing both pore elastic and load elastic deformation mechanisms. elastic,非深层 G represents the deformation response coefficient of non-deep groundwater, mainly representing elastic deformation under load. M,深层 and G M,非深层 These represent the response coefficients to gravity changes caused by deep and shallow groundwater, respectively.

[0046] Deep groundwater refers to groundwater within confined aquifers, while non-deep groundwater includes shallow groundwater, surface water, and soil water.

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

[0048]

[0049] In the formula, For the weighted residuals, W is the weight matrix, ||DH|| 2 λ is the regularization term, λ is the regularization parameter, and D is the Laplace operator;

[0050] The formula for calculating the weight matrix W is as follows:

[0051]

[0052] Where, Σ I and Σ T These are the inverse matrices of the error covariance matrices of the InSAR observations and the gravity satellite observations, respectively; a and b are relative weighting 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 as follows:

[0054]

[0055] In the formula, G T Let G be the transpose of the joint Green's function matrix, W be the weight matrix, and D be the weight matrix. T Let D be the transpose of the Laplace operator matrix, (·) -1 This represents matrix inversion. The regularization parameter λ is determined through L-curve analysis or cross-validation to balance the roughness of the data fitting residuals and solutions.

[0056] In one embodiment of the present invention, the water balance relationship on which the vertical separation of groundwater storage in step S4 is based is:

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

[0058] In the formula, ΔTWS represents the change in total water storage in the study area, ΔSWS represents the change in surface water and soil water storage, and ΔGWS represents the change in total water storage in the study area. 浅层 For changes in shallow groundwater storage, ΔGWS 深层 This represents changes in deep groundwater reserves.

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

[0060] In one embodiment of the present invention, the gravity satellite data adopts GRACE or GRACE-FO satellite gravity field data, the InSAR data adopts SAR image pairs of ascending and descending orbit satellites, and the method further includes the step of modifying the 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 physical model-based vertical signal separation method for groundwater storage variables from gravity satellites provided by this invention has significant advantages in vertically separating changes in groundwater storage. Compared with traditional methods, this invention avoids the problem of difficulty in vertical separation when directly using GRACE data. Through a physical model-driven joint inversion framework, it achieves 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. Attached Figure Description

[0062] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

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

[0064] Figure 2 This is a schematic diagram of the physical response model of surface deformation and groundwater storage variables in an embodiment of the present invention.

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

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

[0067] Figure 1This is a schematic diagram of the overall process of the technical solution of the present invention, which shows the complete process of the vertical signal separation method of gravity satellite groundwater storage variable based on physical model, including data input module, physical model construction module, joint inversion module and 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 influence mechanism of groundwater changes at different depths on surface deformation. Figure 3 This diagram illustrates the framework for the joint inversion of gravity satellite and InSAR data in this invention, demonstrating how gravity change observations and deformation observations jointly constrain the vertical separation process of groundwater storage variables, and showcasing the mechanism of spatial continuity and temporal consistency constraints. Figures 1-3 As shown, this invention provides a physical model-based method for separating vertical signals of groundwater storage variables from gravity satellite data. This method is a joint inversion method using gravity satellite and InSAR data based on a physical model. Addressing the limitations of existing groundwater monitoring technologies in vertical signal separation, this invention aims to achieve high-precision monitoring of regional groundwater storage variables, particularly those at deep and shallow depths. The core innovations of this invention are: 1) Constructing differentiated physical models of deep and shallow groundwater responses, revealing the differences in deformation response mechanisms at different depths; 2) Proposing a joint inversion framework for gravity satellite and InSAR data based on differences in physical response mechanisms, enabling effective separation of vertical groundwater signals using only two remote sensing data sources; 3) Developing an adaptive regularized solution algorithm, overcoming the ill-conditioned solution problem in traditional methods and improving the stability and reliability of the solution.

[0068] In establishing the physical model, this invention clearly distinguishes two main types of groundwater reserves: deep groundwater specifically refers to groundwater located within the main confined aquifer system in the region, surrounded by relatively impermeable layers above and below, with strong confinement, and changes in water reserves cause surface deformation through both pore elastic deformation and load elastic deformation mechanisms; shallow groundwater refers to groundwater without a stable impermeable layer above or with weak impermeability, with weak confinement, and changes in water reserves mainly cause surface deformation through load elastic deformation mechanisms.

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

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

[0071] This step aims to establish the 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. This physical response model consists of two parts: a surface deformation response model and a gravity change response model, specifically including:

[0072] Step S101: Based on the dual response characteristics of deep groundwater, which simultaneously possesses pore elastic deformation and load elastic deformation, a surface deformation response model for deep groundwater is constructed.

[0073] Step S102: Based on the load elastic deformation response characteristics of non-deep groundwater, construct a surface deformation response model for non-deep groundwater;

[0074] Surface deformation response model: This model takes into account two deformation mechanisms related to groundwater: (1) pore elastic 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] In the formula, L InSAR G is the InSAR surface deformation observation vector (unit of length). poro G is the Green's function for pore elastic deformation, used to describe the pore compression / expansion deformation caused by changes in groundwater pressure. elastic dH is the Green's function for elastic deformation under load, used to describe the elastic deformation of the surface caused by groundwater mass loading / unloading, and dH is the vector of groundwater storage change (mm).

[0078] To facilitate understanding, the detailed principle of the pore elastic deformation mechanism will be explained below. First, according to Terzaghi's effective stress principle, the effective stress σ' is the difference between the total stress σ and the pore pressure P, specifically calculated as follows:

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

[0080] When pore pressure decreases (e.g., during groundwater extraction), effective stress increases, leading to compressive deformation of the aquifer framework, manifested as surface subsidence, where the volumetric strain dε v The following relationships exist between the effective stress change dσ', the total stress change dσ, and the pore pressure change dP:

[0081]

[0082] In the formula, dε v Let d be the volumetric strain (dimensionless), db be the aquifer thickness variation (unit of length), b be the initial aquifer thickness (unit of length), dσ' be the effective stress variation (unit of pressure), K be the bulk modulus of the aquifer skeleton (unit of pressure), dσ be the total stress variation (unit of pressure), α be the Biot-Willis coefficient (dimensionless), and dP be the pore pressure variation (unit of pressure). The Biot-Willis coefficient α characterizes the magnitude of deformation caused by pore pressure variation, ranging from 0 to 1. Its value is related to the degree of pressure bearing of the aquifer; the higher the degree of pressure bearing, 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 actual hydrogeological conditions of the study area, rather than simply using a fixed value. Therefore, this invention does not limit its value.

[0083] The total stress change dσ is related to the pore pressure P as follows:

[0084]

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

[0086] Substituting equation (4) into equation (3), we obtain the volumetric strain dε. v The relationship between dP and the pore pressure change is as follows:

[0087]

[0088] Equation (5) shows that the volumetric strain dε v It is linearly positively correlated with the change in pore pressure dP, with a proportionality coefficient of This coefficient is related to the pore elastic parameters B, α, and K of the aquifer, all of which need to be determined based on the specific hydrogeological conditions of the study area. When the pore pressure decreases (dP < 0), the volumetric strain dε... v <0 indicates volume compression.

[0089] Further assuming that the pore pressure change dP and the groundwater storage change dH approximately satisfy a hydrostatic pressure relationship:

[0090] dP=ρ w gdH (6)

[0091] In the formula, ρ wLet ρ be the density of water, g be the acceleration due to gravity, and dH be the change in groundwater storage. Substituting the hydrostatic pressure relationship into formula (5), we obtain the volumetric strain dε. v The relationship between dH and the change in groundwater storage is as follows:

[0092]

[0093] In the formula, ρ w Let dε be the density of water (in units of mass / volume), g be the acceleration due to gravity (in units of length / square of time), and dH be the change in groundwater storage (in units of length). Formula (7) establishes a quantitative relationship between volumetric strain and change in groundwater storage. When the water storage decreases (dH<0), the volumetric strain dε v <0 indicates volume compression.

[0094] The vertical deformation of the Earth's surface, db, can be approximated as the volumetric strain, dε. v The integral over the deformation depth b0 is as follows:

[0095]

[0096] In the formula, db is the vertical deformation of the land surface (where subsidence is negative and uplift is positive, and the unit is length), and b0 is the deformation depth, which is approximately equal to the aquifer thickness (in length). Formula (8) establishes a linear relationship between the vertical deformation of the land surface and the change in groundwater storage.

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

[0098]

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

[0100] Load-elastic deformation response model: Changes in water mass impose or alleviate loads on the Earth's surface, leading to elastic deformation of the Earth's crust. Based on elastic load theory, the relationship between the vertical surface deformation δ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] In the formula, r and r' are the planar coordinates of the observation point and the load location, respectively, and G e It is the Green's function of elastic deformation under load. In the spherical harmonic domain, the vertical deformation of the land surface caused by surface load can be expressed as:

[0103]

[0104] In the formula, a is the Earth's radius, m E For Earth's mass, h' l For the load Love number, ΔC lm Y is the spherical harmonic coefficient for varying mass load. lm (θ,φ) are spherical harmonic functions. The load Love number h' l It describes the Earth's elastic response characteristics to loads, and its values ​​are related to the Earth's internal elastic structure (such as the PREM model).

[0105] In the discretized model mesh, the vertical deformation at the i-th observation point caused by the change in mass load of the j-th mesh element can be expressed as:

[0106]

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

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

[0109] Substituting formula (13) into formula (12), we obtain the load elastic deformation caused by changes in water storage:

[0110]

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

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

[0113] In the formula, G e (i,j) is the Lovegreen function describing the surface deformation caused by a unit mass load, A cell This represents the area of ​​the model's mesh cells.

[0114] The elastic deformation under load 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, a significant difference exists between the deformation response mechanisms of deep and shallow groundwater—deep groundwater (confined aquifers) exhibits both pore elastic deformation and load elastic deformation responses simultaneously. Therefore, the total deformation is:

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

[0118] Non-deep groundwater (unconfined aquifers, soil water, etc.) mainly exhibits a load-elastic deformation response, producing almost no pore elastic 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 should be emphasized that the determination of whether the groundwater in the study area is deep or shallow should be based on the actual hydrogeological conditions and quantitatively assessed through aquifer characteristic parameters (such as confined aquifer, Biot-Willis coefficient, etc.), rather than simply classifying it according to a fixed depth.

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

[0122] Gravity Variation Response Model: This model, based on elastic loading theory, describes the changes in the Earth's gravity field caused by variations in groundwater mass. These changes are caused by various factors, including variations in groundwater mass, soil moisture, glacial melting, and tectonic activity. This invention primarily focuses on gravity variations related to groundwater. Under the assumption of an elastic Earth model, changes in the Earth's surface mass load cause Earth deformation and changes in the gravity field. These changes in the Earth's gravity field caused by surface mass load can be described using elastic loading theory. In the spherical harmonic domain, the gravitational potential U caused by the change in Earth's surface mass load dm is...g It can be represented as:

[0123]

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

[0125] To connect the gravity variation model with the discretized model grid cells and gravity satellite observation points, the aforementioned continuous spherical harmonic function expression needs to be discretized. The study area is discretized into N grid cells, and the mass change of the j-th grid cell is assumed to be Δm. j Then the gravity change L at the i-th gravity satellite observation point (such as the mascon center) is... T (i) can be represented as a linear superposition of the contributions of all mesh element mass variations:

[0126]

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

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

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

[0130] In the formula, L T Let ΔM be the gravity change observation vector (N-dimensional) of a gravity satellite, where ΔM = [Δm1, Δm2, ..., Δm N ] T Let A be the mass variation vector (N-dimensional) of the model mesh elements. G Let L be the Green's function matrix for gravity variation (M x N dimensions, where M is the number of gravity satellite observation points). Equation (21) is the gravity variation response model, which uses the gravity satellite gravity variation observation L... T It is also an important component of the joint inversion framework and is related to the change in the mass of the model grid cells, ΔM.

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

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

[0133] In the formula, ΔM G ρ represents the change in groundwater mass (unit: mass). w A is the density of water (in units of mass / volume). cell Let A be the area of ​​the model grid cell (in area units), and dH be the change in groundwater storage in the model grid cell (in length units). Formula (22) establishes the relationship between changes in groundwater quality and changes in water storage. It should be noted that formula (22) assumes that the change in groundwater quality is uniformly distributed across the entire area A of the model grid cell. cell Above, and in the vertical direction, is the unit depth. More accurate calculations of groundwater quality changes may require consideration of parameters such as aquifer thickness and specific storage ratio. To simplify the model, this invention uses formula (22) to represent the above relationship.

[0134] If we assume that the gravity change is mainly caused by the change in groundwater mass, then substituting formula (22) into formula (21) yields the gravity change response model:

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

[0136] In the formula, G M =A G ·ρ w ·A cell The Green's function matrix for gravity variation (gravitational acceleration / length unit) characterizes the change in gravity caused by a change in unit water storage. Formula (23) uses gravity satellite gravity variation observations L T It is directly linked to changes in groundwater storage dH, providing another observation equation for joint inversion.

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

[0138] L T =G M ·dH

[0139] In the formula, L T dH is the gravity change observation vector from the gravity satellite, dH is the groundwater storage change vector, and G is the groundwater storage change vector. M This is the Green's function matrix for gravity variation, used to characterize the change in gravity caused by a change in unit water storage.

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

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

[0142] In the formula, A G Let A be the Green's function matrix representing the change in gravitational force of a single unit mass, and let its element A be the Green's function matrix. G (i,j) represents the contribution of the unit mass change of the j-th grid cell to the gravity change at the i-th gravity satellite observation point, ρ w Let A be the density of water. cell This represents the area of ​​the model's mesh cells.

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

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

[0145] Specifically, it includes:

[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] In the formula, H 深层 Represents changes in deep groundwater reserves, typically referring to changes in groundwater within confined aquifers, H 非深层 It represents changes in shallow groundwater reserves, including changes in shallow groundwater, soil water, and surface water.

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

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

[0151] In this embodiment, the established joint inversion observation equation set is as follows:

[0152]

[0153] In the formula, L InSAR For InSAR surface deformation observation, L T For gravity satellite observation of gravity changes, H 深层 and H 非深层 Let G represent the changes in deep and shallow groundwater reserves, respectively, and G be the joint Green's function matrix, containing the deformation and gravity response coefficients of deep and shallow groundwater. poro +G elastic ) 深层 G represents the complete deformation response coefficient of deep groundwater, encompassing both pore elastic and load elastic deformation mechanisms. elastic,非深层 G represents the deformation response coefficient of non-deep groundwater, mainly representing elastic deformation under load. M,深层 and G M,非深层 These represent the response coefficients to gravity changes caused by deep and shallow groundwater, respectively.

[0154] Deep groundwater refers to groundwater within confined aquifers, while non-deep groundwater includes shallow groundwater, surface water, and soil water.

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

[0156] Step S3: Introduce regularization constraints and solve using the weighted least squares method; specifically including:

[0157] Step S301: Set the weight matrix to balance the contributions of different observation data;

[0158] Step S302: Introduce a Laplace regularization term to apply a spatial smoothing constraint;

[0159] Step S303: Solve for the changes in deep and non-deep groundwater storage by minimizing the objective function;

[0160] Since inversion problems are typically ill-conditioned, regularization constraints are needed to obtain stable and reliable solutions. Furthermore, considering the potentially different error characteristics of InSAR and gravity satellite data, a weighted least squares method is used for solving the problem. To improve the stability and physical plausibility of the inversion results, regularization constraints, such as spatial smoothing constraints, are introduced. This invention employs Laplacian regularization to apply spatial smoothing constraints, and the regularized objective function Φ(H) is the weighted sum of the weighted residual terms and the regularization term. In this embodiment, the objective function introduced by the regularization constraint in step S3 is:

[0161]

[0162] In the formula, For the weighted residual term, W is the weight matrix, reflecting the weights of observation data errors. It is used to balance the contributions of InSAR and gravity satellite observation data and reduce the impact of data with larger errors on the inversion results. ||DH|| 2Let λ be the regularization term, λ be the regularization parameter controlling the regularization strength, and D be the Laplace operator used to apply spatial smoothness constraints. Laplace regularization assumes that the change in groundwater storage is spatially smooth. By minimizing the sum of squares of the second derivatives (curvature) of the solution, high-frequency oscillations of the solution can be suppressed, and the stability of the solution can be improved. The choice 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 over-smoothing and loss of solution resolution. This invention uses methods such as L-curve analysis or cross-validation to select the optimal regularization parameter λ in order to balance the data fitting residuals 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 (minimizing the weighted residual term).

[0163] The weight matrix W can be expressed in detail as follows:

[0164]

[0165] Where, Σ I and Σ T These are the inverse matrices of the error covariance matrices for InSAR and gravity satellite observations (results from GRACE mascon calculations), respectively. The diagonal elements of these matrices are the reciprocals of the variance of each observation, representing the uncertainty of each observation. The larger the variance (the higher the uncertainty), the smaller its reciprocal, 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 optimization methods to achieve optimal data fit and solution smoothness. The weighting factors a and b are primarily determined using the L-curve method, achieved by running a series of inversions with different values ​​of a, b, and the smoothing factor λ. The goal is to find the optimal weight combination that allows the inversion results to 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 relationship between the residual norm and the roughness norm of the solution. The optimal solution is usually located at the "corner" of the L-curve, indicating an optimal balance 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, the optimal estimate of the groundwater storage change H can be obtained by minimizing the objective function Φ(H) in step S3. That is, the least squares solution, the optimal estimated solution. It can be represented as:

[0169]

[0170] Formula (28) gives the analytical solution to the regularized least squares problem, where G T Let G be the transpose of the joint Green's function matrix, W be the weight matrix, and D be the weight matrix. T Let D be the transpose of the Laplace operator matrix, (·) -1 The matrix inversion is represented by the regularization parameter λ, which is determined by L-curve analysis or cross-validation to balance the roughness of the data fitting residuals and the solution. In actual calculations, when the number of model grid cells is large, the computational cost of directly solving formula (28) is large. Iterative optimization algorithms, such as the Gauss-Newton method or the conjugate gradient method, can be used to iteratively solve the minimum value of formula (23) to obtain the optimal estimated solution of the water storage change H.

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

[0172] Step S401: Obtain information on changes in total water storage in the study area from gravity satellite data;

[0173] Step S402: Extract changes in surface water and soil storage by combining the land surface hydrological model;

[0174] Step S403: Calculate the changes in deep groundwater storage using the joint inversion framework;

[0175] Step S404: Calculate the change in shallow groundwater storage based on the water balance equation;

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

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

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

[0179] In the formula, ΔTWS represents the change in total water storage in the study area, ΔSWS represents the change in surface water and soil water storage, and ΔGWS represents the change in total water storage in the study area. 浅层 For changes in shallow groundwater storage, ΔGWS 深层 This represents changes in deep groundwater reserves.

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

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

[0182] In practice, the method first obtains information on total water storage changes in the study area from GRACE / GRACE-FO gravity satellite data. Then, it extracts surface water and soil water storage changes using a land surface hydrological model (such as GLDAS). Simultaneously, it directly calculates deep groundwater storage changes using the joint inversion framework of this invention. Finally, it calculates shallow groundwater storage changes based on the water balance equation. This method fully utilizes the sensitivity of gravity satellite data to total water storage and the ability of InSAR data to monitor groundwater-induced deformation, achieving effective separation of deep and shallow groundwater storage changes.

[0183] To ensure the scientific validity and reliability of the separation results, this invention designs a systematic verification scheme, specifically including:

[0184] (1) Quantitative index 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. The basic verification criteria are r>0.7 and NSE>0.5.

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

[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 observation data to verify whether the inversion results can accurately capture the time dynamic characteristics of groundwater changes.

[0187] Deep groundwater validation primarily utilizes time-series data of water levels from monitoring wells within the study area. Measured water level changes are converted into storage changes and then systematically compared with the inversion results. Similarly, shallow groundwater validation employs data from monitoring wells at the corresponding strata. The validation process focuses not only on numerical consistency but also on analyzing the degree of agreement on spatial distribution characteristics. Furthermore, time-series analysis assesses whether the temporal evolution characteristics of the inversion results match actual observations from two dimensions: seasonal variation patterns and long-term trends.

[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 orbit satellites, and the method may further include the step of modifying the parameters in the physical response model to make the method applicable to different hydrogeological conditions, including karst areas, plains and mountainous areas.

[0189] Please refer to the table below for a summary description of the parameters, variables, etc. involved in this invention:

[0190]

[0191]

[0192]

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

[0194] The following example demonstrates a scheme for vertical separation of groundwater storage variables, using a region in the North China Plain as an example. This example selects a typical over-extraction area in the North China Plain as the study area. The aquifer system in this area is mainly composed of Quaternary loose sediments, exhibiting a multi-layered structure. The shallow aquifer lithology is dominated by sand and gravel, while the deep aquifer is dominated by medium to fine sand. Over-extraction of both shallow and deep groundwater coexists in this area, and the groundwater level continues to decline, making it one of the most typical groundwater over-extraction areas in my country.

[0195] In the data preparation phase, GRACE monthly gravity field data was selected as the primary data source. Using the RL06 series data products, the influence of non-hydrological signals such as atmospheric and oceanic signals was first removed, followed by Gaussian filtering. Finally, gravity field variations were converted into water equivalent altitude. Simultaneously, SAR data from the Sentinel-1 satellite was acquired, selecting image pairs from both ascending and descending orbits. In deformation field processing, the SBAS method was used to extract surface deformation information, and orbital error and atmospheric delay corrections were applied. Additionally, 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 are set according to the hydrogeological conditions of the study area. The formation compressibility coefficient is determined by layering based on the lithological characteristics of the aquifer, and the elastic storage coefficient is determined based on pumping test data. The calculation of the Green's function needs to consider the study area and data resolution requirements.

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

[0198] Regarding the setting of constraints, spatial continuity constraints consider the integrity of hydrogeological units in the study area. Temporal consistency constraints focus on seasonal variation patterns, and the weights of constraint terms 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] Those skilled in the art can achieve vertical separation of groundwater storage variables by implementing the above technical solution. Specific results may vary depending on the conditions of the study area and the selection of parameters. Those skilled in the art can make equivalent substitutions or appropriate adjustments to the above implementation scheme according to specific application scenarios; as long as these do not violate the core technical concept of this invention, they should all fall within the protection scope of this invention.

[0200] The vertical signal separation method for groundwater storage variables based on a physical model provided by this invention has the following significant advantages compared with existing technologies:

[0201] 1) This invention 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 the vertical separation of groundwater signals from gravity satellites. By introducing 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 dependence of existing methods on dense water well observation data. By constraining spatial continuity and temporal consistency, the reliability and stability of the separation results are improved, making this method applicable to sparse or blank areas of water well sites.

[0203] 3) The technical solution proposed in this invention has strong adaptability. By introducing different parameters and correction terms, it can be applied to different hydrogeological conditions such as karst areas, plains, and mountainous areas. This invention can provide important technical support for groundwater resource management at the national and regional scales, and has significant application value, especially in the dynamic monitoring, water resource assessment, and development and utilization planning of groundwater over-extraction areas.

[0204] Compared to existing technologies, this invention not only solves the problem of separating vertical signals in gravity satellite groundwater monitoring, 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 water resource utilization.

[0205] Those skilled in the art will understand that the accompanying drawings are merely schematic diagrams of one embodiment, and the modules or processes shown in the drawings are not necessarily essential for implementing the present invention.

[0206] Those skilled in the art will understand that the modules in the apparatus of the embodiments can be distributed in the apparatus of the embodiments as described in the embodiments, or they can be located in one or more devices different from this embodiment with corresponding changes. The modules of the above embodiments can be combined into one module, or they can be 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, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate 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, Includes the following steps: Step S1: Construct a physical response model of surface deformation, gravity change, and groundwater storage change. This physical response model includes a surface deformation response model and a gravity change response model; specifically: Step S101: Based on the dual response characteristics of deep groundwater, which simultaneously possesses pore elastic deformation and load elastic deformation, a surface deformation response model for deep groundwater is constructed. Step S102: Based on the load elastic deformation response characteristics of non-deep groundwater, construct a surface deformation response model for non-deep groundwater; Step S103: Based on the elastic load theory, construct a gravity change response model caused by deep and non-deep groundwater; Step S2: Construct a joint inversion framework for gravity satellite and InSAR data, and establish a set of observation equations; specifically including: Step S201: Decompose the unknown parameter vector into two parts: changes in deep groundwater and non-deep groundwater reserves; Step S202: Construct a joint observation vector, which includes InSAR surface deformation observations and gravity satellite gravity change observations; Step S203: Establish a joint Green's function matrix, which includes the surface deformation response coefficients caused by deep and non-deep groundwater and the gravity change response coefficients; Step S3: Introduce regularization constraints and solve using the weighted least squares method; specifically including: Step S301: Set the weight matrix to balance the contributions of different observation data; Step S302: Introduce a Laplace regularization term to apply a spatial smoothing constraint; Step S303: Solve for the changes in deep and non-deep groundwater storage by minimizing the objective function; Step S4: Vertical separation and verification of groundwater storage based on the principle of mass conservation, specifically including: Step S401: Obtain information on changes in total water storage in the study area from gravity satellite data; Step S402: Extract changes in surface water and soil water storage by combining the land surface hydrological model; Step S403: Calculate the changes in deep groundwater storage using the joint inversion framework; Step S404: Calculate the change in shallow groundwater storage based on the water balance equation; Step S405: Verify the separation results using monitoring well water level data; Deep groundwater refers to groundwater within confined aquifers, while non-deep groundwater includes shallow groundwater, surface water, and soil water.

2. The method according to claim 1, characterized in that, The surface deformation response model in step S1 is expressed as follows: ； In the formula, L InSAR G is the InSAR surface deformation observation vector. poro G is the Green's function for pore elastic deformation. elastic For the load elastic deformation Green's function, The deep groundwater storage is expressed in terms of equivalent thickness.

3. The method according to claim 2, characterized in that, Define the porosity elastic deformation Green's function G poro Then the deformation coefficient G of the j-th model grid cell caused by the change in unit water storage on the surface directly above it is... poro (j) is: ； In the formula, B j Skempton coefficient, α j K represents the Biot-Willis coefficient. j ρ is the bulk modulus of the aquifer. w Let g be the density of water, g be the acceleration due to gravity, and b be the acceleration due to gravity. 0,j Deformation depth; Define the Green's function G for elastic deformation under load. elastic (i,j) is: ； In the formula, G e (i,j) is the Lovegreen function describing the surface deformation caused by a unit mass load, A cell This represents the area of ​​the model's mesh cells.

4. The method according to claim 1, characterized in that, The calculation formula for the gravity change response model in step S1 is as follows: ； In the formula, L T dH is the gravity change observation vector from the gravity satellite, dH is the groundwater storage change vector, and G is the groundwater storage change vector. M This is the Green's function matrix for gravity variation, used to characterize the change in gravity caused by a change in unit water storage. Among them, the Green's function moment G of the unit equivalent water gravity change M The calculation formula is: ； In the formula, A G Let A be the Green's function matrix representing the change in gravitational force of a single unit mass, and let its element A be the Green's function matrix. G (i,j) represents the contribution of the unit mass change of the j-th grid cell to the gravity change at the i-th gravity satellite observation point, ρ w Let A be the density of water. cell This represents the area of ​​the model's mesh cells.

5. The method according to claim 1, characterized in that, The observation equations for step S2 are: ； In the formula, L InSAR For InSAR surface deformation observation, L T For gravity satellite observation of gravity changes, H 深层 and H 非深层 Let G represent the changes in deep and shallow groundwater reserves, respectively, and G be the joint Green's function matrix, containing the deformation and gravity response coefficients of deep and shallow groundwater. poro +G elastic ) 深层 G represents the complete deformation response coefficient of deep groundwater, encompassing both pore elastic and load elastic deformation mechanisms. elastic,非深层 G represents the deformation response coefficient of non-deep groundwater, mainly representing elastic deformation under load. M,深层 and G M,非深层 These represent the response coefficients to gravity changes caused by deep and non-deep groundwater, respectively.

6. The method according to claim 1, characterized in that, The objective function introduced by the regularization constraint in step S3 is: ； In the formula, For the weighted residuals, W is the weight matrix. λ is the regularization term, λ is the regularization parameter, and D is the Laplace operator; The formula for calculating the weight matrix W is as follows: ； in, and These are the inverse matrices of the error covariance matrices of the InSAR observations and the gravity satellite observations, respectively; a and b are relative weighting factors, and satisfy a+b=1.

7. The method according to claim 1, characterized in that, The formula for calculating the solution obtained by the weighted least squares method in step S3 is as follows: ； In the formula, G T Let G be the transpose of the joint Green's function matrix, W be the weight matrix, and D be the weight matrix. T Let D be the transpose of the Laplace operator matrix, (⋅) -1 This represents matrix inversion. The regularization parameter λ is determined through L-curve analysis or cross-validation to balance the roughness of the data fitting residuals and solutions.

8. The method according to claim 1, characterized in that, The water balance relationship on which the vertical separation of groundwater storage in step S4 is based is: ； In the formula, ΔTWS represents the change in total water storage in the study area, ΔSWS represents the change in surface water and soil water storage, and ΔGWS represents the change in total water storage in the study area. 浅层 For changes in shallow groundwater storage, ΔGWS 深层 This represents changes in deep groundwater reserves.

9. The method according to claim 1, characterized in that, The verification in step S405 is specifically carried out by calculating the correlation coefficient, root mean square error, Nash-Sutcliffe efficiency coefficient, and the degree of agreement between the inversion results and the measured monitoring well water level data, as well as 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 orbit satellites. The method also includes a step of modifying the parameters in the physical response model to make the method applicable to different hydrogeological conditions, including karst areas, plains and mountainous areas.

Images