Sea Surface Height Estimation with Multi-GNSS and Wavelet De-noising

This paper presents a new sea surface height (SSH) estimation using GNSS reflectometry (GNSS-R). It is a cost-effective remote sensing technique and owns long-term stability besides high temporal and spatial resolution. Initial in-situ SSH estimates are first produced by using the SNR data of BDS (L1, L5, L7), GPS (L1, L2, L5), and GLONASS (L1, L2), of MAYG station, which is located in Mayotte, France near the Indian Ocean. The results of observation data over a period of seven days showed that the root mean square error (RMSE) of SSH estimation is about 32 cm and the correlation coefficient is about 0.83. The tidal waveform is reconstructed based on the initial SSH estimates by utilizing the wavelet de-noising technique. By comparing the tide gauge measurements with the reconstructed tidal waveform at SSH estimation instants, the SSH estimation errors can be obtained. The results demonstrate that the correlation coefficient and RMSE of the wavelet de-noising based SSH estimation is 0.95 and 19 cm, respectively. Compared with the initial estimation results, the correlation coefficient is improved by about 14.5%, while the RMSE is reduced by 40.6%.

Obtaining accurate Sea Surface Height (SSH) is significant for human living, especially for those live along ocean coasts. Furthermore, it is also important to coastal ecosystems and coastal morphology 1,2 . SSH has been obtained primarily with tide gauges during the last centuries. While for the last three decades, satellite altimetry has been the dominating technique [3][4][5] . The satellite altimetry has unique advantages in monitoring sea surface topography in open ocean research such as ocean circulation 6 . However, due to tidal aliasing effect, it's spatial and temporal resolution are insufficient to observe complex and rapidly changing dynamics which makes it difficult to use near the coast 7 . Meanwhile, tide gauges are affected by both sea level and land surface changes since the measurements are related to a benchmark on the land, where they were established. For these reasons, it is difficult to use traditional tide gauges for sea-level studies in tectonically active areas or applications related to changes in the global ocean volume, e.g., the global sea level budget 8 . These applications need absolute observations of sea level, in other words, measurements with respect to a terrestrial reference frame 9 .
In recent years, the reflected signals of Global Navigation Satellite System (GNSS) are used to retrieve a range of geophysical parameters such as SSH, soil moisture, ocean wind speed, etc. [10][11][12][13][14][15][16][17][18] . This technique is also termed GNSS-Reflectometry (GNSS-R), which exploits conventional receivers or low-cost, low-mass, and low-power characteristics of GNSS-R receivers. Different GNSS-R based methods have been proposed in the literature to measure a range of geophysical parameters. Basically, either a single antenna or two separate antennas are used for the direct and reflected signal reception. A number of researchers have investigated GNSS-based ocean surface altimetry by using the measurement of signal arrival time obtained by either mounting the receiver on an aircraft or fixing it on the ground such as on a bridge 19 23 . They successfully demonstrated its ability to estimate local SSH and improved accuracy. Nevertheless, many of their studies used single frequency or single system. With constantly developing and perfecting of other navigation satellite systems such as Galileo and BDS, can we improve the accuracy of SSH estimation using multi frequencies or multi-systems? There are still many challenging problems to be resolved 24 .
In this paper, we propose an improvement of the method based on the SNR analysis of a single antenna. And the focus is on GNSS-R ocean surface altimetry based on the SNR data of multiple satellite constellations and multiple frequencies, namely BDS (L1, L5, L7), GPS (L1, L2, L5) and GLONASS (L1, L2). The purpose is to exploit measurements diversity to achieve a performance gain. The structure of the next part of this paper is organized as follows. Firstly, the fundamental theory of SSH GNSS-R and contains a wavelet analysis in detail are given. Thereafter the data introduction is presented. Then SSH estimation results by using individual GNSS-R and multi-GNSS-R are firstly shown. The SNR data collected by an IGS station are used for evaluation and analysis. Tide gauge observations are used as ground-truth data for performance comparison. Furthermore, the initial discrete SSH estimation results are used to reconstruct the tidal waveform by using wavelet de-noising technique. Finally, the concluding remarks are declared.

