GNSS-IR Model of Sea Level Height Estimation Combining Variational Mode Decomposition

—The global navigation satellite system-reﬂections (GNSS-R) signal has been conﬁrmed to be useful for retrieving sea level height. At present, the GNSS-interferometric reﬂectometry (GNSS-IR) technology based on the least square method to process signal-to-noise ratio (SNR) data is restricted by the satellite elevation angle in terms of accuracy and stability. This article proposes a new GNSS-IR model combining variational mode decomposition (VMD) for sea level height estimation. VMD is used to decompose the SNR data into intrinsic mode functions (IMF) of layers with different frequencies, remove the IMF representing the trend item of the SNR data, and reconstruct the remaining IMF components to obtain the SNR oscillation item. In order to verify the validity of the new GNSS-IR model, the measurement data provided by the Onsala Space Observatory in Sweden is used to evaluate the performance of the algorithm and its stability in high-elevation range. The experimental results show that the VMD method has good results in terms of accuracy and stability, and has advantages compared to other methods. For the half-year GNSS SNR data, the root mean square error and correlation coefﬁcient of the new model based on the VMD method are 4.86 cm and 0.97, respectively.


GNSS-IR Model of Sea Level Height Estimation Combining Variational Mode Decomposition
Yuan Hu, Xintai Yuan, Wei Liu , Jens Wickert , Zhihao Jiang, and Rüdiger Haas Abstract-The global navigation satellite system-reflections (GNSS-R) signal has been confirmed to be useful for retrieving sea level height.At present, the GNSS-interferometric reflectometry (GNSS-IR) technology based on the least square method to process signal-to-noise ratio (SNR) data is restricted by the satellite elevation angle in terms of accuracy and stability.This article proposes a new GNSS-IR model combining variational mode decomposition (VMD) for sea level height estimation.VMD is used to decompose the SNR data into intrinsic mode functions (IMF) of layers with different frequencies, remove the IMF representing the trend item of the SNR data, and reconstruct the remaining IMF components to obtain the SNR oscillation item.In order to verify the validity of the new GNSS-IR model, the measurement data provided by the Onsala Space Observatory in Sweden is used to evaluate the performance of the algorithm and its stability in high-elevation range.The experimental results show that the VMD method has good results in terms of accuracy and stability, and has advantages compared to other methods.For the half-year GNSS SNR data, the root mean square error and correlation coefficient of the new model based on the VMD method are 4.86 cm and 0.97, respectively.

