A GNSS-IR inverse modeling method and system constrained by random walk and time-varying physical model

Abstract

The application discloses a GNSS-IR inverse modeling method and system constrained by random walk and time-varying physical models, and the method comprises the following steps: acquiring a signal-to-noise ratio observation value and preprocessing; adopting a random walk model to describe antenna height variation; splitting an initial phase of a signal-to-noise ratio oscillation term into a frequency deviation constant term and a common time-varying random walk term; constructing a nonlinear least square loss function containing data fitting residuals and double random walk constraints; adopting an iterative solution strategy and adaptively updating a standard deviation factor based on a median absolute deviation, and outputting an antenna height inversion result. In a gentle water level scene, the RMSE is 0.94 cm (the accuracy is improved by 13.0%); in a severe water level scene, the traditional method fails, the RMSE of the application is only 2.14 cm, and the correlation coefficient is 0.9992. The application does not need fine priori, is self-adaptive, and is suitable for dynamic monitoring of water levels in various shorelines.

Description

Technical Field

[0001] This invention relates to a GNSS-IR inverse modeling method and system with random walk and time-varying physical model constraints. Background Technology

[0002] Water level is a core monitoring element in water resource management, natural disaster early warning, and shipping control. Accurate, continuous, and dynamic water level monitoring is of great significance for understanding hydrological processes and ensuring the safe operation of water-related businesses. As the core transition zone between land and water, shoreline waters are affected by multiple factors such as runoff, tides, and storm surges, resulting in complex dynamic processes that place high demands on the temporal continuity, automation, and real-time performance of monitoring technologies.

[0003] Currently widely used water level monitoring technologies mainly include water level gauge measurement, satellite altimetry, and satellite gravity inversion. Water level gauge measurement has high temporal resolution, but it can only provide point-like information with limited coverage. Furthermore, it is a relative measurement and is susceptible to the effects of crustal movement and subsidence. Contact-based measurements also face problems such as siltation and corrosion. Satellite altimetry can obtain absolute water levels, but it suffers from low spatiotemporal resolution and waveform interference from land in small-scale nearshore waters. Satellite gravity inversion can be used to monitor large-scale sea-level changes, but the leakage of hydrological signals from land within a 300 km radius near the coast makes it difficult to reflect true nearshore water level changes. Therefore, existing technologies still have significant shortcomings in continuous, dynamic, and high-precision water level monitoring in coastal waters.

[0004] GNSS Interferometric Reflectometry (GNSS-IR) technology offers a new technical approach to solving the aforementioned problems. This technology utilizes the interference information between direct and reflected signals received by a shore-based GNSS receiver to invert the vertical distance from the antenna to the water surface from signal-to-noise ratio (SNR) observation data, thereby obtaining water level changes. GNSS-IR offers advantages such as all-weather, all-time operation, non-contact operation, low cost, and the ability to provide absolute water level readings, making it a research hotspot in the fields of hydrological remote sensing and marine mapping.

[0005] Existing GNSS-IR water level inversion methods are mainly divided into two categories: spectrum analysis method and simulation method.

[0006] Spectral analysis methods (such as the Lomb-Scargle periodogram method) calculate reflection height by extracting the dominant frequency of the SNR oscillation sequence. These methods are simple and computationally efficient. However, they assume a constant water level within each satellite observation arc, a static assumption that conflicts significantly with dynamic water level scenarios such as tides and floods. Although existing research has improved these methods through dynamic height correction and wavelet analysis, it remains difficult to fundamentally eliminate the dependence on the static assumption within the arc. Furthermore, in highly dynamic scenarios, spectral multi-peak interference is prone to occur, resulting in low data utilization and difficulty in meeting real-time inversion requirements.

[0007] The simulation method constructs a complete physical model of SNR interferometry and simultaneously solves for water level and other model parameters using nonlinear least squares and other parameter estimation methods. The inverse modeling method uses reflection height as a time-varying parameter and employs B-spline functions for modeling, enabling the output of continuous water level sequences. Its inversion accuracy is superior to spectral analysis and it is more suitable for highly dynamic water level scenarios. However, existing inverse modeling methods based on B-spline functions still have the following shortcomings:

[0008] Parameter settings are sensitive and lack universally applicable criteria: the number and distribution (node ​​spacing) of B-spline control points significantly affect the inversion results. Too few control points lead to underfitting, failing to capture rapid water level fluctuations; too many control points lead to overfitting, resulting in unstable solutions or even failures. The distribution of control points (such as quasi-uniform distribution) is difficult to adapt to complex scenarios where rapid and gradual water level changes coexist.

[0009] Strong dependence on prior information: Inverse modeling methods require reasonable initial water level values, B-spline node parameters, and other prior information. In areas with drastic water level changes (such as tidal ranges exceeding 3m), a single fixed initial value deviates significantly from the actual water level at other times, easily leading to nonlinear least squares solutions getting trapped in local optima or diverging directly.

[0010] The physical model is too simplified: Existing models usually set the initial phase of the interference signal oscillation term as a constant, ignoring the time-varying influence of environmental factors such as water temperature, salinity, wind and waves on the reflection phase. This leads to deviations between the model and the real physical process, especially in scenarios with drastic water level changes, where the inversion error increases sharply or even fails.

[0011] Therefore, existing inverse modeling methods based on B-spline functions still have shortcomings in terms of parameter adaptation capability, dependence on prior information, and physical model's ability to represent time-varying characteristics, making it difficult to meet the engineering and automated monitoring needs of various water level change scenarios, such as gentle and drastic changes. Summary of the Invention

[0012] To address the technical problems of existing GNSS-IR inverse modeling methods, such as the sensitivity of B-spline function parameter settings, poor universality, difficulty in adapting to scenarios with drastic water level changes, and low inversion accuracy or even failure under high dynamic conditions due to neglecting the time-varying characteristics of reflection phase, this invention provides a GNSS-IR inverse modeling method based on random walk and time-varying physical model constraints. This method uses a random walk model to characterize the continuous dynamic changes in antenna height and decomposes the initial phase of the signal-to-noise ratio oscillation term into an inter-frequency deviation constant term and a common time-varying random walk term. Simultaneously, an iterative solution strategy is introduced to achieve adaptive estimation of model parameters. This enables continuous high-precision water level inversion without the need for detailed prior information or human intervention, and can operate stably under various hydrological scenarios, including gentle changes and drastic fluctuations.

