GNSS land water storage joint inversion method and system constrained by measured water level data

Abstract

The present application aims at the problems of strong ill-posedness, significant multi-solution, excessive smoothing and boundary leakage which are prone to occur when only relying on GNSS vertical displacement inversion, and provides a GNSS land water storage joint inversion method and system constrained by actual water level data, through proposing a "observation observability evaluation-lake / constraint source automatic screening-water level to grid equivalent water height (EWH) physical mapping-uncertainty propagation to form dynamic box constraint-constraint adaptive update-projection quasi-Newton fast solution" joint inversion framework, realizing TWS change inversion for GNSS stations with sparse, non-uniform spatial distribution or existing observation hollow area.

Description

Technical Field

[0001] This invention relates to a GNSS method for joint inversion of terrestrial water storage constrained by measured water level data, and more particularly to a method for joint inversion of terrestrial water storage (TWS) by fusing hydrological data and geodetic data at the observation level. This method can be effectively applied to high-precision TWS inversion in areas with many lakes and belongs to the field of GNSS hydrogeometry. Background Technology

[0002] Terrestrial water storage, including surface water, groundwater, soil water, glaciers, snow cover, and water stored in vegetation, directly reflects the water surplus or deficit caused by extreme hydrological events. It serves as an important indicator of global or regional TWS distribution and a marker of climate change. Therefore, developing a high spatiotemporal resolution TWS retrieval method is one of the urgent problems to be solved in the field of GNSS hydrogeometry.

[0003] Quantitative assessment of surface water and soil water (TWS) at the regional or global scale is extremely difficult due to the inability to conduct observations of each component at the same scale and with the same precision. Traditional hydrological and meteorological observations and manual surveys, based on discrete stations, can obtain high-precision hydrological information, but they are costly to build and have limited coverage, failing to provide observational results with sufficient spatiotemporal resolution for large areas. Optical or microwave remote sensing satellites can identify the boundaries and area changes of surface water bodies and estimate shallow soil moisture content; however, they cannot fully reflect TWS due to cloud cover or limitations in penetration depth and accuracy. Hydrological models simplify and abstract complex hydrological phenomena and processes through simulation methods, and can estimate surface water and soil water changes at the regional scale to some extent. However, limited by assimilation mechanisms, model parameters, and input data quality, the results still have significant uncertainties in spatial distribution and accuracy, making it difficult to truly reflect the dynamic changes of different types of water bodies.

[0004] The migration and redistribution of TWS (Transient Waves) cause changes in the gravitational field and elastic deformation of the solid Earth. Satellite gravity technologies (such as GRACE and GRACE-FO) provide high-precision global time-varying gravity fields, offering a new perspective for the study of global and large-scale TWS changes. However, GRACE / GRACE-FO observation technologies still face insurmountable technical obstacles, such as limited temporal resolution (10 days or one month) and spatial resolution (approximately 300 km), large uncertainties (instrument measurement errors, signal leakage, background model bias, etc.), and a high data loss rate. Global Navigation Satellite Systems (GNSS) measure surface displacement with all-weather, all-time, and millimeter-level accuracy. By accurately extracting displacements caused by hydrological loads, the dynamic changes of global or regional TWS can be accurately inverted, far exceeding the spatiotemporal resolution of GRACE / GRACE-FO. Therefore, GNSS provides a TWSSC monitoring method independent of satellite gravity technologies and traditional hydro-meteorological observations. However, existing inversion methods have certain requirements on the distribution and density of GNSS station networks, and are only suitable for small areas with dense and uniformly distributed reference stations. In particular, the TWS inversion effect is poor in areas with sparse station networks.

[0005] Existing joint GNSS and GRACE / GRACE-FO or hydrological models can obtain more homogeneous and higher spatiotemporal resolution TWS variations. However, the weights determined by different methods vary significantly, and existing Earth stratification models using GRACE or hydrological models to calculate virtual station displacements introduce new errors. Furthermore, current joint inversion methods have not inherently overcome the inherent accuracy limitations of either technique. Summary of the Invention

[0006] To address the problems of strong ill-posedness, significant multiple solutions, excessive smoothing, and boundary leakage that easily occur when relying solely on GNSS vertical displacement inversion, this invention provides a joint inversion method for GNSS land water storage constrained by measured water level data. By proposing a joint inversion framework of "observation observability assessment - automatic lake / constraint source screening - water level to grid equivalent water height (EWH) physical mapping - uncertainty propagation to form dynamic box constraints - constraint adaptive update - projection quasi-Newton fast solution", it can realize TWS change inversion for sparse, non-uniform GNSS station spatial distribution or areas with observation gaps.

