Land water storage change inversion method considering scale factor and intelligent terminal

Abstract

The application discloses a land water storage change inversion method considering a proportional factor, and the method comprises the following steps: preprocessing a GNSS three-dimensional coordinate time sequence; obtaining GNSS observation displacement caused by hydrological load according to the processed GNSS three-dimensional coordinate time sequence; obtaining GRACE model prediction displacement caused by hydrological load according to equivalent water height data of a GRACE-Mascon product; obtaining a proportional factor of a GNSS real observation station position according to the GNSS observation displacement and the GRACE model prediction displacement; interpolating to obtain a proportional factor of a virtual station according to the proportional factor of the GNSS real observation station position; and establishing a joint inversion ΔTWS model according to observation values of the GNSS real observation station and the virtual station and the proportional factor of the virtual station, and then performing inversion. The application relies on the constraint of the proportional factor from the virtual station, so that the system difference between observation means can be reduced to the greatest extent, a more accurate ΔTWS is provided, and the application has important significance for hydrological application and water resource management of a sparse GNSS station coverage area.

Description

Technical Field

[0001] This invention relates to the field of hydrological analysis, specifically to a method for inverting changes in terrestrial water storage that takes into account proportional factors, and an intelligent terminal thereof. Background Technology

[0002] ΔTWS (ΔTWS) refers to the sum of changes in surface water, groundwater, soil moisture, rainfall, snow, and canopy water. Accurate estimation of the spatiotemporal evolution of ΔTWS is crucial for studying the hydrological cycle and water resource management. In particular, against the backdrop of current global warming, the risk of extreme hydrological events (floods and droughts) has significantly increased, not only causing serious water and food security problems but also potentially leading to economic risks and regional conflicts. Monitoring ΔTWS helps analyze the frequency and intensity changes of extreme meteorological and hydrological events, thereby providing possible response and prevention strategies. Traditional point-based observation methods (such as hydrological stations and well logging) can independently measure different water components, but due to their poor spatiotemporal coverage, they cannot obtain accurate ΔTWS on a global or regional scale. Hydrological models are widely used to estimate global or regional ΔTWS to address and mitigate water resource pressures. However, hydrological models do not incorporate observations of all hydrological components, thus failing to distinguish the impacts of climate change and human activities on water resources.

[0003] Currently, with the development of various geodetic methods, dense Global Navigation Satellite System (GNSS) networks and the Gravity Recovery and Climate Experiment (GRACE) provide reliable independent observations of ΔTWS. GRACE monitors the redistribution of ΔTWS mass by observing changes in the Earth's gravitational field, and can display hydrological cycle processes in coarse space (>300 km), providing important information for large-scale drought and flood monitoring and long-term groundwater changes. However, due to the coarse spatial resolution of GRACE itself, there is significant uncertainty in inferring ΔTWS using GRACE over shorter spatial distances. GNSS can accurately monitor load deformation caused by the redistribution of surface mass in local areas (0-200 km) and infer ΔTWS based on the elastic load theory of the solid Earth. With the improvement of GNSS positioning algorithm accuracy and signal processing algorithm performance, GNSS is increasingly used in research on extreme hydrological drought events, extreme precipitation events, glacier mass loss estimation, and snowmelt equivalent estimation. However, the reliability of GNSS inversion of ΔTWS depends on the spatial distribution and density of GNSS stations. Inversion performance may be affected in areas with sparse GNSS station coverage, which greatly limits the applicability and accuracy of the method.

[0004] Recent research has begun to attempt to combine GNSS and GRACE methods to study ΔTWS, thereby generating more reliable ΔTWS with higher spatiotemporal resolution than using either method alone. However, due to differences in modeling methods and observation techniques between the two geodetic instruments, how to reasonably address the systematic biases related to spatial location and distance remains an unresolved issue. Failure to properly handle these systematic differences in the joint inversion model means introducing unrealistic observations into the model, which significantly impacts the accuracy and reliability of the inversion results.

[0005] Therefore, existing technologies still need to be improved and developed. Summary of the Invention

[0006] The technical problem to be solved by the present invention is to provide a method for inverting changes in land water storage that takes into account the scaling factor, in order to address the above-mentioned deficiencies of the prior art. This method aims to solve the problem that the existing joint inversion model cannot properly handle system differences, which affects the accuracy and reliability of the inversion results.

[0007] The technical solution adopted by this invention to solve the technical problem is as follows:

[0008] In a first aspect, the present invention provides a method for inverting changes in terrestrial water storage that takes into account a scaling factor, wherein the method includes:

[0009] Acquire GNSS three-dimensional coordinate time series, actual GNSS observation station locations, and equivalent water height data from GRACE-Mascon products;

[0010] The GNSS three-dimensional coordinate time series is preprocessed to obtain the processed GNSS three-dimensional coordinate time series;

[0011] Based on the processed GNSS three-dimensional coordinate time series, the GNSS observation displacement caused by hydrological load is obtained;

[0012] Based on the equivalent water height data of the GRACE-Mascon product, the displacement predicted by the GRACE model caused by hydrological load is obtained;

[0013] Based on the GNSS observed displacement and the displacement predicted by the GRACE model, the scaling factor of the actual GNSS observation station position is obtained.

[0014] Based on the scaling factor of the actual GNSS observation station location, the scaling factor of the virtual station is interpolated.

[0015] A joint inversion ΔTWS model was established based on the observations from the actual GNSS observation station and the virtual station, as well as the scaling factor of the virtual station.

[0016] The changes in terrestrial water storage are obtained based on the inversion ΔTWS model.

[0017] In one implementation, the preprocessing of the GNSS three-dimensional coordinate time series to obtain the processed GNSS three-dimensional coordinate time series includes:

[0018] The rectangular coordinate system of the GNSS three-dimensional coordinate time series is converted to the station-centered coordinate system to obtain the NEU coordinate time series of the GNSS.

[0019] The NEU coordinate time series of the GNSS was subjected to gross error removal using the interquartile range method to obtain the NEU coordinate time series after removal.

[0020] A GNSS three-dimensional coordinate time series model is established, and the NEU coordinate time series after elimination is processed according to the GNSS three-dimensional coordinate time series model to obtain the processed GNSS three-dimensional coordinate time series.

[0021] In one implementation, obtaining the GNSS observation displacement caused by hydrological load based on the processed GNSS three-dimensional coordinate time series includes:

[0022] Solve the posterior probability distribution of the unknown parameters in the variational Bayesian independent component analysis framework to obtain the source signal probability distribution function of the GNSS three-dimensional coordinate time series;

[0023] The posterior probability of the hidden variables is calculated based on the probability distribution function of the source signal to obtain the estimated value of the unknown parameters;

[0024] Based on the estimated values ​​of the unknown parameters, an independent component analysis model is obtained;

[0025] The GNSS observation displacement caused by hydrological load is obtained through the independent component analysis model.

[0026] In one implementation, obtaining the GRACE model-predicted displacement caused by hydrological load based on the equivalent water height data of the GRACE-Mascon product includes:

[0027] Calculate the radial and horizontal displacements caused by a unit hydrological mass load on the Earth's surface;

[0028] The Green's functions for vertical and horizontal loads are obtained based on the radial and horizontal displacements.

[0029] The displacement predicted by the GRACE model is obtained based on the Green's functions of the vertical and horizontal loads and the equivalent water height data of the GRACE-Mascon product.

