Terrestrial water storage collaborative inversion method fusing hydrological model and sparse GNSS data

Abstract

The application discloses a kind of land water storage collaborative inversion methods of fusing hydrological model and sparse GNSS data, including obtaining and weighting pre-processing GNSS vertical displacement sequence, decomposing site time series, obtaining hydrological dominant residual, then spatial interpolation is carried out to obtain continuous residual field;The land water storage change given by hydrological model is converted into model deformation by load Green function convolution, and is weighted and fused with residual field to construct initial field;With land water storage change grid as state quantity, establish observation equation, construct joint objective function, iteratively solve until convergence;Linear approximation or Monte Carlo is given to pixel level uncertainty and quality assessment.The application can realize high-fidelity reconstruction and uncertainty quantification of TWS and deformation field under the actual conditions of sparse station, significant observation noise and prior existing deviation in the unified physical-statistical framework.

Description

Technical Field

[0001] This invention belongs to the field of GNSS hydrogeometry and geophysical inversion technology, specifically involving a method for the collaborative inversion of land water storage by integrating hydrological models and sparse GNSS data. Background Technology

[0002] Terrestrial water storage (TWS) is a key state quantity in the global water cycle, and its spatiotemporal variations are closely related to climate variability, droughts and floods, water resource regulation, cryosphere evolution, and the load response of the solid Earth. In addition to the monthly mass variations provided by the GRACE / GRACE-FO satellite gravity data, the Global Navigation Satellite System (GNSS), by measuring three-dimensional displacement of the Earth's surface, especially vertical displacement, can record elastic deformation caused by hydrological loads from another observational dimension, providing a powerful supplement to regional to watershed-scale hydrological monitoring.

[0003] However, due to the sparse and uneven distribution of GNSS networks, incomplete data acquisition, or insufficient openness, the extraction of hydrological signals relying solely on GNSS faces multiple difficulties: (1) Insufficient spatial sampling and interpolation errors: The sparse and uneven distribution of stations leads to unstable spatial reconstruction, and conventional inverse distance, spline, or kriging interpolation significantly amplifies errors at boundaries, in the window area, and in high gradient regions. (2) Mixing of non-hydrological components: GNSS displacement sequences contain multiple sources of interference, such as tectonic deformation, thermoelastic effects of equipment and base, atmospheric and oceanic loads, jumps, and drifts. If only simple trend / seasonal decomposition and low-order filtering are used, it is difficult to cleanly separate and maintain the physical consistency of hydrological signals. (3) Window periods and missing measurements: Intermittent missing measurements at stations, extreme weather / equipment interruptions, and the alternation of old and new stations lead to discontinuities in the time axis, making traditional batch processing methods prone to breakpoints and boundary effects.

[0004] To address this, researchers considered combining prior models with GNSS station data for more accurate hydrological monitoring. However, commonly used hydrological models, such as the Water-Global Analysis and Prognosis Model (WGHM) and the Global Land Data Assimilation System (GLDAS), are susceptible to strong uncertainties caused by precipitation, evapotranspiration, soil parameters, and human activities, leading to prior model bias. Furthermore, the model resolution differs from the representative scale of GNSS stations, and direct difference between the GNSS and the model easily introduces prior bias into the hydrological estimation results. Additionally, TWS exhibits significant seasonal-interannual-multi-interannual amplitude and phase drift, demonstrating spatial heterogeneity and temporal non-stationarity. Filtering / interpolation with fixed windows or parameters can easily result in underfitting or oversmoothing.

[0005] Furthermore, existing hydrological monitoring methods often present a "single optimal field," lacking pixel-level error bands and sensitivity analysis to prior / weight selection, making it difficult to quantify risks in operational applications and decision-making. Direct global optimization of the inverse problem combining GNSS and hydrological priors faces bottlenecks in computation and convergence stability; a robust engineering framework for rapid iteration and local updates is also lacking. Summary of the Invention

[0006] The present invention aims to at least partially solve one of the technical problems existing in the related art.

[0007] The purpose of this invention is to address the problems of sparse GNSS stations, mixed non-hydrological components, biases in prior models, and difficulty in quantifying uncertainties. It provides a collaborative inversion method for land water storage that integrates hydrological models and sparse GNSS data. The method uses GNSS vertical displacement denoising residuals as observations, establishes a "mass-deformation" mapping using the load Green's function, and uses TWS and raster as state variables. It adaptively balances data consistency, spatial smoothness, and prior shrinkage, obtaining high-fidelity TWS and deformation fields through residual backfeedback and iterative solution, and providing pixel-level uncertainties.

[0008] To achieve the above objectives, the present invention provides a method for the collaborative inversion of land water storage by integrating hydrological models and sparse GNSS data, comprising the following steps:

[0009] S1. Obtain the coordinate time series of each GNSS station in the study area, perform quality control, and extract the vertical component coordinate time series of the GNSS stations; simultaneously obtain the TWS prior field provided by the hydrological model, and generate the model-predicted deformation field through convolution with the Earth's elastic load Green's function.

[0010] S2. Decompose the time series of the vertical component coordinates of each GNSS station, and obtain the GNSS residual dominated by hydrological load after deducting non-hydrological terms; perform spatial interpolation on the GNSS residual to form a continuous residual deformation field.

[0011] S3. The continuous residual deformation field and the model-predicted deformation field are weighted and fused to obtain the initial fused deformation field, and the initial land water storage is obtained by reverse calculation.

[0012] S4. Using the rasterized terrestrial water storage as the state variable, establish a linear observation equation and construct a joint objective function that includes a data consistency term, a spatial smoothing term, and a priori shrinkage term.

[0013] S5. Using the initial land water storage as the initial field, the joint objective function is solved by residual recharge and gradient-type iterative algorithm. In each iteration, deformation prediction, residual calculation, residual space interpolation and state variable update are performed sequentially until the convergence condition is met, and the final land water storage is obtained.

[0014] S6. Based on linear covariance propagation or Monte Carlo perturbation, perform uncertainty assessment on the final terrestrial water storage; output the final terrestrial water storage and the corresponding uncertainty assessment results.

[0015] A further preferred technical solution of the present invention is that, in step S1, the coordinate time series of each GNSS station in the study area are obtained, and quality control is performed to extract the vertical component coordinate time series of the GNSS stations; specifically:

[0016] The coordinate time series of each GNSS station in the study area were obtained, and missing values ​​and outliers were filled, outliers were removed, and the reference frame and benchmark were made consistent.

[0017] Extracting the time series of vertical component coordinates of GNSS stations ,in For station indexing, For the number of stations, For time epoch index, The total number of epochs;

[0018] Based on the estimated observation noise level, adaptive weighting is applied according to the observation quality of each station at that epoch, forming the epoch weights of the diagonal weight matrix. .

[0019] Preferably, in step S2, the time series of the vertical component coordinates of each GNSS station is decomposed, and the GNSS residual dominated by hydrological load is obtained after deducting non-hydrological terms; specifically:

[0020] The coordinate time series of each GNSS station is decomposed, and a trajectory model is established, represented as follows:

[0021] ;

[0022] ;

[0023] ;

[0024] ;

[0025] The GNSS residuals, dominated by hydrological load, are obtained by subtracting non-hydrological terms from the time series of vertical component coordinates of GNSS stations, and are expressed as follows:

[0026] ;

[0027] in, This represents the trend term indicating long-term deformation and baseline drift; For periods representing six months or a year; A step term caused by changes in equipment or environment; The residuals after deducting non-hydrological items are dominated by deformation caused by hydrological loads; To observe random noise, its parameters are used to construct epoch weights. ; It is the Heaviside step function; The constant term represents the initial position; Expressed as the rate of long-term linear deformation; The first The sine and cosine coefficients of each periodic term; Angular frequency; Indicates the first The amplitude of a step jump; The moment when the step jump occurs; Indicates the periodic component index. Indicates its total number; Represented as the index of the transition event. The set of its values.

[0028] Preferably, step S2 involves spatial interpolation of the GNSS residuals to form a continuous residual deformation field; specifically:

[0029] Ordinary Kriging interpolation is used for epochs. GNSS residuals Spatial interpolation is performed to obtain the continuous residual deformation field. , represented as:

[0030] ;

[0031] in, Representing geographic coordinates; using an exponential / spherical harmonic Gaussian model to analyze the sample semivariogram. By fitting the data, the corresponding covariance function can be obtained. Based on this, a Kriging covariance matrix is ​​constructed. Covariance vector of the interpolation point ; For Kriging weights, for A set; These are Lagrange multipliers used to ensure that the sum of the Kriging weights is 1.

[0032] Preferably, step S1 involves synchronously acquiring the TWS prior field provided by the hydrological model and generating the model-predicted deformation field through convolution with the Earth's elastic load Green's function; specifically:

[0033] Constructing a vertical load Green's function based on the load Love number. Pre-calculation of design matrix ;

[0034] Obtain changes in terrestrial water storage synthesized from hydrological models And through the Earth's elastic load Green's function Convolution is converted into the model's predicted deformation field, expressed as:

[0035] ;

[0036] in, The geographic coordinates of the integration point. The distance is spherical. For the density of water, It is the acceleration due to gravity. The area of ​​the face unit;

[0037] Using design matrix Discretize the model-predicted deformation field and represent it as follows:

[0038] ;

[0039] in, For the first Each grid area.

[0040] Preferably, in step S3, the continuous residual deformation field and the model-predicted deformation field are weighted and fused to obtain the initial fused deformation field, and the initial terrestrial water storage is then calculated in reverse; specifically:

[0041] continuous residual deformation field With model prediction of deformation field Weighted fusion yields the initial fused deformation field, expressed as:

[0042]

[0043] in, To integrate the weights, an adaptive selection is made using L-curve cross-validation, generalized cross-validation, or K-fold cross-validation. The larger the value, the more trust is placed in the model's prior knowledge;

[0044] Then, the initial terrestrial water storage is determined by inverse calculation based on the initial fusion deformation field or by combining it with the prior TWS convexity. .

[0045] Preferably, step S4 involves establishing a linear observation equation using rasterized terrestrial water storage as the state variable, and constructing a joint objective function that includes a data consistency term, a spatial smoothing term, and a priori shrinkage term; specifically:

[0046] The gridded terrestrial water storage to be sought For the unknown state variables, establish the observation equation, expressed as:

[0047] ;

[0048] in, For the residuals after deducting non-hydrological items, To observe random noise; To design a matrix, design a matrix elements for:

[0049] ;

[0050] in, The distance is spherical. For the density of water, It is the acceleration due to gravity. For the first The area of ​​each grid cell The Green's function for vertical load is constructed based on the load Love number;

[0051] Define the joint objective function:

[0052] ;

[0053] in, This is a data consistency term that minimizes the weighted residuals. For spatial smoothness terms, use the spatial smoothness regularization operator. Suppressing non-physical high frequencies; For the prior reduction term, with the prior weight matrix The constraint solution does not deviate excessively Hyperparameters The method is adaptively determined using L-curve, generalized cross-validation, or K-fold cross-validation.

[0054] The gradient of the joint objective function is:

[0055] ;

[0056] The iteration uses Landweber gradient descent, expressed as:

[0057] ;

[0058] Step length satisfy .

[0059] As a preferred option, step S5 specifically involves:

[0060] Based on initial terrestrial water storage For the initial field, in the first... In the next iteration, according to Computation site predicts deformation ;

[0061] Calculate the observation residuals ;

[0062] For observation residuals Spatial interpolation is performed to obtain the residual field. ;

[0063] Based on the spatial smoothness coefficient of the residual field superposition With prior shrinkage coefficient The updated hydrological signal estimate is expressed as:

[0064] ;

[0065] in To update the step size;

[0066] When the relative changes in the objective function and the residual RMS are both below the threshold, it is expressed as:

[0067] ;

[0068] ;

[0069] If the preset maximum number of iterations is reached, the iteration terminates, and the final terrestrial water storage is obtained. .

[0070] In another aspect, the present invention provides a non-transitory computer-readable storage medium having computer instructions stored thereon, which cause a computer to execute the above-described method for the joint inversion of land water storage by fusing hydrological models and sparse GNSS data.

[0071] In another aspect, the present invention provides an electronic device, comprising: a processor, a communication interface, a memory, and a communication bus, wherein the processor, the communication interface, and the memory communicate with each other through the communication bus, and the processor calls logical instructions in the memory to execute the above-mentioned method for the coordinated inversion of land water storage by fusing hydrological models and sparse GNSS data.

[0072] In another aspect, the present invention provides a computer program product, the computer program product comprising a computer program stored on a non-transitory computer-readable storage medium, wherein when the computer program is executed by a processor, the computer executes the above-described method for the collaborative inversion of land water storage by integrating hydrological models and sparse GNSS data.

[0073] Beneficial effects: The land water storage co-inversion method of the present invention, which integrates hydrological models and sparse GNSS data, achieves high-fidelity reconstruction and uncertainty quantification of TWS and deformation fields under a unified physical-statistical framework for sparse GNSS station networks. It can achieve high-fidelity reconstruction and uncertainty quantification of TWS and deformation fields under actual conditions of sparse GNSS stations, significant observation noise, and prior bias, providing reliable technical support for watershed water resources assessment, drought monitoring, and load response research.

[0074] Compared with the prior art, the present invention achieves the following effects:

[0075] 1) A unified "physical-statistical" framework: GNSS residual observations, hydrological priors and load Green kernels are coupled in the same inverse problem, avoiding the propagation of biases caused by direct differential GNSS and significantly reducing spatial interpolation errors under sparse station networks.

[0076] 2) Balancing fidelity and noise reduction: Through the triple constraints of "data consistency, spatial smoothness, and prior shrinkage", the physical authenticity of TWS amplitude and phase is maintained as much as possible while suppressing non-physical high-frequency and observation noise.

[0077] 3) Adaptation to sparsity and gaps: The residual recharge and prior shrinkage mechanism are robust to station network inhomogeneity and missing values, and can maintain spatiotemporal continuity and stability during gaps.

[0078] 4) Pixel-level uncertainty: Provides covariance, confidence interval and various quality control indicators to facilitate reliability assessment and risk quantification in business applications.