[0013] According to one aspect of the present invention, a GNSS-IR inverse modeling method constrained by random walk and time-varying physical model is provided, comprising: acquiring signal-to-noise ratio (SNR) observation data collected by a GNSS receiver, preprocessing the SNR observation data, and extracting the oscillation component of the reflected signal; describing the change of antenna height over time using a random walk model, representing the antenna height of the current epoch as the antenna height of the previous epoch plus random walk noise following a zero-mean Gaussian distribution; decomposing the initial phase in the oscillation component of the reflected signal into an inter-frequency deviation term and a common time-varying term, wherein the inter-frequency deviation term is a constant term independent of different frequency bands, and the common time-varying term is described by a random walk model to describe its change over time; constructing a nonlinear least squares loss function for the SNR observation based on the antenna height sequence constructed by the random walk model, the inter-frequency deviation term, and the common time-varying term; and adopting an iterative solution strategy to update the antenna height parameters, inter-frequency deviation parameters, and common time-varying parameters in each iteration until the convergence condition is met, and outputting the antenna height inversion results for each epoch.

[0014] This invention makes the following two improvements to the observation model: First, in order to effectively express the continuous change of water level, a random walk is used to represent the change of antenna height; second, considering the time-varying characteristics of GNSS signals and the external environment, the initial phase of the oscillation term of the SNR observation is modeled as a time-varying function.

[0015] As a further technical solution, the iterative solution strategy includes: in the first iteration, the initial value of the antenna height is set to a preset coarse value, and the standard deviation factor of the antenna height and the standard deviation factor of the random walk of the common time-varying term are set to preset initial values, while the initial values ​​of the inter-frequency deviation parameter and the reflected signal amplitude parameter are set to zero; starting from the second iteration, the antenna height results of each epoch output from the previous iteration are used as the initial value of the antenna height in the current iteration, and the standard deviation factor of the antenna height in the current iteration is calculated based on the median absolute deviation of the antenna height results in the previous iteration; the iteration is repeated until the root mean square error of the antenna height inversion results of the two iterations is less than a preset threshold or the preset maximum number of iterations is reached.

[0016] As a further technical solution, the common time-varying term is shared across all frequency bands and is used to characterize the global time-varying characteristics of the reflection phase driven by environmental factors; the inter-frequency deviation term is used to characterize the fixed phase shift caused by differences in wavelength and reflection characteristics in different GNSS frequency bands.

[0017] As a further technical solution, the variance of the random walk noise is proportional to the time difference between the current epoch and the previous epoch, and the variance is expressed as: ,in The standard deviation factor of the antenna height. The variance is adjusted based on the actual time interval when there are missing observation data.

[0018] As a further technical solution, the random walk standard deviation factor of the common time-varying term is calculated based on the median absolute deviation of the common time-varying term sequence obtained in the previous iteration, which is consistent with the calculation method of the standard deviation factor of the antenna height.

[0019] As a further technical solution, the nonlinear least squares loss function consists of three parts: the sum of squared residuals between the signal-to-noise ratio observations and the model predictions, the constraint term of the random walk of the antenna height, and the constraint term of the random walk of the common time-varying term; the loss function is solved iteratively using the Gauss-Newton method.

[0020] According to one aspect of the present invention, a GNSS-IR inverse modeling system constrained by random walk and time-varying physical model is provided, comprising: a data acquisition module for acquiring signal-to-noise ratio (SNR) observation data collected by a GNSS receiver, and preprocessing the SNR observation data to extract the oscillation component of the reflected signal; an antenna height dynamic modeling module for describing the change process of antenna height over time using a random walk model, representing the antenna height of the current epoch as the antenna height of the previous epoch plus random walk noise following a zero-mean Gaussian distribution; and a time-varying phase modeling module for converting the reflected signal oscillation component into a random walk model. The initial phase is decomposed into an inter-frequency deviation term and a common time-varying term, wherein the inter-frequency deviation term is a constant term independent of different frequency bands, and the common time-varying term is described by a random walk model to describe its change over time; the loss function construction module is used to construct a nonlinear least squares loss function for the signal-to-noise ratio observation based on the antenna height sequence constructed by the random walk model, the inter-frequency deviation term, and the common time-varying term; the iterative solution module is used to update the antenna height parameters, inter-frequency deviation parameters, and common time-varying parameters in each iteration using an iterative solution strategy until the convergence condition is met, and output the antenna height inversion results for each epoch.

[0021] As a further technical solution, the iterative calculation module includes: an initialization unit, used to set the initial value of the antenna height to a preset coarse value in the first iteration, set the standard deviation factor of the antenna height and the standard deviation factor of the random walk of the common time-varying term to preset initial values, and set the initial values ​​of the inter-frequency deviation parameter and the reflected signal amplitude parameter to zero; a parameter update unit, used to take the antenna height results of each epoch output by the previous iteration as the initial value of the antenna height in the current iteration starting from the second iteration, and calculate the standard deviation factor of the antenna height in the current iteration based on the median absolute deviation of the antenna height results in the previous iteration; and a convergence judgment unit, used to judge whether the root mean square error of the antenna height inversion results of the two iterations is less than a preset threshold or whether the preset maximum number of iterations has been reached.

[0022] According to one aspect of the present invention, an electronic device is provided, including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor, when executing the program, implements the GNSS-IR inverse modeling method with random walk and time-varying physical model constraints.

[0023] According to one aspect of the present invention, a computer-readable storage medium is provided having a computer program stored thereon that, when executed by a processor, implements the GNSS-IR inverse modeling method with random walk and time-varying physical model constraints.

