Global Positioning System interferometric reflectometry for accurate tide gauge measurement: Insights from South Beach, Oregon, United States

Global Positioning System (GPS) stations located along coastal areas have the ability to measure tide gauge (TG) records by reflected signal reception from the sea water surface. In this study we used the GPS signal-to-noise ratio (SNR) data from the SEPT station (44.63 0N, 124.05 0W) located at South Beach, Oregon, USA, to estimate the TG records using only a few measurements. First, we derived the TG record from a GPS station (GPS-TG) and used traditional TG data from the National Water Level Observation Network (NWLON) sentinel station (Station ID: 9435380) located in Oregon for validation purposes because it was closest to the SEPT station. Our results show that the GPS-TG and NWLON-TG correlate well with the correlation coefficient (CC) of 0.942 and the root mean square (RMS) of their residuals was about 12.90 cm. The corresponding TG prediction by autoregressive moving average (ARMA-TG) and singular spectrum analysis (SSA-TG) models are evaluated for their effectiveness over the station. The comparative analysis demonstrates that the GPS-TG has improved correlation with ARMA-TG (CC of ~0.981 CC, RMS of ~4.80 cm), and SSA-TG (CC of ~0.998 CC, RMS of ~ 0.88 cm) compared to the NWLON-TG (CC of ~0.942 CC, RMS of ~12.90 cm) values. We believe the outcomes from this study contribute to a better understanding of the numerical modeling of TG records as well as other measurements based on reflectometry techniques.