[0079] 5) Computationally efficient and easy to implement: Mapping kernel pre-computation, sparse linear algebra, and optional fast update strategies ensure engineering availability over a wide range and long time series.

[0080] 6) High scalability: It can be extended to joint inversion of E / N / U three components, introduce more priors (such as watershed water conservation), or be jointly assimilated with GRACE / GRACE-FO data to form a more complete "hydrological-gravity" multi-source fusion framework. Attached Figure Description

[0081] Figure 1 This is a flowchart of the method for synergistic inversion of land water storage that integrates hydrological models and sparse GNSS data, as described in this invention. Detailed Implementation

[0082] To make the objectives, technical solutions, and advantages of this invention clearer, the technical solutions of this invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, embodiments of this invention, and should not be construed as limiting the invention. All other embodiments obtained by those skilled in the art based on the embodiments of this invention without creative effort are within the scope of protection of this invention. In the description of this invention, it should be understood that the terminology used is for descriptive purposes only and should not be construed as indicating or implying relative importance.

[0083] The following is combined Figure 1 This invention describes a method for the collaborative inversion of land water storage by integrating a hydrological model and sparse GNSS data.

[0084] Example 1: This example provides a method for the collaborative inversion of land water storage by integrating hydrological models and sparse GNSS data, such as... Figure 1 As shown, the study unfolds within a unified framework of "observation-physical kernel-regularization". Using the GNSS vertical displacement residual as the observed quantity, a physical mapping between TWS and displacement is achieved through a load Green's function. With the TWS grid field as the unknown state, an objective function that balances data consistency, spatial smoothness, and prior shrinkage is constructed. Monthly TWS and its uncertainty are obtained through "residual backfeeding + iterative solution".

[0085] Specifically, the following steps are included:

[0086] S1. Obtain the coordinate time series of each GNSS station in the study area, perform quality control, and extract the vertical component coordinate time series of the GNSS stations; specifically:

[0087] The coordinate time series of each GNSS station in the study area were obtained, and missing values ​​and outliers were filled, outliers were removed, and the reference frame and benchmark were made consistent.

[0088] Extracting the time series of vertical component coordinates of GNSS stations ,in For station indexing, For the number of stations, For time epoch index, The total number of epochs;

[0089] Based on the estimated observation noise level, adaptive weighting is applied according to the observation quality of each station at that epoch, forming the epoch weights of the diagonal weight matrix. .

[0090] The TWS prior field provided by the hydrological model is acquired synchronously, and the model-predicted deformation field is generated through convolution with the Green's function of Earth's elastic load; specifically:

[0091] Constructing a vertical load Green's function based on the load Love number. Pre-calculation of design matrix ;

[0092] Obtain changes in terrestrial water storage synthesized from hydrological models And through the Earth's elastic load Green's function Convolution is converted into the model's predicted deformation field, expressed as:

[0093] ;

[0094] in, The geographic coordinates of the integration point. The distance is spherical. For the density of water, It is the acceleration due to gravity. The area of ​​the face unit;

[0095] Using design matrix Discretize the model-predicted deformation field and represent it as follows:

[0096] ;

[0097] in, For the first Each grid area.

[0098] The Green kernel is constructed from the load Love number and takes into account the land / sea mask / cell area (equal area or area correction); the kernel is pre-computed as a sparse mapping matrix to improve efficiency.

[0099] S2. Decompose the time series of the vertical component coordinates of each GNSS station, and obtain the GNSS residual dominated by hydrological load after deducting non-hydrological terms; specifically:

[0100] The coordinate time series of each GNSS station is decomposed, and a trajectory model is established, represented as follows:

[0101] ;

[0102] ;

[0103] ;

[0104] ;

[0105] The GNSS residuals, dominated by hydrological load, are obtained by subtracting non-hydrological terms from the time series of vertical component coordinates of GNSS stations, and are expressed as follows:

[0106] ;

[0107] in, This represents the trend term indicating long-term deformation and baseline drift; For periods representing six months or a year; A step term caused by changes in equipment or environment; The residuals after deducting non-hydrological items are dominated by deformation caused by hydrological loads; To observe random noise, its parameters are used to construct epoch weights. ; It is the Heaviside step function; The constant term represents the initial position; Expressed as the rate of long-term linear deformation; The first The sine and cosine coefficients of each periodic term; Angular frequency; Indicates the first The amplitude of a step jump; The moment when the step jump occurs; Indicates the periodic component index. Indicates its total number; Represented as the index of the transition event. The set of its values.

[0108] Spatial interpolation is performed on the GNSS residuals to form a continuous residual deformation field; specifically:

[0109] Ordinary Kriging interpolation is used for epochs. GNSS residuals Spatial interpolation is performed to obtain the continuous residual deformation field. , represented as:

[0110] ;

