Fusion of PSO-SVM and ICEEMDAN for high stability GNSS-MR sea level height estimation

ABSTRACT

The Global Navigation Satellite System (GNSS) Multipath Reflection (MR) technique utilises the multipath effects of GNSS signals on the sea surface to retrieve tidal variations, playing a crucial role in tidal monitoring. However, traditional GNSS-MR techniques have certain limitations in terms of accuracy and stability due to restrictions in satellite elevation angles and antenna heights. This study proposes a new GNSS-MR sea surface height retrieval method that combines Particle Swarm Optimization (PSO) optimised Support Vector Machine (SVM) with improved complete ensemble empirical mode decomposition with adaptive noise (ICEEMDAN). The method utilises the GNSS multipath interfering signal frequencies (signal-to-noise oscillation term) extracted by ICEEMDAN as input features to the PSO-SVM model to retrieve the sea level height. Using one year of signal-to-noise ratio data from GNSS stations SC02 and TPW2, the stability and accuracy of the proposed method are evaluated under conditions of high satellite elevation angles and without precise GNSS antenna height information. Experimental results demonstrate that the PSO-SVM-ICEEMDAN method is superior to other existing methods in aspects of retrieval stabilities, root mean square errors, consistent tidal patterns, and so on.

KEYWORDS:

• Sea level height estimation
• GNSS-MR
• ICEEMDAN
• PSO
• SVM

1. Introduction

As a crucial component of the Earth’s ecosystem, the ocean profoundly impacts human production and livelihoods. Therefore, monitoring and studying sea surface height can effectively reduce the harm caused by natural disasters, which is particularly significant for coastal regions vulnerable to marine disasters (Bu et al., 2023; Feng et al., 2013; George et al., 2021). Traditional tide gauges cannot conduct large-scale tidal monitoring and are susceptible to external environmental influences. The emergence of GNSS-R technology and multipath GNSS-R technology has provided a novel avenue for tidal monitoring (Camps et al., 2016; Clarizia et al., 2014; Park et al., 2012, 2014; Southwell et al., 2020). GNSS-MR sea surface height retrieval technology analyzes GNSS interference signals reflected from the sea surface to retrieve changes in sea surface height. Compared to traditional tide gauges, it has the advantages of lower cost, ease of application, and scalability. GNSS-R was first proposed in the 1990s, gradually leading to the development of GNSS-MR remote sensing technology (Hall & Cordey, 1988). As GNSS technology research and applications have advanced, multipath errors that originally affected GNSS positioning accuracy have become a novel source of remote sensing signals for retrieving surface reflection characteristics. Consequently, GNSS-MR remote sensing technology has found widespread applications in various fields, including snow depth retrieval (Hu et al., 2023; Liu et al., 2023; Nievinski & Larson, 2014b, 2014c; X. Wang et al., 2020), soil moisture retrieval (Ding et al., 2023; Ha et al., 2022; Liang et al., 2022; Ran et al., 2022; Roussel et al., 2016), sea surface height retrieval (Hu et al., 2021; N. Wang et al., 2022; Zheng et al., 2022), storm surge (Vu et al., 2019) retrieval and sea ice retrieval (Padullés et al., 2023; Schiavulli et al., 2017). Nikolaidou et al. established and validated comprehensive closed formulas for atmospheric delay and altimetry correction in ground-based GNSS-R (Nikolaidou et al., 2023). Motte et al. conducted surface and vegetation monitoring using the GLORI instrument and achieved promising results (Motte et al., 2016). Gonga et al. utilised GNSS-R technology to detect marine plastic debris, revealing the potential to detect substantial accumulations of certain oceanic garbage under specific conditions (Gonga et al., 2023). Martin-Neira et al. introduced an innovative concept for the TriHex mission, aiming to enable high-resolution remote sensing observations (Martin-Neira et al., 2023). Gholamrezaee et al. inverted tidal variations data for four stations over a three-month period using GPS, GLONASS, Galileo, and BDS, with results indicating a better alignment between GPS L1 data and tide gauge data (Gholamrezaee et al., 2023).