[0024] Compared with the prior art, the beneficial effects of the present invention are as follows: 1. Adaptive water level inversion without requiring detailed prior information was achieved. This invention uses a random walk model to describe antenna height changes and splits the initial phase of the signal-to-noise ratio oscillation term into a frequency band-independent constant term (inter-frequency deviation) and a common time-varying random walk term. Combined with an iterative solution strategy, it gets rid of the strong dependence of traditional inverse modeling methods on the number, distribution and precise initial values ​​of B-spline control points, and achieves adaptive high-precision inversion without human intervention.

[0025] 2. Significantly improved inversion accuracy in scenarios with gradual water level changes. In a verification experiment at one of the monitoring stations (daily variation amplitude less than 0.5 m, average variation rate 0.1 m / h), the root mean square error (RMSE) between the inversion results of the method of this invention and the measured values ​​of the water level gauge was 0.94 cm, with a correlation coefficient of 0.9953. Compared with the traditional inverse modeling method, the RMSE decreased from 1.08 cm to 0.94 cm, improving the inversion accuracy by 13.0%; the median absolute deviation (MAD) decreased from 0.73 cm to 0.68 cm.

[0026] 3. Overcomes the limitations of traditional methods in scenarios with drastic water level fluctuations. In a verification experiment at another monitoring station (daily variation exceeding 3 m, average rate of change 0.2 m / h), the traditional inverse modeling method completely failed (root mean square error reached 48.58 cm, correlation coefficient only 0.4878), while the method of this invention maintained stable and high-precision inversion: root mean square error only 2.14 cm, correlation coefficient as high as 0.9992, and absolute median deviation of 1.32 cm. This indicates that the method of this invention effectively solves the core technical problem of the easy failure of traditional inverse modeling methods under high dynamic water level scenarios.

[0027] 4. The inversion error structure has been optimized, improving long-term stability. By using STL time series decomposition to retrieve residuals, the method of this invention effectively suppresses trend errors and periodic residues in scenarios with gentle water levels; in scenarios with severe water levels, it completely eliminates systematic trends and periodic residues, with residuals stabilizing near zero. The long-term stability of the retrieval results is significantly better than that of traditional methods.

[0028] 5. Possesses broad applicability across various scenarios. The method of this invention has achieved high-precision inversion results in two typical hydrological scenarios: gradual change (GTGU) and violent fluctuation (SC02). It does not require adjustment of the core parameters of the model for different scenarios and can be directly applied to the dynamic monitoring of water levels in various shoreline waters such as oceans, lakes, rivers, and reservoirs. Attached Figure Description

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

[0030] Figure 1 A flowchart illustrating a GNSS-IR inverse modeling method with random walk and time-varying physical model constraints provided in an embodiment of the present invention; Figure 2 The second-order basis function in the embodiments of the present invention A schematic diagram; Figure 3 This is a comparison chart of inverse modeling and inversion results with different numbers of B-spline control points in embodiments of the present invention; Figure 4 The figure shows the experimental results of the inverse modeling method with different numbers of control points for B-spline functions in this embodiment of the invention. Figure 5 This is a schematic diagram of the first Fresnel reflection areas of sites A and B in an embodiment of the present invention; Figure 6 This is a diagram showing the water level inversion results of the original inverse modeling method for stations A and B under the optimal number of control points in this embodiment of the invention. Figure 7 This is a comparison chart of the water level inversion results of the original inverse modeling method and the improved inverse modeling method for sites A and B in this embodiment of the invention; Figure 8 This is a graph showing the STL time series decomposition results of the water level inversion residuals at stations A and B in an embodiment of the present invention. Detailed Implementation

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

[0032] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention. In addition, the technical features of the various embodiments or individual embodiments provided by the present invention can be arbitrarily combined to form new technical solutions. Such combinations are not bound by the order of steps and / or structural composition patterns, but must be based on the ability of those skilled in the art to implement them. When the combination of technical solutions is contradictory or cannot be implemented, it should be considered that such a combination of technical solutions does not exist and is not within the scope of protection claimed by the present invention.

[0033] This embodiment provides a GNSS-IR inverse modeling method constrained by random walk and time-varying physical model. This method can adaptively invert continuous antenna height sequences using multi-frequency and multi-system GNSS signal-to-noise ratio observation data, thereby obtaining water level changes.

[0034] like Figure 1 As shown, the method in this embodiment mainly includes two parts: data preprocessing and nonlinear least squares iterative solution. Data preprocessing provides the oscillation components of the reflected signal for subsequent calculations; the nonlinear least squares iterative solution includes constructing a random walk antenna height model, constructing a time-varying phase model, constructing a loss function, iteratively solving, and outputting the results.

[0035] Step 1: Data Acquisition and Preprocessing.

[0036] First, obtain the signal-to-noise ratio (SNR) observations from the GNSS receiver. To ensure the reflected signal originates from the water surface, determine the azimuth and elevation mask ranges based on the first Fresnel reflection area around the station, and extract the reflected signal arc segment from the water surface. Figure 5 The first Fresnel reflection area of ​​the GTGU and SC02 sites is shown as an example.

[0037] The extracted signal-to-noise ratio (SNR) sequence is fitted with a second-order polynomial to remove the trend term, resulting in an SNR oscillation term that includes interferometric oscillation components. Simultaneously, bending error corrections are applied to the elevation angle (e.g., using the GPT3 model and the Bennett model), and the relative tropospheric delay is calculated (e.g., using the GPT3 model and the VMF3 mapping function), which is then introduced into the observation model as a known correction. These preprocessing methods are well-known techniques in the field and will not be detailed here.

[0038] Step 2: Construct a random walk antenna height model.

[0039] To effectively represent the continuous change in antenna height (the vertical distance from the antenna phase center to the water surface), this embodiment uses a random walk model to describe the change in antenna height over time.

[0040] The random walk model is represented as:

[0041]

[0042] In the formula, This is the antenna height at the current epoch; It is the antenna height of the previous epoch; Follows the variance The zero-mean Gaussian distribution, i.e. ,in,

[0043]