[0111] in, Representing geographic coordinates; using an exponential / spherical harmonic Gaussian model to analyze the sample semivariogram. By fitting the data, the corresponding covariance function can be obtained. Based on this, a Kriging covariance matrix is ​​constructed. Covariance vector of the interpolation point ; For Kriging weights, for A set; These are Lagrange multipliers used to ensure that the sum of the Kriging weights is 1.

[0112] S3. The continuous residual deformation field and the model-predicted deformation field are weighted and fused to obtain the initial fused deformation field, and the initial terrestrial water storage is then calculated in reverse; specifically:

[0113] continuous residual deformation field With model prediction of deformation field Weighted fusion yields the initial fused deformation field, expressed as:

[0114]

[0115] in, To integrate the weights, an adaptive selection is made using L-curve cross-validation, generalized cross-validation, or K-fold cross-validation. The larger the value, the more trust is placed in the model's prior knowledge;

[0116] Then, the initial terrestrial water storage is determined by inverse calculation based on the initial fusion deformation field or by combining it with the prior TWS convexity. .

[0117] S4. Using rasterized terrestrial water storage as the state variable, establish a linear observation equation and construct a joint objective function including a data consistency term, a spatial smoothing term, and a priori shrinkage term; specifically:

[0118] The gridded terrestrial water storage to be sought For the unknown state variables, establish the observation equation, expressed as:

[0119] ;

[0120] in, For the residuals after deducting non-hydrological items, To observe random noise; To design a matrix, design a matrix elements for:

[0121] ;

[0122] in, The distance is spherical. For the density of water, It is the acceleration due to gravity. For the first The area of ​​each grid cell The Green's function for vertical load is constructed based on the load Love number;

[0123] Define the joint objective function:

[0124] ;

[0125] in, This is a data consistency term that minimizes the weighted residuals. For spatial smoothness terms, use the spatial smoothness regularization operator. Suppressing non-physical high frequencies; For the prior reduction term, with the prior weight matrix The constraint solution does not deviate excessively ;

[0126] hyperparameters Adaptive selection is achieved through L-curve, generalized cross-validation (GCV), or K-fold cross-validation; when there are too few available GNSS stations at a given epoch, the efficiency is improved. and increase To enhance prior and smoothness constraints; alternatively, time regularization can be added to the objective function. To constrain temporal continuity; the method supports joint inversion of the E / N / U components and expansion using an elastic Green's tensor kernel.

[0127] The gradient of the joint objective function is:

[0128] ;

[0129] The iteration uses Landweber gradient descent, expressed as:

[0130] ;

[0131] Step length satisfy .

[0132] S5. Using the initial terrestrial water storage as the initial field, the joint objective function is solved using a residual recharge and gradient-based iterative algorithm. Each iteration sequentially performs deformation prediction, residual calculation, residual space interpolation, and state variable update until the convergence condition is met, yielding the final terrestrial water storage. Specifically:

[0133] Based on initial terrestrial water storage For the initial field, in the first... In the next iteration, according to Computation site predicts deformation ;

[0134] Calculate the observation residuals ;

[0135] For observation residuals Spatial interpolation is performed to obtain the residual field. ;

[0136] Based on the spatial smoothness coefficient of the residual field superposition With prior shrinkage coefficient The updated hydrological signal estimate is expressed as:

[0137] ;

[0138] in To update the step size;

[0139] When the relative changes in the objective function and the residual RMS are both below the threshold, it is expressed as:

[0140] ;

[0141] ;

[0142] If the preset maximum number of iterations is reached, the iteration terminates, and the final terrestrial water storage is obtained. .

[0143] S6. Based on linear covariance propagation or Monte Carlo perturbation, perform uncertainty assessment on the final terrestrial water storage; output the final terrestrial water storage and the corresponding uncertainty assessment results.

[0144] For the uncertainty assessment of the final terrestrial water storage, this embodiment provides two complementary solutions:

[0145] (1) Linear approximation: The pixel-level covariance is estimated near the final solution using a second-order approximation (equivalent to the inverse of the Hessian or Fisher information matrix of the objective function) to obtain the standard deviation and confidence band.

[0146] (2) Monte Carlo: Perturbation recalculation of observation noise and hyperparameters, and statistical mean / variance and sensitivity.

[0147] It also outputs quality metrics such as: fitting residual RMS, cross-validation error, residual power spectrum / whitening test, and iterative convergence curve.

[0148] The output includes the following parts:

[0149] 1) Main Results: Monthly TWS Grid Fields and pixel-level standard deviation / confidence interval;

[0150] 2) Derivative products: Integral curves (trend, annual / semi-annual amplitude and phase) of watershed or polygonal regions; observation-prediction comparison and residual statistics; spatial power spectrum and frequency band energy; window period reconstruction comparison;

[0151] 3) Metadata and Quality Control: Design Matrix (F), Grid / Projection, The system includes information on iteration count and threshold, weight construction, time span, and missing test strategies, which facilitates reproduction and comparison.

