Elastic load earth mass inversion method and system based on spatial structure constraint

Abstract

The application provides a space structure constraint-based elastic load ground mass inversion method and system, and belongs to the field of ground mass change monitoring. The method comprises the following steps: obtaining vertical displacement observation data and ground deformation observation data and performing pretreatment to obtain vertical displacement observation and a vertical displacement field; establishing a GNSS observation equation based on the vertical displacement observation and in combination with an elastic load theory; extracting a space structure feature of the vertical displacement field, constructing an InSAR space structure constraint according to the space structure feature, and introducing a Laplace smoothing constraint; comprehensively constructing a target function by combining the GNSS observation equation, the InSAR space structure constraint and the Laplace smoothing constraint, optimizing and determining constraint parameters in the target function, and obtaining a ground mass change inversion result by solving the target function. The application significantly improves the spatial resolution of ground mass change inversion in the GNSS sparse area while ensuring the reliability of inversion.

Description

Technical Field

[0001] This invention belongs to the field of surface quality change monitoring technology, and particularly relates to a method and system for inverting surface quality under elastic load based on spatial structure constraints. Background Technology

[0002] The statements in this section are merely background information related to the present invention and do not necessarily constitute prior art.

[0003] Surface mass migration alters surface load distribution and induces elastic vertical deformation of the Earth's crust. Based on elastic load theory, changes in surface mass such as terrestrial water storage and glacier mass can be quantitatively inverted through surface deformation, making it a core technical means for water cycle evolution, regional water resource assessment, and cryosphere monitoring. Currently, mainstream observation methods are divided into two categories: Global Navigation Satellite Systems (GNSS) and Interferometric Synthetic Aperture Radar (InSAR). Both single observation methods and existing fusion schemes have significant shortcomings.

[0004] GNSS can output long-term stable vertical displacement time series at the millimeter level. Based on the load Green's function, an elastic load model can be constructed, and the amplitude accuracy of the inversion results is reliable. However, the spatial distribution of GNSS stations is relatively sparse, with very few stations in small and medium-sized watersheds, plateaus, and mountainous areas. The discrete points are difficult to form continuous spatial constraints, resulting in low spatial resolution of the inversion and an inability to characterize the spatial heterogeneity of water storage within the region, thus lacking refined monitoring capabilities.

[0005] Interferometric Synthetic Aperture Radar (InSAR) technology has advantages such as wide coverage, high spatial resolution, and continuous area observation. However, InSAR observations are susceptible to atmospheric delay, orbital errors, and temporal decoherence, resulting in insufficient long-term stability. Observation errors are large in low-coherence regions, and the response capability of LOS deformation acquired by InSAR to vertical elastic load deformation is limited by satellite observation geometry, making it difficult to directly use for quantitative inversion of regional-scale surface mass changes.

[0006] To fully leverage the advantages of GNSS's high precision and InSAR's high spatial resolution, existing technologies have conducted joint deformation inversion studies using GNSS and InSAR. The basic process includes: acquiring GNSS and InSAR observation data; using GNSS to correct InSAR orbit and long-wave errors; constructing unified observation equations and solving them using least squares or weighted fusion methods to reconstruct the surface deformation field; and then analyzing seismic deformation or ground subsidence processes. However, these methods primarily focus on deformation field reconstruction, using GNSS to improve the long-term stability of InSAR and InSAR to improve the spatial continuity of the deformation field. They have not yet established a joint inversion framework for elastic load response oriented towards surface mass migration, and typically employ fixed weights for data fusion, making it difficult to dynamically adjust observational contributions, thus affecting the stability and reliability of the inversion results. Summary of the Invention

[0007] To overcome the shortcomings of the existing technologies, this invention provides a method and system for inverting surface quality based on spatial structure constraints and elastic loads, aiming to solve the technical defects of insufficient GNSS spatial resolution, poor InSAR absolute accuracy, model mismatch in existing joint methods, easy error propagation, and inability to adaptively optimize weights.

[0008] To achieve the above objectives, one or more embodiments of the present invention provide the following technical solutions: The first aspect of this invention provides a method for inverting surface mass under elastic loads based on spatial structural constraints; The elastic load surface quality inversion method based on spatial structure constraints includes: GNSS vertical displacement observation data and InSAR surface deformation observation data are acquired and preprocessed to obtain GNSS vertical displacement observation data and InSAR vertical displacement field. Based on GNSS vertical displacement observations, GNSS observation equations are established in conjunction with elastic load theory; The spatial structure features of the InSAR vertical displacement field are extracted, and the InSAR spatial structure constraint equations are constructed based on the spatial structure features and Laplace smoothing constraints are introduced. By integrating GNSS observation equations, InSAR spatial structure constraints, and Laplace smoothing constraints, an objective function is constructed. The constraint parameters in the objective function are optimized and determined. The surface quality change inversion results are obtained by solving the objective function.

