Ship-borne wave gauge using GNSS interferometric reflectometry

ABSTRACT

To obtain wide coverage and frequent wave observations, low-cost methods to observe significant wave heights and wave periods from moving vessels were developed by using Global Navigation Satellite System (GNSS) interferometric reflectometry. Because GNSS signals reflected from the sea surface are always delayed respect to those received directly from the satellite, and this delay depends on the distance between the antenna and the sea surface, the amplitude of the received GNSS signals, in which the direct and reflected signals are mixed, shows interferometric variations in the presence of ocean waves. High sampling rate observations of the GNSS signal amplitude can be used to determine the significant wave period and wave height by referring to a look-up table estimated from simulated received GNSS signals based on a realistic wave spectrum. The methods were applied to actual 20-Hz GNSS observations on a ferryboat and qualitatively good estimations of both significant wave heights and wave periods were obtained. Because the look-up tables can be further modified using in situ wave observations, these estimations can be quantitatively improved in the future.

KEYWORDS:

• GNSS reflectometry
• wave gauge
• interferometer
• ship-borne observations

1. Introduction

Ocean waves are one of the fundamental variables in oceanography and marine engineering, and predicting their characteristics is vital for various activities, including ship routing, fisheries, and marine leisure. Due to the high demand for this information, ocean waves have been observed for a long time by coastal wave gauges or offshore buoys, but the observation locations are significantly limited. Therefore, since the mid-1980s, satellites have been used to observe global wave fields. Satellite altimeters have provided a long-term time series of global significant wave height (SWH) and wave period observations (Challenor and Srokosz 1991; Gommenginger et al. 2003; Ribal and Young 2019; Wang and Ichikawa 2016). In addition, both synthetic aperture radar images from various satellites and the recently commissioned surface waves investigation and monitoring (known as SWIM) sensor on China France oceanography satellite (i.e. CFOSAT) can determine additional wave parameters such as wave propagation directions and wavelengths (Collard, Ardhuin, and Chapron 2005; Hauser et al. 2021). However, the spatial and temporal densities of these satellite observations are too sparse to monitor real-time wave fields, although these observations are still useful for obtaining wave statistics, validating wave models, and driving assimilation wave models.

Ship-borne wave observations would be another solution for monitoring real-time wave fields because data from numerous vessels could be accumulated. Although their distributions are not as uniform as observations from satellites, they can cover areas with especially high data demands, including marine traffic routes and coastal areas where waves are often affected by complex topography or interact with ocean currents. It should also be noted that conventional satellite altimeters are not able to observe SWHs accurately in coastal areas (Ichikawa, Wang, and Tamura 2020). Ship-borne wave observations are, however, mostly based on visual inspections, making them fundamentally qualitative. Quantitative wave observations require the installation of microwave or supersonic wave gauges or microwave radars, which would practically limit the number of vessels available. In addition, because wave gauges are subject to ship wakes around the vessel, while vessels are at sail, locations on a vessel to which wave gauges can be attached are limited. In other words, low-cost and easily attachable sensors are necessary to expand ship-borne quantitative wave observations.

Therefore, a new low-cost method is proposed in this study that measures both significant wave heights and wave periods from moving vessels using Global Navigation Satellite System (GNSS) reflectometry. Reflected GNSS signals are generally avoided because they can be a source of multipath error in GNSS positioning, but GNSS reflectometry (GNSS-R) dares to use them because they include information about the reflection surface. Various parameters can be determined by GNSS-R, including sea surface wind speeds, sea surface heights, and inland water distributions (Asgarimehr, Wickert, and Reich 2018; Garrison et al. 2002; Ichikawa et al. 2019; Jin and Komjathy 2010; Qin and Li 2022; C. S. Ruf, Gleason, and McKague 2019; C. Ruf et al. 2021; Ruffini et al. 2004). The present study introduces the use of GNSS-R as a ship-borne wave gauge.

This paper is organized as follows. Section 2 explains the methods in the present study; Section 2.1 first describes the concept of GNSS-R, and Section 2.2 shows simulation results for simple single-frequency waves. Then, more practical multi-frequency waves are examined in Section 2.3. In Section 3, preliminary results are described for applying the method to actual GNSS observations on a ferryboat. These results are discussed in Section 4, followed by a summary in Section 5.

2. Methods

2.1. Concept

As shown in Figure 1, signals from a GNSS satellite can be reflected by the sea surface and received by an antenna on a vessel, together with the signals directly received from the GNSS satellite. The signal received by the antenna on the vessel can be written as:

Figure 1. Schematic diagram of GNSS reflectometry.

(1) $\mathit{S}=S_D+S_R,$(1)