Research on sea surface height retrieval based on reflected signals falls into two categories: the phase delay analysis method between dual antennas (Bai et al., 2012; Kucwaj et al., 2016; Lestarquit et al., 2016; Wu et al., 2020) and the SNR analysis method based on single antennas (Fagundes et al., 2021; Liu et al., 2023; Zheng et al., 2021). The SNR analysis method has lower requirements for GNSS receivers. Moreover, the multipath interference signals formed by the direct and reflected signals from GNSS exhibit good robustness to wind and waves. Therefore, the SNR analysis method has advantages in GNSS-MR sea surface height retrieval. The retrieval of sea surface height using the GNSS-MR technique can achieve good inversion results (Larson, Löfgren, et al., 2013). Using signal-to-noise ratio analysis and phase delay analysis to retrieve sea surface height from GPS and GLONASS observations, respectively, can effectively solve the problem of insufficient temporal resolution of sea surface height inversion for a single system (Löfgren & Haas, 2014) and the accuracy of the GLONASS-based sea surface height retrieval method is comparable to or even better than that of the traditional GPS-based sea surface tide retrieval method (Hobiger et al., 2014). In addition to GPS, GLONASS, the BeiDou satellite navigation system can also retrieve the sea surface height (Jin et al., 2017). Cahyaningtyas et al. performed tide level inversion using polynomial fitting and wavelet decomposition at two stations, CBRS and CMOR, in Indonesia, and the results showed that polynomial fitting and wavelet have similar performance (Cahyaningtyas et al., 2023). Empirical mode decomposition (EMD) is considered to be an effective method for extracting more advantageous multipath interference signals and for further improving the accuracy of GNSS-MR-based sea level height retrieval (Zhang et al., 2019). EMD has recognised drawbacks such as mode mixing (Colominas et al., 2014; Ghimire et al., 2022; S. Huang et al., 2019), the use of multiple EMD-kind algorithms (EEMD, CEEMD, CEEMDAN, ICEEMDAN), leading to proposals to extract multipath interference signals. The results show the accuracy of sea surface height retrieval from GNSS-MR is higher than that of the EMD algorithm when using the ICEEMDAN algorithm, which effectively overcomes the end-point effect and mode mixing problems associated with the EMD algorithm (Jian et al., 2023). Since the retrieval accuracy of GNSS-MR mainly depends on the extraction accuracy of the multipath interference signal and the accuracy of the GNSS antenna height, while the sea level height retrieval method used by Jian et al. is able to extract the multipath interference signal more efficiently, but it relies on the high-precision GNSS antenna height information. Machine learning can effectively mine the hidden feature information and autonomously learn external information such as antenna height during the training process, thus effectively solving the dependence on a priori knowledge such as antenna height. Based on the foregoing and considering the advantages of ICEEMDAN in extracting multipath interference signals, this paper proposes a sea level height retrieval method that integrates Machine Learning PSO-SVM and ICEEMDAN to realise sea level height retrieval with high stability and high accuracy.

The organisation of the remaining parts of this paper is as follows: Section 2 introduces the basic principles of GNSS-MR sea surface height retrieval. Section 3 presents the high-stability sea surface height retrieval method using the proposed PSO-SVM-ICEEMDAN. Section 4 experimentally validated and discussed the proposed method. Finally, Section 5 summarises the main conclusions of this paper.

2. Principles of SNR-based GNSS-MR sea level height retrieval

The principle of the GNSS-MR sea level height retrieval method is shown in Figure 1. Here $\mathit{h}$ is the vertical reflection distance (distance from the phase center of the GNSS antenna to the sea level), $h_i$ is the sea level tide height and $\mathit{H}$ is the GNSS antenna height. Since the sea level reflects the GNSS signal, the SNR observation received by the GNSS receiver is a composite signal generated by the interference of the reflected and direct signals. It can be expressed as (Nievinski & Larson, 2014a):

(1) $\mathit{SNR}=A_d^2+A_r^2+A_d{A}_{r}cos\mathit{\alpha }$(1)

Figure 1. Schematic diagram of the principle of GNSS-MR sea level height retrieval.

where $A_d$ is the direct signal power and $A_r$ is the reflected signal power. $\mathit{\alpha }$ is the phase delay between the direct and reflected signals, $A_d^2+A_r^2$ is the trend term of the composite signal, and $A_d{A}_{r}cos\mathit{\alpha }$ is the oscillatory term of the composite signal caused by the interference of the direct and reflected signals (interference signal, SNR oscillatory term).

Since only the multipath interference signal carries the sea level height information, and the multipath interference signal will lead to oscillatory changes in the SNR observation signal. Therefore, removing the trend component from the original SNR signal and then extracting the SNR oscillatory term as the multipath interference signal is necessary to realise the sea level tide retrieval. Traditional methods generally use polynomial fitting to extract the SNR oscillatory term. The amplitude of the obtained SNR oscillatory term can be expressed as (Larson et al., 2013) follows:

(2) $A_m=\mathit{A}cos\left(\frac{4\mathit{\pi h}}{\lambda }sin\mathit{E}+\mathit{\phi }\right)$(2)

where $\mathit{\lambda }$ is the wavelength of the carrier signal, $\mathit{A}$ and $\mathit{\phi }$ are the magnitude and phase of the SNR oscillatory term, respectively, and $\mathit{E}$ is the satellite elevation angle. We set the following to more accurately characterise the SNR oscillatory term:

(3) $\mathit{t}=sin\mathit{E},\begin{array}{cc}\mathit{\hspace{0.17em}}& \mathit{f}=\frac{2\mathit{h}}{\lambda }\end{array}$(3)

Equation (2)(2) $A_m=\mathit{A}cos\left(\frac{4\mathit{\pi h}}{\lambda }sin\mathit{E}+\mathit{\phi }\right)$(2) can be simplified to the standard cosine function expression as

(4) $A_m=\mathit{A}cos(2\mathit{\pi ft}+\mathit{\phi })$(4)

The peak frequency $\mathit{f}$ is obtained by performing Lomb-Scargle (LSP) spectral analysis spectral analysis of the SNR oscillation term. Then the vertical reflection distance $\mathit{h}$ is obtained from the following equation:

(5) $\mathit{h}=\frac{\mathit{\lambda f}}{2}$(5)

In conjunction with Figure 1, the sea level height can be obtained from the following equation:

(6) $h_i=\mathit{H}-\mathit{h}$(6)

3. High-stability sea level height retrieval method using PSO-SVM-ICEEMDAN

Most of the classical GNSS-MR techniques described in Section 2 use polynomial fitting to extract SNR oscillation terms for low satellite altitude angles (5°-15°). However, LSP spectrum analysis usually requires a sufficient amount of SNR data. Therefore, increasing the satellite elevation angle to some extent will give better retrieval results. Traditional polynomial fitting cannot adapt to the increasing satellite elevation angles. As the satellite elevation angle increases, the extracted SNR oscillatory term includes trend components, leading to distorted inversion results that reduce retrieval stability. Thus, this paper proposes a GNSS-MR sea surface height retrieval method that integrates PSO-SVM and ICEEMDAN. This method does not employ traditional polynomial fitting to extract SNR oscillatory terms; instead, it utilises the ICEEMDAN algorithm for SNR oscillatory term extraction.

3.1. Basic principles of ICEEMDAN

The improved complete ensemble empirical mode decomposition with adaptive noise (ICEEMDAN) was proposed by (Colominas et al., 2014). ICEEMDAN is an enhancement of the empirical mode decomposition (EMD) that effectively addresses the mode mixing and endpoint effects associated with EMD (Colominas et al., 2014). Similar to other EMD-kind algorithms, ICEEMDAN is an adaptive data processing method suitable for handling nonlinear and non-stationary time series data. The main principles of ICEEMDAN can be summarised as follows:

(7) $\mathit{x}(\mathit{t})=\sum _{\mathit{i}=1}^N{c}_i+\mathit{r},$(7)

where $\mathit{x}(\mathit{t})$, $\mathit{N}$, $c_i$, and $\mathit{r}$ are the primitive sequence, number of IMF, IMF, and residual term, respectively. Because ICEEMDAN approaches decomposition from the perspective of adaptive noise, its decomposition method is best suited to resolving the issues of pseudo modes and residual noise in other EMD-kind algorithms.

3.2. Basic principles of PSO-SVM

Support vector machine (SVM) was originally conceptualised by Corinna et al. at the end of the 20th century (Cortes & Vapnik, 1995). It has shown significant advantages in regression prediction, time series forecasting, and comprehensive evaluation fields. SVM regression prediction methods typically divide sample data into training and testing sets and map the training data into a high-dimensional space for linear regression. Therefore, the main principle of SVM can be summarised as solving the optimal decision function:

(8) $\mathit{F}(\mathit{x})=\sum _{\mathit{i}=1}^m(u_i^{\ast }-u_i)\mathit{K}(x_i,x_j)+\mathit{b}$(8)

where $\mathit{m}$, $\mathit{K}(x_i,x_j)$,$\mathit{u}$, and $\mathit{b}$ are the kernel function, number of support vectors, Lagrange multiplier, and offset, respectively. Commonly used kernel functions include the radial basis function (RBF), linear kernel, polynomial kernel, and sigmoid kernel. RBF kernel function is widely used, and this paper employs the RBF kernel function for SVM regression prediction.

Particle swarm optimisation (PSO) is a relatively recent evolutionary algorithm developed in recent years (Parsopoulos & Vrahatis, 2002). The essence of the PSO algorithm lies in the ability of individuals in the population to share information, allowing the entire population to transition from local to global optima in the solution space. The PSO algorithm treats each potential solution as a particle. After initialising a random particle swarm, it explores the solution space. Once it finds local and global optima, the particles update their velocities and new positions based on certain mathematical formulas through iterations to reach the global optimum. Therefore, PSO holds significant importance in the context of parameter optimisation. Considering the substantial impact of the penalty factor and kernel function parameters on the predictive accuracy of SVM, the PSO-SVM method aims to enhance the predictive capabilities of the method by employing PSO to optimise these two parameters (C. L. Huang & Dun, 2008). In this paper, the PSO algorithm has a population size of 5, and the maximum number of iterations is set to 100. The optimisation range for the kernel function parameter $\mathit{c}$ and the penalty factor $\mathit{g}$ is set to 0.1 to 100.