[0009] As a further technical solution, based on GNSS vertical displacement observations and combined with elastic load theory, GNSS observation equations are established, including: The study area was divided into multiple quality cell grids; The vertical deformation response of each grid mass change to the GNSS station is calculated based on the load Green's function. The GNSS observation equations are then established by constructing a design matrix.

[0010] in, This is the GNSS vertical displacement observation vector; Design the matrix for the load Green's function; These are estimated values ​​for the surface mass variation parameters to be solved. This represents the observation error.

[0011] As a further technical solution, the spatial structure features of the InSAR vertical displacement field are extracted, and the InSAR spatial structure constraint equations are constructed based on these features, with Laplace smoothing constraints introduced, including: The InSAR vertical displacement field is subjected to mean-removal processing, the spatial gradient of the mean-removed InSAR vertical displacement field is calculated, and the magnitude of the spatial gradient is normalized. Construct an InSAR spatial structure constraint matrix based on the spatial relationships between adjacent grids, and establish spatial structure constraint equations. A Laplace smoothing constraint is introduced, and a Laplace smoothing constraint equation is established to ensure a smooth transition of surface quality changes between adjacent grids.

[0012] As a further technical solution, the objective function is:

[0013] in, This is the GNSS vertical displacement observation vector; Design the matrix for the load Green's function; These are estimated values ​​for the surface mass variation parameters to be solved. The spatial gradient of the InSAR vertical displacement field; These are the InSAR structural constraint parameters. For smoothing constraint parameters; Corresponding to the two-dimensional Laplace difference operator; This is the InSAR spatial structure constraint matrix.

[0014] As a further technical solution, the constraint parameters in the objective function are optimized and determined, and the surface quality change inversion results are obtained by solving the objective function, including: The optimal weight parameters of the InSAR spatial structure constraint term are determined using the ABIC method, and the optimal weight parameters of the Laplace smoothing constraint term are determined using the GCV method. Substitute the determined optimal weight parameters into the objective function to establish a regularized inversion equation, and solve the regularized inversion equation to obtain the inversion results of surface quality changes.

[0015] As a further technical solution, the ABIC method is used to determine the optimal weight parameters of the InSAR spatial structure constraint term, specifically:

[0016] in, The optimal weighting parameters for the InSAR spatial structure constraint term; The optimal weight parameters of the Laplace smoothing constraint term are determined using the GCV method, specifically as follows:

[0017] in, The optimal weight parameters are the Laplace smoothing constraint terms.

[0018] As a further technical solution, the determined optimal weight parameters are substituted into the objective function to establish a regularized inversion equation, and the regularized inversion equation is solved to obtain the surface quality change inversion results, including: Substituting the optimal weight parameters of the InSAR spatial structure constraint term and the optimal weight parameters of the Laplace smoothing constraint term into the objective function, we obtain the regularized inversion equation:

[0019] Solving the regularized inversion equation yields estimated values ​​of the surface mass change parameters, thus providing the spatiotemporal distribution results of surface mass change in the study area.

[0020] in, These are estimated values ​​for the surface mass variation parameters to be solved. Design the transpose of the matrix for the load Green's function; Design the matrix for the load Green's function; This is the transpose of the InSAR spatial structure constraint matrix; This is the InSAR spatial structure constraint matrix; This is the transpose of the corresponding two-dimensional Laplace difference operator; Corresponding to the two-dimensional Laplace difference operator; This is the GNSS vertical displacement observation vector; This represents the spatial gradient of the InSAR vertical displacement field.

[0021] A second aspect of the present invention provides an elastic load surface quality inversion system based on spatial structure constraints.