[0030] In one implementation, obtaining the scaling factor of the actual GNSS observation station based on the GNSS observed displacement and the displacement predicted by the GRACE model includes:

[0031] Calculate the ratio between the GNSS observed displacement and the GRACE model predicted displacement, and use the ratio as the scaling factor for the actual GNSS station.

[0032] In one implementation, the step of interpolating the virtual station's scaling factor based on the actual GNSS observation station location includes:

[0033] Generate Thiessen polygons based on the actual locations of the GNSS observation stations;

[0034] The virtual station is obtained based on the Thiessen polygon and the buffer.

[0035] Based on the Pearson correlation coefficient between the GNSS observed displacement and the displacement predicted by the GRACE model, the scaling factor of the virtual station is obtained as follows:

[0036]

[0037] Where i is the number of GNSS stations adjacent to the virtual station; ρ i PCC represents the displacement predicted by GNSS observations and the GRACE model for the i-th station; β i This represents the scaling factor for the i-th station.

[0038] In one implementation, the joint inversion ΔTWS model established based on the observations from the actual GNSS observation station and the virtual station, as well as the scaling factor of the virtual station, is as follows:

[0039]

[0040] Among them, A v and A u b represents the vertical and horizontal Green's function matrices generated based on actual GNSS observation stations; n ,b e ,b u These are the processed GNSS three-dimensional coordinate time series in the GNSS station-centered coordinate system (ENU); b G represents the observations from the virtual station; W is the weighting matrix of the GNSS observations, defined based on the standard error of the GNSS observations; U is the weighting matrix of the GRACE observations; the elements in the weighting matrix W... σ e , σ n , σ u These represent the uncertainties in the station-centered coordinate system ENU; and the U elements of the weighting matrix. i represents the number of actual GNSS observation stations; β j is the scaling factor for virtual stations; j represents the number of virtual stations.

[0041] Secondly, embodiments of the present invention also provide a land water storage change inversion device that takes into account a scaling factor, wherein the device includes:

[0042] The data acquisition module is used to acquire GNSS three-dimensional coordinate time series, GNSS actual observation station locations, and equivalent water height data of GRACE-Mascon products;

[0043] The GNSS three-dimensional coordinate time series acquisition module is used to preprocess the GNSS three-dimensional coordinate time series to obtain the processed GNSS three-dimensional coordinate time series.

[0044] The GNSS observation displacement acquisition module is used to obtain the GNSS observation displacement caused by hydrological load based on the processed GNSS three-dimensional coordinate time series.

[0045] The GRACE model predicted displacement acquisition module is used to obtain the GRACE model predicted displacement caused by hydrological load based on the equivalent water height data of the GRACE-Mascon product.

[0046] The module for obtaining the scale factor of the real GNSS observation station is used to obtain the scale factor of the real GNSS observation station based on the GNSS observation displacement and the displacement predicted by the GRACE model.

[0047] The virtual station scaling factor acquisition module is used to interpolate the scaling factor of the virtual station based on the location of the actual GNSS observation station.

[0048] The model building module is used to build a joint inversion ΔTWS model based on the observations of the real GNSS observation station and the virtual station, as well as the scaling factor of the virtual station.

[0049] The deduction module is used to obtain the change in terrestrial water storage based on the inversion ΔTWS model.

[0050] Thirdly, embodiments of the present invention also provide a smart terminal, wherein the smart terminal includes a memory, a processor, and a land water storage change inversion program considering the scaling factor stored in the memory and executable on the processor. When the processor executes the land water storage change inversion program considering the scaling factor, it implements the steps of the land water storage change inversion method considering the scaling factor as described in any of the above.

[0051] Fourthly, embodiments of the present invention also provide a computer-readable storage medium, wherein the computer-readable storage medium stores a land water storage change inversion program taking into account a scaling factor, and when the land water storage change inversion program taking into account a scaling factor is executed by a processor, it implements the steps of the land water storage change inversion method taking into account a scaling factor as described in any of the above claims.

[0052] Beneficial Effects: Compared with existing technologies, this invention provides a method for inverting land water storage changes that takes into account scaling factors. First, it acquires GNSS three-dimensional coordinate time series, actual GNSS station locations, and equivalent water height data from GRACE-Mascon products. This invention provides comprehensive observations of ΔTWS (ΔTWS), rather than a single hydrological component. Compared to traditional point-based and hydrological model-based methods, it more comprehensively reflects regional ΔTWS information, helping to distinguish the impacts of climate change and human activities on total water resources. Then, based on the processed GNSS three-dimensional coordinate time series, it obtains the GNSS observation displacement caused by hydrological load to recover the vertical crustal deformation caused by hydrological movement, minimizing pollution from non-hydrological deformation sources in the GNSS time series. Secondly, it extracts model-predicted displacements from GRACE products using the load Green's function and determines the scaling factor between the GNSS observation displacements and GRACE model-predicted displacements using least squares. This introduces new GRACE observation data without changing the form of traditional methods, offering higher spatial resolution and reliability compared to single data sources, and capturing detailed features that cannot be captured using single observations. By correcting for the systematic bias caused by the different spatial sensitivities between GNSS and GRACE observations, the differences between the two observation methods are minimized. Finally, by incorporating the displacement predicted by the GRACE model with a scaled-factor adjusted into the GNSS inversion model to obtain ΔTWS, a more reliable ΔTWS is generated. This solves the problem of unrealistic ΔTWS estimation caused by low GNSS station density in local areas. The addition of GRACE observations provides important data constraints for areas with sparse GNSS station distribution, enabling accurate estimation of ΔTWS without adding new ground observations. Attached Figure Description

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

[0054] Figure 1This is a schematic diagram of the process for inverting changes in terrestrial water storage considering the scaling factor, provided in an embodiment of the present invention.

[0055] Figure 2 This is a flowchart illustrating the framework of the method for inverting changes in terrestrial water storage that takes into account the scaling factor, as provided in this embodiment of the invention.

[0056] Figure 3 This is a schematic diagram of the time function (V) and spatial response (U) derived from vbICA according to an embodiment of the present invention.

[0057] Figure 4 This is a schematic diagram of displacement and scaling factor provided in an embodiment of the present invention.

[0058] Figure 5 This is a schematic diagram showing the virtual station site selection scheme provided in an embodiment of the present invention.

[0059] Figure 6 This is a schematic diagram provided by an embodiment of the present invention, showing the spatial distribution of PCC and the reciprocal of the scaling factor between GNSS vertical crustal displacement and GRACE model predicted displacement.

[0060] Figure 7 This is a schematic diagram of the annual amplitude of the equivalent water level provided in an embodiment of the present invention.

[0061] Figure 8 This is a schematic diagram showing the equivalent water height changes in the Jinsha River, Lancang River, and Nanpan River basins in Yunnan Province, provided by an embodiment of the present invention.

[0062] Figure 9 This is a schematic diagram of the principle of the land water storage change inversion device that takes into account the scaling factor provided in the embodiment of the present invention.

[0063] Figure 10 This is a block diagram illustrating the internal structure of a smart terminal provided in an embodiment of the present invention. Detailed Implementation

[0064] To make the objectives, technical solutions, and effects of this invention clearer and more explicit, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and are not intended to limit the invention.