3.3. Basic principles of PSO-SVM-ICEEMDAN sea-level height retrieval method

As mentioned in the previous section, the accuracy of the GNSS-MR sea level height retrieval method based on the traditional method is limited by the antenna height. In contrast, PSO-SVM can build a regression prediction method by learning the external information autonomously without considering the external information, such as the antenna height. ICEEMDAN can decompose the signals more efficiently compared to EMD. Considering the advantages of the two methods, we propose a new sea surface height retrieval method that integrates PSO-SVM and ICEEMDAN as illustrated in Figure 2.

Figure 2. Flowchart of sea level height retrieval method based on PSO-SVM-ICEEMDAN.

Figure 2 shows the basic principle of the proposed method is to use PSO-SVM to establish the regression relationship between the peak frequency (obtained by ICEEMDAN decomposition) and the sea level height and then realise the retrieval of sea level height. The main process is as follows. First, the original SNR sequence is decomposed by ICEEMDAN to obtain some IMF components, and then the appropriate IMF component or combination is selected as the SNR oscillatory term, and then the peak frequency is obtained from the LSP spectrum analysis of the SNR oscillatory term. Then, a regression method between the frequency and sea level height is established by the PSO-SVM, which can be achieved by predicting the sea level height from the peak frequency to predict sea level height.

4. Experiments and discussion on PSO-SVM-ICEEMDAN-based sea level retrieval

Data from two GNSS station s will be carried out to validate the stability and accuracy of the proposed PSO-SVM-ICEEMDAN sea level retrieval model.

4.1. Sea level height retrieval experiment at SC02 station

We utilised SNR observations in the GPS L1 frequency band with a 15-second interval from the continuous operation of the GNSS tracking station SC02 (48.546°N, 123.007°W). The SC02 station is situated along the Washington Friday Harbor coast in the United States, with the antenna height elevated at least 3 metres above the distance from the ground to the sea surface. It is capable of receiving GNSS multipath interference signals originating from the sea surface. The receiver used at the SC02 station is the Trimble NETR9 equipped with a SCIT radome, and the antenna method is TRM59800.80. We utilised observed actual sea level height data from the Friday Harbor tide gauge station (48.545°N, 123.013°W), located at a distance of 359 metres from the SC02 station, as a reference to validate the accuracy and effectiveness of the proposed PSO-SVM-ICEEMDAN sea level retrieval method. Within the satellite elevation angle range of 5° to 25° and satellite azimuth range of 90° to 150°, we selected SNR observation data from the SC02 station for the entire year of 2021 to investigate the stability of GNSS-MR sea level retrieval at high satellite elevation angles. Figure 3 illustrates the environment and location of the SC02 station.

Figure 3. The left map shows the environment of station SC02 (facing east) and the right map shows the location of station SC02 (scale: 1:10km).

To comprehensively demonstrate the advantages of the PSO-SVM-ICEEMDAN sea level retrieval method, the ICEEMDAN decomposition results of the raw SNR observations for satellite G28 on the second day of 2021 at SC02 station are shown in Figure 4. Here the ICEEMDAN effectively decomposes the original SNR signal into IMF components ranging from high to low frequency, along with a residual component. This decomposition process highlights the capability of ICEEMDAN in extracting valuable information from the SNR data, which is crucial for sea level retrieval using the proposed method.

Figure 4. Graph of the original SNR decomposition results based on the ICEEMDAN algorithm.

The spectral diagrams of the oscillatory component of SNR obtained by polynomial fitting and wavelet decomposition (the wavelet function is set to db4, and the number of decomposition layers is 8), as well as the spectral diagrams of the IMF components obtained by ICEEMDAN decomposition, are compared in Figure 5. As discussed in Section 2, the antenna height is the sum of the vertical reflection distance and the tidal level. Therefore, the approximate antenna height can be obtained as 6.8 metres using polynomial fitting. Considering the observed sea level variations near the 2021 SC02 station, ranging from −1.04 metres to −3.10 metres, the peak frequency range of the sea surface multipath interference signals received at SC02 station can be derived as 39 to 82 Hz using formulas 5 and 6.