I. INTRODUCTION
A S THE global temperature rises, the continuous melting of glaciers has caused the average sea level height to rise, which has brought many adverse effects to coastal countries.Therefore, effective and accurate monitoring of sea level height has important practical significance.Traditional sea level height monitoring based on tide gauges has certain limitations, such as limited monitoring area, high cost, and susceptibility to external environmental influences.The global navigation satellite system-reflectometry (GNSS-R) signal proposed by Hall and Cordey is a new branch of GNSS that has been gradually developed since the 1990s and is one of the research hot spots in the field of remote sensing detection and navigation technology [1].Since the discovery of the surface reflection signal, which was originally a source of multipath error, can be used as a new type of remote sensing signal source for the inversion of the characteristic parameters of the surface reflection surface, many scholars launched the research on GNSS-R signal [2]- [5].At present, the research of sea level height measurement based on reflected signal mainly includes phase delay analysis method [6]- [9] and SNR analysis method.SNR data only need to be obtained by a single antenna [10], which has lower requirements for geodetic receivers.At the same time, the SNR signal formed by the interference of the direct signal and the reflected signal has better robustness to wind and waves.Larson et al. proposed GNSS-interferometric reflectometry (GNSS-IR) technology and used it for surface monitoring, including soil moisture [11] and height measurement [12].GNSS-IR technology can be used to detect the earth's surface environment, and has the characteristics of wide coverage, low cost, and multiple signal sources, which provides a new possibility for monitoring the sea level height [13]- [17].Therefore, the sea level height inversion method based on the GNSS-IR model has received widespread attention [18]- [20].
In the study of the GNSS-IR model, Larson et al. [21] combined the least square method and Lomb-Scargle Periodogram (LSP) spectrum method to measure the sea level height of two stations.They achieved the sea level heights with root-meansquare error (RMSE) of 5-10 cm and the correlation coefficients were greater than 0.97 [22].Löfgren et al. [23] used SNR data to measure the sea level height at five stations around the world, and found that the inversion results were not systematically different from the tide gauge measurements.Löfgren et al. [24] used SNR analysis and phase delay analysis to measure the height of the rough sea level height, and the results showed that SNR analysis performed better.Larson et al. [25] pointed out the term dynamic sea level height correction and applied it to Kachemak Bay where the daily sea level height changes greater than 7 m, and achieved the sea level heights with RMSE of 2.3 cm.Strandberg et al. [26] proposed the B-spline method to measure sea level height through GNSS combination to determine the parameters, which further improved the measurement accuracy, reaching a standard deviation of 1.4 cm at Onsala, Sweden, and 3.1 cm at Spring Bay, Tasmania.Santamaría-Gómez et al. [27] improved the measurement accuracy by using an extended Kalman filter/smoothing algorithm and tropospheric refraction correction processing, and obtained an average level difference of about 5 mm and an RMSE of about 3 cm.Zhang et al. [28] found that the empirical mode decomposition (EMD) can be used to detrend the SNR data, which improves the utilization of SNR data, and can still obtain inversion results with higher accuracy in high-elevation range [29].Wang et al. [30] used a combination of four-constellation multi-GNSS multipath reflection method to estimate the sea level height of the BRST station, which improved the accuracy by about 40%-75%.Wang et al. [31] used wavelet decomposition to extract the frequency of SNR data for sea level height measurement, and obtained inversion results that are not much different from LSP spectrum analysis.Jin et al. [32] first applied BDS-reflectometry to estimate sea level height changes based on SNR data and triple-frequency phase and code combinations, and they achieved good agreement with the observation results of the tide gauge.The SNR-based GNSS-IR technology mainly uses the least square method and LSP spectrum analysis to invert the sea level height altimetry, and the research focuses on improving the inversion accuracy and stability.In addition, there are methods for processing SNR data such as the B-spline method [26], the EMD method [29], wavelet decomposition [33], and singular spectrum analysis (SSA) [34].Although SNR-based GNSS-IR technology has gained a lot of progress, the existing GNSS-IR sea level height estimation model still has insufficient accuracy, stability, and utilization of GNSS data.
To solve this problem, this article proposes a new GNSS-IR model of sea level height estimation based on variational mode decomposition (VMD) [35] method.This new model uses the VMD method to decompose the SNR data, and combines the LSP spectrum analysis to extract the oscillation frequency of the SNR.VMD can decompose SNR data into residual sequence and intrinsic mode function (IMF) components representing different local features.Through observation and comparison, we can find the IMF component representing the trend of the SNR data, and then reconstruct the remaining IMF components to obtain the oscillation term of the SNR data, and then LSP spectrum analysis extracts the oscillation frequency of the SNR sequence to obtain the sea level height.Compared with the traditional model, the new model improves the utilization of SNR data and ensures that the algorithm always maintains high stability and accuracy at high-elevation angles.
The rest of this article is organized as follows.Section II presents the basic principle of the GNSS-IR model.Section III introduces the VMD method.Following, in Sections IV and V, the model verification test and discussion are presented.Finally, Section VI concludes this article.