[0007] According to one aspect of the present invention, a joint inversion method for GNSS land water storage constrained by measured water level data is provided, comprising:

[0008] The vertical displacement time series of GNSS stations in the study area were obtained and preprocessed to extract the GNSS vertical displacement signal mainly based on hydrological load.

[0009] The study area was discretized into multiple grids, and a set of linear observation equations between the equivalent water height change of the grids and the GNSS vertical displacement signal was established based on the elastic load theory.

[0010] Identify grids with weak direct GNSS constraints, and automatically select lake level observation data that meet preset conditions as additional constraint sources within the weak grids and their neighborhoods;

[0011] Based on the proportion of water surface coverage area of ​​the selected lakes in each target constraint grid, a mapping relationship is established between lake water level changes and grid equivalent water height changes; and combined with the water level observation error, area estimation error and the uncertainty of non-lake components, a dynamic box constraint interval that varies with time is constructed for each target constraint grid.

[0012] The linear observation equations and the dynamic box constraint interval are combined to form a constrained optimization problem and solved to output a gridded time series of changes in terrestrial water storage.

[0013] As a further technical solution, weak grids directly constrained by GNSS are identified, and lake level observation data that meet preset conditions are automatically selected within the weak grids and their neighborhoods as additional constraint sources, including:

[0014] A spatial analysis method combining Thiessen polygons and a preset buffer zone is used to define the effective area of ​​GNSS direct constraints, and grids located outside the effective area are identified as weak grids.

[0015] The lake constraint credibility score is calculated based on the lake's coverage area ratio within the grid, water level data integrity, and noise level. Based on the score results, the lake and its corresponding target constraint grid are selected.

[0016] As a further technical solution, the construction of the dynamic box constraint interval includes:

[0017] Calculate the comprehensive lake constraint index based on the cumulative coverage ratio coefficient and lake reliability;

[0018] Based on the comprehensive lake constraint index and the equivalent water height change of the grid, the box constraint strength of the grid is dynamically determined.

[0019] Based on the box constraint strength and the equivalent water height change of the grid, the upper and lower boundaries of the box constraint interval are constructed.

[0020] As a further technical solution, the box constraint strength of the grid is dynamically determined as follows:

[0021] ,

[0022] in, Let the box constraint strength of grid j be , For the equivalent water height change of grid j, The constraint strength coefficient, The comprehensive index of lake constraints for grid j.

[0023] As a further technical solution, the upper and lower boundaries of the box constraint interval are:

[0024] ,

[0025] in, These are the upper and lower boundaries of the box constraints of grid j, respectively.

[0026] As a further technical solution, solving the constrained optimization problem includes:

[0027] A quasi-Newton-like algorithm capable of handling box constraints is used for iterative calculation, and after each iteration, the solution is projected onto the dynamic box constraint interval through a projection operator to ensure that the solution always satisfies the constraints.

[0028] As a further technical solution, during the solution process, the width or constraint strength of the dynamic box constraint interval is adaptively adjusted based on the consistency between the GNSS observation residual and the lake constraint residual.

[0029] According to one aspect of the present invention, a GNSS land water storage joint inversion system constrained by measured water level data is provided, comprising:

[0030] The first main module is used to acquire the vertical displacement time series of GNSS stations in the study area and perform preprocessing to extract the GNSS vertical displacement signal mainly based on hydrological load.

[0031] The second main module is used to discretize the study area into multiple grids and establish a set of linear observation equations between the equivalent water height change of the grid and the GNSS vertical displacement signal based on the elastic load theory.

[0032] The third main module is used to identify grids with weak direct GNSS constraints and automatically select lake water level observation data that meet preset conditions as additional constraint sources within the weak grids and their neighborhoods.

[0033] The fourth main module is used to establish a mapping relationship between lake water level changes and grid equivalent water height changes based on the proportion of water surface coverage area of ​​the selected lakes in each target constraint grid; and to construct a dynamic box constraint interval that varies with time for each target constraint grid by combining water level observation error, area estimation error and uncertainty of non-lake components.

[0034] The fifth main module is used to construct a constrained optimization problem by combining the linear observation equations and the dynamic box constraint interval, and to solve it, outputting a gridded time series of changes in terrestrial water storage.