[0044] In the formula, It is the standard deviation factor of the antenna height, and the unit is _____. ; This represents the time difference between the current epoch and the previous epoch. The reason for designing the variance to be related to the time interval is to account for the sampling interval of the observation data and the case where some epochs of the observation data are missing. Compared to B-spline functions, random walks only need to estimate the standard deviation factor. That's fine. Too big or too small. Both can lead to parameter fitting failure, inversion failure, and inaccurate standard deviation factor. This is a prerequisite for a successful random walk.

[0045] Step 3: Construct a time-varying phase model.

[0046] Those skilled in the art know that the reflection characteristics of GNSS signals are closely related to the electromagnetic properties (dielectric constant and conductivity) of the reflecting surface. The dielectric constant of water is extremely sensitive to temperature changes; as water temperature increases, the dielectric constant decreases linearly (for example, the dielectric constant of freshwater is approximately 88 at 0°C, decreasing to around 78 at 25°C). The conductivity of water is also temperature-dependent; higher temperatures result in higher conductivity. Furthermore, increased salinity significantly increases the conductivity of water (seawater conductivity is much higher than that of freshwater). In addition to the electromagnetic properties of the reflecting surface, the Fresnel reflection coefficient is also related to the signal wavelength and the satellite elevation angle, which itself exhibits significant temporal dynamic variations.

[0047] In summary, the Fresnel reflection coefficient of GNSS signals dynamically changes with external environmental conditions and time. This dynamic change in the reflection coefficient directly affects the amplitude and phase of the reflected signal, thus influencing the amplitude and phase characteristics of the interference signal. Therefore, treating the amplitude and phase of the interference signal as constants does not conform to objective physical laws. On the other hand, the influence of environmental factors on the signal is mainly reflected in the phase of the reflected signal; the amplitude change is almost negligible and can still be considered constant. Therefore, to reflect the time-varying characteristics of the initial phase of the interference signal, this embodiment models the initial phase in the SNR oscillation term as a constant term independent of the common random walk and frequency band, with the following expression:

[0048]

[0049] in, It is the inter-frequency deviation term (i.e., a constant term independent of frequency bands). These are different frequency bands, reflecting the inherent differences in reflection characteristics of different frequency bands; This is a common time-varying term, shared by all frequency bands, used to characterize the global time-varying characteristics of the reflection phase driven by environmental factors. In this embodiment, the common time-varying term is also modeled as a random walk process:

[0050]

[0051] In the formula, It is the standard deviation factor of the common time-varying term.

[0052] Step 4: Construct a nonlinear least squares loss function.

[0053] Compared to existing inverse modeling methods, the improved inverse modeling method in this embodiment uses antenna height... The common time-varying part of the initial phase is modeled as a random walk. Therefore, the loss function of the improved inverse modeling method includes not only the residual term of the data fitting but also the state constraint term of the random walk.

[0054] The data fitting residual term of the SNR oscillation term is:

[0055]

[0056] In the formula, yes Number of observations These are observed values. These are model predictions.

[0057] high Common time variables of the first phase random walk constraint , They are respectively

[0058]

[0059]

[0060] Adding the two loss functions mentioned above, the total loss function to be minimized is:

[0061]

[0062] This is an overdetermined equation, and the Gauss-Newton method is also used to solve this nonlinear least squares problem.

[0063] Step 5, iterative solution strategy.

[0064] In nonlinear least squares calculations, the selection of initial parameter values ​​has a crucial impact on the solution performance: improper initial values ​​or values ​​deviating from the true values ​​can easily lead to the solution process converging to a local optimum. Furthermore, the higher the nonlinearity of the loss function, the more significant the negative impact of initial value deviations, making it not only more prone to getting trapped in local optima but also potentially causing solution failure. In inverse modeling methods, a common solution to the initial value problem is to obtain the initial water level using LSP spectral analysis. However, this approach has significant limitations in regions with drastic water level changes: a single, fixed initial height may deviate significantly from the actual water level at other times, leading to solution errors or even convergence failure.

[0065] To address the issues of inappropriate initial value selection and difficulty in accurately setting the standard deviation factor for random walks, this embodiment introduces an iterative solution strategy. Through multiple iterations, accurate initial values ​​and standard deviation factors for water level height are obtained. The specific solution process of the improved inverse modeling method is as follows: Figure 1 As shown.

[0066] The specific iteration steps are described below:

[0067] The initial parameter values ​​for the first iteration can be approximate. Taking observation data with an antenna height variation range of 0-10m, a variation rate of 0-10cm / s, and a sampling interval of 30s as an example, the initial antenna height can be set to 5m, and the standard deviation factor of the antenna height random walk can be taken as 1e-4. (Corresponding to a rate of change of 1 cm / s). For the standard deviation factor of the global time-varying term reflecting environmental impact, the initial value is generally taken as 1e-3. That's fine. It's recommended to choose a slightly larger initial value, as it will converge to a suitable value as the iterations progress. For the amplitude, inter-frequency deviation in the initial phase, and roughness parameters, set them to 0.

[0068] The second iteration is based on the solution results of the first iteration. The initial antenna height values ​​for each epoch are directly taken from the antenna height results output in the first iteration, while the standard deviation factor is calculated based on the antenna height results from the first iteration. To avoid outliers caused by noise or model bias, the Median Absolute Deviation (MAD) algorithm is used to calculate the standard deviation factor. The standard deviation factor of the global time-varying term is also calculated using the same method.

[0069] The third iteration takes the solution results from the second iteration as input and further optimizes the solution by using more reliable intermediate results. The method for obtaining the initial values ​​of the parameters remains the same as in the second iteration.

[0070] Subsequent iterations are the same as in the third round, until the iteration convergence condition is met.

[0071] The convergence condition is set as follows: the root mean square error of the antenna height inversion results in the two iterations is within 1 mm, or the number of iterations reaches the maximum number of iterations (e.g., 10). Typically, it takes five iterations to converge to the optimal solution; the result of the fifth iteration is almost identical to that of the fourth, showing no significant improvement.

[0072] Step 6, output the results.

[0073] Output the antenna height sequence obtained from the final iteration.

[0074] The final antenna height sequence obtained from the iteration is output. This antenna height is defined as the vertical distance from the phase center of the GNSS antenna to the water surface.