II. BASIC PRINCIPLE OF GNSS-IR MODEL
The basic GNSS-IR model of sea level height estimation is shown in Fig. 1.The SNR data received by the GNSS antenna is formed by the interference of the direct signal and the reflected signal is a measurement of signal strength.In Fig. 1, h is the vertical reflection distance, which represents the vertical distance from the center of the antenna to the sea surface, Δd is the path delay between the reflected signal and the direct signal, and θ is the angle between the direct signal and the sea surface.
According to a near-ground direct reflection combined signal power expression proposed by Nevinski and Larson, SNR can be expressed as [36] where A 2 d and A 2 r are the direct signal power and the reflected signal power, respectively, and φ γ is the phase delay between the direct signal and the reflected signal.According to Fig. 1, φ γ can be expressed as where λ is the carrier wavelength, h r is the vertical reflection distance, and θ is the satellite elevation angle.In the expression of the SNR, A 2 d + A 2 r is the direct and multipath signal power trend terms, and 2A d A r cos φ γ is the SNR oscillation term caused by direct and multipath signal interference.
In the coastal sea level height measurements, GNSS-IR techniques based on SNR analysis can eliminate the Doppler effect caused by receiver dynamics.At the same time, the polarization component of the reflected signal changes with the elevation angle.When the elevation of the navigation satellite is low, the proportion of the right-handed polarization component contained in the reflected signal is high, and the lower elevation, the more coherent components of the reflected signal, and the smoother the sea surface represented by the signal.Therefore, the SNR sequence with elevation angles of 5°-12°is often processed to invert surface parameters.However, sometimes enough time series are needed for analysis.The SNR sequence with satellite elevation angles between 5°and 30°have better inversion results, and the algorithm must comprehensively consider the trajectory of the navigation satellite and the location of the station.The SNR data are composed of power trend terms and SNR oscillation terms [37], and the oscillation amplitude, frequency and phase of the SNR oscillation term have a direct correspondence with the surface parameters.Therefore, in order to obtain the GNSS multipath information caused by the surface reflection in the SNR, some methods can be used to remove the trend item of the SNR sequence in the low altitude range, thereby providing the detrended SNR data.This phenomenon can be easily observed in Fig. 2. As the elevation increases, the amplitude of the oscillation will gradually decrease.The frequency f of the SNR oscillation term can be obtained by performing LSP spectrum analysis of the detrended SNR data, and then the vertical reflection distance h reflecting the change of sea surface height can be obtained by f = 2h/λ.

III. INTRODUCTION TO VMD
The innovation of the new GNSS-IR model proposed in this article is to use the VMD method to decompose the SNR data to provide high-quality SNR data for LSP spectrum analysis.In the traditional GNSS-IR model, the detrend processing of SNR data is very critical.The traditional method is to process the SNR sequence through the least square method to obtain the SNR oscillation term.In addition, common SNR data processing methods include wavelet decomposition, EMD and SSA methods.In this article, we select an original SNR sequence on the first day of 2016 of the Onsala Space Observatory, and extracted the SNR observation data at an elevation angle of 5°-30°.It can be seen that the SNR sequence in low-elevation angle range is an oscillating sequence with a trend term.In Fig. 3, the horizontal axis represents the satellite elevation angle, and the vertical axis represents the corresponding SNR amplitude.Obviously, the SNR oscillation term obtained by the least square method is quite complicated, and contains a large number of high-frequency and low-frequency signals, which will introduce errors in the spectrum analysis.Compared with the least square method, the SNR oscillation term obtained by processing the SNR data by the VMD method is quite smooth, which is beneficial to the subsequent LSP spectrum analysis.Therefore, this article proposes to use VMD to decompose the original SNR data to obtain the SNR oscillation term corresponding to the frequency of the coherent signal, and then compare the inversion accuracy of different SNR data processing methods.
The VMD method is a completely nonrecursive decomposition model.In theory, any type of time series can be decomposed into modes with different center frequencies by the VMD, including nonstationary time and nonlinear time [38].Therefore, VMD can be summarized as where X(t) is the original time series; u k (t) is the IMF; r n (t) is residual term.The IMF of VMD is defined as an amplitude-modulatedfrequency-modulated signal [35].The kth mode u k (t) is written as where A k (t) is the instantaneous amplitude; φ k (t) is the instantaneous phase, and its derivative ω k (t) = φ k (t) is the instantaneous frequency.
For each mode u k (t), VMD constructs the analytic signal by means of Hilbert transform and calculate the unilateral frequency spectrum.Then, the frequency spectrum of the modal function is corrected to the estimated center frequency by Fourier transform.In the following step, the bandwidth of the modal component can be calculated by Gaussian smoothing.Therefore, the variational constraint model that minimizes the sum of the spectral widths of all IMFs can be obtained as where {u k } is the set of all modes; {ω k } is the set of corresponding center frequencies; K is the mode number; the constraint is the sum of all modes, which is also equal to the original signal.
In order to solve the optimal solution of the abovementioned variational constrained problem, VMD introduces the Lagrangian multiplication operator and the quadratic penalty factor to convert the problem into a nonconstrained variational problem as follows: where α is the penalty parameter; λ is the Lagrangian multiplier; is the vector inner product.Before satisfying the iterative stop condition, VMD solves the abovementioned nonconstrained variational problem iteratively by introducing alternating direction method of multipliers, and K IMFs components decomposed from the original signal finally can be obtained.
The VMD method can decompose the original signal into a limited number of IMFs that contain local characteristic signals of the original signal.Therefore, by permuting and combining the IMF components, the SNR oscillation term required for subsequent spectrum analysis can be obtained.Generally speaking the SNR oscillation term can be obtained by removing IMF components with the same characteristics as the trend term of the original SNR data, and then reconstructing the remaining IMF components.Considering that the sea level height changes greatly, the frequency change of the coherent signal may be more complicated.In this case, using a combination of multiple IMF components can ensure that more accurate inversion results can be obtained through the VMD method.
In order to better illustrate the advantages of the VMD method in the process of sea level height inversion, the satellite observation signal from the Onsala Space Observatory on the first day Fig. 4. Environment map of GTGU station.We only used the data from the zenith-looking antenna (GTGU) [41]. of 2016 is selected for analysis, and the SNR data at elevations of 5°-30°is decomposed by VMD, as shown in Fig. 3.
In Fig. 3, the horizontal axis represents the elevation angle of the satellite, and the vertical axis represents the signal amplitude.From Fig. 3, we can clearly see that the VMD decomposes SNR data into four layers of signals with different frequencies and characteristics.The VMD method is improved based on the EMD method, and can choose the number of decomposition layers independently.Comparing the amplitude and frequency of each IMF component with the original SNR data, it can be found that the IMF 4 component is a low-frequency component and the amplitude is basically the same as the original SNR data.Therefore, the detrended SNR signal can be obtained by reconstructing the remaining IMF components.In Fig. 4, it is clear shown that the VMD method can effectively remove the trend term of the SNR data.In addition, the VMD method is based on each IMF component obtained by frequency decomposition.The SNR oscillation term corresponding to the frequency of the coherent signal can be obtained by reconstructing the IMF component.Compared with the complex SNR sequence trend item obtained by the least squares fitting method, the SNR oscillation item obtained by the VMD method has less influence on the spectrum analysis.
In previous studies, the traditional GNSS-IR models based on the least squares method usually only achieve good results at low-elevation angles.However, as the elevation increases, the accuracy of the algorithm gradually reduction or even fail.The VMD algorithm proposed in this article still maintains good accuracy with high elevation.In order to verify the improvement effect of the VMD algorithm in different elevation intervals, this article conducted some experiments at the OSO of Chalmers University of Technology in Sweden.

