GNSS-IR real-time water level inversion method and system based on RB particle filtering

Abstract

The application discloses a GNSS-IR real-time water level inversion method and system based on RB particle filtering. The application constructs a Rao-Blackwellized particle filtering state space model, decomposes a water level inversion state vector into a nonlinear substate and a linear substate, and estimates the nonlinear substate and the linear substate by using particle filtering and Kalman filtering respectively. In each signal-to-noise ratio sampling epoch, the nonlinear substate is sequentially estimated by particle filtering, and the linear substate is analytically estimated by Kalman filtering. The particle weight is updated by a Gaussian likelihood function constructed by an observation value of a current epoch. The filtering recursion only depends on the observation value of the current epoch and a posterior estimation result of a previous epoch, and real-time water level inversion is realized epoch by epoch. While maintaining the robustness of inversion of a strong nonlinear system, the application significantly reduces the calculation cost, and meets the real-time water level monitoring demand of a shore-based GNSS receiver with a high sampling frequency.

Description

Technical Field

[0001] This invention relates to the field of global navigation satellite system reflection measurement technology, and in particular to a GNSS-IR real-time water level inversion method and system based on RB particle filtering, which is applicable to continuous, automatic, and real-time water level monitoring of water bodies such as rivers, lakes, reservoirs, and nearshore sea areas. Background Technology

[0002] Global Navigation Satellite System (GNSS) interferometric reflectance measurement technology utilizes the interference effect formed at the receiver antenna between the direct GNSS signal and the signal reflected from the water surface. It extracts the oscillation frequency of the reflected signal from the signal-to-noise ratio (SNR) observation data, and inversely calculates the vertical distance from the antenna phase center to the water surface. Combined with a GNSS antenna elevation reference, it achieves non-contact, high-precision monitoring of absolute water levels. Compared to traditional water level gauges, GNSS-IR technology can monitor water levels without contact with the water body, making it suitable for scenarios where traditional water level gauges are difficult to deploy, such as shallow waters and muddy areas. Compared to satellite altimetry, GNSS-IR can achieve long-term, continuous, and uninterrupted observations using shore-based receivers. Compared to satellite gravity inversion technology, GNSS-IR can accurately capture local dynamic changes in water levels along the shoreline, filling the gap in near-shore area monitoring capabilities. The four major global navigation satellite systems—GPS, BDS, GLONASS, and Galileo—along with multiple regional navigation satellite systems, constitute a signal network containing hundreds of satellites and dozens of frequency bands, providing abundant redundant observation data for GNSS-IR.

[0003] Currently, GNSS-IR water level inversion methods are mainly divided into two categories. The first category is based on spectrum analysis, which uses Lomb-Scargle spectral analysis to perform spectral transformation on the detrended signal-to-noise ratio (SNR) reflection signal sequence, extracts the dominant frequency component, and then calculates the reflector height. To improve the robustness of spectrum analysis, some studies introduce empirical mode decomposition to preprocess the SNR sequence to separate the effective reflection signal from the noise component. Other studies use wavelet time-frequency analysis to replace Lomb-Scargle spectral analysis to improve time-frequency resolution. However, spectrum analysis methods are essentially batch post-processing methods, requiring a certain length of satellite arc observation data to complete an effective inversion. The observation time for a single arc is usually 2 to 3 hours, resulting in low temporal resolution and discontinuous output time in the inversion results. To address the aforementioned issues, some studies have introduced a sliding window mechanism to compress the effective data span from the traditional 2 to 3 hours to as short as 15 minutes, achieving near real-time water level inversion with centimeter-level accuracy. However, the sliding window is essentially still a batch processing mechanism, and the output frequency is constrained by the window length, making it impossible to truly achieve real-time calculation per epoch. Furthermore, the spectral resolution decreases when the window length is shortened, making it difficult to balance real-time performance and inversion accuracy.

[0004] The second category is based on modeling and fitting methods. Some studies model water level changes as time-varying parameter B-spline curves, constructing an inverse modeling optimization problem using the signal-to-noise ratio reflected signal observation equation. The B-spline control point parameters are then solved using nonlinear least squares estimation to obtain a smooth and continuous water level change sequence. Building on this, some studies further introduce unscented Kalman filtering to replace nonlinear least squares estimation, expanding the B-spline control point parameters into a state vector, thus achieving epoch-by-epoch recursive updates of the model parameters.

[0005] However, existing real-time solutions still have significant limitations: near-real-time methods based on spectrum analysis are constrained by the sliding window, making it difficult to balance real-time performance and inversion accuracy; real-time methods based on modeling and fitting, while capable of epoch-by-epoch solutions, suffer from sensitivity to initial parameter values ​​and are prone to filter divergence in high-dynamic water level scenarios. Gaussian filtering methods using unscented Kalman filtering are only suitable for mildly nonlinear, Gaussian noise scenarios, and are prone to filter divergence and a sharp drop in inversion accuracy in extreme environments with sudden water level changes and strong multipath interference. While standard particle filtering can adapt to strongly nonlinear non-Gaussian systems, it suffers from severe particle degradation, requires a large number of particles to ensure accuracy, and incurs high computational costs, failing to meet the real-time inversion requirements of shore-based GNSS stations at frequencies of 1Hz or even higher.

[0006] In summary, existing GNSS-IR water level inversion methods struggle to simultaneously balance the inversion robustness of strongly nonlinear non-Gaussian systems with the computational real-time requirements for epoch-by-epoch output. Summary of the Invention

[0007] To address the shortcomings of existing GNSS-IR water level inversion methods in balancing the robustness of inversion for strongly nonlinear non-Gaussian systems with the real-time computational requirements for epoch-by-epoch output, this invention provides a GNSS-IR real-time water level inversion method and system based on Rao-Blackwellized particle filtering. By constructing a Rao-Blackwellized particle filter state-space model, the water level inversion state vector is decomposed into nonlinear and linear sub-states, which are then estimated in parallel using particle filtering and Kalman filtering, respectively. This achieves high computational efficiency for epoch-by-epoch real-time water level inversion while maintaining the robustness of inversion for strongly nonlinear systems.

