Improved Global Navigation Satellite System–Multipath Reflectometry (GNSS-MR) Tide Variation Monitoring Using Variational Mode Decomposition Enhancement

Abstract

Accuracy and resolution are the two primary challenges that impose limitations on the practical implementation of classical tide-level remote sensing. To improve the accuracy and applicability and increase the temporal resolution of the inversion point near the shore area, the influence of coastal reflection signals in the signal-to-noise ratio (SNR) residual sequence should be weakened significantly. This contribution proposes an anti-interference GNSS Multipath Reflectometry (GNSS-MR) algorithm called VMD_SNR, which is enhanced using variational mode decomposition (VMD). Compared with wavelet decomposition and empirical mode decomposition (EMD) methods, VMD_SNR exhibits superior capabilities in reducing the interference caused by noisy signals. The measurements of ground-based GNSS stations are used to verify the performance improvement in the VMD_SNR algorithm. The results show that the proposed algorithm is better than the wavelet decomposition method and EMD method in terms of accuracy and stability in the shore area, where the effective number is higher than 99% of the total number, and the accuracy is better than 13.80 cm. Moreover, the accuracy improvement is more significant in the high-elevation range, which is 30.16% higher than the wavelet decomposition method and 38.34% higher than the EMD method.

1. Introduction

With global warming and the deterioration of the marine environment, coastal regions face significant threats to their security. Therefore, effective long-term tide monitoring has become a matter of great concern. In recent years, Global Navigation Satellite System reflectometry (GNSS-R), as a new remote sensing monitoring method, has been widely used to detect sea level changes. This is a technique that utilizes the characteristics of satellite signals reflecting off the sea surface and returning to the receiver to gather information about the ocean surface. Compared with traditional tide-level measurement methods, GNSS-R technology can not only compensate for the spatial resolution limitations of tide stations but also improve the low precision of satellite measurement in offshore areas.

Martin-Neira [1] proposed the concept of ocean altimetry and discussed the potential applications of GPS-R technology in Earth observation and climate change research. This laid the foundation for the development and application of GNSS-R technology. Subsequently, the National Aeronautics and Space Administration (NASA) proposed the potential application of using GPS-R for remote sensing sea state [2,3].

With the gradual development and application of GNSS-R technology, many researchers have demonstrated that the technology can act as an “alternative tide gauge”. GNSS-R technology can be broadly categorized into two groups: path delay analysis and signal-to-noise ratio (SNR) analysis. The first technique uses a dual antenna, which is used to receive the direct signal from the satellite and the reflected signal. By measuring the phase difference between these two signals, it is possible to calculate the total path delay of the signal from the satellite to the surface reflection point and back to the receiver. This enables the acquisition of information about the reflection point. This approach can be further divided into phase delay analysis [4] and code delay analysis [5]. Phase delay analysis requires a high-precision frequency counter for measuring the phase difference, enabling the acquisition of accurate reflection point information. On the other hand, code delay analysis, due to its sensitivity to code variations, may have slightly lower accuracy compared to phase delay analysis. The latter uses a single antenna to determine the height information of the reflection point by analyzing the characteristics of the SNR sequence [6,7]. While this method is relatively simple and practical, its accuracy is lower than that of the path delay method. A benefit of using the SNR method is greater robustness in handling wind and wave conditions [8]. Furthermore, it has been demonstrated that the method is useful for determining other crucial sea state parameters, such as significant wave height [9].

In the classical GNSS-MR method, Lomb–Scargle periodogram (LSP) analysis is commonly employed to extract the multipath frequency. However, it is important to note that this method displays the power of all frequencies present in an SNR series, including noise [10,11]. Consequently, the power of the reflected signal from the sea surface can be diminished or overshadowed, potentially leading to a shift or error in the results [12]. Therefore, quality control measures are very important in the process of LSP analysis, and increasing the anti-interference of the algorithm has become the focus of research. Strandberg, et al. [13] relied on a B-spline representation of the temporal sea-level variations to describe the reflection and propagation process of GNSS signals and then used a reverse modeling algorithm to fit SNR data and extract sea-level height information from it to reduce noise interference. Wang, et al. [12], starting from signal decomposition, proposed the wavelet decomposition (WD) method to decompose the SNR residual sequence after removing direct signal into multiple bands and then estimated the water-level frequency for each band using the LSP method to obtain the sea surface height. Zhang, et al. [14] introduced an SNR signal decomposition method based on empirical mode decomposition (EMD). This method sidesteps the loss of vital information while enhancing the accuracy of GNSS-MR tidal-level monitoring. Although the above signal processing methods can effectively alleviate the interference of noisy signals, they also have some limitations. Wavelet decomposition, for instance, is limited by a fixed basis function and constant multi-resolution, which limits its adaptability. On the other hand, EMD uses interpolation to separate the “modes”, which greatly reduces the robustness of the algorithm due to the temporary nature of the method and the lack of mathematical theory; in addition, the use of recursion to separate the remaining components in EMD can amplify the transmission of errors [15,16].