IV. MODEL VERIFICATION TEST IN GTGU, SWEDEN
To verify the feasibility of the new model, we conduct the experiments at the GTGU station, and the GNSS-R data were provided by the Onsala Space Observatory [39].Many scholars have conducted experiments at this GNSS station earlier, which indirectly explained the appropriateness of the location selected in this experiment [40].But this article also gives some explanations on this station.The GNSS station is located at the OSO on the west coast of Sweden (57.4°N, 11.9°E).The equipment at the site includes two Leica GRX1200 receivers and two Leica AR 25 antennas, as shown in Fig. 4. The antenna is a part of a two-antenna experimental installation for measuring sea level height, with one zenith and one nadir-looking antenna.In this study, we only used the data from the zenith-looking antenna (GTGU).The GTGU is installed about 4 m above the MSL, and the receiver is set to a sampling frequency of 1 s.In order to maximize the number of satellite orbits, GTGU is installed in the south direction.Fig. 5 shows the GPS L1 first Fresnel zone that appeared in GTGU in the first week of 2016.Each color represents a GPS satellite and the location of the GNSS station is in the white circle in Fig. 5 [41].
The tide gauge is located about 300 m from the GNSS station and records sea level height every minute.As this tide gauge reports sea level height every minute and the location is sheltered from breaking waves, it can be used as a references to compare different algorithms.According to the sea level height data obtained from the selected tide gauges, it can be found that the sea level height fluctuations from January 2016 to June 2016 are relatively stable.However, due to the proximity of the GTGU station to the Arctic Circle, its tidal changes are irregular, the sea level height curve is not smooth, and there are many burrs.In order to better analyze the experiment, this article obtained the rainfall, wind speed, and temperature of the experimental site in 2016 from the Swedish Meteorological and Hydro logical Institute (SMHI).
In order to verify the feasibility of the VMD method to process SNR data in the sea level height inversion model, this article combines the VMD method and the LSP spectrum analysis to invert the sea level height in the first 6 months of 2016.The satellite elevation range is 5°-30°.The VMD can choose a fixed number of decomposition layers to decompose the original SNR signal.The number of layers used in this article is all four layers unless otherwise specified.As shown in Fig. 3, the VMD decomposes the original SNR data into four layers of IMF components.By observing and comparing the consistency between the IMF components and the SNR sequence, we can find that IMF 4 is a local feature component that characterizes the trend of SNR data.Therefore, accord to reconstruct the remaining three-layer IMF components, the oscillation term of the SNR data can be obtained, and then the LSP spectrum analysis can be performed.
Simultaneously, we also used the least square method, EMD method, wavelet decomposition, and SSA for SNR data processing.According to experience, the EMD method selects the maximum number of decomposition layers as six layers, and the wavelet decomposition method chooses the eight-layer db5 wavelet function.For SSA, we chose to use reconstructed component 1 of the original times series as the trend item signal.In addition, we also conduct a comparative analysis with the results of the first IAG intercomparison campaign.From the data uploaded by the four participating teams, we selected sea level height data consistent with our experimental time period [42].Table I presents the statistics of sea level height inversion result obtained by different methods from DOY 1, 2016 to DOY 181, 2016 (DOY 85-88 is missing).The satellite elevation range and retrieval method selected by each solution are different, resulting in different performance of each solution.For group a and group b, they established a model that relied on sea level to determine model parameters to improve the accuracy of the results and used a moving window method to ensure that the retrieval was updated every minute.The VMD method only needs to process the SNR data and then obtain the measured reflector height through the LSP to achieve the purpose of quality control.From the statistical results, the SNR data processing method based on the VMD method has achieved good accuracy results.The RMSE and Fig. 6.Top is the time series of retrieval results for DOY 1-181, 2016; bottom is the sea level bias with respect to the tide gauge.correlation coefficient are 4.86 cm and 0.97, respectively.The VMD method does not have a specific target expression when separating signals, so it can flexibly fit complex SNR data.Fig. 6 shows the sea level height retrieval results and sea level biases obtained by different methods.
On the basis of observing the sea level inversion results of the long-term series, we continue to study the monthly sea level height inversion results of the VMD method and other commonly used SNR data processing methods.The statistical results are recorded in Table II.It can be clearly seen that the VMD method has advantages in accuracy, correlation coefficient, and GNSS data utilization.The inversion result of the sea level height of the VMD method in January is shown separately, as shown in Fig. 7. Tables III and IV present the comparison of the inversion accuracy of each method at a shorter period.