Theory and Methods
Methods of GNSS-R. Multipath propagation is one of the main GNSS error sources, which constrains high-precision positioning. The multipath caused by the difference of the phase between the direct and reflected signals will affect GNSS observations and give rise to oscillations in the observations. As a result, the SNR data recorded by a GNSS receiver with a single antenna will contain information about the interference, antenna height and hence SSH variation. Figure 1 shows the diagram of SSH estimation based on the SNR method. Here h(i) is the antenna height relative to reflection surface, and θ(i) is the angle between the direct signal and the instantaneous sea surface of a specific scattering point. In Fig. 1, the reflected signal has one excess phase delay when compared to the direct signal, which is related to the antenna height h(i). Figure 2 shows an example of the SNR variations recorded by a BDS receiver. The thin blue fast-varying curve is the observed data, while the red thick smooth curve is the fitted SNR of the direct signal. As we can see, the multipath SNR is proportional to the difference between the observed curve and the fitted curve. It can be seen that the multipath SNR during the rising and setting of the satellite is significantly higher than that elsewhere, as shown in the two black ellipses.
SNR is one of the main observables of the GNSS receiver, which is mainly related to the satellite signal transmission power, antenna gain, and multipath effect. As the elevation angle increases, the antenna gain raises, subsequently, the direct signal SNR increases. In this paper, we use the observed SNR data from the GNSS receiver to estimate SSH. As described in 22 , the SNR is given by  www.nature.com/scientificreports www.nature.com/scientificreports/ where P n is the noise power, A d indicates amplitude of direct signal, A r refers amplitude of the reflected signal, δφ(t) is the reflection excess phase with respect to direct phase (also known as the interferometric phase), given by where h stands for antenna height from the specular scattering point on the reflected surface, λ is the carrier wavelength, θ(t) represents the elevation angle. Since (i) GNSS antennas are designed to filter reflected signals, (ii) the reflected signal is attenuated upon reflection, and (iii) the azimuth angles we choose to make reflected signals come from the sea surface, we can assume that Therefore, the direct signal determines the overall trend of the SNR observations (see Fig. 2). So, we can use a low-order polynomial to fit the observed SNR data and then remove the trend. As a result, a detrended SNR time series is produced, which can be used to estimate the surface parameters such as SSH. As described in 22 , after removing the SNR trend, the remaining SNR can be approximated as where ϕ indicates a phase offset, A denotes the amplitude which is given by By defining ω = sin θ and f = 2h/λ, Eq. (4) can be rewritten as Due to the tidal effect, the antenna height relative to the surface may change smoothly by around 0.72 m during half an hour at the selected observation station. Thus, the frequency of the detrended SNR signal varies gently with time. However, the mean frequency of the signal would be dominant, which basically corresponds to the antenna height at the middle of the observation time period.
From Eq. (6) and the definition of f, the antenna height is associated with the frequency of the detrended SNR signal by Spectral analysis is applied to obtain the spectral peak frequency of the time series. In this paper, the Lomb-Scargle periodogram (LSP) method is used since it can handle unevenly spaced samples 22 . The recorded SNR data are typically evenly sampled in time, but the sine of elevation angle is not evenly distributed. The LSP is a well-known algorithm for detecting and characterizing periodicity in unevenly-sampled time-series, even more, has seen particularly wide use within the astronomy community 25 . As described in 25 , the LSP method calculates the power spectral density by where X and δ 2 are the mean and the variance of the observed sequence respectively, ω denotes the angular frequency, and τ symbolizes the phase which can be inferred by As an example of obtaining the peak frequency using LSP method, Fig. 3 shows the results of LSP analysis of the detrended SNR time series from the SNR data observed by BDS satellite PRN09 on two different days. The elevation angles change from 5° to 20° and from 20° to 5° for day of year (DOY) 191 and DOY 193, respectively. The peak spectral frequency can be clearly observed. The peak frequency difference is caused by the variation in the SSH over the two various observation periods. The estimated peak frequency, besides, Eq. www.nature.com/scientificreports www.nature.com/scientificreports/ where a indicates the scale parameter, b refers to the time parameter, ψ(t) is an analyzing wavelet, and ψ • ( ) symbolizes the complex conjugate of ψ(•). The inverse transform of CWT is given by 1/2 2 where a, b and t are continuous. However, they must be discretized in analyzing real data. The discrete method is called discrete wavelet transform (DWT). In practice, the scale a and the time b are discretized as following j j 0 00 where j and k are integers. And the continuous wavelet function W ψ f(a, b) become the DWT, given by The inverse transform of DWT is computed by Mallat proposed a fast DWT algorithm in 1989 28 , which decompose and reconstruct the signal using the wavelet filters. The decomposition algorithm is given by where t = 1, 2 …, N is the discrete time serial number, N denotes the signal length, f(t) symbolizes the original signal, j = 1, 2 …, J is the decomposition level, J is the maximum decomposition level; H and G are the wavelet low-frequency and high-frequency pass decomposition filters, respectively; A j and D j are the wavelet coefficients of f(t) in the low-frequency and the high-frequency part of the Jth layer, respectively. Figure 4 shows the wavelet decomposition filter model.  www.nature.com/scientificreports www.nature.com/scientificreports/ Each signal has the following relationship: S = A1 + D1; A1 = A2 + D2; A2 = A3 + D3; A3 = A4 + D4. As an example of decomposition, Fig. 5 shows the results of wavelet decomposition under Db6 wavelet with decomposition level 3, and the original signal S is the in-situ SSH estimation of this paper. In addition, the relationship of each signal is: The reconstruction is given by where h and g represents the wavelet low pass and high pass reconstruction filter, respectively. Other symbols are the same as Eq. (15). Figure 6 shows the wavelet reconstruction filter model. One important step of wavelet de-nosing is choosing the appropriate wavelet function W ψ f(a, b). Wavelet functions such as Daubechies (dbN), Haar, and Symlets (SymN) are commonly used for wavelet analysis. To choose a suitable wavelet function, the performance of db6, Haar and Sym2 at different levels are analyzed in section "Improved SSH estimates with Wavelet De-noising".  Tide gauge observations. The observation data produced by the Dzaoudzi tide gauge, which is about ten meters from the MAYG station, were used as the ground-truth data for evaluating the performance of the proposed SSH estimation method. The sampling rate of the Dzaoudzi tide gauge data is 1 min, which can be downloaded from the Intergovernmental Oceanographic Commission (IOC) website (http://www.ioc-unesco.org/). Figure 7 shows the in-situ SSH estimation results using the L1 signals of GPS, GLONASS, and BDS. The horizontal axis represents the DOY, and the tide gauge observations (black solid line) are also shown for comparison. The GNSS-R based SSH estimates are represented by solid circle, solid square, and solid inverted triangle, respectively for signals of three different satellites.

SSH estimation results with Multi-GNSS.
As can be seen from Fig. 7, there is a significant daily periodicity of in-situ SSH changes observed by tide gauge, which is mainly caused by tide and wind variability, as well as sea level pressure variability. There is a good agreement between tide gauge SSH observations and SSH estimates obtained by the GNSS-R methods over a period of seven days. Only a limited number of GNSS satellites can be seen over a certain period of time. Thus, SSH estimations are not uniformly distributed. In some cases, no in-situ SSH estimates are available over a quite long time interval. In general, the number of GPS satellites is the largest among the numbers of GPS, BDS, and GLONASS satellites. So, GPS-R has the highest temporal resolution for SSH estimation on average. Figure 8 shows the combined results of SSH estimates with all available GNSS satellites and frequencies. It is interesting to note that the time resolution of SSH is greatly improved. Table 1 presents the in-situ SSH estimation performance in terms of RMSE and correlation coefficients with GPS, BDS, and GLONASS satellite signals of eight different frequency bands. It can be seen from Table 1 that GNSS-R with GPS L1 signal achieves the best performance, i.e., maximum correlation coefficient (0.91) and highest accuracy (RMSE = 0.23 m). Figure 9 shows the histogram of the in-situ SSH estimation errors of the three different navigation systems. It can be seen that the errors approximately have a normal distribution. The mean  Figure 10 shows the scatter plot between the tide gauge observations and the GNSS-R SSH estimates. Clearly, the GNSS-R estimates are closely around the tide gauge observations. The RMSEs and correlation coefficients of the combined GNSS-R based estimates are 0.32 m and 0.83, respectively. Although the precision and correlation coefficient are slightly degraded compared with the best results using single frequency band GPS L1, the temporal resolution is increased significantly. Improved SSH estimates with Wavelet De-noising. As observed from Fig. 8, compared with the tide gauge observations, the GNSS-R based estimates are rather noisy. In order to improve the GNSS-R based in-situ SSH estimation precision, we exploit the wavelet de-noising method to reconstruct the tidal waveforms. To produce the best possible wavelet de-noising results, a number of parameters should be selected appropriately 29,30 . Three of the wavelets i.e. Db6, Haar, and Sym2 are used to reconstruct the tidal waveform. The number of levels for decomposition and reconstruction is set to be 1-8. Table 2 presents the RMSEs and correlation coefficients by using the three different wavelets. It can be seen that when the decomposition level is less than 5, the three different wavelets produce very similar performance. Db6 slightly outperforms the other two wavelets and Db6 with decomposition level of 3 achieves the minimum RMSE. Meanwhile, the decomposition level is greater than 4, the performance of RMSEs and correlation coefficients are degraded, as well as Sym2 with decomposition level of 8 achieves the maximal RMSE. Therefore, only results produced by Db6 wavelet with decomposition level 3 are presented in the following. Figure 11 shows the reconstructed tidal waveforms after wavelet de-noising. For comparison, the tide gauge observations, as well as the initial GNSS-R in-situ SSH estimates are also presented. Basically, there is a nice match between the observed and the reconstructed tidal waveforms. The histogram of SSH estimation errors is shown in Fig. 12. By comparing with the results have shown in Fig. 9, the in-situ SSH estimation error is significantly reduced and their distribution is more reasonable. The correlation coefficient and RMSE of the wavelet de-noising method are respectively 0.95 and 0.19 m, which are improved by 14.5% and 40.6% over the initial multi-GNSS-R in-situ SSH estimation results.  www.nature.com/scientificreports www.nature.com/scientificreports/    Table 2. The precisions and correlation coefficients after wavelet de-noising with different levels. Figure 11. Time series of in-situ SSH estimates from initial results and de-noising results.