[0022] The elastic load surface quality inversion system based on spatial structure constraints includes: The data acquisition and preprocessing module is configured to: acquire GNSS vertical displacement observation data and InSAR surface deformation observation data and perform preprocessing to obtain GNSS vertical displacement observation data and InSAR vertical displacement field. The GNSS observation equation construction module is configured to: establish GNSS observation equations based on GNSS vertical displacement observations and combined with elastic load theory; The InSAR spatial structure constraint construction module is configured to: extract the spatial structure features of the InSAR vertical displacement field, construct the InSAR spatial structure constraint equation based on the spatial structure features, and introduce Laplace smoothing constraints. The joint inversion solution module is configured to: integrate GNSS observation equations, InSAR spatial structure constraints, and Laplace smoothing constraints to construct an objective function, optimize and determine the constraint parameters in the objective function, and obtain the inversion results of surface quality changes by solving the objective function.

[0023] A third aspect of the present invention provides a computer-readable storage medium having a program stored thereon, which, when executed by a processor, implements the steps in the elastic load surface quality inversion method based on spatial structure constraints as described in the first aspect of the present invention.

[0024] A fourth aspect of the present invention provides an electronic device, including a memory, a processor, and a program stored in the memory and executable on the processor, wherein the processor executes the program to implement the steps in the elastic load surface quality inversion method based on spatial structure constraints as described in the first aspect of the present invention.

[0025] The above one or more technical solutions have the following beneficial effects: (1) This invention uses GNSS vertical displacement as the main observation and extracts only the spatial structure features of the InSAR deformation field as constraints. It does not directly introduce the absolute deformation amplitude of InSAR into the inversion, which effectively avoids the influence of insufficient long-term stability of InSAR and systematic errors such as atmospheric delay and orbital error on the inversion of the amplitude of surface quality change. It makes full use of the high-resolution spatial information of InSAR while ensuring the reliability of the inversion results.

[0026] (2) By constructing the InSAR spatial structure constraint equation and introducing the Laplace smoothing constraint, this invention effectively transmits the spatial structure information provided by the InSAR high spatial resolution deformation field to the surface quality change inversion process. This can significantly improve the spatial continuity and spatial resolution of the surface quality change inversion results in sparse GNSS stations, and obtain the spatiotemporal distribution results of land water storage changes with high spatial resolution. This effectively makes up for the shortcomings of the existing GNSS vertical displacement inversion method, which is difficult to characterize the spatial heterogeneity of the region due to the sparse distribution of stations.

[0027] (3) The present invention adopts the ABIC method to adaptively determine the InSAR spatial structure constraint parameters and the GCV method to determine the Laplace smoothing parameters, thereby realizing the objective optimization selection of constraint parameters and overcoming the problem that the fixed weight or empirical weight fusion method in the prior art is difficult to dynamically adjust the contribution degree according to the quality of the observation data. At the same time, the present invention can realize the continuous time series reconstruction of surface quality changes based on continuous GNSS observation, which can effectively make up for the time discontinuity problem of gravity satellite observation data and provide reliable technical support for watershed-scale hydrological change monitoring and regional water resources assessment.

[0028] Advantages of additional aspects of the invention will be set forth in part in the description which follows, and in part will be obvious from the description, or may be learned by practice of the invention. Attached Figure Description

[0029] The accompanying drawings, which form part of this invention, are used to provide a further understanding of the invention. The illustrative embodiments of the invention and their descriptions are used to explain the invention and do not constitute an improper limitation of the invention.

[0030] Figure 1 This is a flowchart of the method in the first embodiment.

[0031] Figure 2 This is a schematic diagram of GNSS stations and InSAR coverage in a certain study area of ​​the first embodiment.

[0032] Figure 3 The diagram illustrates the spatial distribution of the long-term trend, annual amplitude, and annual phase of land surface mass change obtained by the GRACE / GRACE-FO product of the first embodiment and the method of the present invention. (a) shows the spatial distribution of the long-term trend of land surface mass change obtained by GRACE / GRACE-FO; (b) shows the spatial distribution of the annual amplitude of land surface mass change obtained by GRACE / GRACE-FO; (c) shows the spatial distribution of the annual phase of land surface mass change obtained by GRACE / GRACE-FO; (d) shows the spatial distribution of the long-term trend of land surface mass change obtained by the method of the present invention; (e) shows the spatial distribution of the annual amplitude of land surface mass change obtained by the method of the present invention; and (f) shows the spatial distribution of the annual phase of land surface mass change obtained by the method of the present invention.

[0033] Figure 4 This is a schematic diagram illustrating the GRACE / GRACE-FO product of the first embodiment and the time series of surface quality changes retrieved by the present invention.

[0034] Figure 5 This is a system structure diagram of the second embodiment. Detailed Implementation

[0035] It should be noted that the following detailed descriptions are exemplary and intended to provide further illustration of the invention. Unless otherwise specified, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains.