[0152] Example 2: This example provides a non-transitory computer-readable storage medium storing computer instructions that cause a computer to execute a method for co-retrieval of land water storage by fusing hydrological models and sparse GNSS data. The method includes the following steps:

[0153] S1. Obtain the coordinate time series of each GNSS station in the study area, perform quality control, and extract the vertical component coordinate time series of the GNSS stations; simultaneously obtain the TWS prior field provided by the hydrological model, and generate the model-predicted deformation field through convolution with the Earth's elastic load Green's function.

[0154] S2. Decompose the time series of the vertical component coordinates of each GNSS station, and obtain the GNSS residual dominated by hydrological load after deducting non-hydrological terms; perform spatial interpolation on the GNSS residual to form a continuous residual deformation field.

[0155] S3. The continuous residual deformation field and the model-predicted deformation field are weighted and fused to obtain the initial fused deformation field, and the initial land water storage is obtained by reverse calculation.

[0156] S4. Using the rasterized terrestrial water storage as the state variable, establish a linear observation equation and construct a joint objective function that includes a data consistency term, a spatial smoothing term, and a priori shrinkage term.

[0157] S5. Using the initial land water storage as the initial field, the joint objective function is solved by residual recharge and gradient-type iterative algorithm. In each iteration, deformation prediction, residual calculation, residual space interpolation and state variable update are performed sequentially until the convergence condition is met, and the final land water storage is obtained.

[0158] S6. Based on linear covariance propagation or Monte Carlo perturbation, perform uncertainty assessment on the final terrestrial water storage; output the final terrestrial water storage and the corresponding uncertainty assessment results.

[0159] Example 3: This example provides an electronic device that may include a processor, a communication interface, a memory, and a communication bus. The processor, communication interface, and memory communicate with each other via the communication bus. The processor can call logical instructions from the memory to execute a method for the collaborative inversion of land water storage by fusing hydrological models and sparse GNSS data. This method includes the following steps:

[0160] S1. Obtain the coordinate time series of each GNSS station in the study area, perform quality control, and extract the vertical component coordinate time series of the GNSS stations; simultaneously obtain the TWS prior field provided by the hydrological model, and generate the model-predicted deformation field through convolution with the Earth's elastic load Green's function.

[0161] S2. Decompose the time series of the vertical component coordinates of each GNSS station, and obtain the GNSS residual dominated by hydrological load after deducting non-hydrological terms; perform spatial interpolation on the GNSS residual to form a continuous residual deformation field.

[0162] S3. The continuous residual deformation field and the model-predicted deformation field are weighted and fused to obtain the initial fused deformation field, and the initial land water storage is obtained by reverse calculation.

[0163] S4. Using the rasterized terrestrial water storage as the state variable, establish a linear observation equation and construct a joint objective function that includes a data consistency term, a spatial smoothing term, and a priori shrinkage term.

[0164] S5. Using the initial land water storage as the initial field, the joint objective function is solved by residual recharge and gradient-type iterative algorithm. In each iteration, deformation prediction, residual calculation, residual space interpolation and state variable update are performed sequentially until the convergence condition is met, and the final land water storage is obtained.

[0165] S6. Based on linear covariance propagation or Monte Carlo perturbation, perform uncertainty assessment on the final terrestrial water storage; output the final terrestrial water storage and the corresponding uncertainty assessment results.

[0166] Furthermore, the logical instructions in the aforementioned memory can be implemented as software functional units and sold or used as independent products, and can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of the present invention, or the part that contributes to the prior art, or a part of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the steps of the methods described in the various embodiments of the present invention. The aforementioned storage medium includes various media capable of storing program code, such as USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks.

[0167] Example 4: This example provides a computer program product, which includes a computer program that can be stored on a non-transitory computer-readable storage medium. When the computer program is executed by a processor, the computer can perform a method for the collaborative inversion of land water storage by fusing hydrological models and sparse GNSS data. This method includes the following steps:

[0168] S1. Obtain the coordinate time series of each GNSS station in the study area, perform quality control, and extract the vertical component coordinate time series of the GNSS stations; simultaneously obtain the TWS prior field provided by the hydrological model, and generate the model-predicted deformation field through convolution with the Earth's elastic load Green's function.

[0169] S2. Decompose the time series of the vertical component coordinates of each GNSS station, and obtain the GNSS residual dominated by hydrological load after deducting non-hydrological terms; perform spatial interpolation on the GNSS residual to form a continuous residual deformation field.

[0170] S3. The continuous residual deformation field and the model-predicted deformation field are weighted and fused to obtain the initial fused deformation field, and the initial land water storage is obtained by reverse calculation.

[0171] S4. Using the rasterized terrestrial water storage as the state variable, establish a linear observation equation and construct a joint objective function that includes a data consistency term, a spatial smoothing term, and a priori shrinkage term.