[0075] If a normal water level (relative to a certain reference surface) is required, the following conversion method can be used: The absolute elevation of the antenna phase center can be obtained through static precise single-point positioning or known control point coordinates. Subtract the antenna height obtained from the inversion of the current epoch from the absolute elevation of the antenna phase center to obtain the absolute elevation of the water surface, i.e., the normal water level.

[0076] If we only need to focus on the change in water level, we can calculate it relative to a certain initial moment: the change in water level is the antenna height at the initial moment minus the antenna height at the current moment.

[0077] In the example verification of this embodiment, in order to compare with the measured value of the water level gauge, the antenna height obtained by inversion is negative and the average value of each time period is subtracted to obtain the water level change sequence, thereby eliminating the difference between different benchmarks.

[0078] To verify the effectiveness of the improved inverse modeling method described in this embodiment, two GNSS monitoring stations (station A and station B) with significantly different water level change characteristics were selected. Based on the measured data of the two stations, a comparative experiment on water level dynamic inversion between the original inverse modeling method and the improved inverse modeling method was carried out simultaneously.

[0079] Station A is equipped with a Leica AR 25 antenna and a Leica GRX1200 receiver, with a sampling interval of 30 seconds. Detailed data is shown in Table 1. A Campbell CS476 radar level gauge is located 300 meters from this station, with a sampling interval of 1 minute; its observation data is collected as a reference for the true water level. The daily water level fluctuation at station A is less than 0.5 meters, the standard deviation of the water level series is approximately 0.06 meters, and the average rate of change is 0.1 meters per hour, exhibiting overall characteristics of small fluctuations and a gentle trend.

[0080] Station B is equipped with a Trimble NETR9 geodetic receiver and a TRM59800.80 antenna, with a sampling interval of 15 seconds. Detailed data is shown in Table 1. A tidal observation station is located 350 meters away, equipped with an Aquatrak acoustic tide gauge, with a sampling interval of 6 minutes. Its observation data is also collected as a reference for the true water level at station SC02. The daily water level fluctuation at station B exceeds 3 meters, with a standard deviation of 0.5 meters and an average rate of change of 0.2 meters per hour, exhibiting characteristics of large fluctuations and drastic changes. Since signals in the same frequency band in GNSS have different channel characteristics, the relevant information of the SNR signal used in this embodiment is shown in Table 2.

[0081] Table 1. Detailed Description of GNSS Data from Station A and Station B

[0082]

[0083] Table 2 Brief characteristics of the GNSS signals used in this embodiment

[0084]

[0085] GNSS observation data preprocessing is fundamental to GNSS-IR water level inversion, providing water surface reflection signals with oscillating characteristics for subsequent inversion. To ensure the reflection signals originate from the water surface, SNR data within a specific angular range is used to eliminate reflection signals from outside the water surface (outside the angular range); this step is called angular masking. The elevation and azimuth ranges for each station are confirmed by plotting the first Fresnel zone of the station and combining it with satellite imagery, such as... Figure 5 As shown.

[0086] Based on the angle mask, a second-order polynomial was used to remove the trend term of the satellite arc segment. Then, the global pressure-temperature model GPT3 was used for elevation angle curvature correction, and the relative tropospheric delay was calculated using the GPT3 and VMF3 mapping model. The relevant data preprocessing strategies are shown in Table 3.

[0087] Table 3 GNSS Data Preprocessing Strategies

[0088]

[0089] Furthermore, since the water level measured by the water level gauge and the antenna height retrieved by GNSS-IR have different references, the average values ​​of their respective time periods were subtracted when calculating the relevant accuracy indicators, and the negative value of the antenna height retrieval result was taken.

[0090] The inversion results of the improved method in this embodiment and the original method are compared and analyzed below.

[0091] To facilitate subsequent comparison and analysis of experimental results, it is necessary to determine the optimal nodes for the B-spline function of the inverse modeling method. Therefore, experiments were first conducted using the inverse modeling method with different numbers of control points, and three accuracy indicators—root mean square error (RMSE), median absolute deviation (MAD), and correlation coefficient (Corr)—were statistically analyzed to find the optimal number of control points. The node spacing adopted a quasi-uniform distribution, with the initial height of station A set at 4.0m and the initial height of station B set at 5.0m. The experimental results are as follows: Figure 4 As shown.

[0092] Figure 4 The results show that, for 3 days of observation data at site A, the optimal number of control points for the inverse modeling method is 44 (RMSE=1.08cm, MAD=0.73cm, Corr=0.9955), with a corresponding node time interval of approximately 1.5 hours. For this site, as the number of control points increases, the water level inversion error of the inverse modeling method first decreases and then increases; both excessively small and excessively large control point numbers reduce the accuracy of water level inversion. In contrast, for 3 days of observation data at site B, the correlation between the number of control points and water level inversion accuracy shows no clear pattern. Only a relatively optimal number of control points can be determined as 22 (RMSE=48.58cm, MAD=5.18cm, Corr=0.4878), with a corresponding node time interval of approximately 3 hours. This result indicates that the application effect of the inverse modeling method at site B is significantly reduced, resulting in poor inversion accuracy. One of the main reasons for these differences may be related to the different ranges and amplitudes of water level fluctuations between sites A and B. To further explore the underlying reasons, this embodiment plots the water level inversion results of the inverse modeling method with the optimal number of control points, as shown in the following figure. Figure 6 As shown.

[0093] Figure 6 (a) shows that the water level at station A changes gradually, making the use of a constant parameter model reasonable and effectively adapting to the gradual fluctuations in water level at this station. However, Figure 6(b) shows that the water level at station B changes drastically. Using constant parameters will cause the estimation results to fail to keep up with the rapid dynamic changes in water level, resulting in increased inversion bias and reduced correlation. This reveals the limitations of the original inverse modeling method in adaptability to scenarios with drastic water level fluctuations.