V. DISCUSSION
Experiments conducted at the GTGU station show that compared with the traditional model, the new GNSS-IR model combined with the VMD method has a considerable improvement in the accuracy and the utilization of GNSS data of sea level height  inversion.The superiority of the VMD method also is reflected in the processing of SNR data with high-elevation angles.
In VMD, the value of the decomposition layer k is a custom variable, the decomposition result will get different results as the value of k changes when the value of k is selected.The value directly affects the accuracy of the result.If the value of k is too large or too small, it will affect the result.The number of decomposed layers K and is selected empirically based on experience and observation which impacts its adaptability.In the experiment carried out in this article, the number of decomposition layers selected by the VMD method is four layers.According to the principle of the VMD method to process data, different decomposition levels will affect the SNR oscillation term obtained after reconstruction and thus affect the inversion result [43].In order to study the effect of the number of decomposition layers of the VMD on the inversion results,    the best decomposition layer is beneficial to the improvement of inversion accuracy.We use the VMD method to decompose a piece of SNR data, and calculated the correlation coefficient and kurtosis of each IMF component corresponding to different K values.The results are recorded in Table VI and Fig. 9, respectively.
The SNR of the received signal has an important influence on the measurement accuracy.Earlier, Löfgren [45], Larson [22], and others talked about the influence of wind on the inversion results.When the sea breeze produces obvious roughness due to changes in the wind field, there are a large number of diffuse reflection points around the points that satisfy the geometric relationship of specular reflection.The existence of diffuse reflection points will affect the power curve of the reflected signal, thereby affecting the SNR.This article uses the SMHI to query the climate data such as wind speed, rainfall, and air temperature near the OSO, and analyze whether it has an impact on the inversion results.In Fig. 10, we found that the error distribution of the inversion results has no correlation with the wind speed, rainfall, and air temperature in the monitored area, which is consistent with the research conclusions of previous scholars.