[0036] It should be noted that the terminology used herein is for the purpose of describing particular implementations only and is not intended to limit the exemplary implementations of the present invention.

[0037] Where there is no conflict, the embodiments and features in the embodiments of the present invention can be combined with each other.

[0038] Example 1 This embodiment discloses an elastic load surface quality inversion method based on spatial structure constraints. An elastic load inversion model is established using GNSS vertical displacement observation data. Spatial gradient features of the InSAR vertical displacement field are extracted to construct a spatial structure constraint equation and a Laplace smoothing constraint is introduced. The objective function is constructed by combining the GNSS observation equation and the spatial structure constraint. The constraint parameters are optimized and determined using the ABIC and GCV methods, and the spatiotemporal distribution results of surface quality changes with high spatial resolution are obtained.

[0039] like Figure 1 As shown, the elastic load surface quality inversion method based on spatial structure constraints includes: Step S1: Acquire GNSS vertical displacement observation data and InSAR surface deformation observation data and perform preprocessing to obtain GNSS vertical displacement observation data and InSAR vertical displacement field.

[0040] First, time series data of vertical displacement from continuous GNSS observation stations within the study area were acquired. In this embodiment, the GNSS observation data came from the GNSS displacement product released by the China Crustal Movement Observation Network (CMONOC). This product was processed using GAMIT / GLOBK and QOCA software to perform GNSS observation data calculation and regional network adjustment.

[0041] Gross error removal and outlier detection were performed on the acquired GNSS vertical displacement observation data to remove abnormal observations caused by instrument malfunctions, environmental interference, and other factors. Then, non-hydrological signal corrections were applied to the GNSS vertical displacement time series, including removing the influence of surface deformation caused by non-hydrological factors such as glacier isostatic adjustment (GIA), non-tidal atmospheric load (NTAL), and non-tidal ocean load (NTOL), to obtain GNSS vertical displacement observations that primarily reflect elastic deformation caused by changes in terrestrial water storage. .

[0042] Secondly, InSAR surface deformation observation data covering the study area are acquired. In this embodiment, Sentinel-1 SAR image data is used, and the Small Baseline Set (SBAS) temporal InSAR processing method is employed to obtain the surface deformation field of the study area.

[0043] The acquired InSAR surface deformation observation data were then processed. The specific processing steps included: selecting interferometric pairs based on spatiotemporal baseline thresholds to generate differential interferograms; correcting orbital errors using precise orbital ephemeris data; applying the Goldstein filtering method to filter the interferograms and suppress phase noise; performing phase unwrapping using the minimum cost flow (MCF) method; estimating and removing atmospheric delay phase using spatiotemporal filtering methods; and removing residual orbital errors. Finally, the LOS-direction surface deformation rate and temporal deformation field were obtained. Finally, based on the satellite's incident angle and azimuth parameters, the LOS-direction deformation was converted into a vertical deformation field. Furthermore, bilinear interpolation or kriging interpolation methods are used to sample the InSAR vertical deformation field onto a unified spatial grid, such as a 0.2°×0.2° grid, to achieve spatial reference unification with GNSS observations.

[0044] Step S2: Based on GNSS vertical displacement observations, establish GNSS observation equations in conjunction with elastic load theory.

[0045] The study area is spatially divided into several regular grids. In this embodiment, based on the study area's extent and InSAR spatial resolution, the study area is divided into 0.1°×0.1° equally spaced mass cell grids. The surface mass variation within each grid is considered uniformly distributed, and the grid center point is used as the integration point for subsequent calculations. Assuming there are M grid cells in total, the surface mass variation parameter to be solved is an M×1 dimensional column vector, where each element corresponds to the equivalent water height of the mass variation in one grid cell.

[0046] Based on the elastic load theory, the vertical deformation response caused by a unit mass load change at the i-th GNSS station location is calculated. In this embodiment, the load Green's function is calculated using the Farrell formula, taking into account the Earth's elastic structure and radial stratification characteristics. For a unit mass load within the j-th grid cell, the vertical deformation response coefficient generated at the i-th GNSS station is:

[0047] in, The load Green's function value, Let be the spherical distance from the i-th GNSS station to the center of the j-th grid. Let j be the area of ​​the j-th grid cell. The density of water, Let be the mass change of the j-th grid cell.

[0048] By calculating the vertical deformation response of each grid mass change to all GNSS stations, an M×N load Green's function design matrix H is constructed, where M is the number of grid cells and N is the number of GNSS stations. The element in the i-th row and j-th column represents the vertical deformation response caused by a unit mass change of the j-th grid cell at the i-th GNSS station.