[0094] Based on the optimal number of B-spline control points, experiments were conducted on stations A and B using both the original inverse modeling method and the improved inverse modeling method to verify the effectiveness of the improved method. GNSS data was preprocessed according to Table 2, including angle masking, detrending, and elevation angle curvature correction. Furthermore, the corresponding relative tropospheric delay was calculated and added to the observation model. The specific data solution strategy for the first iteration is shown in Table 4, and the initial parameters are automatically updated with each iteration. In the experiment, the initial parameters for stations A and B were set to be consistent to verify the universality of the improved method.

[0095] Table 4. First-round iteration parameters of the improved inverse modeling method

[0096]

[0097] Figure 7 The results of water level inversion using different methods at different monitoring stations are shown. Figure 7 (a) in the figure represents the water level inversion experiment results at site A; Figure 7 (b) in the figure shows the water level inversion experiment results at site B; Figure 7 (c) in the figure is a 1:1 comparison scatter plot of the experimental results at site A; Figure 7 (d) in the figure is a 1:1 comparison scatter plot of the experimental results at site B. Table 5 summarizes the three precision indicators for the two methods at the two sites: root mean square error (RMSE), median absolute deviation (MAD), and correlation coefficient (Corr).

[0098] Table 5. Statistics on water level inversion accuracy of improved and original inverse modeling methods at different stations.

[0099]

[0100] Figure 7 The time series curve in (a) shows that the improved inverse modeling method (red curve) has a higher fit with the water level gauge observations (black curve) than the original method (blue curve), and is better able to capture small fluctuations in water level. Figure 7The scatter plot in (c) also shows that the red clusters are more densely clustered near the 1:1 reference line. Statistical indicators of the inversion results at site A (Table 5) indicate that the RMSE of the improved inverse modeling method decreased from 1.08 cm to 0.94 cm, improving accuracy by 13.0%, and the MAD decreased from 0.73 cm to 0.68 cm, while the correlation coefficient remained at an extremely high level above 0.99. This demonstrates that in scenarios with gentle water level fluctuations, the improved method further optimizes the inversion accuracy based on the high adaptability of the original method, controlling the error to the millimeter level.

[0101] Figure 7 The time series in (b) clearly shows that the original inverse modeling method (blue curve) has a huge deviation from the water level gauge observation (black curve) and is completely unable to track the rapid dynamic changes in water level; while the improved method (red curve) almost coincides with the water level observation and accurately captures the violent fluctuations in water level. Figure 7 The scatter plot (d) in the comparison further confirms this: the blue clusters in the original method are extremely scattered and have very weak correlations; while the red clusters in the improved method closely adhere to the 1:1 reference line. Combined with the accuracy indicators in Table 5 (RMSE drops sharply from 48.58 cm to 2.14 cm, MAD drops from 5.18 cm to 1.32 cm, and the correlation coefficient jumps from 0.4878 to 0.9992), this result directly proves that the improved method solves the failure problem of the original method in scenarios with drastic fluctuations and significantly improves the model's adaptability to rapid water level changes.

[0102] To deeply analyze the optimization mechanism of the improved inverse modeling method on GNSS-IR water level inversion error, STL time series decomposition was performed on the inversion residuals of stations A (gentle water level fluctuation) and B (severe water level fluctuation). The inversion residuals were decomposed into trend, seasonal, and residual residual terms. The differences in error composition between the traditional and improved methods were compared, and the results are as follows: Figure 8 As shown, in a scenario with gentle water levels (A), the traditional method exhibits cumulative drift and periodic residues in the residuals, resulting in drastic residual fluctuations. The improved method effectively suppresses trend errors and significantly reduces periodic residues, leading to more stable residuals. In a scenario with drastic water level fluctuations (B), the traditional method shows drastic fluctuations in the trend and seasonal terms, resulting in catastrophic divergence of the residuals. The improved method completely eliminates systematic trends and periodic residues, stabilizing the residuals near zero. These results, from the perspective of error decomposition, confirm that random walk and time-varying physical constraints can significantly optimize the GNSS-IR inversion error structure, improving inversion accuracy and long-term stability in both gentle and drastic water level scenarios.

[0103] Combining the water level inversion results from both stations, it can be concluded that the improved inverse modeling method, by optimizing the observation model and parameter update mechanism (including highly random walk modeling, adaptive time-varying terms, and dynamic iteration strategies), not only further improves the inversion accuracy at station A where the water level fluctuates gently, but more importantly, it overcomes the adaptability limitations of the original method at station B where the water level fluctuates dramatically. This means that the improved method is no longer limited by the amplitude of water level fluctuations, possesses stronger scenario universality, and provides a more reliable solution for the application of GNSS-IR water level inversion technology in complex hydrological scenarios.

[0104] To illustrate the technical advantages of the method of this invention, a brief analysis of the traditional inverse modeling method based on B-spline functions is presented below. The principle of this method can be found in existing literature; its core is to represent the antenna height as a linear combination of B-spline basis functions, and obtain a continuous height sequence by estimating the control point weights.

[0105] However, this method has an inherent parameter sensitivity problem. Figure 2 This is a schematic diagram of second-order B-spline basis functions. Figure 3 The comparison of inversion results under different numbers of B-spline control points is shown. Figure 4 The figures show the RMSE, MAD, and correlation coefficient curves for different numbers of control points. It can be seen that: too few control points (e.g., 5) lead to underfitting, failing to capture rapid water level fluctuations; too many control points (e.g., 50, 100) lead to overfitting, causing violent oscillations in the inversion results or even solution failure; the number of control points and node spacing lack universally applicable selection criteria, requiring parameter readjustment for different scenarios, and the distribution pattern (e.g., quasi-uniform distribution) is difficult to adapt to complex scenarios where rapid and gradual water level changes coexist.

[0106] In addition, the traditional inverse modeling method sets the initial phase of the signal-to-noise ratio oscillation term to a constant, ignoring the time-varying influence of environmental factors on the reflection phase, which leads to a sharp increase in error or even failure in scenarios with drastic water level changes (such as station B) (see the results of the traditional method in Table 1).

[0107] In contrast, the method of this invention uses a random walk model instead of a B-spline function, eliminating the need to select the number of control points and node spacing. Simultaneously, it decomposes the initial phase into an inter-frequency deviation constant term and a common time-varying random walk term, and introduces an iterative adaptive strategy based on MAD (Multi-Aspect Ratio). Therefore, the method of this invention fundamentally overcomes the shortcomings of traditional methods, possessing advantages such as parameter adaptation, independence from refined priors, and strong universality.