[0008] According to one aspect of this invention, a GNSS-IR real-time water level inversion method based on Rao-Blackwellized particle filtering is provided, comprising: acquiring signal-to-noise ratio (SNR) observation data and satellite ephemeris data from a GNSS receiver; preprocessing the SNR observation data and satellite ephemeris data to obtain a processed SNR sequence of reflected signals and satellite elevation angle data; constructing a state-space model based on Rao-Blackwellized particle filtering, decomposing the state vector of GNSS-IR water level inversion into nonlinear sub-states and linear sub-states, wherein the nonlinear sub-states include reflector height and reflector surface roughness parameters, and the linear sub-states include in-phase component coefficients and quadrature component coefficients of reflected signals in each frequency band; at each SNR sampling epoch, sequentially estimating the nonlinear sub-states using particle filtering, and analytically estimating the linear sub-states using Kalman filtering, with the two filtering operations performed in parallel and the particle weights updated using a Gaussian likelihood function calculated from the current epoch observation values; weighted averaging of the state estimates of all particles according to their respective particle weights, outputting the optimal estimate of the reflector height for the current epoch, and outputting water level values ​​epoch by epoch based on a known GNSS antenna elevation reference.

[0009] As a further technical solution, a state-space model is constructed based on Rao-Blackwellized particle filtering. The state vector of GNSS-IR water level inversion is decomposed into nonlinear sub-states and linear sub-states, including: based on the Rao-Blackwellization theorem, the reflector height and reflector surface roughness parameters are taken as nonlinear sub-states, and their state transition adopts a random walk model. The prediction particles for the next epoch are generated by superimposing process noise; the in-phase component coefficients and orthogonal component coefficients of the reflected signals of each frequency band are taken as linear sub-states. Under the given nonlinear sub-state conditions, the posterior distribution of the linear sub-states follows a Gaussian distribution, and the in-phase component coefficients and orthogonal component coefficients of each frequency band are estimated using Kalman filters.

[0010] As a further technical solution, the filtering recursion for each sampling epoch is executed in the following order: based on the posterior state estimate of the previous epoch, the state prediction of the current epoch is completed through the state transition model; using the preprocessed signal-to-noise ratio observation of the reflected signal in the current epoch, sequential Kalman filtering update and particle weight update are executed in sequence; the weighted average of all updated particle states is calculated to output the optimal estimate of the reflector height in the current epoch; it is determined whether to trigger resampling, and a particle set is prepared for the filtering recursion of the next epoch; the data input for the filtering recursion is limited to the signal-to-noise ratio observation of the current epoch and the posterior state estimate result of the previous epoch, excluding the participation of observation data from future epochs, forming a causal recursive structure for each epoch.

[0011] As a further technical solution, the particle weights are updated using the Gaussian likelihood function calculated from the current epoch observations. This includes: calculating the prediction residual for each reflected signal-to-noise ratio observation in the current epoch based on the updated linear substate posterior estimate; constructing a Gaussian observation likelihood function based on the prediction residuals of all observations, and calculating the likelihood value of each nonlinear particle in the current epoch; multiplying the particle weights of the previous epoch by the likelihood value of the current epoch to obtain the unnormalized weights of the current epoch; and normalizing the unnormalized weights of all particles.

[0012] As a further technical solution, in the observation equation of the signal-to-noise ratio sequence of the reflected signal, the reflector height and the sine value of the satellite elevation angle form a nonlinear functional relationship, the in-phase component coefficients and quadrature component coefficients of each frequency band form a linear combination relationship, and the surface roughness parameter is introduced into the observation equation in the form of exponential decay.

[0013] As a further technical solution, the construction of the state space model also includes: setting an initial particle distribution for the reflector height in the nonlinear substate based on the prior height range of the GNSS receiver antenna phase center distance from the water surface, wherein the initial distribution follows a Gaussian distribution; setting the initial value of the reflector surface roughness parameter to zero, setting the initial mean of the in-phase component coefficient and the quadrature component coefficient of each frequency band to zero, and setting the initial covariance to a preset value; and determining the process noise variance of each nonlinear substate and linear substate respectively through the first-order difference statistical results of historical observation data.

[0014] As a further technical solution, the signal-to-noise ratio observation data of the GNSS receiver comes from multiple satellite systems and multiple frequency bands. The linear sub-state includes the in-phase component coefficient and quadrature component coefficient corresponding to each frequency band. The number of components in each frequency band is dynamically determined according to the number of effective satellite-frequency band combinations received in the current epoch. The observation noise variance of each frequency band is independently estimated based on the statistical characteristics of the signal-to-noise ratio observation residuals of different frequency bands.

[0015] As a further technical solution, after updating the particle weights for the current epoch, the number of effective particles is calculated: , Where N is the total number of particles. The normalized weight of the i-th particle;

[0016] When the number of effective particles is less than half of the total number of particles, a resampling operation is triggered. The resampling adopts a low-variance resampling method, which eliminates low-weight particles and replicates high-weight particles. After resampling, a small noise is superimposed on the new particles, and the weights of all particles are reset to equal values.

[0017] As a further technical solution, the number of nonlinear particles in the Rao-Blackwellized particle filter is 1500, and the average calculation time per epoch in the parallel filter is in the millisecond range and is less than the sampling interval of the signal-to-noise ratio data.