[0035] According to one aspect of the present invention, a GNSS land water storage joint inversion device constrained by measured water level data is provided, comprising a memory and a processor, wherein the memory stores program instructions that are executed by the processor, and the processor invokes the program instructions to execute the GNSS land water storage joint inversion method constrained by measured water level data.

[0036] According to one aspect of the present invention, a non-transitory computer-readable storage medium is provided, the non-transitory computer-readable storage medium storing computer instructions that cause the computer to execute the GNSS land water storage joint inversion method constrained by measured water level data.

[0037] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0038] 1. Addressing the problem of strong ambiguity caused by sparse station network: Automatically locate weakly constrained GNSS grids through observability assessment, and introduce physically interpretable dynamic water level boundaries into these grids to weaken ill-posedness from a mechanism perspective and reduce the dependence of inversion on a single regularized smoothing.

[0039] 2. The source of constraints can be quantitatively mapped and traced: This invention transforms water level observations into dynamic box constraints of grid EWH by using "lake coverage ratio coefficient / mapping matrix + uncertainty propagation", making the source of constraints, contribution range and credibility explainable and reproducible, rather than empirically imposed boundaries.

[0040] 3. Enhanced robustness: Through a constraint-adaptive update mechanism, when lake regulation, human engineering, or data anomalies cause inconsistencies between water levels and regional TWS, the algorithm can automatically relax constraints or reduce their weights, avoiding distortion of inversion results caused by hard constraints.

[0041] 4. Higher computational efficiency and scalability: The projection quasi-Newton method is used to handle large-scale grid and long-term series problems. Compared with directly constructing augmented Lagrange large matrices for solution, it is more suitable for engineering implementation and large-area applications.

[0042] 5. Wide applicability: It can be used in areas with a high density of lakes and can also be extended to areas with a high density of reservoirs; it can be used on both daily and monthly scales; and it can be extended to other hydrological observations (such as deep well water levels) in the same "observed constraint interval" manner. Attached Figure Description

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

[0044] Figure 1 A schematic diagram of the process for the joint inversion method of GNSS land water storage constrained by measured water level data provided in an embodiment of the present invention.

[0045] Figure 2 This is a schematic diagram of constraint grid determination provided in an embodiment of the present invention.

[0046] Figure 3 A schematic flowchart of the constraint inversion process for the water level-EWH physical mapping provided in an embodiment of the present invention.

[0047] Figure 4 The flowchart for automatic constraint source selection is provided for embodiments of the present invention.

[0048] Figure 5 This is a system structure block diagram provided for an embodiment of the present invention. Detailed Implementation

[0049] Currently, numerous hydrological stations can monitor lake or reservoir water level changes with millimeter-level accuracy, effectively characterizing regional TWS (Terrain Water Storage) variations. Especially in areas with sparse GNSS reference station networks, lake water level changes can effectively compensate for the shortcomings of GNSS-based TWS inversion. Therefore, addressing the problems of high ambiguity and insufficient accuracy in TWS inversion in areas with sparse GNSS station networks, this invention proposes a method for joint GNSS inversion of terrestrial water storage changes based on measured lake water level data constraints. This method directly integrates high-precision lake water level observations and GNSS deformation observations at the observational level. Without relying on the GRACE / hydrological model to construct virtual station displacements as external observations, it introduces lake water level changes to physically constrain the GNSS inversion results, effectively reducing inversion ambiguity and improving the stability and accuracy of terrestrial water storage estimation in areas with sparse station networks. This provides a new technical approach for high-precision TWS inversion in densely populated lake areas and large-scale regions.

[0050] The terms “comprising” and “having”, and any variations thereof, in the specification, claims, and accompanying drawings of this invention are intended to cover a non-exclusive inclusion, such as a process, method, system, product, or apparatus that includes a series of steps or units, not necessarily limited to those explicitly listed, but may include other steps or units not explicitly listed or inherent to such processes, methods, products, or apparatus.

[0051] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, 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, 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. In addition, the technical features of the various embodiments or individual embodiments provided by the present invention can be arbitrarily combined to form new technical solutions. Such combinations are not bound by the order of steps and / or structural composition patterns, but must be based on the ability of those skilled in the art to implement them. When the combination of technical solutions is contradictory or cannot be implemented, it should be considered that such a combination of technical solutions does not exist and is not within the scope of protection claimed by the present invention.