[0172] S5. Using the initial land water storage as the initial field, the joint objective function is solved by residual recharge and gradient-type iterative algorithm. In each iteration, deformation prediction, residual calculation, residual space interpolation and state variable update are performed sequentially until the convergence condition is met, and the final land water storage is obtained.

[0173] S6. Based on linear covariance propagation or Monte Carlo perturbation, perform uncertainty assessment on the final terrestrial water storage; output the final terrestrial water storage and the corresponding uncertainty assessment results.

[0174] The device embodiments described above are merely illustrative. The units described as separate components may or may not be physically separate. The components shown as units may or may not be physical units; that is, they may be located in one place or distributed across multiple network units. Some or all of the modules can be selected to achieve the purpose of this embodiment according to actual needs. Those skilled in the art can understand and implement this without any creative effort.

[0175] Through the above description of the embodiments, those skilled in the art can clearly understand that each embodiment can be implemented by means of software plus necessary general-purpose hardware platforms, and of course, it can also be implemented by hardware. Based on this understanding, the above technical solutions, in essence or the part that contributes to the prior art, can be embodied in the form of a software product. This computer software product can be stored in a computer-readable storage medium, such as ROM / RAM, magnetic disk, optical disk, etc., and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute the methods described in the various embodiments or some parts of the embodiments.

[0176] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for collaborative inversion of land water storage by integrating hydrological models and sparse GNSS data, characterized in that, Includes the following steps: S1. Obtain the coordinate time series of each GNSS station in the study area, perform quality control, and extract the vertical component coordinate time series of the GNSS stations; simultaneously obtain the TWS prior field provided by the hydrological model, and generate the model-predicted deformation field through convolution with the Earth's elastic load Green's function. S2. Decompose the time series of the vertical component coordinates of each GNSS station, and obtain the GNSS residual dominated by hydrological load after deducting non-hydrological terms. Spatial interpolation is performed on the GNSS residuals to form a continuous residual deformation field; S3. The continuous residual deformation field and the model-predicted deformation field are weighted and fused to obtain the initial fused deformation field, and the initial land water storage is obtained by reverse calculation. S4. Using the rasterized terrestrial water storage as the state variable, establish a linear observation equation and construct a joint objective function that includes a data consistency term, a spatial smoothing term, and a priori shrinkage term. S5. Using the initial land water storage as the initial field, the joint objective function is solved by residual recharge and gradient-type iterative algorithm. In each iteration, deformation prediction, residual calculation, residual space interpolation and state variable update are performed sequentially until the convergence condition is met, and the final land water storage is obtained. S6. Based on linear covariance propagation or Monte Carlo perturbation, perform uncertainty assessment on the final terrestrial water storage; output the final terrestrial water storage and the corresponding uncertainty assessment results.

2. The method for synergistic inversion of land water storage by integrating hydrological models and sparse GNSS data according to claim 1, characterized in that, In step S1, the coordinate time series of each GNSS station in the study area are obtained, and quality control is performed to extract the vertical component coordinate time series of the GNSS stations; specifically: The coordinate time series of each GNSS station in the study area were obtained, and missing values ​​and outliers were filled, outliers were removed, and the reference frame and benchmark were made consistent. Extracting the time series of vertical component coordinates of GNSS stations ,in For station indexing, For the number of stations, For time epoch index, The total number of epochs; Based on the estimated observation noise level, adaptive weighting is applied according to the observation quality of each station at that epoch, forming the epoch weights of the diagonal weight matrix. .

3. The method for synergistic inversion of land water storage by integrating hydrological models and sparse GNSS data according to claim 2, characterized in that, In step S2, the time series of the vertical component coordinates of each GNSS station are decomposed, and after deducting non-hydrological terms, the GNSS residuals dominated by hydrological load are obtained; specifically: The coordinate time series of each GNSS station is decomposed, and a trajectory model is established, represented as follows: ； ； ； ； The GNSS residuals, dominated by hydrological load, are obtained by subtracting non-hydrological terms from the time series of vertical component coordinates of GNSS stations, and are expressed as follows: ； in, This represents the trend term indicating long-term deformation and baseline drift; For periods representing six months or a year; A step term caused by changes in equipment or environment; The residuals after deducting non-hydrological items are dominated by deformation caused by hydrological loads; To observe random noise, its parameters are used to construct epoch weights. ; It is the Heaviside step function; The constant term represents the initial position; Expressed as the rate of long-term linear deformation; The first The sine and cosine coefficients of each periodic term; Angular frequency; Indicates the first The amplitude of a step jump; The moment when the step jump occurs; Indicates the periodic component index. Indicates its total number; Represented as the index of the transition event. The set of its values.