VI. CONCLUSION
This article proposes a GNSS-IR model combining VMD for sea level height altimetry.Compared with the traditional GNSS-IR model, the new model improves the accuracy and stability of the inversion results.In this article, we used the VMD method to remove the trend items of SNR data, and verified its feasibility by comparing with other SNR data processing methods.Experimental results show that the new GNSS-IR model has improved both its accuracy and stability under high satellite elevation angle ranges.In addition, the new GNSS-IR model is also suitable for sea areas with extremely complicated sea level height changes.

Fig. 2 .
Fig. 2. Comparison of oscillation term obtained by least squares and VMD (The black line is the result obtained by the least square method and the thick red line is the result obtained by the VMD.).

Fig. 3 .
Fig. 3. Composition of the IMF component of the SNR sequence, high frequency to low frequency from top to bottom.

Fig. 5 .
Fig. 5. Satellite remote sensing images of GTGU station.We can see the GPS L1 first Fresnel zone of GTGU station [41].

Fig. 7 .
Fig. 7. Comparison of the sea level height change obtained by the tide gauge station and GNSS-IR technology based on VMD method in January.

Fig. 8 .Fig. 9 .
Fig. 8. Process of the VMD k value determination method.TABLE VI CORRELATION COEFFICIENTS OF IMF COMPONENTS UNDER DIFFERENT K VALUES

Fig. 10 .
Fig. 10.Comparison of the climate and the sea level bias in GTGU.(a) Comparison of the rainfall and the sea level bias.(b) Comparison of the wind speed and the sea level bias.(c) Comparison of the temperature and the sea level bias.

TABLE I STATISTICS
OF SEA LEVEL HEIGHT INVERSION RESULT OBTAINED BY DIFFERENT METHODS FROM DOY 1, 2016 TO DOY 181, 2016 (DOY 85-88 IS MISSING)

TABLE II COMPARISON
OF THE ACCURACY OF THE INVERSION RESULT OF FIVE METHODS

TABLE III COMPARISON
OF THE ACCURACY OF THE INVERSION RESULT OF FIVE METHODS DURING THE FIRST SEVEN DAYS OF 2016 WITH THE ELEVATION CHANGED FROM 5°-15°(LEFT) TO 5°-35°(RIGHT), FROM TOP TO BOTTOM ARE VMD, SSA, LEAST SQUARES, EMD, AND WAVELET DECOMPOSITION

TABLE IV COMPARISON
OF THE ACCURACY OF THE INVERSION RESULT OF FIVE METHODS WITH THE ELEVATION CHANGED FROM 5°-15°(LEFT) TO 5°-35°( RIGHT) ON DOY 61-67, 2016

TABLE V COMPARISON
[44]MSE OF THE INVERSION RESULT OF DIFFERENT IMF COMPONENT LAYERS, FROM 3 TO 7we analyzed the effect of processing the SNR data of January 2016 with the IMF components of the three to seven layers.In theory, the more decomposition layers, the more accurate the local feature extraction of SNR data.However, due to the influence of the surrounding environment of the receiver, the SNR data is quite complicated.From the results in TableV, it can be seen that as the number of decomposition layers increases, the RMSE of the inversion result first becomes larger and then smaller, but the correlation coefficient still maintains 0.98 and the number of inversion points increases from 1736 to 2269.However, we found that the best decomposition layers would be different when we selected the different experimental intervals.Therefore, different decomposition layers can be selected in different experimental environments.Wu et al.[44]pointed out a novel method based on kurtosis to select K, which may solve this problem.The specific process is shown in Fig.8.A positive integer with k value of 2 − n is selected, starting from k being 2, the kurtosis of the component with the largest correlation coefficient of the original signal under each k value is calculated.If the kurtosis increases monotonously, calculate the k kurtosis as n+1, and repeat the above steps.The maximum kurtosis is used as the optimization criterion and k is the best when the kurtosis is the maximum.The advantage of this model is that it only needs to provide GNSS SNR data without resorting to TG data to calculate the optimal number of decomposition layers.Choosing