Variational mode decomposition (VMD), as a novel signal decomposition method, boasts both strong adaptability and rigorous mathematical derivation. This characteristic aids in mitigating the noise sensitivity problem encountered in the EMD method [17]. Recently, VMD has been successfully tested in the fields of structural health diagnosis [18], wind speed forecasting [19] and denoising processing of seismic exploration data [20]. This paper introduces a VMD_SNR algorithm for GNSS-MR tide monitoring combined with VMD enhancement. The VMD method is combined with the GNSS-MR principle to effectively separate the noise signal from the reflected signal of different reflectors in the original SNR sequence. Through the use of signal extraction criteria, the algorithm extracts the reflected signals of the target, enabling the assessment of tidal changes. The aim is to enhance the stability of GNSS-MR tidal monitoring results in coastal areas and improve the temporal resolution of the inversion point. This method can also be applied in areas, such as GNSS-MR snow depth detection and dynamic sea surface correction, to enhance the accuracy of snow depth detection and tidal monitoring, as well as to increase the temporal resolution of GNSS-MR algorithm inversion points.

2. Improved GNSS-MR Tide Monitoring Algorithm

2.1. Variational Mode Decomposition

Variational mode decomposition (VMD) is a novel decomposition method that can decompose the original signal into bandwidth-limited intrinsic mode functions (IMFs) and compute their respective center frequencies. Unlike the recursive mode extraction of EMD, VMD extracts modes simultaneously, thus appropriately balancing the errors between them. These modes collectively reproduce the input signal, and each mode, when demodulated to the baseband, becomes a smooth AM-FM signal. Its mathematical expression is as follows:

$$u_k\left(t\right)=A_k\left(t\right)\mathrm{c}\mathrm{o}\mathrm{s}\left({\Phi }_k\left(t\right)\right)$$

(1)

where ${\Phi }_k\left(t\right)$ is the phase, a non-decreasing function, ${\Phi }_k\left(t\right)>0$; $A_k\left(t\right)$ is the envelope, a non-negative value.

The corresponding constrained variational expression for the VMD algorithm is:

$$\left\{\begin{array}{c}\begin{array}{c}\mathrm{m}\mathrm{i}\mathrm{n}\\ \left\{u_k\right\},\left\{{\omega }_k\right\}\end{array}\left\{{\sum _K‖{\partial }_t\left[\left(\delta \left(t\right)+\frac{j}{\pi t}\right)\ast u_k\left(t\right)\right]e^{-j{\omega }_{k}t}‖}_2^2\right\}\\ s. t. \sum _K{u}_k=S\end{array}\right.$$

(2)

where $K$ is the number of modes to be decomposed, and its determination methods mainly include Particle Swarm Optimization, Center Frequency Observation, and EMD pre-decomposition. In this paper, the commonly used Center Frequency Observation method is adopted to determine the value of $K$. $\left\{u_k\right\}$ and $\left\{{\omega }_k\right\}$ represent the $k$-th mode component and the central frequency, respectively, while $\delta \left(t\right)$ stands for the Dirac delta function, and $\ast$ denotes the convolution operator.

Transforming the constrained variational problem into an unconstrained variational problem, the augmented Lagrange expression is obtained as follows:

$$\begin{gathered} {L\left(\left\{u_k\left(t\right)\right\},\left\{{\omega }_k\right\},\lambda \left(t\right)\right)=\alpha {\sum }_K‖{\partial }_t\left[\left(\delta \left(t\right)+\frac{j}{\pi t}\right)\times u_k\left(t\right)\right]e^{-j{\omega }_{k}t}‖}_2^2 \\ +{‖S\left(t\right)-{\sum }_K{u}_k\left(t\right)‖}_2^2+〈\lambda \left(t\right),S\left(t\right)-{\sum }_K{u}_k\left(t\right)〉 \end{gathered}$$

(3)

where $\alpha$ represents the penalty factor, which is typically set to a default value of $\alpha =2000$. $\lambda$ denotes the Lagrange multiplier. Each mode component $u_k$ can be obtained using the following formula:

$${\widehat{u}}_k^{n+1}\left(\omega \right)\leftarrow \frac{\widehat{f}\left(\omega \right)-\sum _{i<k}{\widehat{u}}_i^{n+1}\left(\omega \right)-\sum _{i>k}{\widehat{u}}_i^n\left(\omega \right)+\frac{{\widehat{\lambda }}^n\left(\omega \right)}{2}}{1+2\alpha {(\omega -{\omega }_k^n)}^2}$$

(4)

$${\omega }_k^{n+1}\leftarrow \frac{{\int }_0^{\infty }\omega {\left|{\widehat{u}}_k^{n+1}\left(\omega \right)\right|}^{2}d\omega }{{\int }_0^{\infty }{\left|{\widehat{u}}_k^{n+1}\left(\omega \right)\right|}^{2}d\omega }$$

(5)

$${\widehat{\lambda }}^{n+1}\left(\omega \right)\leftarrow {\widehat{\lambda }}^n\left(\omega \right)+\tau (\widehat{f}\left(\omega \right)-{\sum }_K{\widehat{u}}_k^{n+1}\left(\omega \right))$$

(6)

where $\tau$ represents the noise tolerance. Under conditions of higher noise levels, the Lagrange multiplier may compromise the accuracy of the modes. Therefore, by setting $\lambda =0$, the Lagrange multiplier is discarded. ${\widehat{u}}_k^{n+1}\left(\omega \right)$, ${\widehat{u}}_i^n\left(\omega \right)$, $\widehat{f}\left(\omega \right)$ and ${\widehat{\lambda }}^n\left(\omega \right)$ correspond to the Fourier transforms of $u_k^{n+1}\left(\omega \right)$, $u_i^n\left(\omega \right)$, $f\left(\omega \right)$ and ${\lambda }^n\left(\omega \right)$, respectively. Here, $f\left(\omega \right)$ represents the frequency domain representation of the signal to be decomposed.

Compared with previous studies [12,15], the algorithm presented in this research avoids modeling a single mode as a signal with an explicit IMF formula. Instead, it uses the associated narrowband properties, provides mathematical formulas, clarifies the theoretical basis, establishes relationships between the parameters of the IMF model and the estimated signal bandwidth and enhances the algorithm’s robustness. Moreover, each mode is constrained around its central frequency in a band-limited manner. Consequently, this approach effectively addresses the issue of noise sensitivity inherent in the EMD method.

2.2. VMD_SNR

The SNR, which is a measure of the signal strength received by a GPS antenna, can be expressed in the classical GNSS-MR:

$$SNR=\frac{1}{P_n}\left(P_d+P_r+2\sqrt{P_d{P}_r}cos{\phi }_i\right)=\frac{\left(P_d+P_r\right)}{P_n}+\frac{2\sqrt{P_d{P}_r}}{P_n}cos{\phi }_i$$

(7)

where $P_d$ is the power received directly from the satellite, $P_r$ is the power received indirectly from the satellite, $P_n$ is the noise power and ${\phi }_i$ is the phase difference between the direct signal and the $i$th reflected signal. SNR is commonly divided into two components through low-order polynomial fitting [21]: the part $\overline{s}$, which mainly contains direct signals ($\frac{\left(P_d+P_r\right)}{N}$), and the part $\delta S$, which contains interference signals ($\frac{2\sqrt{P_d{P}_r}}{N}cos{\phi }_i$). Since most of the object information exists in the reflected signal, we pay more attention to the segment $\delta S$ that contains the interference signal for the segmented result.

Under the static assumption of a constant reflector surface, when there is a one-ray reflection condition, the phase difference between direct signal and reflected signal can be expressed as [22] (Figure 1):

$$\phi =2\pi ∆S/\lambda =4\pi hsin\theta /\lambda$$

(8)

where $h$ is the vertical distance between the antenna and the reflecting surface (RH), $\theta$ is the elevation angle of the GPS satellite and $\lambda$ is the wavelength of the GPS signal. Therefore,

$${\omega }_{\phi }=2\pi f=\frac{d\phi }{dsin\theta }=4\pi h/\lambda$$

(9)

From Equation (9), the sinusoid frequency $f$, which can be acquired through LSP analysis applied to $\delta S$, can be expressed by $h=\lambda f/2$. All the LSPs in this study are normalized for the ordinate and converted from frequency to RH along the abscissa. However, due to the reception of complex mixed signals by GNSS, interference from other signals can introduce bias into the results (Figure 2). Figure 2 shows the spectral results of the SNR residual sequence LSP at different elevation angles. The horizontal coordinate represents the vertical reflection height (RH), while the vertical coordinate represents the spectrum amplitude. Larson, Löfgren and Haas [6] believed that reflecting surfaces with RH values of 4–7 m tend to be associated with the sea. Therefore, frequencies associated with reflecting surfaces with RH values of <4 or >7 m tend to be noisy frequencies. The reflected signal at RH < 3 m is more likely to come from the ground and objects on it. Figure 2 shows that, as the elevation angles increase, the number of extreme points in the LSP also increases, with a majority of them occurring at RH < 2 m. While extreme points can be observed in the previous ocean values across all four diagrams, only the LSP for the SNR with elevation angles of 5°–13° and 5°–20° exhibit peaks corresponding to the sea surface RH. If higher elevation angles, such as 5°–25° or 5°–30°, were used, these RH peaks would be disregarded due to being considered nonphysical. Consequently, employing higher elevation angles would result in skewed results.

Therefore, this paper proposes the VMD_SNR algorithm. The specific process is as follows:

• Initially, the SNRs were selected based on azimuths and elevation angles corresponding to the sea.
• Separate these SNRs to $\overline{s}$ and $\delta S$ by fitting a polynomial.
• Configure VMD parameters, using the center frequency observation method to determine the value of $K$.
• Using the configured VMD method, calculate each IMF component of the $\delta S$ sequence according to Equation (2).
• Determine the reflected signal of the target based on the dual constraints of the mutual relation number and the prior value of the tide level.

  $$R\left(\delta S, {imf}_i\right)=\frac{\sum _{i=1}^N\left[\delta S\left(t\right)-\overline{\delta S\left(t\right)}\right]\left[{imf}_i\left(t\right)-\overline{{imf}_i}\right]}{\sqrt{{\sum }_{i=1}^N{\left[\delta S\left(t\right)-\overline{\delta S\left(t\right)}\right]}^2}\sqrt{{\sum }_{i=1}^K{\left[{imf}_i\left(t\right)-\overline{{imf}_i}\right]}^2}}$$

  (10)

  where $\overline{\delta S\left(t\right)}$ and $\overline{{imf}_i}$, respectively, represent the average values of SNR residuals and the average values of each IMF component.
• Calculate the reflected signal of the target using the LSP analysis method, and then obtain the corresponding frequency value.

2.3. Tidal Monitoring Using VMD_SNR Algorithm

To illustrate the tide monitoring process based on the VMD_SNR algorithm, we selected a sequence of SNR data from the SC02 station of PRN11. The data were obtained on day of year (DOY) 232 in 2011, with an azimuth angle range of 50°–240° and different elevation angle ranges, as shown in Figure 3. In the figure, the horizontal axis represents the sine of the elevation angle, and the vertical axis represents the signal amplitude. The last line represents the original $\delta S$ signal, the penultimate line represents the trend term of the signal and IMF 1–5 represents the five eigenmode components obtained via decomposition. It is evident that VMD can effectively decompose the original signal into components with different frequencies, allowing for extraction and analysis of the target signal through the L-S spectrum while minimizing interference from noisy signals. This paper employs a dual-constraint approach involving cross-correlation coefficients and prior tidal values. By calculating the correlation between the intrinsic mode function (IMF) components, from which the LSP results within the tidal prior values, and $\delta S$, the IMF component with the highest correlation is selected as the target signal for SNR analysis. Since the ocean RH value can be obtained only by knowing the ocean signal frequency, only the IMF component containing ocean information needs to be extracted; it is not necessary to extract all components containing ocean signals. From

Table 1. The IMF components’ probability density for the SC02 station.

IMF Component (%)	1	2	3	4	5
Probability density	0.17	0.16	0.22	0.27	0.74

, it can be seen that imf5 has the highest correlation, and its frequency range is more consistent with the tidal variation of the SC02 station. Therefore, imf5 is selected as the target signal for further analysis. It is worth noting that the IMF number $K$ resulting from VMD decomposition is not fixed and will change with the complexity of the signal. Consequently, the choice of the target component is also not fixed; it is not limited only to imf5. This correlation-based selection method depends on the signal itself and is not artificially constrained or defined.

Figure 4a shows the target signal IMF 5 and $\delta S$ sequence, and Figure 4b shows the LSP results of the two sequences. Obviously, the sequence extracted using the VMD_SNR algorithm is smoother, and its oscillation is more consistent with the characteristics of reflected signals, indicating that the signal mixing degree is reduced and the single signal feature is enhanced. This point is confirmed in Figure 4b. Compared with the LSP result of the original $\delta S$ sequence, the LSP result of the sequence extracted via the VMD_SNR algorithm greatly reduces the interference of the ground-reflected signal and enhances the energy of the reflected signal from the sea surface, thus avoiding the influence of high-frequency noise and other reflected signals.

Compared with the classical GNSS-MR tidal monitoring method, the VMD_SNR algorithm effectively separates the direct signal from other reflected signals, and the algorithm is not constrained for low elevations (less than 15°), thereby enhancing the applicability of the GNSS-MR tide monitoring algorithm in high-elevation areas. This improvement contributes to enhancing the time resolution of sampling points.

To validate the performance of the VMD method in signal extraction, a set of signal simulation experiments was designed to compare it with the two commonly used signal decomposition methods, wavelet decomposition and EMD. Furthermore, this method was combined with the GNSS-MR algorithm, and two sets of experimental results were designed to verify the enhanced capability of VMD in tide monitoring.

3. GNSS Aliasing Reflection Signal Extraction

In order to evaluate the performance of signal extraction using the VMD method, three sets of sinusoidal signals with different frequencies ($f_1=10,{ f}_2=43, f_3=100$) are simulated in the GNSS-MR process, in which the signal with the set frequency $f_2=43$ is the target signal. The sinusoidal signal expression is:

$$s_i\left(t\right)=A_{mi}sin\left(2\pi \times f_i\times t\right)(i=1,2,3)$$

(11)

where $A_{mi}$ represents the amplitude of the $i$th signal. The experiment set the sampling points to 1000.

Mixture matrix $A$ is a randomly generated 3-by-3 matrix:

$$A=\left(\begin{array}{ccc}0.4227& -0.3017& -0.1624\\ -0.5040& 0.2433& -0.1649\\ -0.1173& 0.7599& -0.2275\end{array}\right)$$

(12)

Mixed signals are represented as:

$$S={\sum }_{i=1}^{3}A{s}_i\left(t\right)+r\left(t\right)$$

(13)

where $r\left(t\right)$ is the white Gaussian noise signal.

To verify the universality of the method, four different levels of white noise are added to the simulation experiment, corresponding to SNRs of 3.7 dB, 8.6 dB, 10 dB and 13.6 dB, respectively. The wavelet decomposition method, EMD method and VMD method were employed to extract the original signals with three frequencies. The parameter configuration for wavelet decomposition was kept consistent with that in Wang, Zhang and Zhang [12], involving Daubechies4 wavelet and Mallat algorithm for sequence decomposition with eight layers. Figure 5 illustrates the partial simulation waveforms obtained using the three methods at an SNR of 10 dB. It is evident from Figure 5 that the waveforms obtained via wavelet decomposition and EMD decomposition exhibit multiple signals, with severe mode aliasing occurring in the separated $s_2$ and $s_3$ signals obtained through the EMD method. Conversely, the signal waveforms obtained via the VMD method appear regular, with relatively singular signal frequencies. It is worth noting that the amplitude and phase of the signals obtained through the separation of the three methods are different from the source signal, which is affected by the characteristics of the mixed signal, signal decomposition method and parameter settings. To gain a more intuitive understanding of the decomposition capabilities of the three methods, the similarity coefficient approach is employed to assess the accuracy of the separation results:

$$\rho =\frac{\sum _{i=1}^N\left[A_i-\overline{A}\right]\left[B_i-\overline{B}\right]}{\sqrt{{\sum }_{i=1}^n{\left[A_i-\overline{A}\right]}^2}\sqrt{{\sum }_{i=1}^n{\left[B_i-\overline{B}\right]}^2}}$$

(14)

where $A$ represents the source signal, and $B$ represents the decomposed signal. Because signal decomposition generates a series of components, it is common to calculate the similarity coefficient between all components and the source signal, excluding noise and trend components. This forms a similarity coefficient matrix, and the component with the maximum similarity coefficient is typically selected as the corresponding component separated.

Table 2. Maximum similarity coefficients of separated signals and source signals using three methods.

SNR/dB	Separation Method
	Wavelet			EMD			VMD
	${\mathit{s}}_1$	${\mathit{s}}_2$	${\mathit{s}}_3$	${\mathit{s}}_1$	${\mathit{s}}_2$	${\mathit{s}}_3$	${\mathit{s}}_1$	${\mathit{s}}_2$	${\mathit{s}}_3$
3.7	−0.8791	0.9306	−0.8486	−0.8761	0.4327	−0.8371	−0.9949	0.9596	−0.9680
8.6	−0.8810	0.9308	−0.8432	−0.6233	0.7127	−0.6920	−0.9955	0.9744	−0.9675
10	−0.8805	0.9312	−0.8444	−0.6885	0.4936	−0.8146	−0.9957	0.9720	−0.9720
13.6	−0.8783	0.9306	−0.8421	−0.6497	0.6578	−0.7611	−0.9941	0.9542	−0.9723

calculates the maximum similarity coefficient between the signals separated by the three methods and the source signal across four different signal-to-noise ratios. As can be seen from the table, the maximum similarity coefficient between the signal components separated by the VMD method and the source signal is more than 95%, the wavelet decomposition method is between 84% and 94% and the EMD method is less than 88%. It is proved that the VMD method is superior to the other two signal decomposition methods in restoring the original signal. In order to highlight the signal extraction capability of the VMD method, the $s_2$ signal was analyzed in spectrum (Figure 6). It is obvious from the figure that the spectrogram of the $s_2$ signal obtained using the VMD method is more single, while the wavelet decomposition method has individual extremes in addition to a maximum value. Due to modal aliasing, the EMD method has more irrelevant frequencies in its spectrogram.

Based on the conducted simulation experiments, it is evident that the VMD method has great advantages in decomposing the signal mixed by a variety of reflected signals, and it can well restore the source signal in the mixed signal. Although the wavelet analysis method can also restore the source signal well, it is affected by the basic function and the number of decomposition layers, and the decomposition result is uncertain. The EMD method has large deviation in the results due to its modal aliasing problem.

4. Tidal Monitoring Experiment

The SC02 station belongs to the EarthScope Plate Boundary Observatory (PBO) network, and its data sampling period is 15 s. It is located on the coast of Friday Harbor in Washington, USA. This station receives GNSS reflection signals from the sea. For comparison purposes, the nearest Friday Harbor tide station is 359 m away and can provide measured water depth data. The station is surrounded by open sea without significant obstructions, allowing it to receive a broad range of reflected signals from the sea. In order to ensure that the experimental results are only affected by the elevation angle, this paper chose an azimuth angle range of 50°–240°. In this work, the period of interest was chosen as DOY 232–239 of 2011. During most of the period, the wind speed of the SC02 was less than 5 m/s.

4.1. Precision Analysis of Tide Monitoring

To fully demonstrate the performance of GNSS-MR tide monitoring based on the VMD_SNR algorithm, this paper compares it with the method of removing the direct signal using only polynomial fitting and the method of extracting sea-surface-reflected signal based on wavelet decomposition and EMD. Figure 7 illustrates the inversion results obtained from the four methods, as well as the measurements from the tide station. In order to show the impact of noise signals on the inversion results, all the results, including outliers, are presented in the figure without any error processing. From the figure, it is evident that when using only the polynomial fitting method, significant outliers are present in the tidal results, reaching approximately 5 m relative to the sea level. Moreover, the number of outliers increases noticeably as the elevation angle expands. The reason is that with the elevation angle rising, the first Fresnel reflection area is closer to the measuring station. When more coastal signals are included, the cleaner ocean signals cannot be obtained only by using polynomial fitting to remove direct signals. After signal decomposition, the outliers of tidal results obtained based on VMD_SNR, wavelet analysis and EMD are significantly reduced, indicating that signal decomposition can effectively improve the influence of coastal reflection signals on inversion results. By comparing the three signal decomposition methods, it is found that the tidal results based on the VMD_SNR algorithm have high similarity and good coincidence with the tidal station results. Despite the presence of individual deviation values in the results based on wavelet analysis, the overall trend is consistent with the results of tide stations. However, the results based on the EMD method differ greatly from those of tide stations. The reason for the analysis is that the EMD method is an adaptive decomposition method, which cannot manually decompose the original signal according to the complexity of the signal, so there are problems such as mode aliasing, which cannot correctly separate the ocean signal.

To provide a clearer comparison among the four methods, an error diagram was generated depicting the difference between the inversion results after removing the two-fold standard deviation and the measured results (Figure 8). The stability and accuracy of the inversion results were analyzed, and four evaluation metrics, including efficiency, average deviation, RMSE and correlation coefficient, were computed and are presented in

Table 3. Comparison of the accuracy of the inversion result of four methods with the elevation changed, from top to bottom, followed by 5°–13°, 5°–15°, 5°–20°.

Elevation	Method	Availability (%)	Mean Bias (cm)	RMSE (cm)	Correlation Coefficient (%)
5°–13°	Quadratic fitting	83.88	11.44	13.71	98.22
	Wavelet	99.77	10.57	13.00	98.36
	EMD	98.60	11.14	13.12	98.41
	VMD	99.77	8.42	10.14	99.05
5°–15°	Quadratic fitting	67.51	12.02	14.80	97.95
	Wavelet	99.77	11.59	13.55	98.24
	EMD	97.24	14.26	17.22	97.43
	VMD	1	6.87	8.56	99.30
5°–20°	Quadratic fitting	43.36	11.87	14.79	97.80
	Wavelet	99.77	10.83	13.07	98.24
	EMD	95.33	13.71	16.88	97.52
	VMD	1	7.51	9.23	99.17

. The availability of the inversion results can be expressed by the formula: Availability = Number of sea-level values greater than 3/Total number of results. Figure 8 illustrates that the inversion points from the VMD_SNR method are evenly distributed near the diagonal line, indicating relatively accurate inversion results obtained using this method. Conversely, the deviation results from polynomial fitting exhibit a notable deviation from the diagonal line for most of the points, indicating substantial deviations in the results obtained using this method. Although the wavelet analysis and EMD methods do not display significant deviations from the diagonal line, the distance from the diagonal line is significantly greater. From Figure 8, it can be observed that the inversion points of the VMD_SNR method are evenly spread around the diagonal line, indicating that the results obtained using this method are relatively accurate. In contrast, the other three methods exhibit a significantly greater distance from the diagonal line.

Combined with the value in Table 3, the results obtained after three signal processing methods consistently maintain a sufficient number of effective results, with a utilization rate of >95% in the three sets of elevation angle ranges. This is significantly higher compared to the polynomial fitting method, especially within an elevation angle range of 5°–20°, where the effectiveness rate is increased by 1.20-times. This result demonstrates the usefulness of employing signal extraction methods for removing interference signals, leading to a significant enhancement in the effectiveness of result sampling points. Comparing the mean bias and RMSE results of the four methods, it can be seen that the accuracy of VMD_SNR is the highest, with the mean bias being 8.42 cm at the maximum and 6.87 cm at the minimum. The maximum size of RMSE was 10.14 cm and the minimum was 8.56 cm. Wavelet analysis was followed by mean biases of 11.59 cm and 10.57 cm, respectively. The maximum RMSE was 13.55 cm and the minimum was 13.00 cm. The maximum mean bias of the EMD method was 14.26 cm and the minimum was 11.14 cm. The maximum RMSE was 17.22 cm and the minimum RMSE was 13.12 cm. The maximum mean deviation of the quadratic fitting method was 12.02 cm and the minimum was 11.44 cm. The maximum size of RMSE was 14.80 cm and the minimum was 13.71 cm. The mean bias of VMD_SNR results was 42.85% higher than that of the quadratic fitting method. Improved by 40.72% compared with wavelet analysis, it was 51.82% higher than the EMD method. The accuracy of RMSE was 42.16% higher than that of the quadratic fitting method. It was 36.83% higher than wavelet analysis. It was 50.29% higher than the EMD method. It is proved that GNSS-MR based on VMD_SNR is feasible in tidal monitoring and has higher performance than the other two signal extraction methods. It is worth noting that the comparative results reveal that the accuracy of the polynomial fitting method is higher than that of the wavelet decomposition and EMD methods. This is because in the processing of the results, the abnormal values with fixed bias in this method are discarded, resulting in higher precision. However, the availability of result sampling points is significantly lower than that of the other two signal processing methods.

4.2. Analysis of Interference Signals

The Fresnel reflection zone refers to the area where electromagnetic waves undergo reflection on a reflective surface [23]. It can serve as a standard for measuring the effective range of ground detection. When electromagnetic waves pass through the first Fresnel zone, the signal at the receiving point is strongest. Therefore, to ensure that the receiving antenna has the strongest capability to receive reflected signals, the first Fresnel reflection zone is often discussed in GNSS-MR. The elevation angle is an important factor that affects the size of the Fresnel reflection zone, and it exhibits an inverse proportional relationship. A smaller elevation angle corresponds to a larger Fresnel reflection zone, while a larger elevation angle results in a smaller Fresnel reflection zone that is closer to the measuring station. Since most stations are situated along the coast, the receiver will inevitably receive reflected signals from the coastline. Additionally, as the elevation angle decreases, the receiver may capture a greater number of signals from the sea surface. Conversely, with a higher elevation angle, the receiver may pick up more signals from the coastal area. To obtain the “pure” seawater reflection signals, the SNR sequence used for inversion is limited to certain elevation and azimuth angle ranges [23]. However, even when the elevation angle is reduced to 10°, some shallow coastal areas remain included within the azimuth angle of 180°–240° (Figure 9). Too small an elevation angle range will result in too short an SNR sequence, which may fail to capture the signal frequency accurately and subsequently affect the accuracy of the results. Moreover, seawater changes dynamically over time, and real-time monitoring of its changes is the goal of tide monitoring. SNR sequences with high elevation angles are also critical to improve the temporal resolution of the results. In order to solve the influence of coastal reflection signal on the accuracy of inversion results, signal decomposition is an effective method. Based on the distribution of the first Fresnel reflection region around the station and the minimum length requirement of the SNR sequence, three sets of elevation angle ranges (5°–10°, 7°–12°, 10°–15°) are designed for comparison, in order to obtain comparable accuracy, and compared with wavelet analysis and EMD signal decomposition methods to verify the performance of VMD_SNR at high elevation angles.

Figure 10 illustrates the GNSS-MR tide monitoring results and the measured results of tide stations using the VMD_SNR method after removing the two-fold standard deviation values of three different elevation angle boundary conditions. It can be clearly observed from the figure that the inversion results obtained using the three groups of experiments have a good correspondence with the measured results, and the correlation coefficient is higher than 98%. However, the deviation values of 7°–12° and 10°–15° are more than those of 5°–10°, and the number of deviation values increases with the increase in elevation angles, indicating that elevation angle does indeed impact the accuracy of the inversion results. It is worth noting that as the elevation angle increases, the number of results gradually decreases. This is because the first Fresnel zone corresponding to some results contains too little ocean area to form an effective frequency, so this part of the results is eliminated. In order to more intuitively compare the accuracy of the three groups of experimental results, the mean deviation, root mean square error and the number of mutual relations were counted to quantitatively analyze the denoising effect (

Table 4. Comparison of the accuracy of the inversion result of three methods with the elevation changed, from top to bottom, followed by 5°–10°, 7°–12°, 10°–15°.

Elevation	Method	Mean Bias (cm)	RMSE (cm)	Correlation Coefficient (%)
5°–10°	Wavelet	11.01	12.75	98.63
	EMD	13.19	15.52	98.14
	VMD	8.58	10.06	99.05
7°–12°	Wavelet	13.27	15.81	97.75
	EMD	13.90	16.02	97.75
	VMD	11.13	12.84	98.46
10°–15°	Wavelet	17.01	19.76	97.15
	EMD	19.00	22.38	95.57
	VMD	11.63	13.80	98.19

). As can be seen from the table, the average deviation of the three groups of experimental results based on the VMD_SNR method is a maximum value of 11.63 cm, a minimum value of 8.58 cm and a difference of 3.05 cm. The maximum value of RMSE was 13.80 cm, the minimum value was 10.06 cm and the difference was 3.74 cm. It shows that the accuracy of the three sets of elevation angle inversion results is similar, which proves that the VMD_SNR method can well extract ocean signals from SNR sequences containing more coastal reflected signals, increasing the temporal resolution of the inversion point near the shore area.

Figure 11 shows the GNSS-MR tide monitoring results based on wavelet analysis, EMD and VMD_SNR measured results of tide stations and statistical values (Table 4) at three different elevation angle ranges. To more intuitively compare the advantages and disadvantages of the three methods, all results are drawn in the graph without eliminating the median error. It can be clearly seen from the figure that in the three elevation angle ranges, the tidal results obtained based on the VMD_SNR method are more consistent with the results of the tide survey station, the RMSE is better than 13.80 cm and the correlation is higher than 98%. The tidal monitoring results based on wavelet analysis have the highest accuracy at an elevation angle range of 5°–10°, with an RMSE of 19.76 cm, but the deviation increases gradually with the rising elevation angle. The deviation of results obtained based on the EMD method is larger than that of the first two methods, and this phenomenon is most obvious at an elevation angle range of 10°–15°, where the RMSE is 22.38 cm. Compared with the three methods, it can be seen that in the high-elevation angle range, the GNSS-MR tidal monitoring performance based on VMD_SNR is the best, which is 30.16% higher than the wavelet analysis method and 38.34% higher than the EMD method.

5. Conclusions

To improve the accuracy and resolution in the practical application of classical GNSS-MR, this paper proposes an anti-interference GNSS-MR algorithm (VMD_SNR) based on VMD enhancement to monitor tide-level changes. The algorithm decomposes the SNR residual sequence using a completely non-recursive VMD method. Compared with the commonly used wavelet decomposition and EMD signal decomposition methods, VMD has a strict mathematical derivation formula and its results are optimized, and it also has better numerical stability. As a result, it can better deal with noise and non-stationarity in signals. Through the design of simulation experiments, the VMD method is compared with wavelet decomposition and the EMD method in the decomposition of complex reflected signals, and the results show that for mixed signals containing multiple reflected signals, VMD shows stronger signal decomposition ability, and the maximum similarity between the separated signal components and the source signal is the highest.

The VMD_SNR algorithm effectively mitigates interference from noise signals and other reflected signals on oceanic reflection signals. It reduces the energy of noise frequencies in the LSP, prevents the occurrence of non-physical peaks and, thereby, enhances the accuracy of inversion results. Additionally, this algorithm extends its applicability to coastal areas, improves the accuracy of results at high altitudes and increases the temporal resolution of inversion points. Using the SC02 station as an example, we evaluate the accuracy and inversion point temporal resolution performance of the VMD_SNR algorithm through two sets of experiments. The experiments are designed with different initial elevation angles but the same elevation angle range, and with the same initial elevation angle but different elevation angle ranges. From the results, it can be observed that compared to the traditional polynomial fitting method, wavelet decomposition and EMD signal processing methods, the proposed approach yields a higher availability of tidal monitoring results, exceeding 99% in validity, over 99% in correlation coefficient and providing accuracy greater than 10.14 cm. It outperforms the polynomial fitting, wavelet decomposition and EMD methods by 42.16%, 36.83% and 50.29%, respectively. Furthermore, the proposed method exhibits excellent applicability, even at high elevation angles, with accuracy comparable to that at low-elevation angles, surpassing 13.80 cm, and improving by 30.16% and 38.34% compared to wavelet decomposition and EMD methods, respectively. It is proved that the proposed method can remove the interference of coastal signals and other noisy signals better than the wavelet decomposition and EMD methods, with stronger applicability and higher reliability.

Monitoring real-time changes in sea level is essential for ocean remote sensing, and in future work, attempts will be made to extract the instantaneous frequencies of SNR sequences using the VMD method, which is of great significance for real-time monitoring of changes in tide gauge stations.