[0052] This invention uses GNSS vertical displacement time series as the main observation. Without relying on the construction of virtual stations, it introduces measured water level information of lakes in the study area, transforming point / area water level observations into quantifiable constraint intervals of grid average mass change. Based on the geometric distribution of the GNSS station network and the sensitivity of the Green's function, it automatically identifies weak grids that need to be constrained, thereby achieving synergistic constraints between hydrological and geodetic observations at the observation level. It is particularly suitable for high spatiotemporal resolution TWS inversion in areas with dense lakes but sparse GNSS station networks.

[0053] Reference Figure 1 As shown, the specific implementation steps of the GNSS land water storage joint inversion method constrained by measured water level data of the present invention are as follows:

[0054] Step 1: Accurate Extraction and Error Removal of GNSS Hydrological Load Signals. This step preprocesses the GNSS vertical displacement time series of the study area, removing the effects of non-tidal oceanic and atmospheric load displacements, thermal expansion effects, and post-ice rebound effects; it then uses least squares fitting to remove long-term linear trends, steps, post-earthquake deformations, and outliers; further, it employs principal component analysis or independent component analysis to identify and eliminate common-mode errors and spatial correlation noise from the reference station network, obtaining the GNSS vertical displacement time series dominated by hydrological load. Details are as follows:

[0055] Daily vertical displacement time series of continuous GNSS reference stations in the study area, as well as water level time series and corresponding water surface area data (which can be measured or remote sensing inversion products) of lake / reservoir water level stations, were collected. GNSS vertical displacement was corrected using non-tidal atmospheric and marine load products; corrections were made for thermal expansion and post-ice rebound effects; and least squares fitting was used to remove long-term linear trend terms, step terms, and post-earthquake exponential decay terms, resulting in a residual series dominated by hydrological load. The fitting model is as follows:

[0056] (1),

[0057] in, for Time and location For reference location of the site, u For linear trend terms, and All are periodic terms of amplitude ( j =1 is the anniversary signal. j =2 is a half-year signal). For step jump, For the step jump moment, w The number of step jumps. v Indicates the first v There are several step jumps, where H is the step function. It is random noise.

[0058] To further remove spatially correlated errors, the residual time series of all stations were constructed into a data matrix. Furthermore, PCA is employed to further refine the hydrological signals. PCA projects relevant variables onto a linearly independent principal component space through orthogonal linear transformation, maximizing data variance. Its eigenvectors are obtained by eigenvalue decomposition of the covariance matrix, describing the main directions of data variation. The principal components are sorted according to the percentage of variance they explain, with the first principal component representing the direction of data variation with the largest variance and contributing the most to the motion of the GNSS network.

[0059] (2),

[0060] in, Let U be the data matrix, U be the principal component matrix (also known as the time function), and V be the spatial weighting function. Principal components are selected based on the cumulative variance contribution rate (e.g., ≥90%) to reconstruct and eliminate common-mode errors, thereby reducing the influence of spatially correlated noise and obtaining a purer hydrological displacement signal.

[0061] Simultaneously, GNSS and water level data are unified to the same sampling time (e.g., daily scale or monthly average), and the observation noise standard deviation is calculated for each GNSS station and water level station (e.g., estimated by residual variance) for subsequent weight matrix construction. Finally, the water level time series of each water level station is uniformly subtracted from the same reference time (e.g., the first day or the average of a specified base period) to obtain the water level change series ΔH, which is then time-aligned with the GNSS hydrological displacement series.

[0062] Step 2, study gridding and forward response modeling. (Refer to...) Figure 2As shown, this step discretizes the study area into a regular grid with a preset spatial resolution, using the TWS variation within the grid (represented by the equivalent water height EWH) as the unknown parameter; each grid is represented by an equivalent disk parameterization, and the Green's function matrix between the grid and the station is calculated using elastic load theory to form a linear observation equation set; and a spatial regularization term is constructed for stable solution. The details are as follows:

[0063] The study area is divided according to a preset spatial resolution. M Each grid cell, with the TWS equivalent water height vector of each grid cell at each time step. As an unknown quantity to be estimated, the grid is represented by an equivalent disk, and the calculation is based on the elastic load theory. M grid N Green's function response matrix of vertical displacement of GNSS stations S , of which elements This represents the contribution of the j-th grid to the vertical displacement of the i-th GNSS station (i.e., the Green's function between the load of the j-th grid and the i-th station):