[0065] Those skilled in the art will understand that, unless specifically stated otherwise, the singular forms “a,” “an,” “the,” and “the” used herein may also include the plural forms. It should be further understood that the term “comprising” as used in this specification means the presence of the stated features, integers, steps, operations, elements, and / or components, but does not exclude the presence or addition of one or more other features, integers, steps, operations, elements, components, and / or groups thereof. It should be understood that when we say an element is “connected” or “coupled” to another element, it can be directly connected or coupled to the other element, or there may be intermediate elements. Furthermore, “connected” or “coupled” as used herein can include wireless connections or wireless coupling. The term “and / or” as used herein includes all or any units and all combinations of one or more associated listed items.

[0066] It will be understood by those skilled in the art that, unless otherwise defined, all terms used herein (including technical and scientific terms) have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains. It should also be understood that terms such as those defined in general dictionaries should be understood to have the same meaning as in the context of the prior art, and should not be interpreted in an idealized or overly formal sense unless specifically defined as herein.

[0067] Monitoring changes in terrestrial water storage (ΔTWS) helps analyze the frequency and intensity of extreme meteorological and hydrological events, thus providing possible response and prevention strategies. Traditional point-based observation methods (such as hydrological stations and well logging) can independently measure different water components, but due to their poor spatiotemporal coverage, they cannot obtain accurate ΔTWS globally or regionally. Traditional GNSS inversion methods for ΔTWS are suitable for areas with dense station coverage, and the inversion results are affected by the existing station distribution and density, therefore, accurate ΔTWS cannot be obtained in scenarios with sparse GNSS station coverage.

[0068] To address the aforementioned technical issues, this invention proposes using the GRACE model with scale factor correction to predict displacement as an additional constraint introduced into the inversion model. This corrects for unrealistic inversion results caused by insufficient ground observations. Relying on the scale factor constraint between the two observation methods, it can minimize the system differences between the observation methods and provide more accurate ΔTWS. This is of great significance for hydrological applications and water resource management in areas with sparse GNSS station coverage.

[0069] Exemplary methods

[0070] This embodiment provides a method for inverting changes in terrestrial water storage that takes into account the scaling factor.

[0071] like Figure 1As shown, the method includes the following steps:

[0072] Step S100: Obtain GNSS three-dimensional coordinate time series, GNSS actual observation station locations, and equivalent water height data from GRACE-Mascon products;

[0073] Specifically, the Global Navigation Satellite System (GNSS), also known as a global satellite navigation system, is a space-based radio navigation and positioning system that provides users with all-weather 3D coordinates, velocity, and time information at any location on the Earth's surface or in near-Earth space. It comprises one or more satellite constellations and augmentation systems required to support specific tasks. Utilizing carrier phase differential technology (RTK) in GNSS, its measurement accuracy can reach the centimeter level. Compared to traditional manual surveying, it offers unparalleled advantages such as high accuracy, ease of operation, portable equipment, all-weather operation, and no line-of-sight between measurement points. It has been widely applied in fields such as geodesy, crustal movement, resource exploration, and cadastral surveying.

[0074] Thanks to the rapid development of satellite gravity observation technology, especially since the successful launch of the Gravity Recovery and Climate Experiment (GRACE) satellite in 2002, Earth's gravity field models have been fundamentally and significantly improved. GRACE provides an unprecedentedly high-precision, time-varying gravity field covering the entire globe, offering entirely new observational data and research methods for scientific research in multiple fields of Earth science.

[0075] In this embodiment, GNSS and GRACE are combined to study ΔTWS, in order to generate ΔTWS that is more reliable and has higher spatiotemporal resolution than using either method alone. This requires acquiring GNSS and GRACE data, specifically including GNSS three-dimensional coordinate time series, actual GNSS observation station locations, and equivalent water height data from GRACE-Mascon products.

[0076] Step S200: Preprocess the GNSS three-dimensional coordinate time series to obtain the processed GNSS three-dimensional coordinate time series;

[0077] In one implementation, step S200 of this embodiment includes the following steps:

[0078] Step S201: Convert the rectangular coordinate system of the GNSS three-dimensional coordinate time series to the station-centered coordinate system to obtain the NEU coordinate time series of the GNSS.

[0079] Specifically, the XYZ coordinates obtained via GNSS are in a spatial rectangular coordinate system (WGS-84). In time series analysis, these need to be converted to the station-centered coordinate system (ENU). The conversion formula is as follows:

[0080]

[0081] Where e, n, and u are the coordinates of the three components in the ENU coordinate system; x, y, and z are the coordinates in the corresponding WGS-84 coordinate system; x0, y0, and z0 are the WGS-84 coordinates of the initial position of the station; and M is the transformation matrix, the specific expression of which is as follows:

[0082]

[0083] Where λ0 and φ0 represent geodetic longitude and geodetic latitude in the NEU coordinate system.

[0084] Step S202: Use the interquartile range method to remove gross errors from the NEU coordinate time series of the GNSS, and obtain the NEU coordinate time series after removal.

[0085] Specifically, since the station position calculation is affected by various factors such as multipath effects, ionospheric delay, and orbit modeling, gross errors may occur. This invention uses the Interquartile Range (IQR) method to detect and eliminate these gross errors. The basic principle of IQR is as follows:

[0086] IQR = Q2 - Q1

[0087] The time series are arranged in ascending order. Q1 and Q2 represent the lower quantile and upper quantile values ​​closest to 1 / 4 and 3 / 4, respectively. The detection interval for gross errors is [Q1-1.5·IQR, Q2+1.5·IQR]. Values ​​outside this interval are considered outliers.

[0088] Step S203: Establish a GNSS three-dimensional coordinate time series model, and organize the removed NEU coordinate time series according to the GNSS three-dimensional coordinate time series model to obtain the processed GNSS three-dimensional coordinate time series.

[0089] Specifically, the purpose of modeling the GNSS three-dimensional coordinate time series is to remove noise, correct for step changes, and extract useful signals. This embodiment uses the following formula to fit the time series:

[0090]

[0091] Among them, t iRepresents the year; a represents the initial coordinates; b represents the linear motion rate; c and d represent the annual motion amplitude; e and f represent the semi-annual motion amplitude; g j This indicates a step jump caused by factors such as instrument replacement or earthquakes; H(t) i -T gj ) represents the step function; v i The residuals are used to fit the model.

[0092] Then, the NEU coordinate time series after removal is sorted out.

[0093] Assuming a regional network consisting of N continuous GNSS stations, the preprocessed GNSS coordinate time series contains three components: ENU, M = 3N (the total number of time series), and T (the total number of epochs). Then, the processed data matrix X is:

[0094]

[0095] Where M is the number of epochs, T is the number of observation stations, and x represents the observed value.

[0096] For example, in this embodiment, the GNSS coordinate time series is preprocessed. Non-tidal atmospheric load (NTAL) and non-tidal ocean load (NTOL) deformation products under the CF reference frame provided by the Potsdam Geoscience Center (GFZ) are used to eliminate motions associated with non-tidal atmospheric and ocean loads in the coordinate time series. Then, the least squares method is used to fit the linear trend, annual and semi-annual cycle motions, and steps caused by earthquakes and instrument replacements in the three-dimensional coordinate time series. Finally, common-mode errors, long-term linear trends, and steps are removed from the original time series, and the remaining GNSS time series is used for further processing.

[0097] Step S300: Based on the processed GNSS three-dimensional coordinate time series, obtain the GNSS observation displacement caused by hydrological load;

[0098] In one implementation, step S300 of this embodiment includes the following steps:

[0099] Step S301: Solve for the posterior probability distribution of the unknown parameters in the variational Bayesian independent component analysis framework to obtain the source signal probability distribution function of the GNSS three-dimensional coordinate time series;