[0018] According to one aspect of this specification, a GNSS-IR real-time water level inversion system based on Rao-Blackwellized particle filtering is provided, comprising: a data acquisition module for acquiring signal-to-noise ratio (SNR) observation data and satellite ephemeris data from a GNSS receiver; a preprocessing module for preprocessing the SNR observation data and satellite ephemeris data to obtain a processed reflection signal SNR sequence and satellite elevation angle data; and a state decomposition module for constructing a state-space model based on Rao-Blackwellized particle filtering, decomposing the state vector of GNSS-IR water level inversion into nonlinear sub-states and linear sub-states, wherein the nonlinear sub-states include reflection signals. The system includes parameters for reflector height and reflector surface roughness. The linear substate includes the in-phase component coefficients and quadrature component coefficients of the reflected signals in each frequency band. A parallel filtering estimation module is used to perform sequential estimation of the nonlinear substate using particle filtering and analytical estimation of the linear substate using Kalman filtering at each signal-to-noise ratio sampling epoch. The two filtering methods are executed in parallel, and the particle weights are updated using the Gaussian likelihood function calculated from the current epoch observations. An epoch-by-epoch water level output module is used to perform a weighted average of the state estimates of all particles based on their respective particle weights, output the optimal estimate of the reflector height for the current epoch, and output the water level value epoch-by-epoch based on the known GNSS antenna elevation reference.

[0019] Compared with the prior art, the beneficial effects of the present invention are as follows: 1. In this invention, by constructing a Rao-Blackwellized particle filter state space model, the water level inversion state vector is decomposed into nonlinear sub-states and linear sub-states, and particle filtering and Kalman filtering are used for parallel estimation respectively. Only a small number of particles are needed to estimate the nonlinear height parameter, and the linear parameter is estimated in parallel through Kalman filtering. This can improve the inversion robustness while reducing the computational load by several orders of magnitude.

[0020] 2. In this invention, by employing Rao-Blackwellized particle filtering, the number of particles required is significantly reduced while ensuring robustness of inversion in extreme scenarios, thereby lowering computational complexity and ultimately achieving high-precision, high-real-time epoch-by-epoch water level inversion. Experimental results show that when the number of particles is 1500, the algorithm achieves RMSEs of 1.26 cm and 5.88 cm between the inversion results and the measured water levels in two types of stations with significantly different dynamic water level characteristics, respectively, with correlation coefficients exceeding 99%. The average computation time per epoch is only 6 ms, far lower than the data sampling interval, and it can stably track the dynamic changes of water level throughout the entire process.

[0021] 3. In this invention, the real-time inversion results are highly consistent with the measured data from the water level gauge and the high-precision inversion results obtained after the event, throughout the entire experimental period. The water level rise and fall trends, the timing of peaks and troughs, and their amplitudes are all precisely matched, with no obvious phase delay, amplitude distortion, or trend deviation. At high-dynamic stations where the daily water level change exceeds 3m, the real-time algorithm completely reproduces the full-cycle characteristics of the water level change. The peak and trough positions, rise and fall rates are consistent with the measured data and the after-the-fact inversion results, with no trend deviation.

[0022] 4. In this invention, the residual sequences retrieved in real-time from both stations exhibit a random fluctuation distribution centered at 0, without any positive or negative trend drift. The residual fluctuation range of the GTGU station is basically controlled within the range of -5cm to +5cm, with over 99% of the residuals within ±3cm. The mean residual is 0.21cm, and the standard deviation is 1.38cm, indicating that the residuals approximately follow a zero-mean normal distribution without systematic bias. The residual fluctuation range of the SC02 station is within the range of -20cm to +20cm, with nearly 90% of the residuals within ±10cm. The mean residual is 0.35cm, and the standard deviation is 5.88cm, which is close to a zero-mean normal distribution, indicating that even under severe water level fluctuations, no systematic bias will occur.

[0023] 5. In this invention, the real-time water level inversion can converge to the vicinity of the reference water level value within only 10 to 15 epochs after the start, indicating that the algorithm can converge quickly. Attached Figure Description

[0024] 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.

[0025] Figure 1 A schematic diagram of the overall process of a GNSS-IR real-time water level inversion method based on RB particle filtering provided in an embodiment of the present invention; Figure 2 This is a graph showing the water level inversion accuracy and average calculation time per epoch under different particle numbers in this embodiment of the invention. Figure 3 The figure shows the experimental comparison results of improved inverse modeling (post-hoc water level inversion) and RB particle filtering (real-time water level inversion) for the GTGU and SC02 stations in this embodiment of the invention. Figure 4 This is a residual distribution diagram of the real-time water level inversion results for stations GTGU and SC02 in this embodiment of the invention. Detailed Implementation

[0026] In existing technologies, both spectral analysis and inverse modeling methods rely on long-term data accumulation, making epoch-by-epoch real-time inversion impossible. Gaussian filtering methods using unscented Kalman filtering are only suitable for mildly nonlinear, Gaussian noise scenarios, and are prone to filter divergence and a sharp drop in inversion accuracy under extreme environments with sudden water level changes and strong multipath interference. While standard particle filtering can adapt to strongly nonlinear non-Gaussian systems, it suffers from severe particle degradation, requires a large number of particles to ensure accuracy, and incurs high computational costs, failing to meet the real-time inversion requirements of shore-based GNSS stations at frequencies of 1Hz or higher. To address these issues, this invention proposes a GNSS-IR real-time water level inversion method based on an RB particle filter. This method performs marginal decomposition of the state vector for GNSS-IR water level inversion, separating the strongly nonlinear sea surface height parameter and surface roughness from the linear signal amplitude and phase parameters, and then uses particle filtering and Kalman filtering in parallel estimation for each. While ensuring robustness of inversion in extreme scenarios, this method can significantly reduce the number of particles required and lower computational complexity, ultimately achieving high-precision, high-real-time epoch-by-epoch water level inversion.

[0027] 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.

[0028] like Figure 1 As shown, this embodiment of the invention provides a GNSS-IR real-time water level inversion method based on RB particle filtering, including the following steps:

[0029] Step 1: Data Acquisition and Preprocessing.

[0030] Acquire signal-to-noise ratio observation data and satellite ephemeris data from the GNSS receiver.