[0064] (3),

[0065] (4),

[0066] in, The Love number represents the vertical load. G It is the gravitational constant; R Where is the Earth's radius; g is the acceleration due to gravity; For Legendre functions; Angular distance; then ,in This represents the angular distance between the grid load center and the station, where n is the Legendre function order. It is a polynomial. The Green's function represents the relationship between the load of a single grid and that of a single site.

[0067] Due to the limited number of observations, an inversion objective function is constructed, consisting of a data fit difference function and a spatial regularization term. The data fit difference is used to constrain... Sh With observed displacement Spatial regularization is used to suppress ill-conditioned and excessive oscillations caused by sparse station networks (one or a combination of minimum and smoothest models, Laplace smoothing, and anisotropic smoothing can be used). Therefore, the following inversion objective function is constructed:

[0068] (5),

[0069] The first term is the data fit difference function, and the second term is the spatial regularization function, which is used to constrain the smoothness between adjacent grid cells. For regularization parameters, h 0 For reference background scene, This is the data error matrix, where each element is the reciprocal of the error for each observation. The weight matrix has the following form:

[0070] ,

[0071] The individual weight matrix is ​​as follows:

[0072] ,

[0073] in, It is a diagonal matrix representing the spatially related weighting function; It is a finite difference operator. It is the angular distance weighted function matrix.

[0074] Step 3: Identification of Weak GNSS Constraint Grids and Automatic Selection of Lake Constraint Sources. This step, based on the spatial distribution of the GNSS station network and Green's function sensitivity, calculates the "observability index" for each grid, automatically identifying weak grids with insufficient direct GNSS constraints. Within these weak grids and their neighborhoods, lake level stations or water level products are automatically retrieved. A lake constraint credibility score is established by combining factors such as lake surface coverage ratio, data completeness, and noise level, thereby determining the set of lakes used for constraints and the corresponding set of target constraint grids. The approach using Thiessen polygons + buffers in existing manuscripts can be considered as an implementation method. (Refer to...) Figure 3 As shown, the details are as follows:

[0075] Based on GNSS station network distribution and Green's function response matrix S Calculate the observability index for each grid cell. :

[0076] (6),

[0077] in, For the first i Standard deviation of noise observed at each GNSS station. Below the threshold At that time, grid j This means that the grid is identified as having weak GNSS constraints and is included in the external constraint candidate set.

[0078] Thiessen polygons are generated using GNSS stations, and the maximum empty circle center characteristic corresponding to each vertex is utilized to quickly locate the void areas covered by the reference station network. Simultaneously, combining the Green's function decay characteristic, a buffer zone (radius can be 50-100km) is constructed for each station to define areas where direct GNSS constraints are insufficient. Then, referring to... Figure 4 As shown, lake / reservoir water level stations or water level products are retrieved within the weak grid and its neighborhood to establish a candidate lake set. L Each candidate lake was evaluated based on data completeness, lake area ratio, and data noise rate. k Calculate reliability That is, to establish a lake constraint credibility score:

[0079] (7),

[0080] in, For data integrity, Scoring the lake area It is the area of ​​the k-th lake. It's an area error. It is the relative area error of the k-th lake. It is a reference error scale. This represents the area of ​​the smallest lake in the candidate lake set. This represents the area of ​​the largest lake in the set of lakes. clip is the cutoff function. clip(x,0,1) means that when x<0, it takes the value 0; when x>1, it takes the value 1. This represents the noise in the k-level sequence of the lake. Indicates a noise reference scale. This indicates the proportion of missing data for lake k. This represents the data noise rate. (Based on...) Select the top few lakes from high to low as the constraint sources and record their confidence levels.

[0081] Step 4: Water level-grid EWH physical mapping and dynamic box constraint construction.

[0082] Based on the spatial coverage area of ​​lakes within the grid, a mapping coefficient matrix (coverage ratio coefficient) is constructed from lakes to grids. Lake water level changes are converted into equivalent water height changes at the grid scale. Uncertainty propagation is then performed, considering uncertainties related to data completeness, lake area estimation error, and data noise rate, to obtain the boundary constraints for each constrained grid at each time step. Stronger constraints (narrow intervals or even approximate equality constraints) are applied to grids with high coverage ratios and high reliability, while weaker constraints (wide intervals or soft constraints) are applied to grids with moderate coverage ratios or average reliability. (Refer to...) Figure 4 As shown, the details are as follows:

[0083] Calculate lakes k In the grid j Coverage area S j Representing grid j area ， And define the coverage ratio coefficient:

[0084] (8)

[0085] When a grid contains multiple lakes, a mapping matrix can be constructed. Changes in lake water level Mapped to grid-scale EWH contribution:

[0086] (9),

[0087] Finally, based on the cumulative coverage ratio coefficient and lake reliability... Calculate the comprehensive constraint index of lakes :

[0088] (10)

[0089] in, For grid j The sum of the coverage ratios, i.e., the cumulative coverage ratio coefficient. To cover weighted reliability.

[0090] Dynamically determine the box constraint strength of the grid :

[0091] (11),

[0092] In the above formula, Commonly used values ​​are 0.1, 0.4, and 0.8.

[0093] Therefore, the following upper and lower bounds can be constructed:

[0094] (12),

[0095] in, The upper and lower boundaries of the box constraint.

[0096] Step 5: Joint Inversion Solution with Dynamic Box Constraints and Adaptive Constraint Update. This step unifies the data fitting term and spatial regularization term from Step 2 with the dynamic box constraints from Step 4 into a constrained optimization problem. Fast algorithms such as Projective Quasi-Newton (L-BFGS-B) are used to solve the problem, ensuring that the inversion results always satisfy the box constraints during the iteration process. Based on GNSS residual consistency and lake constraint residual consistency, the constraint interval width and confidence weights are adaptively adjusted to output the optimal TWS grid time series and its uncertainty assessment. Details are as follows:

[0097] For grid sets that require water level constraints The dynamic upper and lower bound variables are obtained from step 4. For unconstrained grids, their upper and lower bounds are taken as... and Therefore, based on the objective function in step 2, an optimization problem with box constraints is formed:

[0098] (13)

[0099] First, calculate the gradient of the objective function:

[0100] (14)

[0101] In the t-th iteration, based on the current solution and its gradient Construct the search direction:

[0102] (15)

[0103] in, Let be the direction of the search in the t-th iteration. Let be the quasi-Newton approximation matrix of the Hessian inverse of the objective function, representing the curvature information of the objective function near the current iteration point. This search direction, while ensuring descent performance, has a faster convergence speed compared to the steepest descent method.

[0104] To ensure that the updated inversion results always satisfy the upper and lower boundary constraints of the grid, the search results are projected and corrected:

[0105] (16)

[0106] in, It is the first t The step size parameter for the next iteration. To ensure the objective function decreases, Backtracking can be used to determine this, making The projection operator satisfies the sufficient descent condition and guarantees that the updated value remains within the feasible region of the box constraints.

[0107] (17)

[0108] The above formula guarantees that for grids with water level constraints, the update results are limited to... Within the range; for unconstrained grids, At this point, the projection does not change this component. j This represents the inversion result value from the previous step.

[0109] After completing one iteration, update the quasi-Newton matrix based on the results of two adjacent iterations and the gradient change:

[0110] (18)

[0111] in, This represents the gradient change vector. Therefore, the matrix can be... Update the matrix to obtain the required quasi-Newton approximation matrix for the next iteration. :

[0112] (19)

[0113] In the above formula, For scalar scaling factors (curvature normalization); when Skip or use a damped quasi-Newton method to maintain Zhengding.

[0114] Furthermore, step 6 is included: outputting the grid EWH results obtained from each epoch inversion as time-series products and spatial distribution products, and simultaneously outputting: (1) constrained grid mask file; (2) observability index of each grid. Constraint width Constraint credibility score (3) GNSS fitting residuals, constraint residuals, and logs of abnormal constraint sources are used to support subsequent hydrological interpretation, comparative verification, and long-term operational use.

[0115] The implementation of the various embodiments of the present invention is based on programmed processing by a device with processor functionality. Therefore, in practical engineering, the technical solutions and functions of the various embodiments of the present invention are encapsulated into various modules. Based on this reality, and building upon the above embodiments, the embodiments of the present invention provide a GNSS land water storage joint inversion system constrained by measured water level data. This system is used to execute the GNSS land water storage joint inversion method constrained by measured water level data in the above method embodiments.