[0049] Establish a linear observation equation between GNSS vertical displacement observations and surface mass change parameters:

[0050] in, This is the GNSS vertical displacement observation vector; Design the matrix for the load Green's function; These are estimated values ​​for the surface mass variation parameters to be solved. This equation represents the observation error. It establishes a quantitative response relationship between GNSS vertical displacement and changes in surface mass, providing an observational constraint basis for subsequent joint inversion.

[0051] Step S3: Extract the spatial structure features of the InSAR vertical displacement field, construct the InSAR spatial structure constraint equation based on the spatial structure features, and introduce the Laplace smoothing constraint. The InSAR spatial structure constraint equation is used to constrain the spatial distribution structure of the surface mass variation parameters.

[0052] Because GNSS stations are sparsely distributed in space, it is difficult to reconstruct the distribution of surface quality variations with high spatial resolution using GNSS observations alone. Therefore, in this embodiment, InSAR continuous deformation fields are used to construct spatial structure constraints.

[0053] First, calculate the spatial gradient of the InSAR vertical displacement field:

[0054] in, This represents the spatial gradient of the InSAR vertical displacement field; , represents the two-dimensional spatial gradient operator.

[0055] Subsequently, an InSAR spatial structure constraint matrix is ​​constructed based on the spatial relationships between adjacent grids. Before constructing spatial structure constraints, the InSAR vertical deformation field is mean-reduced to weaken the impact of long-wavelength errors such as orbital errors and atmospheric delay residuals on the absolute deformation amplitude. Next, the first-order spatial difference between adjacent grids is calculated to obtain the spatial gradient direction and relative change intensity of the InSAR vertical deformation field. Then, the spatial gradient amplitude is normalized or scale-transformed to represent only the spatial deformation structure without directly constraining the surface mass change amplitude. The spatial structure constraint equations are then established:

[0056] in, This is the InSAR spatial structure constraint matrix, used to describe the spatial gradient of surface quality variations. This reflects the spatial variation structure revealed by InSAR observations. This constraint incorporates the high spatial resolution structural information provided by InSAR observations into the surface quality variation inversion process. Furthermore, to ensure a smooth transition of surface quality variations between adjacent grids, a Laplace smoothing constraint is introduced:

[0057] in, This corresponds to the two-dimensional Laplace difference operator.

[0058] Step S4: By integrating the GNSS observation equations, InSAR spatial structure constraints, and Laplace smoothing constraints, an objective function is constructed. The constraint parameters in the objective function are optimized and determined. The surface quality change inversion results are obtained by solving the objective function.

[0059] By combining the GNSS observation equations, InSAR spatial structure constraints, and Laplace smoothing constraints, the objective function is constructed as follows:

[0060] in, Let be the objective function. These are InSAR structure constraint parameters used to balance the relative contributions of GNSS observations and InSAR structure constraints. These are smoothing constraint parameters used to control the degree of spatial smoothness.

[0061] The optimal weight parameters for the InSAR spatial structure constraint term are determined using the Akaike Bayesian Information Criterion (ABIC) method. The ABIC criterion seeks the best balance between model fit and model complexity, and its calculation formula is as follows:

[0062] in, These are the optimal weight parameters for the InSAR spatial structure constraint term.

[0063] The optimal weight parameters of the Laplace smoothing constraint term are determined using the generalized cross-validation (GCV) method, specifically as follows:

[0064] in, The optimal weight parameters are the Laplace smoothing constraint terms.

[0065] Substituting the optimal weight parameters of the InSAR spatial structure constraint term and the optimal weight parameters of the Laplace smoothing constraint term into the objective function, we obtain the regularized inversion equation:

[0066] Solving the regularized inversion equation yields estimated values ​​of the surface mass change parameters, thus providing the spatiotemporal distribution results of surface mass change in the study area.

[0067] in, These are estimated values ​​for the surface quality variation parameters to be solved. Each element corresponds to the equivalent water height value of the surface quality variation in one grid cell of the study area. By rearranging the equivalent water height values ​​of each grid cell according to their spatial location, the spatial distribution field of surface quality variation in the study area can be obtained.

[0068] Finally, key parameters such as long-term rate of change, annual amplitude, and annual phase are extracted from the time series of surface mass changes obtained for each grid cell. The inversion results are compared and verified with gravity satellite products. The effectiveness and reliability of the method of this invention are further verified by calculating the spatial correlation coefficient and time series consistency index.