[0031] Input GNSS observation data and corresponding satellite ephemeris data. Preprocess the GNSS observation data and ephemeris data input at the current epoch, performing satellite orbit calculation, receiver position calculation, elevation angle calculation, elevation angle curvature correction, relative tropospheric delay correction, signal-to-noise ratio (SNR) data extraction, and direct component removal (detrending). Simultaneously, perform SNR data validity screening (angle masking) and outlier removal, outputting the effective reflected signal SNR observations and satellite elevation angle data for the current epoch across multiple satellites and frequency bands, providing reliable observation input for subsequent filtering calculations. To reduce redundant unknown parameters, the polynomial coefficients of the direct component of the satellite-frequency band are first calculated using historical observation data to directly obtain the reflected signal SNR (dSNR) data from the SNR observations.

[0032] Step 2: Construct a state-space model based on Rao-Blackwellized particle filtering.

[0033] A state-space model is constructed based on Rao-Blackwellized particle filtering, and the state vector of GNSS-IR water level inversion is decomposed into nonlinear sub-states and linear sub-states.

[0034] RB Particle Filter (RBPF) is an algorithm that combines particle filtering and Kalman filtering. Its core theoretical basis is the Rao-Blackwellization statistical theorem: for a decomposable state vector, its estimated variance can be reduced through marginal decomposition, that is, the estimated variance after decomposition is always no greater than the estimated variance of the full-state particle filter.

[0035] RBPF decomposes the state vector into two parts: a nonlinear substate and a linear substate, denoted as ,in It is a nonlinear substate. It is a linear substate. It is the observation time (epoch).

[0036] Based on Bayesian estimation theory, the posterior probability distribution of the state is:

[0037]

[0038] In the formula, It contains all historical measurement information and observation data up to the current time k; It is the posterior distribution of the nonlinear substates, estimated using particle filtering; It is the posterior distribution of the linear substate when the nonlinear substate is known. It is assumed that the distribution is Gaussian and is estimated analytically using Kalman filtering without consuming particle resources.

[0039] In GNSS-IR, reflector height and surface roughness are slowly time-varying, simplifying to a random walk model. Nonlinear substates. Where h is the height of the reflector, For the roughness parameter of the reflecting surface. Linear substate. = ( It is the number of frequency bands, including the in-phase component coefficients and quadrature component coefficients of the reflected signals in each frequency band.

[0040] Oscillation term of signal-to-noise ratio The physical model (i.e., the observation equation for the signal-to-noise ratio of the reflected signal) is as follows:

[0041]

[0042] Where dSNR is the signal-to-noise ratio of the detrended reflected signal, h is the reflector height, and θ is the satellite elevation angle. L=2π / λ represents the wavenumber, and λ is the carrier wavelength. Here, C1 and C2 are the roughness parameters of the reflecting surface, and C1 and C2 are the coefficients of the in-phase component and the orthogonal component, respectively. A represents the amplitude. For phase, the amplitude parameter A and phase parameter corresponding to each frequency band. Linear substates are constructed, and Kalman filtering is used to estimate them given nonlinear substates.

[0043] Compared to standard particle filtering, which requires a large number of particles to estimate the full state variables, RBPF only requires a small number of particles to estimate the nonlinear height parameters. The linear parameters are estimated in parallel through Kalman filtering, which can improve the robustness of the inversion while reducing the computational load by several orders of magnitude, thus meeting the requirements of real-time water level inversion.

[0044] Step 3: RBPF filter initialization.

[0045] Only It is executed once at time 1 to generate the initial nonlinear substate particles and their corresponding weights, as well as the Kalman filter for each linear substate.

[0046] Sampling is required during initialization One particle, in this embodiment, is taken as... The prior distribution of the nonlinear substate is used to collect particles in the nonlinear substate. This embodiment assumes that the initial height follows a Gaussian distribution, the surface roughness follows a uniform distribution, random particle sampling is performed, and all particles have equal weights, denoted as .

[0047]

[0048] In addition, the Kalman filter for the linear substates is initialized. For each particle, the mean and covariance of the linear substates are initialized.

[0049]

[0050] Considering Therefore, set and The initial mean of all is 0, and the initial variance of all is 400.

[0051] Step 4: Epoch-by-epoch filtering recursion.

[0052] At each signal-to-noise ratio sampling epoch, the nonlinear substate is sequentially estimated using particle filtering, and the linear substate is analytically estimated using Kalman filtering. The two filtering methods are executed in parallel, and the particle weights are updated using the Gaussian likelihood function calculated from the current epoch observations.

[0053] The filtering recursion for each sampling epoch is executed in the following order: based on the posterior state estimate of the previous epoch, the state prediction for the current epoch is completed through a state transition model; using the preprocessed signal-to-noise ratio (SNR) observation of the reflected signal at the current epoch, sequential Kalman filtering update and particle weight update are performed sequentially; a weighted average of all updated particle states is calculated to output the optimal estimate of the reflector height for the current epoch; it is determined whether resampling is triggered, and a particle set is prepared for the filtering recursion of the next epoch. The data input for the filtering recursion is limited to the SNR observation of the current epoch and the posterior state estimate of the previous epoch, excluding the participation of observation data from future epochs, forming a causal recursive structure for each epoch.

[0054] Each epoch of RBPF filtering can be divided into four steps: prediction, update, state estimation, and resampling.

[0055] (1) Prediction:

[0056] from epochal state prediction The epoch states are divided into nonlinear substate predictions and linear substate predictions using Kalman comprehension. Process noise is added to each nonlinear particle to generate... Predicted particles for each epoch:

[0057]

[0058] In GNSS-IR, both reflector height and surface roughness are slowly time-varying, simplifying to a random walk model. ); It is nonlinear process noise, derived from a Gaussian distribution. Random sampling is performed. It is the process noise variance of the nonlinear substate.

[0059] For each linear substate particle, the state of the next epoch is predicted using a Kalman filter:

[0060]