where $S_D$ and $S_R$ are the direct and reflected GNSS signals, respectively. Because the direct signals travel the shortest path length, $S_R$ is always delayed with respect to ${\mathrm{S}}_D$ (Ichikawa et al. 2019; Roussel et al. 2014). The excess path length $\mathrm{\Delta }\mathit{L}$ is geometrically expressed as:

(2) $\mathrm{\Delta L}=2\left(\mathit{H}-\mathit{\eta }\right)sin\mathit{\theta }=2\mathit{h}sin\mathit{\theta },$(2)

where $\mathit{H}$ and $\mathit{\eta }$ are the heights of the onboard GNSS antenna and the sea surface at the reflection point, respectively; $\mathit{h}$ is the vertical distance from the antenna and the reflection point, namely, $\mathit{h}=\mathit{H}-\mathit{\eta };$ and $\mathit{\theta }$ is the elevation angle of the given GNSS satellite.

The direct GNSS signal $S_D$ has right-hand circular polarization (RHCP), whereas the reflected signal $S_R$ has both left-hand circular polarization (LHCP) and RHCP components. The ratio of these components depends on the elevation angle $\mathit{\theta }$, and the LHCP component is dominant when $\mathit{\theta }$ is larger than 30° (Löfgren et al. 2011). Meanwhile, for a GNSS satellite with a lower elevation angle, the two RHCP signals $S_D$ and $S_R$ interfere with each other. Using the L-band microwave frequency $\mathrm{\Omega }$ and wavelength $\mathit{\lambda }$ for a given GNSS satellite allows $S_D$ and $S_R$ to be expressed as:

(3) $S_D=\mathit{A}cos\mathrm{\Omega }\mathit{t},$(3)

(4) $S_R=\mathit{A}\mathit{r}cos\left(\mathrm{\Omega }\mathit{t}+2\mathit{\pi }\Delta \mathit{L}/\mathit{\lambda }\right),$(4)

where $\mathit{t}$ is the arbitrary time, $\mathit{A}$ is the amplitude of the received GNSS direct signal, and $\mathit{r}$ is the amplitude ratio of the received GNSS reflected signal. In general, $\mathit{r}$ depends on the mean squared slope (or roughness) of the sea surface and is smaller than unity (Gleason, Gommenginger, and Cromwell 2010). Therefore, the received signal $\mathit{S}$ can be further expressed as:

(5) $\mathit{S}=\sqrt{A^2\left(1+r^2\right)+2A^2\mathit{r}cos\left(2\mathit{\pi }\Delta \mathit{L}/\mathit{\lambda }\right)}cos\left(\mathrm{\Omega }\mathit{t}+\mathrm{\Phi }\right),$(5)

where $tan\mathrm{\Phi }=\mathit{r}sin\left(2\mathit{\pi }\Delta \mathit{L}/\mathit{\lambda }\right)\mathit{\hspace{0.17em}}/\mathit{\hspace{0.17em}}\left[1+\mathrm{r}cos\left(2\mathit{\pi }\Delta \mathit{L}/\mathit{\lambda }\right)\right]$. Equations (2)(2) $\mathrm{\Delta L}=2\left(\mathit{H}-\mathit{\eta }\right)sin\mathit{\theta }=2\mathit{h}sin\mathit{\theta },$(2) and (5(5) $\mathit{S}=\sqrt{A^2\left(1+r^2\right)+2A^2\mathit{r}cos\left(2\mathit{\pi }\Delta \mathit{L}/\mathit{\lambda }\right)}cos\left(\mathrm{\Omega }\mathit{t}+\mathrm{\Phi }\right),$(5) ) can be used to express the power of the received signal $\mathit{S}$, often referred as the signal-to-noise-ratio (SNR), as

(6) $\mathrm{SNR}=A^2\left(1+r^2\right)+2A^2\mathit{r}cos\left[4\mathit{\pi h}sin\mathit{\theta }/\mathit{\lambda }\right].$(6)

All variables $\mathit{A},\mathit{r},\mathit{\theta }$, and $\mathit{h}$ change over time (Ribot et al. 2014). The amplitude $\mathit{A}$ varies depending on the ionospheric and atmospheric conditions, together with the satellite elevation angle $\mathit{\theta }$ and the antenna gain pattern. The reflection ratio $\mathit{r}$ depends on the elevation angle $\mathit{\theta }$ and the sea surface roughness that depends on wind speed. Most remarkable SNR variations, however, are generated by temporal variations of $\mathit{h}sin\mathit{\theta }$, as cyclic sinusoidal changes; every time $\mathit{h}sin\mathit{\theta }$ changes by a length of 0.5 $\mathit{\lambda }$, the SNR cyclically varies. In other words, the number of cycles of periodical SNR variations gives the ratio of $\mathit{h}sin\mathit{\theta }$ to 0.5 $\mathit{\lambda }$, allowing the vertical distance $\mathit{h}$ between the antenna and the sea surface to be estimated from the number of these cycles. This method is therefore called GNSS interferometric reflectometry (GNSS-IR) because it is based on SNR variations caused by interferometric variations between RHCP $S_D$ and $S_R$ signals.