Figure 5. Comparison of spectrograms of SNR oscillatory terms obtained by polynomial fitting and wavelet decomposition with spectrograms of IMF components obtained by ICEEMDAN.

Figure 5 indicates that the peak frequency of the SNR oscillatory term obtained by polynomial fitting is significantly lower than the frequency range of the multipath interference signals. This discrepancy arises because, at high elevation angles, the multipath interference signals decrease, and polynomial fitting fails to adapt to this variation, resulting in the extracted the SNR oscillatory term containing not only interference signals but also certain direct signals. Since the amplitude of the direct signal is higher and the frequency is lower, the peak frequency of the SNR oscillatory term extracted by the polynomial method in Figure 5 is lower than the frequency range of the interference signals, obscuring the true peak frequency of the interference signals. In contrast, wavelet decomposition, capable of analysing non-stationary signals, decomposes the original signal into detail signals (high-frequency components) and approximation signals (low-frequency components) with outstanding signal decomposition capability. It can more accurately eliminate direct signals, and thus, the peak frequency of the extracted the SNR oscillatory term by wavelet decomposition falls within the frequency range of the interference signals.

ICEEMDAN decomposes the original signal into different IMF components and a residual component in an adaptive manner from high frequency to low frequency. Therefore, there are significant differences between different IMF components. The peak frequency of IMF1 is higher than the frequency range of the interference signals, indicating it as high-frequency noise; IMF4 and IMF5 have peak frequencies lower than the interference signal range, suggesting them as partial direct signals. Only the peak frequencies of IMF2 and IMF3 fall within the interference signal range, implying that both IMF2 and IMF3 components contain multipath interference signals. To prevent interference signals from being decomposed into different IMF components, ICEEMDAN reconstructs IMF2 and IMF3 as the SNR oscillatory term for tidal retrieval. The results in Figure 5 also show that the peak frequency of IMF2 + IMF3 lies within the interference signal range, confirming the correctness of reconstructing IMF2 and IMF3. Additionally, the peak frequency of the Sum of IMFs (IMF1 + IMF2 + IMF3 + IMF4 + IMF5) is found to be still lower than the frequency range of the multipath interference signals. This is because, at high elevation angles, interference signals decrease, and the amplitudes of IMF4 and IMF5, representing direct signals, are higher than the amplitude of interference signals. If reconstructed, their amplitudes would overshadow the interference signals, resembling a polynomial-like result, where the frequency of the direct signal masks the true frequency of the interference signals, ultimately leading to retrieval distortion and reduced stability.

Furthermore, Figure 6 compares the differences between the SNR oscillatory terms extracted by polynomial fitting, wavelet decomposition, and ICEEMDAN, in comparison to the clean multipath interference signal. From Figure 6, it can be observed that, in contrast to polynomial fitting and wavelet decomposition, the SNR oscillatory terms extracted by ICEEMDAN are more consistent with the multipath interference signal and exhibit higher correlation.

Figure 6. Comparison of the SNR oscillatory term extracted by polynomial fitting, wavelet decomposition, and ICEEMDAN with the clean multipath interference signal.

To comprehensively compare the advantages of the proposed PSO-SVM-ICEEMDAN method over other approaches, we simultaneously employed polynomial fitting, wavelet decomposition, ICEEMDAN, and PSO-SVM-ICEEMDAN to retrieve the tidal variations at the SC02 station in 2021, as illustrated in Figure 7. From Figure 7, it is evident that polynomial fitting results in numerous anomalous retrieval values, whereas the methods of wavelet decomposition, ICEEMDAN, and PSO-SVM-ICEEMDAN exhibit significantly fewer anomalous retrieval values. This validates the findings from Figures 5 & 6, indicating that the SNR oscillatory terms extracted by polynomial fitting contain components with larger amplitudes of direct signals, leading to inaccurate frequency estimation of the multipath interference signal and consequently generating a large number of anomalous retrieval values. In contrast, the wavelet decomposition and ICEEMDAN methods are effective in separating direct signals from interference signals, resulting in fewer anomalous retrieval values.

Figure 7. Comparison of the tidal retrieval values of polynomial, wavelet, ICEEMDAN and PSO-SVM-ICEEMDAN methods at station SC02 with the tide observations at tide gauge stations (TG denotes the tide level observations at tide gauge stations, and RV denotes the tidal retrieval values of different methods).

To further compare the retrieval accuracy of each method, Figure 8 conducts a regression analysis on the tidal retrieval values of each method and the observed values at the tidal gauge stations. From Figure 8, it is evident that the tidal retrieval values of the proposed PSO-SVM-ICEEMDAN exhibit the best linear regression properties and the lowest error among the four methods. Therefore, the retrieval accuracy of the PSO-SVM-ICEEMDAN model is the highest among the four methods, followed by ICEEMDAN, with polynomial fitting showing the lowest accuracy. The accuracy of wavelet decomposition falls between polynomial fitting and ICEEMDAN.