4. The method for synergistic inversion of land water storage by integrating hydrological models and sparse GNSS data according to claim 2, characterized in that, In step S2, spatial interpolation is performed on the GNSS residuals to form a continuous residual deformation field; specifically: Ordinary Kriging interpolation is used for epochs. GNSS residuals Spatial interpolation is performed to obtain the continuous residual deformation field. , represented as: ； in, Representing geographic coordinates; using an exponential / spherical harmonic Gaussian model to analyze the sample semivariogram. By fitting the data, the corresponding covariance function can be obtained. Based on this, a Kriging covariance matrix is ​​constructed. Covariance vector of the interpolation point ; For Kriging weights, for A set; These are Lagrange multipliers used to ensure that the sum of the Kriging weights is 1.

5. The method for synergistic inversion of land water storage by integrating hydrological models and sparse GNSS data according to claim 1, characterized in that, Step S1 involves synchronously acquiring the TWS prior field provided by the hydrological model and generating the model-predicted deformation field through convolution with the Earth's elastic load Green's function; specifically: Constructing a vertical load Green's function based on the load Love number. Pre-calculation of design matrix ; Obtain changes in terrestrial water storage synthesized from hydrological models And through the Earth's elastic load Green's function Convolution is converted into the model's predicted deformation field, expressed as: ； in, Represents geographic coordinates; The geographic coordinates of the integration point. The distance is spherical. For the density of water, It is the acceleration due to gravity. The area of ​​the face unit; Using design matrix Discretize the model-predicted deformation field and represent it as follows: ； in, For the first Each grid area.

6. The method for synergistic inversion of land water storage by integrating hydrological models and sparse GNSS data according to claim 1, characterized in that, In step S3, the continuous residual deformation field and the model-predicted deformation field are weighted and fused to obtain the initial fused deformation field, and the initial terrestrial water storage is then calculated in reverse; specifically: continuous residual deformation field With model prediction of deformation field Weighted fusion yields the initial fused deformation field, expressed as: ； in, Represents geographic coordinates; To integrate the weights, an adaptive selection is made using L-curve cross-validation, generalized cross-validation, or K-fold cross-validation. The larger the value, the more trust is placed in the model's prior knowledge; Then, the initial terrestrial water storage is determined by inverse calculation based on the initial fusion deformation field or by combining it with the prior TWS convexity. .

7. The method for synergistic inversion of land water storage by integrating hydrological models and sparse GNSS data according to claim 1, characterized in that, Step S4 describes establishing a linear observation equation using rasterized terrestrial water storage as the state variable, and constructing a joint objective function that includes a data consistency term, a spatial smoothing term, and a priori shrinkage term; specifically: The gridded terrestrial water storage to be sought For the unknown state variables, establish the observation equation, expressed as: ； in, For the residuals after deducting non-hydrological items, To observe random noise; To design a matrix, design a matrix elements for: ； in, The distance is spherical. Represents geographic coordinates; The geographic coordinates of the integration point. For the density of water, It is the acceleration due to gravity. For the first The area of ​​each grid cell The Green's function for vertical load is constructed based on the load Love number; Define the joint objective function: ； in, The epoch weights are formed by adaptively weighting the diagonal weight matrix based on the estimated observation noise level and the observation quality of each station at that epoch. This is a data consistency term that minimizes the weighted residuals. For spatial smoothness terms, use the spatial smoothness regularization operator. Suppressing non-physical high frequencies; For the prior reduction term, with the prior weight matrix The constraint solution does not deviate excessively Hyperparameters The method is adaptively determined using L-curve, generalized cross-validation, or K-fold cross-validation. The gradient of the joint objective function is: ； The iteration uses Landweber gradient descent, expressed as: ； Step length satisfy .

8. The method for synergistic inversion of land water storage by integrating hydrological models and sparse GNSS data according to claim 7, characterized in that, Step S5 is as follows: Based on initial terrestrial water storage For the initial field, in the first... In the next iteration, according to Computation site predicts deformation ; Calculate the observation residuals ; For observation residuals Spatial interpolation is performed to obtain the residual field. ; Based on the spatial smoothness coefficient of the residual field superposition With prior shrinkage coefficient The updated hydrological signal estimate is expressed as: ； in To update the step size; When the relative changes of the objective function and the residual RMS are both below the threshold, it is expressed as: ； ； If the preset maximum number of iterations is reached, the iteration terminates, and the final terrestrial water storage is obtained. .

9. A non-transitory computer-readable storage medium, characterized in that, It stores computer instructions that enable the computer to execute the land water storage co-inversion method according to any one of claims 1-8, which integrates hydrological models and sparse GNSS data.

10. An electronic device, characterized in that, include: The system includes a processor, a communication interface, a memory, and a communication bus. The processor, communication interface, and memory communicate with each other via the communication bus. The processor calls logical instructions from the memory to execute the land water storage collaborative inversion method based on any one of claims 1-8, which integrates hydrological models and sparse GNSS data.

Images