[0116] See Figure 5The system includes: a first main module, namely the GNSS data processing module, used to acquire the vertical displacement time series of GNSS stations in the study area and perform preprocessing to extract the GNSS vertical displacement signal mainly based on hydrological load; a second main module, namely the gridding and Green's function forward modeling module, used to discretize the study area into multiple grids and establish a set of linear observation equations between the equivalent water height change of the grid and the GNSS vertical displacement signal based on the elastic load theory; a third main module, including a weak grid identification module, used to identify grids with weak direct GNSS constraints, and a lake candidate and scoring module, used to automatically select lake water level observation data that meet preset conditions as additional constraint sources within the weak grid and its neighborhood; The four main modules include a water level EWH mapping and uncertainty propagation module, used to establish a mapping relationship between lake water level changes and grid equivalent water height changes based on the proportion of water surface coverage area of ​​selected lakes in each target constraint grid; a dynamic box constraint generation module, used to construct a time-varying dynamic box constraint interval for each target constraint grid by combining water level observation errors, area estimation errors, and uncertainties of non-lake components; and a fifth main module, including a projection quasi-Newton solution module, used to construct a constrained optimization problem by combining the linear observation equations and the dynamic box constraint interval and solve it; and a result output and quality assessment module, used to output the grid time series of terrestrial water storage changes and perform quality assessment.

[0117] The GNSS land water storage joint inversion system constrained by measured water level data provided in this invention addresses the problems of strong ill-posedness, significant multiple solutions, excessive smoothing, and boundary leakage that easily occur when relying solely on GNSS vertical displacement inversion. Figure 5 Several modules in the study propose a joint inversion framework consisting of "observation observability assessment - automatic lake / constraint source screening - water level to grid equivalent water height (EWH) physical mapping - uncertainty propagation to form dynamic box constraints - constraint adaptive update - projection quasi-Newton fast solution". This framework enables TWS change inversion for GNSS stations in sparsely distributed, non-uniform, or observation-hole regions.

[0118] It should be noted that the system embodiments provided by the present invention are used not only to implement the methods in the above method embodiments, but also to implement the methods in other method embodiments provided by the present invention. The only difference is that corresponding functional modules are set. The principle is basically the same as that of the above system embodiments provided by the present invention. As long as those skilled in the art can improve the modules in the above system embodiments by referring to the specific technical solutions in other method embodiments and combining technical features to obtain corresponding technical means and technical solutions composed of these technical means, on the basis of the above system embodiments, and on the premise of ensuring the practicality of the technical solutions, they can obtain corresponding system-like embodiments for implementing the methods in other method-like embodiments.

[0119] Based on the same inventive concept as any of the foregoing embodiments, this embodiment of the invention also provides a GNSS land water storage joint inversion device constrained by measured water level data, including a memory and a processor. The memory stores program instructions that are executed by the processor, and the processor calls the program instructions to execute the GNSS land water storage joint inversion method constrained by measured water level data.

[0120] Based on the same inventive concept as any of the foregoing embodiments, this embodiment of the invention also provides a non-transitory computer-readable storage medium storing computer instructions that cause the computer to execute the GNSS land water storage joint inversion method constrained by measured water level data.

[0121] Those skilled in the art will understand that embodiments of the present invention can be provided as methods, systems, or computer program products. Therefore, the present invention can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the present invention can take the form of a computer program product embodied on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0122] This invention is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the invention. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart illustrations and / or block diagrams. Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0123] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0124] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0125] In summary, the present invention describes a novel method for GNSS joint inversion of terrestrial water storage changes based on measured lake water level data constraints. This method uses high-precision measured lake water level data to constrain the TWS changes in sparse areas of GNSS stations, thereby reducing the ambiguity of TWS retrieval by GNSS and improving the accuracy of TWS estimation in sparse areas of the GNSS station network.

[0126] 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 or all of the technical features therein; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the technical solutions of the embodiments of the present invention.

Claims

1. A joint inversion method for GNSS land water storage constrained by measured water level data, characterized in that, include: The vertical displacement time series of GNSS stations in the study area were obtained and preprocessed to extract the GNSS vertical displacement signal mainly based on hydrological load. The study area was discretized into multiple grids, and a set of linear observation equations between the equivalent water height change of the grids and the GNSS vertical displacement signal was established based on the elastic load theory. Identify grids with weak direct GNSS constraints, and automatically select lake level observation data that meet preset conditions as additional constraint sources within the weak grids and their neighborhoods; Based on the proportion of water surface coverage area of ​​the selected lakes within each target constraint grid, a mapping relationship is established between lake water level changes and grid equivalent water height changes. Furthermore, considering water level observation errors, area estimation errors, and uncertainties in non-lake components, a dynamic box constraint interval is constructed for each target constraint grid, varying over time. The construction of the dynamic box constraint interval includes: calculating a comprehensive lake constraint index based on the cumulative coverage ratio coefficient and lake reliability; dynamically determining the box constraint strength of the grid based on the comprehensive lake constraint index and grid equivalent water height changes; and constructing the upper and lower boundaries of the box constraint interval based on the box constraint strength and grid equivalent water height changes. The linear observation equations and the dynamic box constraint interval are combined to form a constrained optimization problem and solved to output a gridded time series of changes in terrestrial water storage.