[0069] Specifically, to verify the effectiveness of the method proposed in this invention, a certain region was selected as the test area. Figure 2 The spatial distribution of InSAR measurement points and GNSS stations within the test area is presented. This includes 795 InSAR observation points and 1 continuous GNSS station. GNSS vertical displacement data are derived from CMONOC observation products, and InSAR deformation data are derived from Sentinel-1 SAR image processing results. This invention selects recent years' GNSS observation data and InSAR deformation field data to conduct a study on the inversion of surface quality changes.

[0070] Using acquired GNSS vertical displacement observation data and spatial structure constraints extracted from the InSAR vertical deformation field, surface mass variation inversion was performed. During the inversion process, the ABIC method was used to determine the InSAR spatial structure constraint parameters. When the structure constraint parameter was set to 3.7, the inversion results achieved the best balance between spatial continuity and observation residuals; therefore, 3.7 was adopted as the optimal structure constraint parameter. The GCV method was used to determine the optimal smoothing factor, yielding a smoothing parameter of 1.9. The spatiotemporal distribution of terrestrial water storage variation with a spatial resolution of 0.1° × 0.1° was obtained.

[0071] The inversion results of this invention are compared and verified with GRACE / GRACE-FO products, such as... Figure 3 As shown. From Figure 3 (a) and Figure 3 As shown in (d), the rate of change in water storage in the study area generally shows an increasing trend. The inversion results indicate that the long-term change rate within the region is mainly distributed between 10 and 30 mm / a, with a relatively larger increase in the southwest region, reaching over 30 mm / a in some areas; the change rate in the northeast region is relatively smaller. The inversion results are basically consistent with the overall spatial distribution characteristics reflected by the GRACE / GRACE-FO products, but the inversion results have higher spatial resolution and can reveal more detailed spatial differences within the region. Figure 3 (b) and Figure 3 As shown in (e), the annual amplitude of water storage in the study area mainly ranges from 20 to 60 mm. The amplitude is larger in the southwest region, with a maximum value approaching 60 mm; the amplitude is relatively smaller in the northeast region. The inversion results are in good agreement with the GRACE / GRACE-FO products in the high-amplitude and low-amplitude regions, while also exhibiting richer local spatial details. Figure 3 (c) and Figure 3 As shown in (f), the annual phase of the study area is relatively stable, mainly concentrated between 270 and 330 days, with the corresponding peak in water storage changes occurring primarily from autumn to early winter each year. The inversion results exhibit the same seasonal variation characteristics as the GRACE / GRACE-FO products, but provide phase distribution information with higher spatial resolution.

[0072] The above results demonstrate that the GNSS observation and InSAR spatial structure-constrained elastic load inversion method for land mass change can obtain high spatial resolution results of land water storage changes in sparse GNSS observation areas, and maintains good consistency with GRACE / GRACE-FO, thus verifying the effectiveness and reliability of the method.

[0073] To further verify the reliability of the method of this invention in terms of time series, the time series of terrestrial water storage changes obtained by inversion was compared with GRACE / GRACE-FO products. Figure 4 As shown in the time series comparison results, the inversion results are in good agreement with GRACE / GRACE-FO in terms of seasonal fluctuations and overall trends. Furthermore, since there is a data gap of approximately 11 months in the latter half of 2017 for GRACE / GRACE-FO, the method of this invention, based on continuous GNSS observations and InSAR spatial structure constraints, can achieve continuous reconstruction of water storage changes during this period, thus effectively compensating for the time discontinuity problem in gravity satellite observations.

[0074] Table 1 presents the fitting results for the long-term rate of change, annual amplitude, and annual phase. Regarding the long-term rate of change, GRACE / GRACE-FO is 34.5 ± 6.4 mm / a, while the inverted result is 18.8 ± 9.6 mm / a. Regarding the annual amplitude, GRACE / GRACE-FO is 25.6 ± 10.3 mm, while the inverted result is 41.2 ± 15.4 mm, indicating a stronger seasonal amplitude characteristic in the inverted result. The annual phases of GRACE / GRACE-FO and the inverted result are 264.7 ± 20.5 days and 271.2 ± 21.6 days, respectively, both corresponding to the peak period from autumn to early winter, indicating a high degree of consistency between the two in their seasonal cycle structure.

[0075] Table 1 Comparison of GRACE / GRACE-FO with the time series of terrestrial water storage retrieved by this invention