Introduction
Hydrometeorological hazards are relatively frequent and constitute a major threat to life and property within the palette of natural hazards [1].The Global Positioning System-Interferometric Reflectometry (GPS-IR) technique is used to monitor sea-level measurements by using a geodetic quality GPS receiver for identifying the nearby environment, including the benefits of climate dependent and independent sea level information [2].Additionally, the GPS-IR has been established and validated for measuring snow depth [3], tsunamis [4], surface soil moisture [5], flooding inundation [6], inland water [7], storm signature tide records [8], and vegetation water content [9].Several studies have been approved that a single geodetic receiver is able to provide comparable accuracy to traditional tide gauges (TG) [e.g.2,10,11].The GPS antenna is designed in such a way that it receives the direct signal that becomes Right-Handed Circularly Polarized (RHCP), similar to the transmitted signal and the reflected signal becomes Left-Handed Circularly Polarized (LHCP) (Fig. 1).The planar surface distance (H0) is measured from the RHCP antenna phase center because both reflected and direct signals are received [2].The GPS signal strength, including code and phase observables is recorded as signal-to-noise ratio (SNR) from the receiver independent exchange format (RINEX) observation data.The SNR measurement unit is decibel-Hertz (dB-Hz) or decibel (dB) and indicates the GPS signal power ratio to the measured noise.The power of signals transmitted by the GPS satellite, tracking algorithm in the receiver, and antenna gain pattern are the main SNR factors [12].The mathematical description of the SNR equation is described in section 2.
The quantity and location of TG do not yield sufficiently accurate observations and require frequent maintenance [13].Sometimes a storm surge causes an abnormal water surface rise as strong winds push water onshore.These storm surges are measured by the water height above the normal predicted astronomical tides.Löfgren et al. [14] proved that sea surface roughness increased because of increased wind speed causing decreased reflected signal power and SNR oscillation amplitude; while the SNR data L1 signal from a single geodetic receiver remains unaffected.However, the weaker reflected signal power still interfered with the direct signal and created SNR oscillations [15].The damaging impacts of these storm surge related TGs can be reduced by proper modeling and forecast.The comparison analysis and validation performance of TG data have been done by several authors in the past decade [e.g., [16][17][18].Some hydrological models were established by performing simulation modeling with and without forecast [e.g.13,19].The problems related to physical based models leads to numerous investigation methods, in which singular spectrum analysis (SSA) and autoregressive moving average (ARMA) methods seem to offer interesting alternative paradigms.We collected SNR data from the SEPT station (44.63 ⁰N, 124.05 ⁰W) located at South Beach, Oregon, USA for the current study and forecasted that with the SSA and ARMA models.The station consists of standard geodetic quality choke-ring zenith pointing antenna with dual-frequency carrier-phase GNSS ASTERX-U receivers operating at 1 sample per second.Peng et al [15] found that the reflected signal power of GLONASS-SNR data is generally weaker compared to GPS satellite SNR data and the generated sea surface height by GLONASS satellites contains more outliers such as noise error and non-physical reflector heights.Hence, in the current study we used only GPS satellite data with the L-Band dual-frequency of GPS signals having wavelengths of 19.05 cm (L1 signal) and 24.45 cm (L2 signal) and estimated the GPS-TG records from the SEPT station.The National Water Level Observation Network (NWLON) sentinel station (Station ID: 9435380) located in Oregon is used for validation purposes because it the closest station to SEPT.The detailed study described here is separated into the following sections.In Section 2, a brief methodology of TG estimation from GPS signals and its numerical model is presented.The study based on GPS-TG forecasting is outlined in Section 3 and the study conclusions are provided in Section 4. .

Estimation of TG from GPS data
The study of GPS-IR is based on the multipath, or interference of direct and reflected GPS signals [20].The reflected signal oscillation waves are a direct influence on SNR and can be given by the most basic signal wave equation form as follows: where A stands for the function deviation from zero known as amplitude, ω is angular frequency in radians per second, t is time and ϕ is phase in radians.If f is the ordinary frequency equivalent to number of oscillations that occur each second of time, then the Equation ( 1) is rewritten as: where the velocity and signal wavelength are v and λ, respectively, we can transfer the equation into the form, As we can see from the Fig. 1, the travel path is transferred in terms of the priori reflector height (H0); or the vertical distance between the horizontal reflecting surface and the GPS antenna phase center, Now, comparing Equation (4) with Equation (1) we get, Equation ( 5) represents the reflected signal phase delay in terms of priori reflector height where the left side varies with time.Axelrad et al. [21] suggested that the time varying multipath effect can be eliminated by direct consideration of the rate of phase delay change with respect to elevation angle sine.Hence differentiating Equation (5) with respect to sin θ, ( ) ( ) According to Axelrad et al. [21] the 2H0/λ is the spatial frequency (fM) of the multipath measurement cycles per full satellite arc from 0 to 90°.Hence, we can write from Equation (6), To remove this influence and only record the interference pattern, the SNR interferogram is detrended with a low-order polynomial (here, second order) and converted to a linear scale (volt/volt) estimation for each satellite track on each day, where SNRdb-Hz is the raw data, and SNRV/V is the data transformed to a linear, volt/volt scale.The SNRV/V obtained value from Equation ( 8) is used to estimate the fM by using Lomb-Scargle periodogram [LSP; 22], which is known as the dominant frequency (fm,t), and simple planar reflector distance will convert to an effective reflector height (Heff) [see e.g., 2, 5, 15],

Forecasting of GPS-TG Records
The applicability of the ARMA model algorithm can be reviewed by taking a time-series of TG records as Z (t) at time (t), then the general mathematical form of the ARMA (p, q) model is given by: where ϕi is auto-regressive (AR) and θj is the moving average (MA) coefficient.The p and q are AR and MA orders, respectively.The e(t) stands for white noise, a normal distribution with zero mean and σ standard deviation.
The mathematical form of the AR operator is given by, and similarly, the mathematical form of the MA operator is given by, These are the ARMA modeled equation mathematical forms for TG record modeling.The working SSA method procedure is divided into four main steps: embedding, singular value decomposition, eigentriple grouping, and diagonal averaging.In TG cases, we consider the given time-series of GPS-TG records as Z (z1, z2, z3, ..., zN) of length N, where the length of N must be greater than 2 and have at least one non-zero term.The trajectory matrix (Z) for the SSA method following a Hankel matrix can be given like this [23], where L is an integer and can be represented as 1<L<N/2.The complete SSA method description can easily be understood by reviewing Ansari et al., [24]; here our plan is to show only the reconstructed time-series after the application of SSA methodology.At this step the original matrix is transferred into a new time-series.Let us suppose that each matrix of Ai with a(i,j), elements is an (L×L) matrix such that, (1,1) .... .... .... : : : : ....
Now matrix Ai is transferred into new time-series a1, a2, ... ak ... aN, ( The diagonal averaging practices to every matrix of Ai will produce a reconstructed series in terms of N k a .The original series T are decomposed into the summation of a reconstructed series such that, This decomposition is the main SSA algorithm outcome.

Background
The interference between the direct and reflected GPS signals is recorded in SNR interferograms.Direct signal SNR plots with respect to elevation angle as a typical example of PRN 09 and PRN 17 on 01 January 2018 are shown in Fig. 2a.We can see from the figure when the satellite rises overhead the SNR value is increasing while the oscillation pattern effect decreasing.This happens because the oscillating interference pattern frequency is a primary function of the path length difference.As the satellites achieved a high angle, the relative proportion of the directly received signal increases because of the antenna gain pattern, which increases the SNR and decreases the oscillation pattern effect [25].The GPS antenna receives both the direct and reflected signals from the neighboring water surface when any satellite travels along its arc.These two separate signals (direct and reflected) interfere with each other and create a characteristic oscillation pattern covering a long-time SNR data trend [15].We used a second order polynomial to remove the trend, retaining only the interference pattern.The polynomial represents the direct signal and has no interest in measuring the sea level variation [2].The remaining SNR signal with respect to the sine function of elevation angle is shown in Fig. 2b.The plotted oscillations are caused by ocean reflections and should be a constant frequency [2].These oscillation frequencies belong to the reflector height H and can be estimated by translating the SNR data from a logarithmic scale unit (dB-Hz unit) to a linear scale unit (volt/volt unit) as shown in Fig. 2c.The remainder of SNR with respect to the elevation angle sine do not have equally spaced points, hence we used LSP [22] to create a power spectrum for detecting the periodic component.The time periods between randomization of the observation and the actual observations lead to a randomization of Fourier peak height and location.The dominant peak corresponds to the multipath frequency in the SNR data [26] for the selected satellites (PRN 09 and PRN 17) as shown in Figure 3a.The dominant frequency is converted to an effective reflector height by using Equation ( 2) and as shown in Fig. 3b.

ARMA-TG Method
We selected three days (04 April, 10 August and 10 October of 2017) of SNR data from the SEPT site (the only data available to us) and validated its performance compared to the NWLON sentinel station.For each day three epochs of TG in different colors are shown in Fig. 4. We estimated the correlation coefficient (CC) by taking the points above zero because the model does not accept negative values.The GPS-TG data with NWLON data shows a 0.942 CC and the root mean square (RMS) of their residuals is 12.90 cm (Table 1).We predicted this observed GPS-TG data by using the ARMA (p, q) method, examining their variation on the selected days by taking the order p=1 and q=1 (note that p and q are the well-known symbols used in ARMA methods) and their CCs are listed in Table 1; showing that the ARMA model results highly correlate to the measured TG with a CC of approximately 0.981.Therefore, it can be concluded that the model's results perfectly represent the observed TG.The model's TG values accuracy is validated by the residuals between the ARMA-TG and observed GPS-TG values and assessed by RMS (Table 1).The RMS values between the observed and ARMA-TG values is 4.80 cm; comparatively better than those of the NWLON-TG model.According to these results, we can state that the proposed ARMA forecasting method operates perfectly as evidenced by its better performance and reliability.

SSA-TG Method
We did the SSA analysis by taking a window length of eight and checked the reconstructed components (RC) from the beginning and compared them with the original time-series.The CC and their respective plots are shown in Fig. 5.The first RC shows that the SSA modeled TG (SSA-TG) data and observed GPS-TG have a CC equal to 0.5102, meaning they have a 51.02% relationship to each other.We added all the RCs (RC1+RC2+…+RC8) and CC continues increasing until finally, at the 8th RC (RC8) CC reached approximately 0.998 meaning the original GPS-TG and predicted SSA TG values have almost 1 to 1 relationship.Hence, supporting that the SSA method is a successful tool in forecasting the TG series and to the best of our knowledge no better model has been tested for the TG variation prediction over the region.The RCs are used to discriminate the TG variation by determining their significant contribution towards overall variance.In brief, the TG variability is highlighted by the eigenvector whereas its temporal characteristics are presented by the RC and are collectively referred to as a mode.The primary mode bearing the highest eigenvalues contains a major portion of the variance followed by subsequent modes with relatively lesser eigenvalues.A similar explanation for variability modes was demonstrated by Vu et al [8].They found that the contribution to the variance of the first mode is 47.2% while the second mode contribution to the variance is 44.2% and concluded that these components correspond to the tidal components of the sea surface height and represent the tides.The contribution to the variance of the third (RC3) and fourth (RC4) mode were much less and accounted for only 2.7% and 1.0%, respectively.Although the sea surface height is affected by storm surge, there are other factors, such as wind, atmospheric pressure, waves, etc. also responsible for surges.Vu et al [8] used the third component, RC3 and correlated their results with atmospheric pressure, significant wave height, and wind related to the 2010 Xynthia Storm and found that they correlated with RC3 at about 70 %, 53% and 53 %, respectively.We compared the RC3 component with wind speed and noticed that in the SEPT case, wind speed shows only a 14% correlation; much lower than the Xynthia Storm, meaning the RC3 and surge are not strongly related and accounts for minimal surge impact on the sea surface height.

Conclusions
We used SNR data from the SEPT GPS station (44.63 ⁰N, 124.05 ⁰W) located at South Beach, Oregon, USA to estimate the TG records from only a few measurements.First, a mathematical description of GPS-IR is presented and prediction methods like ARMA and SSA are outlined.The observed GPS-TG results were then compared with the NWLON sentinel station (Station ID: 9435380) located in Oregon for validation purposes.We found that GPS-TG and NWLON-TG have a 0.942 CC and residual RMS is equivalent to ~12.90 cm.Although the comparison results are satisfactory, we proposed application of some universal techniques for improvement.The techniques work and the CC of GPS-TG versus ARMA-TG and GPS-TG versus SSA-TG improved to 0.981 and 0.998, respectively.The RMS results also show improvement, reaching approximately 4.80 cm and zero cm by ARMA and SSA predicted values, respectively.Hence, in conclusion, we can say these prediction methodologies are superior options for TG record forecasting.Finally, we compared the RC3 component with the wind speed and noticed that in the SEPT case the wind speed shows 14% correlation.This correlation with wind speed is very low, indicating that surge is not strongly TG related and accounts for less sea surface height impact.

Figure 1 .
Figure 1.The GPS interferometric reflectometry (GPS-IR) technique where the land-based receiver receives the GPS signal electromagnetic wave phase.The GPS direct signal reaches the satellite directly while the reflected signal arrives after reflection from the nearby water surface.

Figure 2 .
Figure 2. a.The SNR L1 and L2 frequency measurements for two randomly selected satellites (PRN 09 and PRN 17) recorded by SEPT Station on 01 January 2018 with respect to elevation angle.

Figure 2 .Figure 2 .
Figure 2. b.The detrended SNR L1 and L2 frequency measurements (dB-Hz) for the randomly selected satellites (PRN 09 and PRN 17) recorded by SEPT Station on 01 January 2018 with respect to elevation angle sine.

Figure 3 .Figure 3
Figure 3. a. Lomb Scargle Periodograms (LSP) of the SNR data (from Fig.2c).The LSP peaks determine the reflector frequency used to estimate sea level.

Figure 4 .
Figure 4.The GPS-TG records (dotted points with vertical lines) estimated from GPS satellites (three points per day) and measured NWLON (horizonal line) TG records.

Figure 5 .
Figure 5. Reconstruction of TG time-series plot of SEPT site.