2. The GNSS land water storage joint inversion method constrained by measured water level data according to claim 1, characterized in that, Identify grids with weak direct GNSS constraints, and automatically select lake level observation data that meet preset conditions within the weak grids and their neighborhoods as additional constraint sources, including: A spatial analysis method combining Thiessen polygons and a preset buffer zone is used to define the effective area of ​​GNSS direct constraints, and grids located outside the effective area are identified as weak grids. The lake constraint credibility score is calculated based on the lake's coverage area ratio within the grid, water level data integrity, and noise level. Based on the score results, the lake and its corresponding target constraint grid are selected.

3. The GNSS land water storage joint inversion method constrained by measured water level data according to claim 1, characterized in that, The box constraint strength of the grid is dynamically determined as follows: ， in, Let the box constraint strength of grid j be , For the equivalent water height change of grid j, The constraint strength coefficient, The comprehensive index of lake constraints for grid j.

4. The GNSS land water storage joint inversion method constrained by measured water level data according to claim 3, characterized in that, The upper and lower boundaries of the box constraint interval are: ， in, These are the upper and lower boundaries of the box constraints of grid j, respectively.

5. The GNSS land water storage joint inversion method constrained by measured water level data according to claim 1, characterized in that, Solving the constrained optimization problem includes: A quasi-Newton-like algorithm capable of handling box constraints is used for iterative calculation, and after each iteration, the solution is projected onto the dynamic box constraint interval through a projection operator to ensure that the solution always satisfies the constraints.

6. The GNSS land water storage joint inversion method constrained by measured water level data according to claim 1 or 5, characterized in that, During the solution process, the width or constraint strength of the dynamic box constraint interval is adaptively adjusted based on the consistency between the GNSS observation residual and the lake constraint residual.

7. A GNSS land water storage joint inversion system constrained by measured water level data, characterized in that, include: The first main module is used to acquire the vertical displacement time series of GNSS stations in the study area and perform preprocessing to extract the GNSS vertical displacement signal mainly based on hydrological load. The second main module is used to discretize the study area into multiple grids and establish a set of linear observation equations between the equivalent water height change of the grid and the GNSS vertical displacement signal based on the elastic load theory. The third main module is used to identify grids with weak direct GNSS constraints and automatically select lake water level observation data that meet preset conditions as additional constraint sources within the weak grids and their neighborhoods. The fourth main module is used to establish a mapping relationship between lake water level changes and grid equivalent water height changes based on the proportion of water surface coverage area of ​​selected lakes within each target constraint grid. It also constructs a dynamic box constraint interval for each target constraint grid, taking into account water level observation errors, area estimation errors, and uncertainties in non-lake components. The construction of the dynamic box constraint interval includes: calculating a comprehensive lake constraint index based on the cumulative coverage ratio coefficient and lake reliability; dynamically determining the box constraint strength of the grid based on the comprehensive lake constraint index and the grid equivalent water height changes; and constructing the upper and lower boundaries of the box constraint interval based on the box constraint strength and the grid equivalent water height changes. The fifth main module is used to construct a constrained optimization problem by combining the linear observation equations and the dynamic box constraint interval, and to solve it, outputting a gridded time series of changes in terrestrial water storage.

8. A GNSS land water storage joint inversion device constrained by measured water level data, characterized in that, It includes a memory and a processor, the memory storing program instructions that are executed by the processor, and the processor calling the program instructions to execute the GNSS land water storage joint inversion method constrained by measured water level data as described in any one of claims 1 to 5.

9. A non-transitory computer-readable storage medium, characterized in that, The non-transitory computer-readable storage medium stores computer instructions that cause the computer to execute the GNSS land water storage joint inversion method constrained by measured water level data as described in any one of claims 1 to 5.

Images