[0108] This embodiment describes a GNSS-IR inverse modeling method constrained by random walk and time-varying physical models. Existing primitive inverse modeling methods based on B-spline functions suffer from limitations such as the lack of universality in control point selection criteria, insufficient dynamic adaptability of nodes, difficulty in accurately depicting dramatic water level fluctuations in nearshore waters, and heavy reliance on external prior data. To address these issues, this embodiment employs a random walk model to characterize the continuous dynamic change of antenna height. Simultaneously, a time-varying physical model is constructed, decomposing the initial phase of the signal-to-noise ratio oscillation term into an inter-frequency deviation constant term and a common time-varying random walk term, along with an iterative solution strategy, significantly reducing the dependence on prior data. Finally, comparative experiments are conducted in typical water areas with different water level fluctuation characteristics, verifying the stable and high-precision inversion capability of the proposed method under various hydrological scenarios, demonstrating its scenario universality and engineering application value.

[0109] Based on the same inventive concept as the aforementioned method embodiments, this invention also provides a GNSS-IR inverse modeling system constrained by random walk and time-varying physical model, the system comprising the following modules:

[0110] Data acquisition module: This module acquires the signal-to-noise ratio (SNR) observations collected by the GNSS receiver and preprocesses these SNR observations to extract the oscillation term of the reflected signal. Specifically, this module determines the azimuth and elevation mask ranges based on the first Fresnel reflection region around the station, extracts the reflected signal arc segments from the water surface, performs polynomial fitting on the extracted SNR sequence to remove trend terms, and corrects for elevation curvature error and tropospheric delay error.

[0111] Antenna Height Dynamic Modeling Module: Used to describe the change of antenna height over time using a random walk model, representing the antenna height of the current epoch as the antenna height of the previous epoch plus random walk noise following a zero-mean Gaussian distribution.

[0112] Time-varying phase modeling module: used to split the initial phase in the oscillation term of the reflected signal into an inter-frequency deviation term and a common time-varying term, wherein the inter-frequency deviation term is a constant term independent of different frequency bands, and the common time-varying term is described by a random walk model to describe its change with time.

[0113] Loss function construction module: used to construct a nonlinear least squares loss function for the signal-to-noise ratio observation based on the antenna height sequence constructed by the random walk model, the inter-frequency deviation term, and the common time-varying term.

[0114] Iterative solution module: Used to update antenna height parameters, inter-frequency offset parameters and common time-varying parameters in each iteration using an iterative solution strategy until the convergence condition is met, and output the antenna height inversion results for each epoch.

[0115] Furthermore, the iterative solution module includes:

[0116] Initialization unit: Used to set the initial value of antenna height to a preset coarse value in the first iteration, set the standard deviation factor of antenna height and the standard deviation factor of random walk of common time-varying term to preset initial values, and set the initial values ​​of inter-frequency deviation parameter and reflected signal amplitude parameter to zero.

[0117] Parameter update unit: Starting from the second iteration, it uses the antenna height results of each epoch output from the previous iteration as the initial value of the antenna height for the current iteration, and calculates the standard deviation factor of the antenna height for the current iteration based on the median absolute deviation of the antenna height results from the previous iteration.

[0118] Convergence Judgment Unit: Used to determine whether the root mean square error of the antenna height inversion results in the two rounds is less than a preset threshold or whether the preset maximum number of iterations has been reached.

[0119] The specific workflow of each module in this embodiment corresponds to the steps in the aforementioned method embodiments, and will not be repeated here. This system can achieve adaptive high-precision antenna height (water level) inversion without requiring fine prior information, and is suitable for scenarios with both gentle and drastic water level changes.

[0120] Based on the same inventive concept as the foregoing method embodiments, this invention also provides an electronic device, including a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the program, it implements the GNSS-IR inverse modeling method with random walk and time-varying physical model constraints described in the above method embodiments.

[0121] The electronic device may be a computing unit built into a GNSS receiver, an edge computing device, a server, or a personal computer, etc. When the processor executes the program, it preprocesses the signal-to-noise ratio observations, constructs a random walk model, constructs a time-varying phase model, constructs a loss function, iteratively solves the problem, and outputs the antenna height inversion result, according to the steps in the method embodiment.

[0122] Based on the same inventive concept as the foregoing method embodiments, this embodiment of the invention also provides a computer-readable storage medium storing a computer program thereon. When executed by a processor, this program implements the GNSS-IR inverse modeling method with random walk and time-varying physical model constraints described in the above method embodiments.

[0123] The storage medium includes, but is not limited to, non-volatile storage media such as ROM, RAM, hard disk, solid-state drive, USB flash drive, and optical disk. This storage medium can be installed in GNSS monitoring equipment or data processing terminals to store computer programs that implement the method of this invention.

[0124] In summary, this invention has analyzed in detail the limitations of traditional inverse modeling methods based on B-spline functions, including the lack of universality in selecting the number of B-spline control points, insufficient dynamic adaptation capability of node spacing, and the potential for calculation errors due to improper parameter settings. Piecewise smooth B-spline functions are difficult to accurately characterize irregular and drastic dynamic changes in water level. Therefore, this invention employs a random walk model to characterize the continuous dynamic change process of antenna height. Simultaneously, the initial phase of the signal-to-noise ratio oscillation term is decomposed into an inter-frequency deviation constant term representing the difference between frequency bands and a common time-varying term representing the dynamic changes in the environment. The latter is modeled as a random walk process, and an iterative solution strategy is introduced to achieve accurate solution of model parameters, enabling high-precision antenna height inversion without relying on fine prior parameters. Experimental verification at two typical monitoring stations, Station A (stable water level) and Station B (drastic water level fluctuation), shows that the root mean square error (RMSE) of the inversion results of the method of this invention is 0.94 cm and the correlation coefficient is 0.9953 in the scenario of stable water level, which is 13.0% higher than the accuracy of the traditional method. In the scenario of drastic water level fluctuation, the traditional method has failed, but the method of this invention still achieves excellent results with an RMSE of 2.14 cm and a correlation coefficient of 0.9992, effectively breaking through the limitation of the traditional method in adapting to the amplitude of water level fluctuation, and proving the universality of the method of this invention in various scenarios and its engineering application value.

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