Figure 8. Linear regression analysis of the retrieved tide level values of different methods at station SC02 with the observed tide level values at tide gauge stations, and the heat map in the figure shows the absolute error between the retrieved values and the observed values.

As indicated by the analysis above, the poor stability of polynomial fitting, leading to a large number of anomalous retrieval values, decreases the final retrieval accuracy. To quantitatively describe the stability of different methods, we introduce the concept of a stability index (SI), which represents the percentage of effective retrieval values out of the total retrieval values. Considering that while machine learning can effectively improve retrieval accuracy, it may increase computational costs, we also compare the computational complexity of the four methods (i.e. algorithm runtime), as shown in Table 1. In Table 1, ‘R’ represents the correlation coefficient.

Table 1. Accuracy and stability statistics of different sea level height retrieval methods at SC02 station.

Method	Polynomial	Wavelet	ICEEMDAN	PSO-SVM-ICEEMDAN
RMSE	292.21	29.06	24.08	23.05
R	0.1968	0.929	0.9499	0.953
Time	1.6	513.18	704.83	715.35
SI	56.65%	95.64%	95.42%	96.13%

From Table 1, it is observed that, whether in terms of RMSE, correlation coefficient, or stability, wavelet decomposition, ICEEMDAN, and the proposed PSO-SVM-ICEEMDAN method outperform the traditional polynomial fitting. The stability and accuracy of the PSO-SVM-ICEEMDAN method are superior to wavelet decomposition and ICEEMDAN, with a stability of up to 96.13%, an RMSE of 23.05 cm, and a correlation coefficient of 0.953. Additionally, it can be noted that, due to the lower computational complexity of polynomial fitting, its time complexity is also the lowest. While the time complexity of wavelet decomposition and ICEEMDAN increases due to the multiple operations required for signal decomposition, the PSO-SVM-ICEEMDAN algorithm, despite incorporating machine learning for tidal regression prediction, does not significantly increase computational costs compared to wavelet decomposition and ICEEMDAN.

In order to comprehensively compare the differences between tidal data obtained by different methods and observed tidal data from tide gauge stations, we conducted tidal harmonic analysis on the tidal data obtained from tide gauge stations, polynomial fitting, wavelet decomposition, ICEEMDAN, and PSO-SVM-ICEEMDAN. The amplitudes and phases of their four main constituents (S2, M2, K1, and O1 tidal components) were calculated and are presented in Table 2. Figure 9 further illustrates the absolute errors of the tidal constituent amplitudes obtained by each method from the tidal gauge station’s amplitudes (phases show similar conclusions and are therefore not listed). Combining Table 2 and Figure 9 reveals that polynomial fitting exhibits the largest differences in all four tidal constituents compared to the tide gauge station, with K1 tidal constituent showing a deviation as high as 35.37 cm. Although the other three methods also exhibit some deviations from the tide gauge station results, these deviations are smaller than those of polynomial fitting. In contrast, the constituents obtained by the other three methods exhibit smaller deviations from the tide gauge station results. Specifically, in the comparison of amplitudes for the M2 and O1 constituents, PSO-SVM-ICEEMDAN shows the smallest deviation, while in the comparison of amplitudes for the S2 and K1 constituents, wavelet decomposition exhibits the smallest deviation. Therefore, at the SC02 station, compared to polynomial fitting and ICEEMDAN, the tidal retrieval values obtained by wavelet decomposition and PSO-SVM-ICEEMDAN are more consistent with the observed values at the tide gauge station in terms of tidal patterns. Additionally, considering the accuracy analysis mentioned earlier, we can conclude that the PSO-SVM-ICEEMDAN method achieves higher accuracy and stronger stability in tidal retrieval with high elevation angle and long period without significantly increasing computational costs. It is closer to the tidal variation patterns observed at the tide gauge station.

Figure 9. Absolute errors of tidal constituents from different sea level retrieval methods at station SC02.

Table 2. Comparison of tidal constituents obtained from different sea level height retrieval methods at station SC02.