The GNSS-IR method has been used to estimate tidal variations by a stationary GNSS antenna on land (Larson, Löfgren, and Haas 2013; Löfgren et al. 2011). In this case, the target temporal variation of $\mathit{h}$ is the sea surface height $\mathit{\eta }$ due to tidal elevations. For example, for the M₂ tidal constituent with a 2-m range, $\mathit{h}sin\mathit{\theta }$ changes by 0.68 m during a 6.2-hour flood or ebb tidal period if the elevation angle $\mathit{\theta }$ of a GNSS satellite remains around 20°, which is 6.8 times 0.5 $\mathit{\lambda }$ (or 0.1 m), suggesting 6.8 cycles of SNR variation over 6.2 hours. The maximum temporal change rate of $\mathit{h}sin\mathit{\theta }$ would be 4.4$\times$10⁻⁵ m s⁻¹, meaning that it will take at least 37.8 minutes to change 0.1 m to achieve one periodic SNR change. Therefore, if the time series of SNR variations is long enough to cover a few semiannual periods and its sampling intervals are fine enough to resolve a 37.8-min cyclic change, tidal elevations can be recognized by GNSS-IR.

Note that variation of $\mathit{\eta }$ due to ocean waves is regarded as noise for these estimations of tidal elevations so that SNR variations due to ocean waves are removed by temporal averaging because wave periods are much shorter than tidal periods. For example, due to the above-mentioned M₂ constituent, ocean waves with a 2-m height experience the same 6.8 SNR interferometric cycles in a half-wave period, namely within only a few seconds rather than 6.2 hours. In other words, high-rate SNR sampling is required to observe ocean waves by GNSS-IR. In the case of ocean waves with a 2-m height and a 10 s period, as in the previous example, the shortest SNR interferometric cycle would be 0.5 s; thus, the sampling rate should be higher than 10 Hz to obtain five samples per cycle.

If the antenna is on a vessel and not stationary as on land, variations of $\mathit{h}$ are caused not only by the sea surface elevation $\mathit{\eta }$ but also by the antenna height displacement $\mathit{H}$ due to vessel motions such as heaves. However, if the sailing vessel is large enough, the motion of the antenna $\mathit{H}$ would be significantly smaller than the sea surface height change $\mathit{\eta }$caused by ocean waves. In addition, the sea surface height $\mathit{\eta }$ is defined at the reflection point far from the vessel, which means that $\mathit{\eta }$ is not synchronized with the antenna height change $\mathit{H}$. Furthermore, even in a rare case when $\mathit{H}$ perfectly follows the surface height $\mathit{\eta }$, $\mathit{H}$ itself can be measured by GNSS positioning, making it still possible to perform on-board wave observations from vessels. Hereafter, therefore, we regard the height change $\mathit{h}$ as mainly being caused by surface waves $\mathit{\eta }$ alone.

2.2. Simulation results for single-frequency waves

Equation (6)(6) $\mathrm{SNR}=A^2\left(1+r^2\right)+2A^2\mathit{r}cos\left[4\mathit{\pi h}sin\mathit{\theta }/\mathit{\lambda }\right].$(6) can be used to simulate SNR variations for single-frequency waves $\mathit{\eta }$. Figure 2 shows simulated SNR variations observed at a 20-Hz sampling rate, for waves with a wave period T of 5 s and wave heights of 1.0 m (a) and 2.8 m (b). Other parameters in Equation (6)(6) $\mathrm{SNR}=A^2\left(1+r^2\right)+2A^2\mathit{r}cos\left[4\mathit{\pi h}sin\mathit{\theta }/\mathit{\lambda }\right].$(6) are set as $\mathit{A}=$ 1 dBHz, $\mathit{r}=$0.6,$\mathit{\lambda }=$ 0.19 m, H = 15 m, and $\mathit{\theta }=$ 18°. Here, the 20-Hz sampling rate is selected as the maximum sampling rate for a low-cost GNSS receiver. For the convenience with the later discussions, we expandingly use the term SWH even for the wave height of single-frequency waves.

Figure 2. Simulated 20-Hz SNR variations for single-frequency waves $\mathit{h}$ with wave period T of 5 s and SWHs of 1.0 m (a) and 2.8 m (b). Horizontal lines at 2.25 dBhz (blue) and 1.0 dBhz (red) in panel (a) are drawn for the convenience of discussion in the text.