Claims

1. A GNSS-IR inverse modeling method with random walk and time-varying physical model constraints, characterized in that, include: Acquire signal-to-noise ratio (SNR) observation data collected by a GNSS receiver, and preprocess the SNR observation data to extract the oscillation component of the reflected signal; A random walk model is used to describe the change of antenna height over time. The antenna height in the current epoch is represented as the antenna height in the previous epoch plus the random walk noise that follows a zero-mean Gaussian distribution. The initial phase of the oscillation component of the reflected signal is split into an inter-frequency deviation term and a common time-varying term. The inter-frequency deviation term is a constant term that is independent of different frequency bands, and the common time-varying term is described by a random walk model to describe its change with time. Based on the antenna height sequence constructed by the random walk model, the inter-frequency deviation term, and the common time-varying term, a nonlinear least squares loss function for the signal-to-noise ratio observation is constructed. An iterative solution strategy is adopted, in which the antenna height parameters, inter-frequency deviation parameters and common time-varying parameters are updated in each iteration until the convergence condition is met, and the antenna height inversion results of each epoch are output.

2. The GNSS-IR inverse modeling method with random walk and time-varying physical model constraints according to claim 1, characterized in that, The iterative solution strategy includes: In the first iteration, the initial value of the antenna height is set to a preset coarse value, and the standard deviation factor of the antenna height and the standard deviation factor of the random walk of the common time-varying term are set to preset initial values. The initial values ​​of the inter-frequency deviation parameter and the reflected signal amplitude parameter are set to zero. Starting from the second iteration, the antenna height results of each epoch output from the previous iteration are used as the initial value of the antenna height for the current iteration, and the standard deviation factor of the antenna height for the current iteration is calculated based on the median absolute deviation of the antenna height results from the previous iteration. Repeat the iteration until the root mean square error of the antenna height inversion results in the two rounds is less than the preset threshold or the preset maximum number of iterations is reached.

3. The GNSS-IR inverse modeling method with random walk and time-varying physical model constraints according to claim 1, characterized in that, The common time-varying term is shared across all frequency bands and is used to characterize the global time-varying characteristics of the reflection phase driven by environmental factors; the inter-frequency deviation term is used to characterize the fixed phase shift caused by differences in wavelength and reflection characteristics in different GNSS frequency bands.

4. The GNSS-IR inverse modeling method with random walk and time-varying physical model constraints according to claim 1, characterized in that, The variance of the random walk noise is proportional to the time difference between the current epoch and the previous epoch, and the variance is expressed as: ,in The standard deviation factor of the antenna height. For time difference; When there are missing observation data, the variance is adjusted according to the actual time interval.

5. The GNSS-IR inverse modeling method with random walk and time-varying physical model constraints according to claim 1, characterized in that, The random walk standard deviation factor of the common time-varying term is calculated based on the median absolute deviation of the common time-varying term sequence obtained in the previous iteration, which is consistent with the calculation method of the standard deviation factor of the antenna height.

6. The GNSS-IR inverse modeling method with random walk and time-varying physical model constraints according to claim 1, characterized in that, The nonlinear least squares loss function consists of three parts: the sum of squared residuals between the signal-to-noise ratio observations and the model predictions, the constraint term for the random walk of the antenna height, and the constraint term for the random walk of the common time-varying term. The Gauss-Newton method is used to iteratively solve the loss function.

7. A GNSS-IR inverse modeling system with random walk and time-varying physical model constraints, characterized in that, include: The data acquisition module is used to acquire signal-to-noise ratio (SNR) observation data collected by the GNSS receiver, and to preprocess the SNR observation data to extract the oscillation component of the reflected signal. The antenna height dynamic modeling module is used to describe the change of antenna height over time using a random walk model. It represents the antenna height of the current epoch as the antenna height of the previous epoch plus random walk noise that follows a zero-mean Gaussian distribution. The time-varying phase modeling module is used to decompose the initial phase of the oscillation component of the reflected signal into an inter-frequency deviation term and a common time-varying term, wherein the inter-frequency deviation term is a constant term independent of different frequency bands, and the common time-varying term is described by a random walk model to describe its change with time. The loss function construction module is used to construct a nonlinear least squares loss function for the signal-to-noise ratio observation based on the antenna height sequence constructed by the random walk model, the inter-frequency deviation term, and the common time-varying term; The iterative solution module is used to update the antenna height parameters, inter-frequency deviation parameters, and common time-varying parameters in each iteration using an iterative solution strategy until the convergence condition is met, and outputs the antenna height inversion results for each epoch.

8. The GNSS-IR inverse modeling system with random walk and time-varying physical model constraints according to claim 7, characterized in that, The iterative solution module includes: The initialization unit is used to set the initial value of the antenna height to a preset coarse value in the first iteration, set the standard deviation factor of the antenna height and the standard deviation factor of the random walk of the common time-varying term to preset initial values, and set the initial values ​​of the inter-frequency deviation parameter and the reflected signal amplitude parameter to zero. The parameter update unit is used to take the antenna height results of each epoch of the previous iteration as the initial value of the antenna height in the current iteration, starting from the second iteration, and calculate the standard deviation factor of the antenna height in the current iteration based on the median absolute deviation of the antenna height results in the previous iteration. The convergence judgment unit is used to determine whether the root mean square error of the antenna height inversion results of the two rounds is less than a preset threshold or whether the preset maximum number of iterations has been reached.

9. An electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the processor executes the program, it implements the GNSS-IR inverse modeling method with random walk and time-varying physical model constraints as described in any one of claims 1 to 6.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that, When executed by the processor, the program implements the GNSS-IR inverse modeling method with random walk and time-varying physical model constraints as described in any one of claims 1 to 6.

Images