[0076] Comprehensive analysis shows that the method of the present invention significantly improves the spatial resolution while maintaining consistency with the temporal variation of GRACE / GRACE-FO products, verifying the effectiveness and reliability of the present invention in high-resolution inversion of terrestrial water storage changes.

[0077] Example 2 This embodiment discloses an elastic load surface quality inversion system based on spatial structure constraints; like Figure 5 As shown, the elastic load surface quality inversion system based on spatial structure constraints includes: The data acquisition and preprocessing module is configured to: acquire GNSS vertical displacement observation data and InSAR surface deformation observation data and perform preprocessing to obtain GNSS vertical displacement observation data and InSAR vertical displacement field. The GNSS observation equation construction module is configured to: establish GNSS observation equations based on GNSS vertical displacement observations and combined with elastic load theory; The InSAR spatial structure constraint construction module is configured to: extract the spatial structure features of the InSAR vertical displacement field, construct the InSAR spatial structure constraint equation based on the spatial structure features, and introduce Laplace smoothing constraints. The joint inversion solution module is configured to: integrate GNSS observation equations, InSAR spatial structure constraints, and Laplace smoothing constraints to construct an objective function, optimize and determine the constraint parameters in the objective function, and obtain the inversion results of surface quality changes by solving the objective function.

[0078] Example 3 The purpose of this embodiment is to provide a computer-readable storage medium.

[0079] A computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the steps in the elastic load surface quality inversion method based on spatial structure constraints as described in Example 1.

[0080] Example 4 The purpose of this embodiment is to provide an electronic device.

[0081] An electronic device includes a memory, a processor, and a program stored in the memory and executable on the processor. When the processor executes the program, it implements the steps in the elastic load surface quality inversion method based on spatial structure constraints as described in Embodiment 1.

[0082] The steps and methods involved in the apparatuses of Embodiments 2, 3, and 4 above correspond to those in Embodiment 1. For specific implementation details, please refer to the relevant description section of Embodiment 1. The term "computer-readable storage medium" should be understood as a single medium or multiple media including one or more instruction sets; it should also be understood as including any medium capable of storing, encoding, or carrying an instruction set for execution by a processor and enabling the processor to perform any of the methods in this invention.

[0083] Those skilled in the art will understand that the modules or steps of the present invention described above can be implemented using general-purpose computer devices. Optionally, they can be implemented using computer-executable program code, thereby allowing them to be stored in a storage device for execution by a computer device, or they can be fabricated as separate integrated circuit modules, or multiple modules or steps can be fabricated as a single integrated circuit module. The present invention is not limited to any particular combination of hardware and software.

[0084] While the specific embodiments of the present invention have been described above in conjunction with the accompanying drawings, this is not intended to limit the scope of protection of the present invention. Those skilled in the art should understand that various modifications or variations that can be made by those skilled in the art without creative effort based on the technical solutions of the present invention are still within the scope of protection of the present invention.

Claims

1. A method for inverting surface mass under elastic load based on spatial structural constraints, characterized in that, include: GNSS vertical displacement observation data and InSAR surface deformation observation data are acquired and preprocessed to obtain GNSS vertical displacement observation data and InSAR vertical displacement field. Based on GNSS vertical displacement observations, GNSS observation equations are established in conjunction with elastic load theory; The spatial structure features of the InSAR vertical displacement field are extracted, and the InSAR spatial structure constraint equations are constructed based on the spatial structure features and Laplace smoothing constraints are introduced. By integrating GNSS observation equations, InSAR spatial structure constraints, and Laplace smoothing constraints, an objective function is constructed. The constraint parameters in the objective function are optimized and determined. The surface quality change inversion results are obtained by solving the objective function.

2. The method for inverting surface mass under elastic load based on spatial structural constraints as described in claim 1, characterized in that, Based on GNSS vertical displacement observations and combined with elastic load theory, GNSS observation equations are established, including: The study area was divided into multiple quality cell grids; The vertical deformation response of each grid mass change to the GNSS station is calculated based on the load Green's function. The GNSS observation equations are then established by constructing a design matrix. in, This is the GNSS vertical displacement observation vector; Design the matrix for the load Green's function; These are estimated values ​​for the surface mass variation parameters to be solved. This represents the observation error.