Both panels show repeated modulation patterns of dense and sparse variations. Because the temporal change rates of the surface height $\mathit{\eta }$ are a minimum at the wave troughs and crests, these sparse variations correspond to the samples when $\mathit{\eta }$ is at the wave troughs or crests so that they would appear repeatedly with the half-wave period, namely 2.5 s in this case. Meanwhile, the variations in Figure 2(b) are significantly denser than those in Figure 2(a), although the 2.5-s modulation patterns are the same. This is simply because the maximum vertical speed of $\mathit{\eta }$ is faster when the amplitude of single waves is larger.

To quantitatively measure these discrepancies in SNR high-frequency densities, we introduce a simple index that counts the number of SNR curves crossing a given SNR value in a fixed duration D. For example, the crossing number at 2.25 dBHz (blue horizontal line) in Figure 2(a) is 24 in 10 s, while it becomes 28 at 1.0 dBHz (red horizontal line). Because the crossing numbers obviously depend on the selected SNR values, we first select 100 SNR reference values within the range of SNR variations in the duration D, and then determine the mean crossing number Nc^D.

From comparisons of high-frequency variation densities in Figure 2(a,b), the crossing number Nc^D is larger when SWH is larger over the same wave period. However, the density of the SNR variation modulates over time with a cycle of the half-wave period 0.5 T, making the crossing number Nc^D depend on how many cycles are included in the observation duration D, that is, the ratio between D and the wave period T. Figure 3(a) shows variations of the crossing number in 10 s Nc^10s, against various SWH values. The parameters for the simulations are the same as in Figure 2, but the period T of the single-frequency waves is set as 2, 4, and 6 s for the yellow, green, and red curves, respectively. As expected, Nc^10s for the yellow curve becomes the largest for a given SWH, because more wave cycles are included within 10 s. However, when the duration D is set to the given wave period T, all curves of the crossing number Nc^T show the same increasing trend as the SWH, except for their trailing edges at larger SWH values.

Figure 3. (a) Crossing number of simulated 20-Hz SNR over 10 s, Nc^10s, for single-frequency waves $\mathit{h}$ with different SWHs $H_S$ and wave periods $\mathit{T}$. Yellow, green and red curves represent results for wave periods of 2, 4, and 6 s, respectively. (b) Same as (a), but when counting duration is not fixed at 10 s but limited to the wave period $\mathit{T}$ in seconds, as $\mathit{N}{c}^T=\mathit{N}{c}^{10\mathrm{s}}\ast \mathit{T}/10\mathrm{s}$.

The distinct discrepancy of the trailing edges in Figure 3(b) is due to the limitation of the 20- Hz sampling rate. Because the Nyquist frequency for the 20-Hz sampling rate is 10 Hz, SNR variations are aliased when the change of $\mathit{h}sin\mathit{\theta }$ within 0.1 s exceeds 0.5 $\mathit{\lambda }$, that is, the critical vertical speed for the 20-Hz sampling rate. Because the maximum vertical speed of ocean waves depends on the ratio of the wave height to the wave period, Nc^T is not accurately counted when the SWH becomes larger than the critical value whose maximum vertical velocity exceeds the critical velocity for a given wave period. Note that the trailing edges in Figure 3 start when the SWH becomes around 2, 4, and 6 m for waves with T = 2, 4, and 6 s, respectively. In practice, however, these waves are not present because such steep waves will break under those conditions.

In the previous examples in Figures 2 and 3, the elevation angle $\mathit{\theta }$ was fixed at 18°. In Figure 4, the dependence on the elevation angle is examined. The crossing number Nc^T becomes larger as the elevation angle $\mathit{\theta }$ increases, because variations of $\mathit{h}$ are emphasized by the factor $sin\mathit{\theta }$. When Nc^T is normalized by $sin\mathit{\theta }$as in Figure 4(b), all these curves become the same, showing an increasing trend with the SWH.

Figure 4. (a) Crossing number Nc^T for simulated 20-Hz SNR observations for single-frequency waves $\mathit{h}$ with a 6-s wave period and various SWHs $H_S$ and elevation angles $\mathit{\theta }$. Green, red, and blue curves represent results with elevation angles of 20°, 15° and 10°, respectively. (b) The same as (a), but when the crossing number Nc^T is normalized by $sin\mathit{\theta }$.

2.3. Simulation results for multi-frequency waves

To investigate SNR variations with more realistic surface waves, multi-frequency waves with the ISSC wave spectrum were used rather than single-frequency waves. The ISSC wave spectrum or Bretschneider or modified Pierson-Moskowitz spectrum, is defined as

(7) $\mathit{Ws}\left(\mathit{f}\right)=\frac{5}{16}{H_S}^2{f_p}^4{f}^{-5}\exp \left[-\frac{5}{4}{\left(\frac{f}{f_p}\right)}^{-4}\right],$(7)