[0061] in, It is the transition matrix of the linear substate. In GNSS-IR, the linear substate is also slowly time-varying, therefore ( (It is the identity matrix); It is the process noise matrix of the linear substate, and it is a diagonal matrix.

[0062] (2) Update:

[0063] The weight update consists of three sub-steps: sequential Kalman filtering, likelihood function calculation, and weight update.

[0064] (2.1) Sequential Kalman Filtering: Based on the preprocessed signal-to-noise ratio data and elevation angle data of the reflected signal, for the linear substate of each particle, the preprocessed multi-satellite and multi-band effective observations are used in conjunction with the GNSS-IR reflected signal physical model to calculate the predicted signal-to-noise ratio of the reflected signal for each particle, and the sequential Kalman filter update is performed: the Kalman gain corresponding to each observation is calculated in turn, and the mean and covariance matrix of the linear substate are updated successively to obtain the posterior estimation result of the linear substate corresponding to each particle, which provides the predictive observation basis for the likelihood function calculation.

[0065] Specifically, using Epochal observations The weights of each particle are updated by calculating the residuals. First, the linear sub-states of each particle in the current epoch are estimated using Kalman filtering. For the current... For each observation in an epoch, the corresponding predicted observation can be calculated using the observation equation and rewritten as follows:

[0066]

[0067] in, and Frequency band The corresponding linear substate. Assume the current... Era, frequency band The number of SNR observations is Then the corresponding observation equation for this epoch and frequency band is:

[0068]

[0069] Where H is the observation matrix, V is the residual matrix, and the corresponding observation noise matrix is... for

[0070]

[0071] Since the satellite signals observed by the receiver at each epoch may be different, in order to avoid complex matrix construction and matrix inversion, this embodiment performs sequential Kalman filtering updates for each particle. Each update only needs to select the corresponding linear sub-state based on the frequency band of the observed value. The update process is as follows:

[0072] ① Calculate the Kalman gain:

[0073]

[0074] in, It is the covariance of the linear substates obtained from the prediction step. It is observation noise. It is a scalar (updated only based on a new observation each time).

[0075] ② Calculate the linear sub-states and their corresponding covariances

[0076]

[0077] in, It is an observation vector. It is the covariance of the linear substates obtained from the prediction step. It is the posterior state covariance matrix updated at time k.

[0078] (2.2) Likelihood function calculation: After obtaining the estimated values ​​of the linear sub-states through Kalman filtering, it is necessary to calculate the likelihood values ​​to update the weights of the nonlinear sub-states. Assuming that the observation noise follows a Gaussian distribution, a Gaussian observation likelihood function is constructed, and the likelihood value corresponding to each particle is calculated to quantify the degree of matching between the particle state and the actual observation.

[0079] Assumption Total of epochs If there are 100 observations, the likelihood value is calculated as follows:

[0080]

[0081] in, It is a sequence of observations. It is the first The residuals of each observation.

[0082] (2.3) Based on the likelihood function calculation results, complete the sequential update of particle weights: Based on the particle weight at time step, multiply by the particle likelihood value at the current epoch to obtain... The original weights of the particles at time 1; all weights are normalized to ensure that the sum of the weights of all particles is 1, thus completing the update of the particle weights.

[0083] Specifically, the particle weights are updated and normalized to ensure that the sum of the particle weights is 1.

[0084]

[0085] (3) State estimation:

[0086] Based on the updated particle weights, a weighted average of the nonlinear and linear substates of all particles is performed to obtain the optimal global state estimate at time k, thus completing the output of the real-time water level result:

[0087]

[0088] The state estimates of all particles are weighted and averaged according to their respective weights to output the optimal estimate of the reflector height for the current epoch. The absolute elevation reference of the GNSS antenna is determined by precise single-point positioning or relative positioning, and the water level value for the current epoch is obtained using the formula Hw=Hant−h, where Hant is the elevation of the GNSS antenna phase center and h is the reflector height.

[0089] (4) Resampling:

[0090] After several rounds of calculations, the weights of most particles approach 0, while only a small number of particles have non-zero weights. Direct weighted summation will cause a sharp drop in the accuracy of the estimation results, so resampling is required.

[0091] Calculate the number of effective particles

[0092]

[0093] When the number of effective particles When the number of effective particles in the current epoch is less than half of the total number of particles, a low-variance resampling operation is triggered. This operation eliminates low-weight particles, replicates high-weight particles, and adds small noise to the new particles to maintain particle diversity. After resampling, the weights of all particles are reset to equal values. Specifically, the resampling operation is as follows:

[0094] (4.1) Calculate the cumulative sum of the weights ;

[0095] (4.2) Generate initial random numbers ;

[0096] (4.3) to ,calculate ;

[0097] (4.4) For each Find the smallest satisfy , replicated particles As new particles (including nonlinear substates and linear substates and their corresponding covariance matrices);

[0098] (4.5) The weights of the new particles are reset to ;

[0099] (5) End judgment:

[0100] Determine if new GNSS observation data exists: If new observation data exists, return to the state prediction step and proceed to the filtering iteration process of the next epoch; if no new observation data exists, end the algorithm and complete the water level inversion for the entire observation cycle.

[0101] This embodiment uses measured data from the GTGU and SC02 stations to experimentally verify the real-time water level inversion algorithm based on the RB particle filter. Detailed information on the two GNSS stations and their corresponding reference water level gauges can be obtained using existing technologies. In the initial stage of the experiment, based on the historical observation data from the previous 7 days at each station, second-order polynomial fitting was used to obtain the polynomial coefficients of the direct signal components of each satellite-frequency band, which were then used to extract dSNR data in the real-time calculation. The strategy used for GNSS data preprocessing remained consistent, with the angle mask employing a more accurate angular boundary of the water surface reflection area.

[0102] Regarding the configuration of the RB particle filter algorithm parameters, this embodiment also uses the post-hoc GNSS-IR water level inversion results from the previous 7 days of historical GNSS data. By statistically analyzing the first-order difference variance of reflector height, water surface roughness, and corresponding parameters for each frequency band, reasonable values ​​for relevant algorithm parameters are determined to ensure the accuracy of real-time water level inversion. Based on the above historical statistical results, Table 1 lists the model parameters of the RB particle filter in real-time dynamic inversion. The process noise variance, as in the improved inverse modeling method, is represented by the standard deviation factor. This is because the data sampling interval is not necessarily 1 second, i.e., the actual process noise variance... ( It is the standard deviation factor. (This refers to the data sampling interval). In the particle initialization step, the initial height is the same as the prior height, the roughness of the reflecting surface is set to 0, and the values ​​for each frequency band are... and The initial values ​​are also set to 0, and the initial state variance is 400. Furthermore, since the actual observation noise of the signal-to-noise ratio (SNR) data is unknown, the observation noise for each satellite system and frequency is estimated based on the average residual from the last hour of SNR observations.

[0103] Table 1 Parameters of the RB Particle Filter Real-Time Water Level Inversion Model

[0104]

[0105] To determine the optimal particle number for the RBPF-based GNSS-IR real-time water level inversion algorithm and achieve the best balance between inversion accuracy and real-time computational efficiency, this embodiment uses data from the SC02 station as the experimental object to conduct a particle number sensitivity experiment. The experiment sets the particle number range to 500–3000 and the step size to 200. Curves showing the inversion RMSE accuracy and average computation time per epoch for different particle numbers are plotted. The results are as follows: Figure 2 As shown. To quantify the optimal balance between accuracy and efficiency, the inversion RMSE and average single-epoch time were normalized to the [0,1] interval and then summed. The number of particles corresponding to the minimum comprehensive index was taken as the optimal parameter. The optimal balance point was finally determined to be 1500 particles, corresponding to an inversion RMSE of 5.88 cm and an average single-epoch computation time of only 5.68 ms. In the particle range of 500-1500, the inversion RMSE decreased significantly with the increase of particle number. In the low particle number range, the inversion accuracy was severely insufficient due to particle degradation. At 1500 particles, the RMSE dropped to the centimeter level, entering the accuracy convergence plateau. After the particle number exceeded 1500, the inversion accuracy did not improve significantly, but the computation time increased linearly. When the particle number was too high, redundant sampling would interfere with the posterior estimation, causing accuracy degradation. In summary, 1500 particles maximized the reduction of computational overhead while ensuring centimeter-level inversion accuracy. Therefore, this embodiment selected the optimal particle number of 1500 for experiments to compare and verify the performance of the real-time algorithm.

[0106] To verify the effectiveness of the RBPF-based real-time water level estimation algorithm, this embodiment selects GNSS observation data from the GTGU and SC02 stations within the same time period as the improved inverse modeling method. The time range for the GTGU station is from day 262 to day 264 of 2023, with a sampling interval of 30 seconds; the time range for the SC02 station is from day 42 to day 44 of 2023, with a sampling interval of 15 seconds. Figure 3 The results of post-hoc water level inversion using the improved inverse modeling method and real-time water level inversion based on RBPF are plotted and compared with a reference water level gauge. The black scatter dots represent the relative water level of the reference water level gauge, the blue curve represents the post-hoc water level inversion results and corresponding residuals of the improved inverse modeling method, and the red curve represents the real-time water level inversion results and corresponding residuals based on RB particle filtering.

[0107] From the water level time series fitting results ( Figure 3As can be seen from (a) and (c) in this embodiment, the real-time inversion algorithm based on RBPF proposed in this embodiment can still accurately capture the dynamic change pattern of water level under the online recursive constraint that only uses the current observation data and cannot obtain future global information.

[0108] For GTGU stations with gentle water level changes, the real-time inverted water level curves, the measured data from the water level gauges, and the high-precision inversion results after the event all showed a high degree of consistency throughout the entire experimental period. Figure 3 In (a) of the data, the water level fluctuation trend, the timing of peaks and troughs, and their amplitudes were all accurately matched, with no significant phase delay, amplitude distortion, or trend deviation. For the SC02 station, where daily water level changes exceed 3m, the real-time algorithm completely reproduced the full-cycle water level variation characteristics. The peak and trough positions, fluctuation rates, and measured data were consistent with the post-analysis results, showing no trend deviation. This verifies the algorithm's universality and tracking stability in highly dynamic water level change scenarios. Furthermore, the residual distribution results ( Figure 3 As can be seen from (b) and (d) in the figures, the residual sequences retrieved in real time from both stations exhibit a random fluctuation distribution centered at 0, without any positive or negative trend drift. Further statistical inversion results and the probability density function of the residual distribution of the water level gauge reference values ​​are shown in... Figure 4 As shown, where Figure 4 In the figure, (a) and (b) represent the residual distributions of the GTGU and SC02 stations, respectively.

[0109] Figure 4 The statistical results in (a) show that the fluctuation range of the real-time inversion residuals of the GTGU site is basically controlled within the range of -5cm to +5cm, and more than 99% of the residuals are within this range. Within 3cm, the residual mean was 0.21cm and the standard deviation was 1.38cm, indicating that the residuals approximately followed a zero-mean normal distribution with no systematic bias. The residual statistics for station SC02 show ( Figure 4 In (b) of the paper, the fluctuation range of the real-time inversion residuals is within -20cm to +20cm, but nearly 90% of the residuals are within this range. Within 10cm, it indicates that during dynamic periods of large water level fluctuations, the real-time inversion residuals do not diverge with the increase of water level amplitude, and the fluctuation range remains stable. The residual mean is 0.35cm and the standard deviation is 5.88cm, which is close to the zero-mean normal distribution, indicating that even in the case of severe water level fluctuations, the method will not have systematic bias.

[0110] Statistical analysis of the elevation inversion sequences at the two stations revealed that the real-time water level inversion converged to near the reference water level within only 10-15 epochs, indicating that the algorithm can converge quickly. After removing the water level inversion results from the first 15 epochs, the statistical results of the core accuracy indicators of the experiment are shown in Table 2, further verifying the accuracy performance of the real-time algorithm from a quantitative perspective.

[0111] Table 2. Statistics on accuracy indicators of post-event and real-time inversion at different sites

[0112]

[0113] Note: The high-precision inversion in the table refers to the improved dynamic water level inversion method based on inverse modeling.

[0114] The high-precision post-hoc inverse modeling method employs a full-time batch processing mode, utilizing global observation data to optimize the inversion results, representing the high-precision results in this experimental scenario. In contrast, the real-time algorithm uses an online recursive framework based on RBPF, which can only use the observation input at the current moment, resulting in accuracy constraints. Under these conditions, the real-time inversion results for the GTGU station show only a 0.40cm increase in RMSE compared to the post-hoc inversion, with an accuracy decay of 46.5% and a correlation coefficient decrease of only 0.54%, maintaining centimeter-level high inversion accuracy. For the SC02 station, which experiences large water level changes, although the RMSE of the real-time inversion results increases by 4.05cm compared to the post-hoc inversion, the correlation coefficient only decreases by 0.42%, indicating that its water level inversion capability has not deteriorated. For the SC02 station, the difference in accuracy between real-time and post-hoc inversion is greater than that for the GTGU station. This is because the amplitude of water level changes at this station is larger, and the noise impact of real-time observation input is more significant, not due to algorithm performance failure. Its correlation coefficient remains above 0.99, proving that the algorithm can consistently and stably track water level changes.

[0115] To address the limitations of traditional GNSS-IR water level inversion methods, which rely on full-time batch processing and cannot achieve real-time output, this invention proposes a real-time GNSS-IR water level inversion method based on RB particle filtering. First, a state-space model adapted for real-time GNSS-IR inversion is constructed, clearly defining the observation equation with dynamic water level changes as the state variable and real-time signal-to-noise ratio data as the observation variable. To solve the problems of particle degradation and high computational cost in standard particle filtering, which make it difficult to meet real-time requirements, RB particle filtering is used to decouple the state space, dividing the state variable into nonlinear and linear sub-states. These sub-states are estimated in parallel and in real-time using particle filtering and Kalman filtering, respectively. This significantly improves the algorithm's computational efficiency while maintaining estimation accuracy, forming an epoch-wise recursive real-time solution framework. Finally, comparative experiments are conducted with typical GNSS stations exhibiting different water level change characteristics, comparing the proposed real-time inversion method with high-precision inverse modeling post-inversion, verifying its effectiveness, stability, and scenario universality in different hydrological scenarios.

[0116] This invention also provides a GNSS-IR real-time water level inversion system based on RB particle filtering, comprising:

[0117] The data acquisition module is used to acquire signal-to-noise ratio observation data and satellite ephemeris data from the GNSS receiver;

[0118] The preprocessing module is used to preprocess the signal-to-noise ratio observation data and satellite ephemeris data to obtain the processed reflection signal signal-to-noise ratio sequence and satellite elevation angle data.

[0119] The state decomposition module is used to construct a state space model based on Rao-Blackwellized particle filtering, and decompose the state vector of GNSS-IR water level inversion into nonlinear substates and linear substates. The nonlinear substates include reflector height and reflector surface roughness parameters, and the linear substates include in-phase component coefficients and quadrature component coefficients of the reflected signals in each frequency band.

[0120] The parallel filtering estimation module is used to perform sequential estimation of the nonlinear substate using particle filtering and analytical estimation of the linear substate using Kalman filtering at each signal-to-noise ratio sampling epoch. The two filtering methods are executed in parallel and the particle weights are updated using the Gaussian likelihood function calculated from the current epoch observation value.

[0121] The water level output module is used to perform a weighted average of the state estimates of all particles based on the weight of each particle, output the optimal estimate of the reflector height for the current epoch, and output the water level value epoch by epoch based on the known GNSS antenna elevation reference.

[0122] The specific implementation methods of the above modules are the same as the specific implementation methods of each step in the method embodiment, and will not be repeated here.

[0123] In summary, this invention proposes a novel method that balances real-time computational efficiency and monitoring accuracy—a real-time water level inversion algorithm based on the RB particle filter (RBPF). Addressing the strong nonlinearity and non-Gaussian observation characteristics of GNSS-IR inversion, and the problems of particle degradation, high computational load, and difficulty in meeting online real-time computation requirements of standard particle filtering, this invention employs RBPF to decouple the state variables into nonlinear and linear sub-states. Parallel estimation is performed using particle filtering and Kalman filtering, significantly reducing the required number of particles while maintaining inversion accuracy, thus significantly improving online computational efficiency and ultimately forming a complete epoch-by-epoch recursive inversion process. Comparative verification at two stations with different water level change characteristics shows that, with 1500 particles, the RMSEs of the inversion results and the water level gauge reference values ​​at the two stations are 1.26 cm and 5.88 cm, respectively, with correlation coefficients exceeding 99%. The residuals follow an approximately zero-mean normal distribution, with no systematic bias or error accumulation. The average computation time per epoch is only about 6 ms, indicating that this real-time algorithm can accurately track dynamic water level changes using only the constraints of current observation data.

[0124] 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.

[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 real-time water level inversion method based on RB particle filtering, characterized in that, include: Acquire signal-to-noise ratio observation data and satellite ephemeris data from the GNSS receiver; The signal-to-noise ratio observation data and satellite ephemeris data are preprocessed to obtain the processed signal-to-noise ratio sequence of the reflected signal and satellite elevation angle data; A state-space model is constructed based on Rao-Blackwellized particle filtering. The state vector of GNSS-IR water level inversion is decomposed into nonlinear substates and linear substates. The nonlinear substates include reflector height and reflector surface roughness parameters, and the linear substates include in-phase component coefficients and quadrature component coefficients of the reflected signals in each frequency band. At each signal-to-noise ratio sampling epoch, the nonlinear substate is sequentially estimated using particle filtering, and the linear substate is analytically estimated using Kalman filtering. The two filtering methods are executed in parallel, and the particle weights are updated using the Gaussian likelihood function calculated from the current epoch observations. The state estimates of all particles are weighted and averaged according to their respective weights to output the optimal estimate of the reflector height for the current epoch. The water level value is then output epoch by epoch based on the known GNSS antenna elevation reference.

2. The GNSS-IR real-time water level inversion method based on RB particle filtering according to claim 1, characterized in that, A state-space model based on Rao-Blackwellized particle filtering is constructed, decomposing the state vector of GNSS-IR water level inversion into nonlinear and linear sub-states, including: Based on the Rao-Blackwellization theorem, the reflector height and reflector surface roughness parameters are taken as nonlinear sub-states, and their state transition adopts a random walk model. The predicted particle for the next epoch is generated by superimposing process noise. The in-phase and quadrature component coefficients of the reflected signals in each frequency band are taken as linear sub-states. Under the condition of a given nonlinear sub-state, the posterior distribution of the linear sub-state follows a Gaussian distribution. The in-phase and quadrature component coefficients of each frequency band are estimated using a Kalman filter.

3. The GNSS-IR real-time water level inversion method based on RB particle filtering according to claim 1, characterized in that, The filtering recursion for each sampling epoch is performed in the following order: Based on the posterior state estimate of the previous epoch, the state prediction of the current epoch is completed through the state transition model. Using the signal-to-noise ratio observation of the preprocessed reflected signal at the current epoch, sequential Kalman filter update and particle weight update are performed sequentially. The updated particle states are weighted and averaged to output the optimal estimate of the reflector height for the current epoch. Determine whether resampling is triggered and prepare the particle set for the filtering recursion of the next epoch; The data input for the filtering recursion is limited to the signal-to-noise ratio observation value of the current epoch and the posterior state estimation result of the previous epoch, excluding the participation of observation data of future epochs, thus forming an epoch-by-epoch causal recursive structure.

4. The GNSS-IR real-time water level inversion method based on RB particle filtering according to claim 1, characterized in that, The particle weights are updated using the Gaussian likelihood function calculated from the current epoch observations, including: For each reflected signal-to-noise ratio observation in the current epoch, the prediction residual is calculated based on the updated linear substate posterior estimate; A Gaussian observation likelihood function is constructed based on the prediction residuals of all observations, and the likelihood value of each nonlinear particle at the current epoch is calculated. Multiply the particle weights of the previous epoch by the likelihood value of the current epoch to obtain the unnormalized weights of the current epoch. Normalize the unnormalized weights of all particles.

5. The GNSS-IR real-time water level inversion method based on RB particle filtering according to claim 1, characterized in that, In the observation equation of the signal-to-noise ratio sequence of the reflected signal, the reflector height and the sine of the satellite elevation angle form a nonlinear functional relationship, the in-phase component coefficients and quadrature component coefficients of each frequency band form a linear combination relationship, and the surface roughness parameter is introduced into the observation equation in the form of exponential decay.

6. The GNSS-IR real-time water level inversion method based on RB particle filtering according to claim 1, characterized in that, The construction of the state-space model also includes: Based on the prior height range of the GNSS receiver antenna phase center distance from the water surface, an initial particle distribution for the reflector height in the nonlinear substate is set, and the initial distribution follows a Gaussian distribution. The initial values ​​of the surface roughness parameters are set to zero, the initial mean values ​​of the in-phase component coefficients and quadrature component coefficients of each frequency band are set to zero, and the initial covariance is set to a preset value. The process noise variances of each nonlinear and linear substate were determined using the first-order difference statistics of historical observation data.

7. The GNSS-IR real-time water level inversion method based on RB particle filtering according to claim 1, characterized in that, The signal-to-noise ratio (SNR) observation data of the GNSS receiver comes from multiple satellite systems and multiple frequency bands. The linear sub-state includes the in-phase component coefficient and quadrature component coefficient corresponding to each frequency band. The number of components in each frequency band is dynamically determined according to the number of effective satellite-frequency band combinations received in the current epoch. The observation noise variance of each frequency band is independently estimated based on the statistical characteristics of the SNR observation residuals of different frequency bands.

8. The GNSS-IR real-time water level inversion method based on RB particle filtering according to claim 1, characterized in that, After updating the particle weights for the current epoch, calculate the number of effective particles: ， Where N is the total number of particles. The normalized weight of the i-th particle; When the number of effective particles is less than half of the total number of particles, a resampling operation is triggered. The resampling adopts a low-variance resampling method, which eliminates low-weight particles and replicates high-weight particles. After resampling, a small noise is superimposed on the new particles, and the weights of all particles are reset to equal values.

9. The GNSS-IR real-time water level inversion method based on RB particle filtering according to claim 1, characterized in that, The number of nonlinear particles in the Rao-Blackwellized particle filter is 1500, and the average calculation time per epoch in the parallel filter is in the millisecond range and is less than the sampling interval of the signal-to-noise ratio data.

10. A GNSS-IR real-time water level inversion system based on RB particle filtering, characterized in that, include: The data acquisition module is used to acquire signal-to-noise ratio observation data and satellite ephemeris data from the GNSS receiver; The preprocessing module is used to preprocess the signal-to-noise ratio observation data and satellite ephemeris data to obtain the processed reflection signal signal-to-noise ratio sequence and satellite elevation angle data. The state decomposition module is used to construct a state space model based on Rao-Blackwellized particle filtering, and decompose the state vector of GNSS-IR water level inversion into nonlinear substates and linear substates. The nonlinear substates include reflector height and reflector surface roughness parameters, and the linear substates include in-phase component coefficients and quadrature component coefficients of the reflected signals in each frequency band. The parallel filtering estimation module is used to perform sequential estimation of the nonlinear substate using particle filtering and analytical estimation of the linear substate using Kalman filtering at each signal-to-noise ratio sampling epoch. The two filtering methods are executed in parallel and the particle weights are updated using the Gaussian likelihood function calculated from the current epoch observation value. The water level output module is used to perform a weighted average of the state estimates of all particles based on the weight of each particle, output the optimal estimate of the reflector height for the current epoch, and output the water level value epoch by epoch based on the known GNSS antenna elevation reference.

Images