Tidal constituents		TG	Polynomial	Wavelet	ICEEMDAN	PSO-SVM-ICEEMDAN
S2	Amplitude(m)	0.0842	0.1326	0.068	0.0489	0.058
	Phase(°）	38.76	207.55	90.95	29.32	38.35
M2	Amplitude(m)	0.3435	0.1064	0.3538	0.372	0.3473
	Phase(°）	13.11	44.51	11.34	14.08	12.09
K1	Amplitude(m)	0.7664	0.4127	0.7658	0.7037	0.676
	Phase(°）	287	310.11	291.92	289.45	289.25
O1	Amplitude(m)	0.355	0.0949	0.2923	0.2973	0.3103
	Phase(°）	250.61	313.46	262.79	244.57	251.65

4.2. Sea level height retrieval experiment at TPW2 station

Identical experiments for sea level height retrieval were conducted to validate the proposed method across different sites. The continuous GNSS tracking station TPW2 (46.207°N, 123.768°W), located in Oregon, U.S.A., served as the research subject. TPW2 station employs a Trimble NETRS receiver, SCIS antenna radome, and TRM29659.00 antenna method, providing SNR observations in the GPS L1 frequency band with a 15-second interval. Measured sea level heights were obtained from the nearby Astoria tide gauge station (46.207°N, 123.768°W). Figure 10 shows the TPW2 station environment and its location. Similar to the SC02 station, SNR observations for sea level height retrieval were selected for the TPW2 station for the entire year of 2022, limited to satellite elevation angles between 5° and 25° and satellite azimuth angles between 170° and 250°.

Figure 10. The left map shows the environment of station TPW2 (facing west) and the right map shows the location of station TPW2 (scale: 1:10km).

Similar to the SC02 station, we applied the traditional polynomial, Wavelet, ICEEMDAN, and PSO-SVM-ICEEMDAN methods to retrieve sea level heights for the TPW2 station in 2022. We compared these results with the measurements from the tide gauge station, as illustrated in Figure 11.

Figure 11. Comparison of the tidal retrieval values of polynomial, wavelet, ICEEMDAN and PSO-SVM-ICEEMDAN methods at station TPW2 with the tide observations at tide gauge station.

From Figure 11, it is evident that, similar to the results at station SC02, the stability of the polynomial fitting method remains poor, with numerous anomalous retrieval values. It is noteworthy that the wavelet decomposition method also exhibits a considerable number of anomalous retrieval values at station TPW2. In contrast, the retrieval values obtained by the ICEEMDAN and PSO-SVM-ICEEMDAN methods show a high degree of consistency with the observed values at the tide gauge station. Similarly, to further explore the tidal retrieval capabilities of different methods, Figure 12 presents regression analyses of the retrieval values and observed values at the tide gauge station for polynomial fitting, wavelet decomposition, ICEEMDAN, and PSO-SVM-ICEEMDAN.

Figure 12. Linear regression analysis of the retrieved tide level values of different methods at station TPW2 with the observed tide level values at tide gauge stations.

The results in Figure 12 indicate that, unlike station SC02, the stability of the wavelet method significantly decreases in tidal retrieval at station TPW2. Compared with SC02, the R-square and slope of its regression model were significantly lower while the intercept and RMSE were significantly higher. This may be related to the GNSS station environment, as seen in Figures 3 and 10. In contrast to the relatively homogeneous environment near station SC02, the environment around TPW2 is more complex, being situated at the mouth of the Columbia River with heavy maritime traffic. This complexity poses significant challenges for extracting interference signals, resulting in poorer stability of the wavelet method at TPW2.

Table 3 summarises stability indicators, accuracy indicators, and computational complexity indicators for the four methods. From Table 3, it is evident that, while wavelet decomposition maintains high stability, its accuracy decreases significantly, falling below that of ICEEMDAN and PSO-SVM-ICEEMDAN. The PSO-SVM-ICEEMDAN method demonstrates good general applicability, maintaining optimal accuracy and stability, with an RMSE of 25.17 cm, a correlation coefficient of 0.9458, and stability of 97.53%. Furthermore, compared to station SC02, the stability of polynomial fitting at TPW2 is significantly reduced, confirming the difficulty of extracting multipath interference signals due to the more complex environment at TPW2. The computational complexity of each model is similar to that at SC02; although PSO-SVM-ICEEMDAN achieves the best stability and accuracy using machine learning algorithms, its computational complexity does not increase significantly compared to ICEEMDAN and wavelet.

Table 3. Accuracy and stability statistics of different sea level height retrieval methods at TPW2 station.

Method	Polynomial	Wavelet	ICEEMDAN	PSO-SVM-ICEEMDAN
RMSE	383.68	148.23	38.07	25.17
R	0.2318	0.4479	0.8916	0.9458
Time	1.33	592.98	737.06	743.2
SI	33.81%	77.22%	97.04%	97.53%

Similar to station SC02, we analysed the four main tidal constituents obtained by the four methods and the tide gauge station at TPW2, as shown in Table 4 and Figure 13. Combining Table 4 and Figure 13, it is evident that, due to the more complex environment at TPW2, the absolute errors in the amplitudes of the four constituents obtained by polynomial fitting are much larger than those at SC02, with the absolute error in the amplitude of the K1 constituent reaching as high as 367.16 cm. Wavelet decomposition exhibits lower errors only in the M2 constituent, while in the comparison of amplitudes for the other three constituents, the PSO-SVM-ICEEMDAN method shows the smallest absolute errors, all below 4 cm. Therefore, it can be concluded that, in the one-year tidal retrieval experiment at TPW2, PSO-SVM-ICEEMDAN demonstrates optimal stability and accuracy, and its tidal retrieval values are more consistent with the observed values at the tide gauge station in terms of tidal variations.

Figure 13. Absolute errors of tidal constituents from different sea level retrieval methods at station TPW2.

Table 4. Comparison of tidal constituents obtained from different sea level height retrieval methods at station TPW2.

Tidal constituents		TG	Polynomial	Wavelet	ICEEMDAN	PSO-SVM-ICEEMDAN
S2	Amplitude(m)	0.1228	0.7373	0.2607	0.1853	0.0979
	Phase(°）	282.87	306.35	298.31	308.77	264.3
M2	Amplitude(m)	0.5801	0.9312	0.5936	0.5251	0.5596
	Phase(°）	263.51	170.58	233.38	255.27	265.85
K1	Amplitude(m)	0.4324	4.104	2.2262	0.5862	0.439
	Phase(°）	258.61	322.93	316.43	291.68	260.8
O1	Amplitude(m)	0.0457	1.1265	0.1738	0.0818	0.0781
	Phase(°）	224.44	156.18	196.64	136.46	274.91

5. Conclusions

This study conducted a one-year sea surface height inversion experiment at the SC02 and TPW2 stations in the United States, successfully validating the effectiveness and superiority of the PSO-SVM-ICEEMDAN method in GNSS-MR sea surface height inversion. A comparative analysis was performed with the retrieval results of the new method against polynomial fitting, wavelet decomposition, and ICEEMDAN. Experimental results indicated that polynomial fitting exhibits lower stability and accuracy in retrieval at high elevation angles. While wavelet decomposition achieved favourable retrieval results at the SC02 station, it is susceptible to station environmental influences. For ICEEMDAN to achieve high-precision tidal retrieval, accurate antenna height information is required. The new method effectively overcomes these limitations. Without precise antenna height information and without significantly increasing computational complexity, the new method maintained good retrieval accuracy and stability at both SC02 and TPW2 stations, with RMSE values of 23.05 cm and 25.17 cm, and stability reaching 96.13% and 97.53%, respectively. These performances surpassed those of polynomial fitting, wavelet decomposition, and ICEEMDAN. Furthermore, the tidal variations obtained by the new method were more consistent with tide gauge station observations. This contribution provides insights for achieving long-period, high-stability tidal retrieval in complex environments, high elevation angles, and without precise antenna height information.

It is noteworthy that the new method still requires the polynomial method to approximately determine GNSS antenna height for selecting suitable IMF components. This may pose limitations on the extensive applicability of the new method. Our future research will focus on realising adaptive IMF component selection to address this limitation.

Author contributions

Conceptualization, X.W. and L.J.; methodology, L. J. and X. W.; software, L. J. and X. W; validation, X. W., L. J. and H.W.; formal analysis, L. J., X. W., H. H. and L. Y.; investigation, L. J. and X. W.; resources, X. W. and H.W.; data curation, X. W. L. J. and H. H.; writing – original draft preparation, L. J. and X. W.; writing – review and editing, L. Y., X. W. and H. W.; visualisation, X. W. and L. J.; supervision, H. W. and X. W.; project administration, H. H. and X. W.; funding acquisition, X. W., H. W. and L. Y. All authors have read and agreed to the published version of the manuscript.

Acknowledgments

The authors will thank all editors and reviewers for their valuable comments and suggestions on this manuscript. The authors will also thank the professors and doctoral students in GNSS Research Center, Wuhan University for advice on essay revising and editing.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The data included in this study are freely available through the following web pages: EarthScope Consortium provided GNSS observations from stations SC02 and TPW2 (https://data.unavco.org/), and the National Oceanic and Atmospheric Administration provided sea surface tide observations near stations SC02 and TPW2 (https://tidesandcurrents.noaa.gov/). Most of the graphs in this paper were plotted by MATLAB software.

Additional information

Funding

This paper is funded by the National Natural ScienceFoundation of China (Grant Nos. 52364008 and 42301440), Key Laboratory of Geo-space Environment and Geodesy, Ministry of Education, Wuhan University (20-01-02), and Introduced Talents Research Fund Project of Guizhou University (2020043), Cultivation Project of Guizhou University, China (No. Gui Da Pei Yu [2020]57).