where $\mathit{f}$ is the frequency, $f_p$ is the peak frequency and $H_S$ is the SWH (Tucker 1991; Tucker, Challenor, and Carter 1984). We simulated the 20-Hz variations of $\mathit{\eta }$ with the ISSC wave spectrum as follows: for each frequency $f_i$ at intervals of 0.05 Hz, sinusoidal variations whose power is given by $\mathit{Ws}\left(f_i\right)$ with a random phase were generated, and then they were accumulated from 0.05 Hz to 10 Hz, or the Nyquist frequency of the 20-Hz sampling. In the present study, waves with periods larger than 20 s were ignored.

From 120-s records of the 20-Hz SNR observations from the simulated $\mathit{\eta }$ variations, the same relationships between the period-averaged crossing number normalized by the elevation angle $\mathit{\theta }$, $\mathit{N}{c}^{T_p}/sin\mathit{\theta },$ and the SWH was determined as in Figure 5. In the figure, the SWH $H_{\mathit{\hspace{0.17em}}1/3}$ is calculated as the highest third of the waves, rather than the parameter $H_S$ in Equation (7)(7) $\mathit{Ws}\left(\mathit{f}\right)=\frac{5}{16}{H_S}^2{f_p}^4{f}^{-5}\exp \left[-\frac{5}{4}{\left(\frac{f}{f_p}\right)}^{-4}\right],$(7) . The wave period $T_p$ is defined using the spectrum peak frequency $f_p$ as $T_p=2\mathit{\pi }/f_p$.

Figure 5. Period-averaged crossing number normalized by elevation angle $\mathit{\theta }$, $\mathit{N}{c}^{T_p}/sin\mathit{\theta }$, plotted for various SWHs and wave periods of multi-frequency waves with ISSC wave spectrum. The elevation angle $\mathit{\theta }$ is fixed at 18° as in Figure 3(b).

Figure 5 shows qualitative tendencies similar to those for the single-frequency waves shown in Figure 3(b). Namely, a common increasing trend with the SWH is found for all periods, although the trend has a slight dependency on the wave period $T_p$ in Figure 5. Meanwhile, the values of $\mathit{N}{c}^{T_p}/sin\mathit{\theta }$ for the multi-frequency waves are quantitatively different from those for the single-frequency waves. For example, $\mathit{N}{c}^T/sin\mathit{\theta }$ for 6-s period and 3-m single-frequency SWH (Figure 4(b)) is approximately 120, but that for multi-frequency waves (Figure 5) is around 90. This would suggest that the relationship between the SWH and the GNSS-IR crossing numbers qualitatively holds even with practical multi-frequency waves, but it quantitatively depends on the wave spectrum of the field to be applied. Note that the $\mathit{N}{c}^T/sin\mathit{\theta }$ values for an 8-s period (blue circles) are closer to those in Figure 4(b). Because the spectrum peak in Equation (7)(7) $\mathit{Ws}\left(\mathit{f}\right)=\frac{5}{16}{H_S}^2{f_p}^4{f}^{-5}\exp \left[-\frac{5}{4}{\left(\frac{f}{f_p}\right)}^{-4}\right],$(7) is narrower when the peak frequency $f_p$ decreases, the ISSC spectrum becomes closer to a line spectrum of single-frequency waves when $T_p$ becomes larger, which confirms the dependency of the present GNSS-IR method on the wave spectrum shapes.

3. Preliminary results with real data

The relationship investigation described in the previous section was applied to actual SNR data. We installed a helix GNSS antenna (TOP107, TOPGNSS) on the international ferryboat New Camellia (Camellia Line Co. Ltd.; 19961 gross tons, 170-m length, and 24-m width), which crosses the Tsushima Strait between Hakata, Japan, and Pusan, Korea (Figure 6) in 6 or 10 h. At a 20-Hz sampling rate, 120-s SNR variations are recorded every hour by a GNSS receiver (F9P, U-blox). Because the duration of high-rate burst samplings should be short to keep data sizes small, we set it as 120 s, considering absence of significant long-period swells in the semi-enclosed Tsushima Strait. Meanwhile, the interval of the burst sampling was set for one hour, accounting for a large spatial scale of the wave field. In practice, the duration and intervals of the burst samplings will be selected depending on characteristics of wave fields and ship speeds of each application case. Note that a few GNSS satellites are always found at low elevation angles, if several satellite systems are combined, including GPS, GLONASS, and BeiDou.

Figure 6. (a) GNSS antenna installed at starboard-side edge of New Camellia compass deck (yellow circle). (b) GNSS antenna (green circle) installed on edge of reflector metal plate overhanging seaward side of edge, which intentionally receives GNSS signals reflected from sea surface.