[0100] Specifically, we solve for the posterior probability distribution of the unknown parameters in the variational Bayesian independent component analysis framework.

[0101] GNSS coordinate time series data is obtained by a linear combination of several source signals (tectonic deformation, non-tectonic deformation caused by geophysical effects, noise, etc.), and can be expressed as:

[0102] X M*T =A M*L S L*T +ε

[0103] Where X is the observation matrix; A is the mixture matrix; S is the source signal matrix; and ε is the noise matrix. The purpose of variational Bayesian independent component analysis is to obtain the source signal matrix S and the mixture matrix A using only the observation matrix X, when both the source signal matrix and the mixture matrix are unknown. In Bayesian networks, Bayesian inference is used to calculate the posterior probability distribution of unknown parameters, while independent component analysis learns the model's hidden variables (i.e., the source signals) under the guidance of the objective function. Therefore, the two are consistent within the Bayesian framework. The source signals are independent, and the probability distribution function of the source matrix S can be written as:

[0104]

[0105] Where P represents the probability distribution function; S i L represents the source signal; L represents the number of source signals.

[0106] Step S302: Calculate the posterior probability of the hidden variable based on the probability distribution function of the source signal to obtain the estimated value of the unknown parameter;

[0107] In independent component analysis models, this can theoretically be achieved by calculating the posterior probabilities of the hidden variables, as shown by the following formula:

[0108]

[0109] Where M represents a specific model chosen; p(W|M) is the signal source model, and p(X|M) is the normalization factor, usually called the marginal probability of model M. Assuming all unknown parameters in model M are represented by the vector W = (ε, A, S, q, θ), where ε, A, S, q, and θ represent the model background noise, mixing matrix, simulated signal source, Gaussian mixture components, and model parameters, respectively, given the observed signal X, the posterior probability of parameter W can be expressed as:

[0110]

[0111] Step S303: Obtain the independent component analysis model based on the estimated values ​​of the unknown parameters;

[0112] Step S304: Obtain the GNSS observation displacement caused by hydrological load through the independent component analysis model.

[0113] Choudrey et al. (2002) transformed equation (9) into a calculation of the negative free energy F of the model, which can be expressed as follows:

[0114] F[X]=[ln(p(X,W)] p′(W) +H[X] (10)

[0116] Here, p′(W) is an approximate estimate of the posterior probability p(W|X,M). It can be further expressed in the form of factorization:

[0117]

[0118] Miskin et al. (2000) proved that when the estimated posterior probability p′(W) in equation (10) is closest to the true posterior probability p(X|M), F[X] reaches an extreme value. Then, by taking the partial derivatives of p′(W) with respect to each parameter, the estimated values ​​of the unknown parameters in model M can be obtained. The source signal can be obtained by taking the estimated values ​​of the position parameters.

[0119] For example, after obtaining the preprocessed coordinate time series, it is necessary to extract the hydrological load deformation from the accurate time series. Figure 3The normalized time response, power spectral density, and spatial response of five ICs are shown, i.e., the time function (V) and spatial response (U) derived from vbICA. The time function and corresponding power spectral density plot for each IC are displayed at the top. The time function is normalized between 0 and 1, and the weight of each IC is shown in the lower left corner. The spatial response is shown below, with red arrows representing the horizontal response and color scales representing the vertical response. The spatial response represents the direction of station movement, and its value is proportional to the importance of a specific component for each station. IC1 represents a common-mode annual signal, and its spatial response shows a consistent movement trend across all stations, primarily representing vertical surface movement caused by the hydrological cycle. Similarly, IC2 and IC4 are also common-mode signals, meaning the stations exhibit periodic movement in the same direction. IC2 is mainly in the NW-SE direction, and IC4 is mainly in the NE-SW direction. The power spectral density of IC2's time response shows a combination of a high-power, long-period signal and lower-power annual and semi-annual periodic signals, and over 90% of its stations exhibit consistent spatial responses. Therefore, IC2 is attributed to local small-scale hydrological signals because the horizontal component is more sensitive to short-wavelength heterogeneity compared to vertical measurements. The power spectral density of the time response of IC4 peaks near the period of the annual intersection error signal (351.5 and 175.9 days), thus attributing this component to the annual intersection error in GNSS measurements. The power spectral densities of IC3 and IC5 both exhibit certain annual and low-frequency signals. IC3 mainly exhibits a high-power low-frequency signal, while IC5 is mainly an annual signal. Therefore, crosstalk may exist between IC3 and IC5 (leakage from one IC to the other), but this does not affect the reconstruction of the hydrological load deformation signal. The spatial responses of IC3 and IC5 both exhibit NE-SW anisotropic patterns, and their combination records small-scale hydrological information of the watersheds in southwestern and northeastern Yunnan.

[0120] Step S400: Based on the equivalent water height data of the GRACE-Mascon product, obtain the displacement predicted by the GRACE model caused by hydrological load;

[0121] In one implementation, step S400 of this embodiment includes the following steps:

[0122] Step S401: Calculate the radial and horizontal displacements caused by a unit hydrological mass load on the Earth's surface;

[0123] Specifically, determine the radial displacement caused by the mass load on the Earth's surface.

[0124] Assuming a uniform mass with a uniform angular radius is placed on the spherically symmetric, self-gravitating, elastic, and isotropic surface of the Earth, the resulting radial displacement on the Earth's surface is given by the following equation:

[0125] ′

[0126] u n =U n (r)P n (cosθ)e iωt =ah n P n (cosθ)e iωt

[0127] Among them, u n It is radial displacement; U n (r) is the radial coefficient of the spherical harmonic expansion; P n (cosθ)

[0128] ′

[0129] For an nth-order Legendre polynomial; h n e is the load Love number; iωt This represents the time evolution of the applied load. 'a' is the Earth's radius; since the boundary conditions are based on the Earth's mass, the above equation can be rewritten as:

[0130]

[0131] Where, m E Let m be the mass of the Earth; m′ be a fixed mass load. Taking this as the load per unit mass, we have:

[0132]

[0133] Step S402: Obtain the Green's functions for the vertical and horizontal loads based on the radial and horizontal displacements;

[0134] Specifically, the Green's functions for vertical and horizontal loads are calculated. The vertical displacement caused by a 1 kg load applied to the Earth's surface can be expressed as:

[0135]

[0136] Similarly, the horizontal displacement of the earth's surface caused by applying a 1 kg load can be expressed as:

[0137]

[0138] It is worth noting that horizontal displacement does not apply to cases where the degree is 0, therefore the degree of v starts from n=1.

[0139] Step S403: Based on the Green's functions of the vertical and horizontal loads and the equivalent water height data of the GRACE-Mascon product, obtain the displacement predicted by the GRACE model.

[0140] Specifically, the mass load deformation in the GRACE-Mascon product is calculated. Based on the load Green's function and the equivalent water height provided by the GRACE-Mascon product, the GRACE model predicted displacement can be calculated using the formula Gm = d, where G is the Green's function matrix generated based on the grid and GNSS stations in the study area; m is the equivalent water height of the study area provided by the GRACE-Mascon product; and d is the GRACE model predicted displacement for each station to be calculated.

[0141] Step S500: Based on the GNSS observed displacement and the displacement predicted by the GRACE model, obtain the scaling factor of the actual GNSS observation station position;

[0142] In one implementation, step S500 of this embodiment includes the following steps:

[0143] Step S501: Calculate the ratio between the GNSS observed displacement and the GRACE model predicted displacement, and use the ratio as the scaling factor for the actual GNSS observation station position.

[0144] Specifically, determining the scaling factor between the displacement d predicted by the GRACE model and the displacement time series x observed by GNSS, i.e., determining the scaling factor between the GNSS observed displacement and the displacement predicted by the GRACE model, can be achieved through a linear model:

[0145] x = β·d + b

[0146] Where β represents the scaling factor or scaling coefficient; b is the translation coefficient.

[0147] For example, the load Green's function is calculated based on the PREM model, and the equivalent water height in the GRACE-Mascon product is converted into model-predicted displacements of GNSS station locations. Then, the vertical surface deformation observed by all GNSS stations is compiled into a monthly average series and aligned with the time series of displacements predicted by the GRACE model. Finally, the scaling factors for both are calculated based on the least squares principle. Figure 4 This paper presents a comparison between GRACE model-predicted displacements and GNSS-observed displacements before and after scaling by a scaling factor for 15 selected stations. The figures show the vertical crustal displacements observed at the selected GNSS stations (black line), the GRACE model-predicted displacements (green line), and the scaled GRACE model-predicted displacements (red dashed line). Labels below each subplot indicate the scaling factor (1 / β) between the GNSS and GRACE-derived displacements and their PCC (Position Controlled Cross-Section). The results demonstrate that the scaling factor-based method effectively reduces the discrepancies between GRACE model-predicted displacements and GNSS-observed displacements in local areas.

[0148] Step S600: Based on the scaling factor of the actual GNSS observation station location, interpolate to obtain the scaling factor of the virtual station;

[0149] In one implementation, step S600 of this embodiment includes the following steps:

[0150] Step S601: Generate Thiessen polygons based on the actual locations of the GNSS observation stations;

[0151] Step S602: Obtain the virtual station based on the Thiessen polygon and the buffer.

[0152] Specifically, due to the limitations of the actual spatial resolution of the GRACE-Mascon product, the spatial resolution of the selected virtual station should not exceed the actual spatial resolution of the GRACE-Mascon product (300km). Furthermore, the number of selected virtual GNSS stations should not be excessive, so that GNSS observations can dominate the inversion process. Based on these two points, this invention proposes using Thiessen polygons and buffer zones to determine the location of virtual stations. Thiessen polygons, also known as von Rooi diagrams, are continuous polygons composed of the perpendicular bisectors of a set of adjacent points, and are currently widely used in station location selection. Thiessen polygons have three characteristics: 1. There is only one discrete point within a Thiessen polygon. 2. Points on the edges of a Thiessen polygon are equidistant from the discrete points on their two sides. 3. Points within a Thiessen polygon are closest to their corresponding discrete points (station locations). Based on these characteristics, the location selection of virtual stations follows these steps:

[0153] 1. Generate Thiessen polygons based on the spatial distribution of existing stations to determine the radiation area of ​​each station;

[0154] 2. To ensure the reliability of the inversion, the distance between the grid of interest and the nearest station should not exceed 100km.

[0155] 3. Vertices of the Thiessen polygons not covered by the buffer are selected first as virtual station locations, and it is observed whether the buffer covers all the grids of interest.

[0156] Step S603: Based on the Pearson correlation coefficient between the GNSS observed displacement and the displacement predicted by the GRACE model, the scaling factor of the virtual station is obtained.

[0157]

[0158] Where i is the number of GNSS stations adjacent to the virtual station; ρ i PCC represents the displacement predicted by GNSS observations and the GRACE model for the i-th station; β i This represents the scaling factor for the i-th station.

[0159] Specifically, the scaling factor of the virtual station is determined. Based on the Pearson correlation coefficient (PCC) between the displacement predicted by the GEACE model and the displacement observed by GNSS, and the scaling factors of several neighboring virtual stations, the scaling factor of the virtual station is interpolated according to the following formula:

[0160]

[0161] Where i represents the number of GNSS stations adjacent to the virtual station; ρ i PCC represents the displacement predicted by GNSS observations and the GRACE model for the i-th station; β i This represents the scaling factor for the i-th station.

[0162] For example, due to the limitations of GRACE's actual spatial resolution, the spatial resolution of the GRACE-VCD selected as a virtual station should not exceed GRACE's actual spatial resolution (300 km). Furthermore, the number of selected virtual GNSS stations should not be excessive, so that GNSS observations can dominate the inversion process. Based on these two points, this invention selects the location of virtual stations based on Thiessen polygons and buffer zones. Figure 5 This demonstrates our proposed site selection scheme. Figure 5 (a) shows the generation of Thiessen polygons (black lines) based on existing GNSS stations (black dots), with gray boxes representing the grid to be estimated. Figure 5 (b) shows a buffer zone (yellow circle) with a radius of 100 km generated based on existing GNSS stations. Figure 5 (c) shows the added virtual stations (black rectangles) generating buffers (green circles). Ultimately, eight virtual stations were selected to constrain the edge regions of the study area, ensuring the minimum number of stations introduced while maintaining the reliability of the inversion results.

[0163] Before determining the scale factor of the virtual stations, it is first necessary to establish the correlation between the measured GNSS time series and the displacement predicted by the GRACE model for all stations in Yunnan Province. A higher PCC indicates that the GNSS is less affected by observation conditions, errors, and other factors, and represents more the motion caused by the hydrological cycle. Figure 6 The spatial distribution of PCC and the reciprocal of the scaling factor between GNSS vertical crustal displacement and GRACE model predicted displacement is shown. Figure 6 In (a), the circles represent GNSS stations, and the colors represent the PCC between GNSS and GRACE displacements. Figure 6 In (b), the circles represent GNSS stations, the rectangles represent virtual stations, the colors represent the reciprocal of the scaling factor, the gray lines represent Thiessen polygons generated based on existing GNSS stations, and the blue lines represent rivers. Figure 6The PCC values ​​for both are between 0.65 and 0.87. The scaling factor between the measured GNSS time series and the displacement predicted by the GRACE model is between 1.49 and 2.85. A larger scaling factor indicates that GRACE underestimates the deformation displacement caused by hydrological activity. Therefore, the scaling factor of the virtual station is determined by combining the correlation and scaling factor between the measured displacements of several nearby GNSS stations and the displacement predicted by the GRACE model.

[0164] Step S700: Establish a joint inversion ΔTWS model based on the observation values ​​of the real GNSS observation station and the virtual station, as well as the scaling factor of the virtual station;

[0165] In this embodiment, GNSS and GRACE jointly invert the regional ΔTWS, that is, by using the obtained GNSS observed displacement and the displacement predicted by the GRACE model after scaling factor correction, a joint inversion model is established, including the following steps:

[0166] First, a GNSS-only inversion ΔTWS model is established.

[0167] GNSS networks can accurately acquire three-dimensional deformation characteristics of the Earth's surface. The Green's function maps ΔTWS to vertical and horizontal displacements to quantify the elastic response of the surface. Therefore, based on the sensitivity of the Green's function to near-field mass changes, effective capture of local fine-grained load signals can be achieved under reasonable station distribution. A weighted least squares inversion algorithm is used to minimize the following mismatch function:

[0168]

[0169] Where W is the weighting matrix; A is the Green's function matrix, which serves as the carrier linking the three-dimensional deformation of each GNSS station to the mass change of each grid; x is the estimated equivalent water height vector; and b is the GNSS observation vector. To avoid unrealistic spatial variations between adjacent grids, we use the Laplace operator. This is used to define model roughness, ensuring a smooth variation in the quality of adjacent grid cells. Smoothing factor Used to balance the relative weights between model roughness and data mismatch, appropriate The value is given by the L curve.

[0170] Secondly, a joint inversion ΔTWS model of GNSS and GRACE (GNSS-GRACE model) is established.

[0171] The sensitivity of GNSS-only models decreases with increasing distance. Therefore, when station coverage is insufficient or the spatial distance between stations is large, the inversion results are affected, leading to unrealistic local mass variations. To address the problem of limited inversion accuracy caused by insufficient or uneven distribution of GNSS stations, this invention uses GRACE as an auxiliary observation method and incorporates the model-predicted displacement derived from the equivalent water height obtained by GRACE into the inversion model, providing constraints for the inversion results. The principle of the new method can be expressed as:

[0172]

[0173] Where W and U are weighted matrices; A and A' are weighted matrices. G b and b represent the Green's function matrices generated based on GNSS and GRACE observations, respectively; G These represent the observations from the GNSS station and the virtual station (the displacement predicted by the GRACE model), respectively. It is a smoothing factor; Let F(x) be the Lappas operator. To minimize F(x), the above equation can be further expressed as:

[0174]

[0175] Among them, A v and A u b represents the vertical and horizontal Green's function matrices generated based on existing GNSS stations. n ,b e ,b u These are the three-dimensional coordinate time series in the GNSS station-centered coordinate system (ENU); b G denoted as the displacement of the virtual station; W is the weighting matrix of GNSS observations, defined based on the standard error of the GNSS observations; U is the weighting matrix of GRACE observations, defined by the scaling factor between GNSS and GRACE.

[0176]

[0177]

[0178] Where, σ e , σ n , σ u These represent the uncertainties in the ENU coordinate system; i represents the number of GNSS stations; β represents the uncertainty. j is the scaling factor for the virtual stations; j represents the number of virtual stations.

[0179] The smoothing factor in the inversion model is then determined using the L-curve method. Inverting ΔTWS is a typical ill-posed problem because the parameters to be estimated are much larger than the number of observations. Therefore, this embodiment uses the Tikhonov regularization principle to address this issue, and derives the relative weights of data mismatch and model roughness from the L-curve.

[0180]

[0181]

[0182] in, The parameters to be estimated; This represents the fitting residual. Therefore, the L-curve is formed by the points... The curve formed by these points has the maximum curvature, and the point on this curve can be chosen as the approximately optimal smoothing factor.

[0183] To evaluate the performance of the GNSS-GRACE inversion method proposed in this invention in practice, we use the inversion results of the GNSS-only model as a comparison, and GRACE-Mascon and GLDAS products as references. Figure 7 The data shows that Yunnan's water resources distribution exhibits a spatial trend of gradually decreasing from southwest to northeast. The EWH amplitude range derived from the GNSS inversion model is approximately 40-300 mm, while the GNSS-GRACE inversion model is 55-320 mm, GRACE is 30-220 mm, and GLDAS is 30-175 mm. The amplitude of ΔTWS inferred by the GNSS-only and GNSS-GRACE inversion models is significantly larger than that of the GRACE and GLDAS data. This is due to the coarse spatial resolution of GRACE and the simplification of the hydrological model by GLDAS. We note that due to insufficient station density in the Nanpanjiang River basin, the maximum deviation of the ΔTWS amplitude inferred by the GNSS-only and GNSS-GRACE inversion models reaches 15-100 mm, indicating that the GNSS-only model significantly underestimates the variation in ΔTWS. Figure 8The data from left to right shows the equivalent water height changes in the Jinsha River, Lancang River, and Nanpan River basins in Yunnan Province. From top to bottom, the data show the equivalent water height time series derived from the GNSS-GRACE model (black line), the GNSS-Only model (red line), the GRACE model (green line), and the GLDAS model (cyan line). The gray dashed line represents the corresponding monthly climatology (multi-year monthly average), and the gray shading represents the monthly equivalent water height deficit (monthly actual value minus monthly climatology). In the Jinsha River and Lancang River basins, the EWH changes inferred by the GNSS-GRACE and GNSS-Only models show strong correlation, with PCCs of 0.88 and 0.86, respectively. Under these scenarios with a reasonable distribution of GNSS stations, both models show high consistency with GRACE and GLDAS. However, in the Nanpan River basin, where the spatial distribution of GNSS stations is relatively sparse, the spatial constraints provided by GRACE become very important. In this scenario, the EWH derived by the GNSS-only and GNSS-GRACE algorithms shows a significant difference, with a PCC of only 0.44. The PCCs of the EWH derived by the GNSS-only algorithm compared to GRACE and GLDAS are -0.13 and -0.17, respectively. Under the same conditions, the inversion results of the GNSS-GRACE model still show high consistency with GRACE and GLDAS, with PCCs of 0.77 and 0.82 for EWH, respectively. Therefore, introducing GRACE observations is essential in cases of sparse station distribution, ensuring that the inversion of local ΔTWS does not undergo unrealistic changes.

[0184] Step S800: Obtain the change in terrestrial water storage based on the inversion ΔTWS model.

[0185] In summary, this embodiment combines GNSS and GRACE observation data, demonstrating advantages in spatial resolution and reliability compared to existing methods, and can capture detailed features that cannot be captured using a single observation method. Compared to existing methods, it improves the reliability and stability of observations in areas with sparse GNSS station coverage, enabling the acquisition of more reliable ΔTWS even when ground observations are insufficient. The method in this embodiment exhibits advantages in simplicity and robustness because it introduces new observation data without altering the form of traditional methods. It has good scalability and can be quickly adapted and applied over a large area.

[0186] Exemplary device

[0187] like Figure 9 As shown in the illustration, this embodiment also provides a land water storage change inversion device that takes into account the scaling factor, the device comprising:

[0188] Data acquisition module 10 is used to acquire GNSS three-dimensional coordinate time series, GNSS real observation station locations, and equivalent water height data of GRACE-Mascon products;

[0189] The GNSS three-dimensional coordinate time series acquisition module 20 is used to preprocess the GNSS three-dimensional coordinate time series to obtain the processed GNSS three-dimensional coordinate time series.

[0190] The GNSS observation displacement acquisition module 30 is used to obtain the GNSS observation displacement caused by hydrological load based on the processed GNSS three-dimensional coordinate time series.

[0191] The GRACE model predicted displacement acquisition module 40 is used to obtain the GRACE model predicted displacement caused by hydrological load based on the equivalent water height data of the GRACE-Mascon product.

[0192] The scaling factor acquisition module 50 for real GNSS observation stations is used to obtain the scaling factor of real GNSS observation stations based on the GNSS observation displacement and the displacement predicted by the GRACE model.

[0193] The virtual station scaling factor acquisition module 60 is used to interpolate the scaling factor of the virtual station based on the scaling factor of the real GNSS observation station.

[0194] The model building module 70 is used to build a joint inversion ΔTWS model based on the observation values ​​of the real GNSS observation station and the virtual station, as well as the scaling factor of the virtual station.

[0195] The inference module 80 is used to obtain the change in terrestrial water storage based on the inversion ΔTWS model.

[0196] In one implementation, the GNSS three-dimensional coordinate time series acquisition module 20 includes:

[0197] The NEU coordinate time series acquisition unit is used to convert the rectangular coordinate system of the GNSS three-dimensional coordinate time series into the station-centered coordinate system to obtain the NEU coordinate time series of the GNSS.

[0198] The NEU coordinate time series acquisition unit is used to perform gross error removal on the NEU coordinate time series of the GNSS using the interquartile range method to obtain the NEU coordinate time series after removal.

[0199] The processed GNSS three-dimensional coordinate time series acquisition unit is used to establish a GNSS three-dimensional coordinate time series model, and to process the removed NEU coordinate time series according to the GNSS three-dimensional coordinate time series model to obtain the processed GNSS three-dimensional coordinate time series.

[0200] In one implementation, the GNSS observation displacement acquisition module 30 of this embodiment includes:

[0201] The source signal probability distribution function acquisition unit is used to solve the posterior probability distribution of unknown parameters in the variational Bayesian independent component analysis framework to obtain the source signal probability distribution function of the GNSS three-dimensional coordinate time series.

[0202] The unknown parameter estimation unit is used to calculate the posterior probability of the hidden variable based on the probability distribution function of the source signal, and obtain the estimated value of the unknown parameter.

[0203] An independent component analysis model acquisition unit is used to obtain an independent component analysis model based on the estimated values ​​of the unknown parameters;

[0204] The GNSS observation displacement acquisition unit is used to obtain the GNSS observation displacement caused by hydrological load through the independent component analysis model.

[0205] In one implementation, the GRACE model prediction displacement acquisition module 40 of this embodiment includes:

[0206] Radial and horizontal displacement acquisition units are used to calculate the radial and horizontal displacements caused by a unit hydrological mass load on the Earth's surface;

[0207] A vertical and horizontal load Green's function acquisition unit is used to obtain the vertical and horizontal load Green's functions based on the radial and horizontal displacements.

[0208] The GRACE model predicted displacement acquisition unit is used to obtain the GRACE model predicted displacement based on the vertical and horizontal load Green's functions and the equivalent water height data of the GRACE-Mascon product.

[0209] In one implementation, the scaling factor acquisition module 50 for the real GNSS observation station described in this embodiment includes:

[0210] The scaling factor acquisition unit for the real GNSS observation station is used to calculate the ratio between the GNSS observed displacement and the GRACE model predicted displacement, and to use the ratio as the scaling factor for the location of the real GNSS observation station.

[0211] In one implementation, the scaling factor acquisition module 60 of the virtual station in this embodiment includes:

[0212] A Thiessen polygon generation unit is used to generate Thiessen polygons based on the actual locations of the GNSS observation stations.

[0213] The virtual station acquisition unit is used to obtain the virtual station based on the Thiessen polygon and the buffer.

[0214] The virtual station scaling factor acquisition unit is used to obtain the virtual station scaling factor based on the Pearson correlation coefficient between the GNSS observed displacement and the displacement predicted by the GRACE model.

[0215]

[0216] Where i is the number of GNSS stations adjacent to the virtual station; ρ i PCC represents the displacement predicted by GNSS observations and the GRACE model for the i-th station; β i This represents the scaling factor for the i-th station.

[0217] In one implementation, the model building module 70 of this embodiment is used to establish a joint inversion ΔTWS model based on the observations of the real GNSS observation station and the virtual station, as well as the scaling factor of the virtual station.

[0218]

[0219] Among them, A v and A u b represents the vertical and horizontal Green's function matrices generated based on actual GNSS observation stations; n ,b e ,b u These are the processed GNSS three-dimensional coordinate time series in the GNSS station-centered coordinate system (ENU); b G represents the observations from the virtual station; W is the weighting matrix of the GNSS observations, defined based on the standard error of the GNSS observations; U is the weighting matrix of the GRACE observations; the elements in the weighting matrix W... σ e , σ n , σ u These represent the uncertainties in the station-centered coordinate system ENU; and the U elements of the weighting matrix. i represents the number of actual GNSS observation stations; β j is the scaling factor for virtual stations; j represents the number of virtual stations.

[0220] Based on the above embodiments, the present invention also provides a smart terminal, the principle block diagram of which can be as follows: Figure 10As shown, the intelligent terminal includes a processor, memory, network interface, display screen, and temperature sensor connected via a system bus. The processor provides computing and control capabilities. The memory includes non-volatile storage media and internal memory. The non-volatile storage media stores the operating system and computer programs. The internal memory provides the environment for the operation of the operating system and computer programs in the non-volatile storage media. The network interface is used to communicate with external terminals via a network connection. When executed by the processor, the computer program implements a method for inverting changes in land water storage considering a scaling factor. The display screen can be an LCD screen or an e-ink screen. The temperature sensor is pre-installed within the intelligent terminal to detect the operating temperature of internal devices.

[0221] Those skilled in the art will understand that Figure 10 The block diagram shown is merely a partial structural diagram related to the present invention and does not constitute a limitation on the smart terminal to which the present invention is applied. A specific smart terminal may include more or fewer components than those shown in the figure, or combine certain components, or have different component arrangements.

[0222] In one embodiment, a smart terminal is provided, comprising a memory, a processor, and a land water storage change inversion program considering a scaling factor stored in the memory and executable on the processor. When the processor executes the land water storage change inversion program considering a scaling factor, it implements the following operation instructions:

[0223] Acquire GNSS three-dimensional coordinate time series, actual GNSS observation station locations, and equivalent water height data from GRACE-Mascon products;

[0224] The GNSS three-dimensional coordinate time series is preprocessed to obtain the processed GNSS three-dimensional coordinate time series;

[0225] Based on the processed GNSS three-dimensional coordinate time series, the GNSS observation displacement caused by hydrological load is obtained;

[0226] Based on the equivalent water height data of the GRACE-Mascon product, the displacement predicted by the GRACE model caused by hydrological load is obtained;

[0227] Based on the GNSS observed displacement and the displacement predicted by the GRACE model, the scaling factor of the actual GNSS observation station position is obtained.

[0228] Based on the scaling factor of the actual GNSS observation station location, the scaling factor of the virtual station is interpolated.

[0229] A joint inversion ΔTWS model is established based on the observations from the actual GNSS observation station and the virtual station, as well as the scaling factor of the virtual station.

[0230] The changes in terrestrial water storage are obtained based on the inversion ΔTWS model.

[0231] Those skilled in the art will understand that all or part of the processes in the methods of the above embodiments can be implemented by a computer program instructing related hardware. The computer program can be stored in a non-volatile computer-readable storage medium. When executed, the computer program can include the processes of the embodiments of the above methods. Any references to memory, storage, operational databases, or other media used in the embodiments provided by this invention can include non-volatile and / or volatile memory. Non-volatile memory may include read-only memory (ROM), programmable ROM (PROM), electrically programmable ROM (EPROM), electrically erasable programmable ROM (EEPROM), or flash memory. Volatile memory may include random access memory (RAM) or external cache memory. By way of illustration and not limitation, RAM is available in a variety of forms, such as static RAM (SRAM), dynamic RAM (DRAM), synchronous DRAM (SDRAM), dual operating data rate SDRAM (DDRSDRAM), enhanced SDRAM (ESDRAM), synchronous link DRAM (SLDRAM), RAMbus direct RAM (RDRAM), direct memory bus dynamic RAM (DRDRAM), and memory bus dynamic RAM (RDRAM), etc.

[0232] In summary, this invention discloses a method for inverting land water storage changes considering the scaling factor. The method includes: preprocessing a GNSS three-dimensional coordinate time series; obtaining the GNSS observation displacement caused by hydrological load based on the processed GNSS three-dimensional coordinate time series; obtaining the GRACE model predicted displacement caused by hydrological load based on the equivalent water height data of the GRACE-Mascon product; obtaining the scaling factor of the actual GNSS observation station based on the GNSS observation displacement and the GRACE model predicted displacement; obtaining the scaling factor of the virtual station based on the location of the actual GNSS observation station; performing a checkerboard test on the actual and virtual GNSS observation stations; establishing an inversion ΔTWS model based on the scaling factor; and performing the inversion. This invention relies on the constraint of the scaling factor of the virtual station, thereby minimizing system differences between observation methods and providing more accurate ΔTWS, which is of great significance for hydrological applications and water resource management in areas with sparse GNSS station coverage.

[0233] 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 land water storage change inversion method considering a scaling factor, characterized in that, The method includes: Acquire GNSS three-dimensional coordinate time series, actual GNSS observation station locations, and equivalent water height data from GRACE-Mascon products; The GNSS three-dimensional coordinate time series is preprocessed to obtain the processed GNSS three-dimensional coordinate time series; Based on the processed GNSS three-dimensional coordinate time series, the GNSS observation displacement caused by hydrological load is obtained; Based on the equivalent water height data of the GRACE-Mascon product, the displacement predicted by the GRACE model caused by hydrological load is obtained; Based on the GNSS observed displacement and the displacement predicted by the GRACE model, the scaling factor of the actual GNSS observation station position is obtained. Based on the scaling factor of the actual GNSS observation station location, the scaling factor of the virtual station is interpolated. The joint inversion is established according to the observation values of the GNSS real observation station and the virtual station, and the scale factor of the virtual station model; According to the inversion Model, the land water storage change is obtained; The step of interpolating the scale factor of the virtual station based on the scale factor of the actual GNSS observation station location includes: Generate Thiessen polygons based on the actual locations of the GNSS observation stations; The virtual station is obtained based on the Thiessen polygon and the buffer. Based on the Pearson correlation coefficient between the GNSS observed displacement and the displacement predicted by the GRACE model, the scaling factor of the virtual station is obtained as follows: ， wherein, is the number of GNSS stations adjacent to the virtual station; represents the PCC of the GNSS observations and the GRACE model predicted displacements for the i-th station; represents the PCC of the GNSS observations and the GRACE model predicted displacements for the i-th station; represents the scaling factor for the i-th station; The inversion was established based on the observations from the actual GNSS observation station and the virtual station, as well as the scaling factor of the virtual station. The model is: ， in, and The vertical and horizontal Green's function matrices are generated based on real GNSS observation stations; , , These are the processed GNSS three-dimensional coordinate time series in the GNSS station-centered coordinate system (ENU); These are observations from a virtual station; This is the weighting matrix for GNSS observations, defined based on the standard error of the GNSS observations; This is the weighted matrix of GRACE observations; weighted matrix elements in , , These represent the uncertainties in the station-centered coordinate system ENU; the weighting matrix... element , This indicates the number of actual GNSS observation stations; The scaling factor for virtual stations; This indicates the number of virtual stations.

2. The method for inverting changes in terrestrial water storage considering scaling factors according to claim 1, characterized in that, The preprocessing of the GNSS three-dimensional coordinate time series to obtain the processed GNSS three-dimensional coordinate time series includes: The rectangular coordinate system of the GNSS three-dimensional coordinate time series is converted to the station-centered coordinate system to obtain the NEU coordinate time series of the GNSS. The NEU coordinate time series of the GNSS was subjected to gross error removal using the interquartile range method to obtain the NEU coordinate time series after removal. A GNSS three-dimensional coordinate time series model is established, and the NEU coordinate time series after elimination is processed according to the GNSS three-dimensional coordinate time series model to obtain the processed GNSS three-dimensional coordinate time series.

3. The method for inverting changes in terrestrial water storage considering proportionality factors according to claim 1, characterized in that, The step of obtaining the GNSS observation displacement caused by hydrological load based on the processed GNSS three-dimensional coordinate time series includes: Solve the posterior probability distribution of the unknown parameters in the variational Bayesian independent component analysis framework to obtain the source signal probability distribution function of the GNSS three-dimensional coordinate time series; The posterior probability of the hidden variables is calculated based on the probability distribution function of the source signal to obtain the estimated value of the unknown parameters; Based on the estimated values ​​of the unknown parameters, an independent component analysis model is obtained; The GNSS observation displacement caused by hydrological load is obtained through the independent component analysis model.

4. The method for inverting changes in terrestrial water storage considering proportionality factors according to claim 1, characterized in that, The step of obtaining the GRACE model-predicted displacement caused by hydrological load based on the equivalent water height data of the GRACE-Mascon product includes: Calculate the radial and horizontal displacements caused by a unit hydrological mass load on the Earth's surface; The Green's functions for vertical and horizontal loads are obtained based on the radial and horizontal displacements. The displacement predicted by the GRACE model is obtained based on the Green's functions of the vertical and horizontal loads and the equivalent water height data of the GRACE-Mascon product.

5. The method for inverting changes in terrestrial water storage considering proportionality factors according to claim 1, characterized in that, The scaling factor for obtaining the actual GNSS observation station position based on the GNSS observed displacement and the displacement predicted by the GRACE model includes: Calculate the ratio between the GNSS observed displacement and the displacement predicted by the GRACE model, and use the ratio as a scaling factor for the actual GNSS observation station locations.

6. An inversion apparatus based on the method for inverting changes in terrestrial water storage considering the scaling factor as described in any one of claims 1-5, characterized in that, The device includes: The data acquisition module is used to acquire GNSS three-dimensional coordinate time series, GNSS actual observation station locations, and equivalent water height data of GRACE-Mascon products; The GNSS three-dimensional coordinate time series acquisition module is used to preprocess the GNSS three-dimensional coordinate time series to obtain the processed GNSS three-dimensional coordinate time series. The GNSS observation displacement acquisition module is used to obtain the GNSS observation displacement caused by hydrological load based on the processed GNSS three-dimensional coordinate time series. The GRACE model predicted displacement acquisition module is used to obtain the GRACE model predicted displacement caused by hydrological load based on the equivalent water height data of the GRACE-Mascon product. The module for obtaining the scale factor of the real GNSS observation station is used to obtain the scale factor of the real GNSS observation station position based on the GNSS observation displacement and the displacement predicted by the GRACE model. The virtual station scaling factor acquisition module is used to interpolate the scaling factor of the virtual station based on the scaling factor of the actual GNSS observation station location. The model building module is used to establish a joint inversion based on the observations from the actual GNSS observation station and the virtual station, as well as the scaling factor of the virtual station. Model; The deduction module is used to perform the inversion based on the inversion. The model yields changes in terrestrial water storage.

7. A smart terminal, characterized in that, The smart terminal includes a memory, a processor, and a land water storage change inversion program that takes into account the scaling factor and is stored in the memory and can run on the processor. When the processor executes the land water storage change inversion program that takes into account the scaling factor, it implements the steps of the land water storage change inversion method that takes into account the scaling factor as described in any one of claims 1-5.

8. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores a land water storage change inversion program that takes into account a scaling factor. When the scaling factor-based land water storage change inversion program is executed by a processor, it implements the steps of the land water storage change inversion method that takes into account a scaling factor as described in any one of claims 1-5.

Images