3. The method for inverting surface mass under elastic load based on spatial structural constraints as described in claim 1, characterized in that, The spatial structure features of the InSAR vertical displacement field are extracted, and the InSAR spatial structure constraint equations are constructed based on these features, with Laplace smoothing constraints introduced, including: The InSAR vertical displacement field is subjected to mean-removal processing, the spatial gradient of the mean-removed InSAR vertical displacement field is calculated, and the magnitude of the spatial gradient is normalized. Construct an InSAR spatial structure constraint matrix based on the spatial relationships between adjacent grids, and establish spatial structure constraint equations. A Laplace smoothing constraint is introduced, and a Laplace smoothing constraint equation is established to ensure a smooth transition of surface quality changes between adjacent grids.

4. The method for inverting surface mass under elastic load based on spatial structural constraints as described in claim 1, characterized in that, The objective function is: in, This is the GNSS vertical displacement observation vector; Design the matrix for the load Green's function; These are estimated values ​​for the surface mass variation parameters to be solved. The spatial gradient of the InSAR vertical displacement field; These are the InSAR structural constraint parameters. For smoothing constraint parameters; Corresponding to the two-dimensional Laplace difference operator; This is the InSAR spatial structure constraint matrix.

5. The method for inverting surface mass under elastic load based on spatial structural constraints as described in claim 1, characterized in that, The constraint parameters in the objective function are optimized and determined. The surface quality change inversion results are obtained by solving the objective function, including: The optimal weight parameters of the InSAR spatial structure constraint term are determined using the ABIC method, and the optimal weight parameters of the Laplace smoothing constraint term are determined using the GCV method. Substitute the determined optimal weight parameters into the objective function to establish a regularized inversion equation, and solve the regularized inversion equation to obtain the inversion results of surface quality changes.

6. The method for inverting surface mass under elastic load based on spatial structural constraints as described in claim 5, characterized in that, The optimal weight parameters for the InSAR spatial structure constraint term are determined using the ABIC method, specifically as follows: in, The optimal weighting parameters for the InSAR spatial structure constraint term; The optimal weight parameters of the Laplace smoothing constraint term are determined using the GCV method, specifically as follows: in, The optimal weight parameters are the Laplace smoothing constraint terms.

7. The method for inverting surface mass under elastic load based on spatial structural constraints as described in claim 5, characterized in that, Substituting the determined optimal weight parameters into the objective function, a regularized inversion equation is established, and the regularized inversion equation is solved to obtain the surface quality change inversion results, including: Substituting the optimal weight parameters of the InSAR spatial structure constraint term and the optimal weight parameters of the Laplace smoothing constraint term into the objective function, we obtain the regularized inversion equation: Solving the regularized inversion equation yields estimated values ​​of the surface mass change parameters, thus providing the spatiotemporal distribution results of surface mass change in the study area. in, These are estimated values ​​for the surface mass variation parameters to be solved. Design the transpose of the matrix for the load Green's function; Design the matrix for the load Green's function; This is the transpose of the InSAR spatial structure constraint matrix; This is the InSAR spatial structure constraint matrix; This is the transpose of the corresponding two-dimensional Laplace difference operator; Corresponding to the two-dimensional Laplace difference operator; This is the GNSS vertical displacement observation vector; This represents the spatial gradient of the InSAR vertical displacement field.

8. A surface mass inversion system for elastic loads based on spatial structural constraints, characterized in that, include: The data acquisition and preprocessing module is configured to: acquire GNSS vertical displacement observation data and InSAR surface deformation observation data and perform preprocessing to obtain GNSS vertical displacement observation data and InSAR vertical displacement field. The GNSS observation equation construction module is configured to: establish GNSS observation equations based on GNSS vertical displacement observations and combined with elastic load theory; The InSAR spatial structure constraint construction module is configured to: extract the spatial structure features of the InSAR vertical displacement field, construct the InSAR spatial structure constraint equation based on the spatial structure features, and introduce Laplace smoothing constraints. The joint inversion solution module is configured to: integrate GNSS observation equations, InSAR spatial structure constraints, and Laplace smoothing constraints to construct an objective function, optimize and determine the constraint parameters in the objective function, and obtain the inversion results of surface quality changes by solving the objective function.

9. A computer-readable storage medium having a program stored thereon, characterized in that, When executed by the processor, the program implements the steps in the elastic load surface quality inversion method based on spatial structure constraints as described in any one of claims 1-7.

10. An electronic device comprising a memory, a processor, and a program stored in the memory and executable on the processor, characterized in that, When the processor executes the program, it implements the steps in the elastic load surface quality inversion method based on spatial structure constraints as described in any one of claims 1-7.

Images