For calibration, we selected three cases: when the sea state was calm, moderate, and rough (Table 1). As for the true sea conditions, wave periods and SWH values were obtained from the Japan Weather Association (JWA) Precise Ocean data Laboratory and Intelligence Service (POLARIS) hindcast data (Sato and Matsuura 2019), whose correlation coefficients with respect to in situ buoy observations are known to be as high as 0.94 for SWHs and 0.83 for wave periods. For each case, two GNSS satellites with lower elevation angles were chosen; these are listed as sub-cases in Table 1.

Table 1. Three cases for calibration of actual SNR observations, with the date in 2022, time in UTC, and locations. The wave periods and SWHs were obtained from JWA POLARIS hindcast data. For each case, two GNSS satellites used in the analysis are listed as sub-cases. For each sub-case, the satellite pseudo random noise code (PRN), elevation, and azimuth angles are shown, where the letters “G,” “R” and “C” stand for the GPS, GLONASS, and BeiDou satellites, respectively. The azimuth angle is measured clockwise from the ship heading (north-northwestward).

Case	Date	Time	Latitude [°N]	Longitude [°E]	SWH [m]	Period [s]	Sub-case	PRN	Elevation angle [°]	Azimuth angle [°]
1	5/18	14:10–14:12	34.196	129.926	0.17	3.6	1G	G08	12.9	103
							1R	R18	8.2	64
2	6/14	16:10–16:12	34.702	129.574	1.51	5.3	2G	G04	12.7	68
							2C	C33	7.6	44
3	6/24	14:10–14:12	34.138	129.951	2.56	6.4	3G	G03	13.8	44
							3R	R24	12.0	104

Figure 7 shows an example of actual 20-Hz observations of SNR variations for sub-case 1 G (Table 1). High-frequency SNR variations are recognized in the actual observations, as seen in simulations (Figure 2). However, Figure 7 also includes long-term variations. Because cycles of interferometric SNR variations are much shorter than the wave period, as long as the wave heights are larger than 0.5 $\mathit{\lambda }/sin\mathit{\theta }$, such long-term SNR variations are not explained by interferometric changes caused by ocean waves. Rather than interferometric phase changes, therefore, they are caused by temporal variations of the received amplitude $\mathit{A}$ and the amplitude ratio $\mathit{r}$ in Equation (6)(6) $\mathrm{SNR}=A^2\left(1+r^2\right)+2A^2\mathit{r}cos\left[4\mathit{\pi h}sin\mathit{\theta }/\mathit{\lambda }\right].$(6) due to changes in the antenna power gain. When the antenna tilts due to rolling or pitching ship motions, the angles from the antenna zenith toward the GNSS satellite and the reflection point (Figure 1) also change, which results in variations of the antenna power gain, leading to fluctuations of $\mathit{A}$ and $\mathit{r}$. Note that the pitching or rolling ship motions also generate changes in the antenna height $\mathit{H}$, but such changes do not affect the received amplitude $\mathit{A}$.

Figure 7. Example of actual 20-Hz observations of SNR variations. Data are for sub-case 1G (Table 1). The horizontal axis shows the elapsed time in seconds from 14:10 on 18 May 2022. The orange curve represents variations longer than 14 s, which could be induced by ship motion.

In the present analysis, therefore, variations longer than 14 s (the orange curve in Figure 7) were estimated by a second-order polynomial fit by the Savitzky-Golay filter (Savitzky and Golay 1964), and then removed from the original SNR data to focus on the high-frequency variations. We selected 14 s as the period of waves with an approximately 300-m wavelength, which is nearly twice the 170-m ship length.

The first step is to estimate the significant wave period $T_g$ from the SNR data (Figure 8). Although temporal modulations of high-frequency SNR variations can be seen in Figure 7, they are not as significant as the modulations for the cases of the single-frequency waves in Figure 2. Therefore, we introduced a new quantitative index to estimate the significant wave period. For each 20-Hz SNR observation, the short-term crossing number in 2 s, $\mathit{N}{c}^{2\mathit{s}},$ is determined (Figure 9), whose time series would efficiently indicate the modulations of the high-frequency SNR variations. We first set common 100 SNR reference values in the rage of SNR variations in the 120-s duration, resulting that each 2-s SNR curve crosses only a part of these SNR reference values due to a smaller variation range in a 2-s duration. Therefore, $\mathit{N}{c}^{2\mathit{s}}$ is determined as the discrete median value of non-zero crossing numbers. Then, the lag for the first peak of the temporal-lagged autocorrelation of its time series is determined as the significant wave period $T_g$. For example, $\mathit{N}{c}^{2\mathit{s}}$ time series in Figure 9 has the first autocorrelation peak at a 4.4-s lag, so $T_g$ for sub-case 1 G is estimated as 4.4 s.

Figure 8. Schematic flow chart of present analysis.

Figure 9. Temporal variations of crossing numbers within 2 s, $\mathit{N}{c}^{2\mathit{s}}$, for SNR records shown in Figure 7.

Long-term crossing number $\mathit{N}{c}^{120\mathit{s}}$ in the 120-s duration is separately determined, and it is normalized by the estimated wave period $T_g$ and the elevation angle $\mathit{\theta }$, as $\mathit{N}{c}^{T_g}/sin\mathit{\theta }=\mathit{N}{c}^{120\mathrm{s}}\ast T_g/120\mathrm{s}/sin\mathit{\theta }$ (Figure 8). Then, the SWH value corresponding to the crossing number $\mathit{N}{c}^{T_g}/sin\mathit{\theta }$ and the wave period $T_g$ is determined by referring to Figure 5. These results are listed in Table 2 and shown in Figure 10.

Figure 10. Results of comparisons for three cases. Red stars represent POLARIS hindcast data, and blue dots are the averages of the estimates of two sub-cases in each case (Table 1). Green numbers indicate cases, and blue letters show sub-cases.

Table 2. Estimated SWH and wave period $T_g$ for each sub-case. Mean SWH and wave period are also listed for each case. The mean period is determined from the mean of the two frequencies. Due to differences in L1-band microwave wavelength $\mathit{\lambda }$ for satellites (0.1903 m for GPS; 0.1868 m for GLONASS; 0.1920 m for BeiDou), $\mathit{N}{c}^{T_g}$ values for GLONASS and BeiDou satellites have been modified by factors (1.02 and 0.99, respectively).

Case	Sub-case	Wave Period $T_g$ [s]		$\mathit{N}{c}^{T_g}/sin\mathit{\theta }$	SWH [m]
			mean			mean
1	1G	4.4	4.1	12.18	0.3	0.4
	1R	4.0		18.15	0.4
2	2G	4.6	5.5	60.72	1.6	1.4
	2C	6.8		49.51	1.2
3	3G	8.5	7.6	97.43	2.3	2.8
	3R	6.9		137.82	3.3

Because each case has two sub-cases, the mean of the two sub-cases was determined. In general, these mean estimates (blue dots in Figure 10) qualitatively agree well with sea state variations indicated by the hindcasted values (stars in Figure 10); namely, both the SWH and wave period become larger when the sea state becomes rougher. Quantitatively, however, both SWH and wave period are overestimated in Cases 1 and 3, although the SWH discrepancies are less than 0.3 m and period discrepancies are less than 20% of the hindcasted periods. Note also that estimations of two sub-cases are close in Case 1, but they become separated in rougher sea state.

4. Discussion

As shown in the flow chart in Figure 8, estimation of the SWH from observed GNSS-IR high-frequency SNR variations is affected by two issues: the estimated wave period $T_g$ and the external look-up table. For the latter look-up table, we used the simulation results calculated from multi-frequency waves with the ISSC wave spectrum (Figure 5) in the present analysis. However, as discussed in Section 2.3, the quantitative relationship between the SWH and the crossing number depends on the wave spectrum. For example, if we use a look-up table based on a single-frequency wave spectrum (Figure 4(b)) instead of the ISSC multi-frequency wave spectrum (Figure 5), all SWH estimations from the given crossing number values can be decreased, especially when the wave period $T_p$ is short. It could be difficult to find the best wave spectrum for a particular field of interest theoretically, but in practice, a look-up table can be established by validation with many in situ observations.

Because the crossing number $\mathit{N}{c}^{T_g}/sin\mathit{\theta }$ is normalized by the wave period $T_g$, uncertainty in the estimation of the wave period $T_g$ also affects the SWH estimation, as shown in Figure 3(a). For example, if the wave period $T_g$ is overestimated by 1.1 times, the crossing number $\mathit{N}{c}^{T_g}$ will be also overestimated by the same factor, and therefore the SWH would be overestimated too. Note that this proportional relationship between $T_g$ and SWH is consistent with overestimated results in Cases 1 and 3 in Figure 10.

In addition, it should be noted that the estimated $T_g$ could be systematically different from the peak period $T_p$ in Figure 5. In other words, if the estimated $T_g$ tends to be larger than the peak period $T_p$, the crossing number $\mathit{N}{c}^{T_g}$ based on $T_g$ is also overestimated when used as replacement for the crossing number $\mathit{N}{c}^{T_p}$ in the look-up table based on Figure 5 that refers to the peak period $T_p$. Also, note that the wave period in the POLARIS dataset is based on the first moment of the calculated wave spectrum, which could be also systematically different from the peak period $T_p$. Because the preliminary comparisons in Table 2 and Figure 10 neglect these differences of the definitions of the wave period, the present discussions would be limited to qualitative results only.

The estimations for the two sub-cases have significant discrepancies, except in Case 1, although the same wave fields were observed. The reflection points are away from the ship by $\mathit{H}/tan\mathit{\theta }$, from 61 m (for sub-case 3 G with a 13.8° elevation angle) to 112 m (for sub-case 2 C with a 7.6° elevation angle), which are far away enough to avoid the ship wake. Among those limited four sub-cases, sub-cases 2 C and 3 G estimate relatively longer significant wave periods $T_g$ and smaller SWHs than the other sub-case pairs, which cannot be explained by the proportional relationship between $T_g$ and SWH discussed above. Note that the azimuth angle in those two sub-cases are the same 44° (Table 1), which could suggest the results of sub-cases depend on the radial directions due to inhomogeneous wave field. For example, the presence of swells propagating in a fixed direction, which are commonly observed in the semi-enclosed Tsushima Strait, will produce inhomogeneous wave field.

In Figure 1, a flat non-tilted sea surface is assumed at the reflection point to establish simple geometric reflection, but this assumption can fail in the presence of large-scale waves (Ichikawa et al. 2019). If the sea surface tilts slightly over a large area, the position of the actual reflection point will shift from the ideal reflection point, and Equation (2)(2) $\mathrm{\Delta L}=2\left(\mathit{H}-\mathit{\eta }\right)sin\mathit{\theta }=2\mathit{h}sin\mathit{\theta },$(2) can no longer be applied to the path length delay $\mathrm{\Delta L}$ for the reflected GNSS signals. Even with such a tilted sea surface, an increase in the SWH would be qualitatively related to an increase of the crossing number Nc^T because the delay $\mathrm{\Delta }\mathit{L}$ and the vertical distance $\mathit{h}$ are still positively correlated. However, their relationship would be quantitatively different from that in Figure 5, so that the estimations may variate depending on the unknown tilting conditions during observations. This suggests that mean estimations of a large numbers of sub-cases with various satellites in different azimuth angles could be more stable to reduce fluctuations in estimations. In the current method, however, the number of sub-cases is limited due to the requirement of low elevation angles for the RHCP components of GNSS reflected signals. But more sub-cases could be available if another GNSS antenna is installed on the portside of the ship.

In the present GNSS-IR method, no interferometric cycles are counted if variations of the path length delay in Figure 1 are shorter than 0.5 $\mathit{\lambda }$, suggesting that the present method may not work well under very calm conditions. Also under calm conditions, slight antenna motion $\mathit{H}$ may not be negligible in the relation $\mathit{h}=\mathit{H}-\mathit{\eta }$ because $\mathit{\eta }$ becomes small too. In addition, hindcasted estimations in POLARIS data might also be too sensitive to noise in the input data in such a low-energy regime. Nevertheless, the estimations in Case 1 agree with the hindcasted values, with discrepancies only 0.23 m and 0.5 s.

5. Conclusions

A method for estimating both wave period and wave height from moving vessels using low-cost GNSS interferometric reflectometry (GNSS-IR) was developed in this study. Reflected GNSS signals are always delayed with respect to the direct signals, and this delay depends on the vertical distance between the reflection point and the GNSS antenna. In the presence of ocean waves, the delay varies over time, making the phase differences between the direct and reflected electromagnetic GNSS signals periodically change every time the delay changes at around 0.1 m. Thus, the power of the mixed direct and reflected GNSS signals (i.e. the SNR) received by the antenna shows high-frequency variations.

In the present study, the density of these high-frequency SNR variations was counted as the mean numbers of occasions when the observed SNR changes across a given SNR value, and they were found to be positively correlated with the SWH. Using a simulated look-up table between the mean crossing number and the SWH allows the SWH to be estimated from observed SNR variations.

In addition, because the density of these high-frequency SNR variations depends on the speed of sea surface height elevation changes, the density decreases corresponding to slower vertical speeds at the wave troughs and crests. Therefore, the modulation of the density enables us to estimate the significant wave period.

The methods were examined with preliminary observations using equipment installed on a ferryboat crossing the Tsushima Strait between Japan and Korea. Calm (when the SWH is less than 0.2 m), moderate (SWH is about 1.5 m), and rough (SWH exceeds 2.5 m) sea state cases were used to compare the hindcasted wave periods and SWHs with those estimated from the present GNSS-IR method. Qualitatively, the estimated results well distinguished sea state changes. Quantitatively, however, some of the estimated values were slightly overestimated compared to the hindcasted values. The present method uses a simulated look-up table that depends on the shape of the wave spectra used, which could explain these quantitative discrepancies, together with naively ignoring the different definitions of wave periods. Practically, however, look-up tables can be established by calibration using many in situ field observations, which will increase the precision of the estimations in the future.

Disclosure statement

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

Additional information

Funding

This work was supported by the Japan Society for the Promotion of Science under KAKENHI Grant Numbers [JP21K19848, JP22H01301, and JP23K22572]. The authors thank members of Camellia Line Co. Ltd. for their helpful cooperation.