Airborne Coherent GNSS Reflectometry and Zenith Total Delay Estimation over Coastal Waters

Abstract

High-precision GNSS (global navigation satellite e system) measurements can be used for remote sensing and nowadays play a significant role in atmospheric sounding (station data, radio occultation observations) and sea surface altimetry based on reflectometry. A limiting factor of high-precision reflectometry is the loss of coherent phase information due to sea-state-induced surface roughness. This work studies airborne reflectometry observations recorded over coastal waters to examine the sea-state influence on Doppler distribution and the coherent residual phase retrieval. From coherent observations, the possibility of zenith total delay inversion is also investigated, considering the hydrostatic mapping factor from the Vienna mapping function and an exponential vertical decay factor depending on height receiver changes. The experiment consists of multiple flights performed along the coast between the cities of Calais and Boulogne-sur-Mer, France, in July 2019. Reflected signals acquired in a right-handed circular polarization are processed through a model-aided software receiver and passed through a retracking module to obtain the Doppler and phase-corrected signal. Results from grazing angle observations (elevation < 15°) show a high sensitivity of Doppler spread with respect to sea state with correlations of 0.75 and 0.88 with significant wave height and wind speed, respectively. An empirical Doppler spread threshold of 0.5 Hz is established for coherent reflections supported by the residual phase observations obtained. Phase coherence occurs in 15% of the observations; however, the estimated zenith total delay for the best event corresponds to 2.44 m, which differs from the typical zenith total delay (2.3 m) of 5%.

1. Introduction

Climate change is one of the foremost topics of study within many fields in the scientific community. Different studies have been carried out in recent decades to determine the possible risks and hazards and their impact on humans and their environments. One of the most important parameters for climate change monitoring is the variability and changes in oceanic waters, particularly in coastal areas. There are two main parameters investigated in coastal sea surface variability. First, the sea-level change is primarily caused by global warming, which generates thermal expansion due to increasing seawater temperatures, the melting of sea ice and glaciers, and the loss of the terrestrial ice mass in Greenland and Antarctica [1]. Second, the sea state is the description of waves generated by the wind, including their height, direction, and period [2], where the wind–wave component can represent considerable changes in sea level along coasts [3,4].

Common techniques for sea-level and sea-state monitoring comprise tide gauges and buoys that generate precise in situ observations but sparse coverage and, therefore, limited spatial resolution. For global scale, satellite altimetric missions with temporal resolution between 6 and 12 days allow high-accuracy sea-level data in open water but reduced performance in coastal zones due to specific wind–wave influence on the sea surface and the combination of sea waters and land at shorelines [5]. In 1993, the European Space Agency (ESA) proposed a multistatic radar concept that relies on the use of global navigation satellite system (GNSS) signals with an interferometric approach, combining the direct signal with the signals reflected off the Earth’s surface to retrieve sea surface height [6]. This technique is called nowadays GNSS reflectometry (GNSS-R). Figure 1 depicts a schematic representation of GNSS-R.

Currently, GNSS-R has shown the capability of retrieving not only sea surface height [7,8,9], but also further applications, such as sea ice [10,11,12] and soil moisture detection [13,14,15]. For sea-state retrievals, represented by wind speed (WS) and significant wave height (SWH), different methods have been proposed in the literature. For example, in [16], the WS is retrieved by relating the waveforms (delay mapping) of the reflected signals to mean square slopes, which are comprised in a theoretical model from a bistatic radar equation (BRE). As the Doppler spreading effect impacts significantly the sea-state retrievals, for more accurate measurements, the paper also suggested the use of one of the main observables in GNSS-R, the delay–Doppler map (DDM), which represents the scattered power as a function of the time delay and Doppler frequency shift. Consequently, multiple studies have investigated the relation of the DDM and the sea surface roughness from ground-based experiments [17,18] or using satellite reflectometry data [19,20].

Water-surface-level monitoring is also a topic widely studied in GNSS-R. Mainly, determining the reflecting surface height requires modeling the interaction of the reflected signal with the water surface, modeling the signal propagation when passing through the atmosphere, and comparison between the observations and the modeled or ancillary data. The signal propagation modeling comprises a signal delay as electromagnetic waves pass through the atmosphere due to changes in medium density. The troposphere accounts for approximately 80% of the atmospheric delay, which is referred to as tropospheric delay. For coastal sea-level studies, from ground-based data, the tropospheric delay has been commonly disregarded apparently because of the lack or insignificant tropospheric effect seen in the data due to the small height of the reflecting surface. However, based on the results in [21], the correction of tropospheric delay is strongly advised regardless of the antenna height with respect to the sea surface, as it represents a scale error in the measurements. The correction of tropospheric delays has been applied in GNSS-R altimetry retrievals in different studies. In [22], tropospheric correction was carried out by ray tracing using atmospheric refractivity profiles. A complementary approach was used in [23] for ray tracing the tropospheric refraction over a spherical Earth. In [24], a combination of zenith delay from GNSS positioning and Neill’s mapping function was implemented for tropospheric delay corrections. From airborne data, an improvement of about 0.7 m was found in the sea surface height inversion in [25] after removing the tropospheric delay based on the model proposed in [26].

This study examines the occurrence of high-precision phase delay information for airborne reflectometry based on the Doppler spreading estimations under the condition of various sea states in coastal waters. In addition, the paper investigates the possibility of inversion of the zenith total delay at the reflecting surface level from coherent observations and using the Vienna mapping function. Multiple flights were conducted along Opal Coast between the cities of Calais and Boulogne-sur-Mer, France, in the North Sea. A software receiver based on [27] processes the direct and reflected signals. The reflected signals are retracked, aided by the signal path difference between the direct and reflected links retrieved from a specular reflection geometrical model [28]. The residual phase observable is computed from the $I$ and $Q$ components of the retracked signal, and the power spectral density (PSD) allows us to obtain the residual Doppler shift to estimate the Doppler spread, which is closely related to the randomness of ocean waves [29]. The inversion of the zenith total delay comprises the linear regression of the coherent residual phase on the hydrostatic mapping factor multiplied by an exponential vertical decay part depending on height receiver changes.

The paper is organized as follows: Section 2 describes the experiment conducted in July 2019. It includes the study area description and the used platform and antenna setup. Section 3 contains the data set collected, the methods implemented, and the processing steps. Section 4 presents the result analyzing the sea-state conditions and the Doppler spread limit considered for coherent observations. Finally, Section 5 discusses the findings, limitations, and possible further applications based on the presented methodology.

2. Experiment

An airborne GNSS-R measurement campaign was carried out in July 2019. Four flights following the same trajectory were performed on the 12th, 15th, 17th, and 19th along the Opal Coast in the North Sea between the cities of Calais and Boulogne-sur-Mer, France. The total length trajectory over the ocean was ~95 km with a duration of ~1 h each flight. To ensure a transmitter–surface–receiver geometry that includes most of the reflection measurements over the sea surface and near the coast, GPS signals coming from the westward direction (azimuths between 180° and 360°) were selected for the analysis. Besides, due to the loss of most of the GPS satellites in the tracking step along the west-to-east flight segment or insufficient reflection events over the sea surface, the north-to-south segment was analyzed in this study (highlighted in orange in Figure 2, right). In the north-to-south flight segment, it was possible to retrieve continuous reflection tracks from the chosen satellites G1, G7, G8, G10, G11, G16, G18, and G30.

The sea-state variations in the study area, given by the mean value of wind speed (WS), wind directions, and significant wave height (SWH), are presented in

Table 1. Sea-state parameters during the measurement campaign.

Date	Wind Speed (m/s)	Wind Direction (deg)	SWH (m)
12 July 2019	5.49	117	0.30
15 July 2019	4.29	67	0.58
17 July 2019	2.92	204	0.26
19 July 2019	6.50	240	0.55

. These parameters are obtained from the ECMWF ERA5 model along the flight trajectory. The lowest sea state can be observed on 17 July 2019.

The used platform was a gyrocopter, an ultralight and very stable aircraft. The flight altitude was ~780 m above sea level at a speed of ~90 km/h. The setup consisted of one antenna to acquire the direct and the reflected signal in a right-handed circular polarization (RHCP) similarly as presented in [30]. The antenna was linked to two receivers: a Syntony Echo-L receiver for the record of raw GPS L1 signals used for reflectometry and a Javad Delta GNSS receiver for the record of RINEX messages for precise trajectory solution in postprocessing. The antenna was mounted on the aircraft tilted ~43° with respect to the zenith direction, as it is shown in Figure 3.

3. GNSS-R Data and Methods

3.1. Data and Processing

The receiver Syntony Echo-L is capable of recording GPS L1 signals sampled at a frequency of 16.368 MHz. Four raw datasets were obtained for each day and classified as data level 0. The duration of each set was ~18 min, which corresponded to the north-to-south segment of the trajectory to obtain reflection points on the sea surface from GPS satellites located in the west-side sky. Data level 1 included the complex direct and reflected signals, that is, the in-phase $\left(I\right)$ and quadrature $\left(Q\right)$ components after the tracking and retracking process at a rate of 50 Hz. Data level 2 comprised the power peak, and relative Doppler shift and Doppler spread $\left({\sigma }_f\right)$ estimates of the residual reflected signal every 10 s from the power spectral density. The processing flowchart is shown in Figure 4.

3.1.1. Geometrical Path Difference Model

The specular point positions and the path difference between the direct and reflected signals are derived according to [23]. A geometrical model is implemented to characterize specular reflections considering surface curvature. This model requires the transmitter $\left(T_x\right)$ and receiver $\left(R_x\right)$ positions in an Earth-centered Earth-fixed (ECEF) frame. The $T_x$ position is calculated from the broadcasted ephemeris. The $R_x$position is calculated by postprocessing using the nearby permanent GNSS antennas of the French National Institute of Geographic and Forest Information (IGN) network as reference stations. The Earth’s curvature is modeled assuming a spherical surface. The latter is associated with the osculation sphere tangent to the ellipsoid World Geodetic System 1984 (WGS-84) in a reference specular point. From this point, an iterative process using different concentric sea levels (osculating spheres radii) is set until finding the sphere that best fits with the ellipsoid and satisfies that the incident angle equals the reflected angle (specular reflection). Once the final specular point (SP) is defined, a geoid undulation correction is applied using the EIGEN-6C2 model [31]. The link $T_x-SP-R_x$ is established as the reflected path $L^R$. The direct path is modeled from the $T_x-R_x$ link and is defined as $L^D$.

The path difference is given by the residual between the reflected and direct path, ${\mathsf{\Delta }}_p\left(t\right)=L^R\left(t\right)-L^D\left(t\right)$, with changes over time mainly influenced by the passage of the GNSS satellite changing its elevation and the aircraft trajectory. Tropospheric and ionospheric delay corrections are disregarded for the geometrical path model. The optical path considering tropospheric delay using the International Standard Atmosphere parameters is implemented for the residual phase analysis in Section 4.2. Figure 5 shows the tracks of specular points retrieved from the geometrical difference path model on July 17 for the 8 satellites analyzed.

Considering a smooth reflecting surface, the spatial resolution of GNSS-R measurement can be linked to the first Fresnel zone [32]. The footprint associated with the first Fresnel zone is the active scattering region where most of the energy is reflected. The size of this ellipsoidal footprint is given by its major and minor axes, depending on the GPS L1 wavelength, satellite elevation, and receiver height [33]. In this experiment, the major and minor axes correspond to 950 and 55 m, respectively, for satellite elevation angles of 4°. For elevation angles of 70°, the axes correspond to 14 and 13 m, respectively.

3.1.2. Tracking and Retracking

The software receiver is constructed according to [27] to derive the in-phase and quadrature components of the reflected signal $, I_0^R$ and $Q_0^R$, after the tracking stage (Equations (1)–(8) may be found in [27] with changes in notation). The direct and reflected signal, $S^D$and $S^R,$ respectively, are defined by:

$$S^D=A^{D}CA\left(t-{\tau }^D\left(t\right)\right)\sin \left(2\pi f^{D}t-{\varphi }^D\left(t\right)\right)+{\eta }^D\left(t\right)$$

(1)

$$S^R=A^{R}CA\left(t-{\tau }^D\left(t\right)-{\mathsf{\Delta }}_p\left(t\right)/c\right)\sin \left(2\pi f^{D}t-{\varphi }^D\left(t\right)-\psi \left(t\right)\right)+{\eta }^R\left(t\right)$$

(2)

where the amplitudes of the direct and reflected signal are given by $A^D$ and $A^R$. $CA$ is the code-division multiple access (CDMA) broadcasted by the GPS satellites. The code delay of the direct signal is represented by ${\tau }^D$, and $f^D$ and ${\varphi }^D$ are the frequency and phase delay. ${\eta }^D$ and ${\eta }^R$ are zero-mean Gaussian noises. The path difference between the direct and reflected signal is represented by ${\mathsf{\Delta }}_p$, and $c$ is the speed of light. The difference in the phase between the direct and reflected signal is denoted by $\psi$.

The tracking module processes the signals by means of delay, phase, and frequency locked tracking loops (DLL, PLL, and FLL) with an integration time of 20 ms. Refined estimates of ${\tau }^D$, $f^D$, and ${\varphi }^D$ are provided after the tracking of the direct signal. With the estimates, the local replica of the direct signal is constructed and expressed by the in-phase and quadrature components as follows:

$$p^{I,D}\left(t\right)=CA\left(t-{\tau }^D\left(t\right)\right)\sin \left(2\pi f^{D}t-{\varphi }^D\left(t\right)\right)$$

(3)

$$p^{Q,D}\left(t\right)=CA\left(t-{\tau }^D\left(t\right)\right)\cos \left(2\pi f^{D}t-{\varphi }^D\left(t\right)\right)$$

(4)

After the demodulation of the reflected signal with the local replica of the direct signal, its components, $I_0^R$ and $Q_0^R$, can be obtained by:

$$I_0^R={{\displaystyle \int }}_{k T_c}^{\left(k+1\right)T_c}{S}^R\left(t\right)p^{I,D}\left(t\right)dt$$

(5)

$$Q_0^R={{\displaystyle \int }}_{k T_c}^{\left(k+1\right)T_c}{S}^R\left(t\right)p^{Q,D}\left(t\right)dt$$

(6)

where $k$is the measurement index, and $T_c$ corresponds to the coherent integration time used in the direct signal tracking stage (20 ms). Making use of the geometrical path difference model, the $I_0^R$ and $Q_0^R$ components are finally given in Equations (7) and (8), where ${\eta }_k^I$ and ${\eta }_k^Q$ are two independent zero-mean Gaussian noises and $\mathsf{\Lambda }( )$ models the normalized correlation function of the CDMA code.

$$I_0^R=\frac{A^R}{2}\mathsf{\Lambda }\left(-{\mathsf{\Delta }}_{p,k}/c\right)\cos \left(-{\psi }_k\right)+{\eta }_k^I$$

(7)

$$Q_0^R=\frac{A^R}{2}\mathsf{\Lambda }\left(-{\mathsf{\Delta }}_{p,k}/c\right)\sin \left(-{\psi }_k\right)+{\eta }_k^Q$$

(8)

The removal of the data bits on $I_0^R$ and $Q_0^R$ components is performed by using the $\mathrm{sign}$ function of the in-phase component of the direct signal [28] as follows:

$$\begin{gathered} D=\mathrm{sign} \left(I_0^D\right), \\ {I}_b^R= I_0^{R}D, \\ {Q}_b^R=Q_0^{R}D \end{gathered}$$

(9)

The complex representation of the resulting reflected signal is given by its phasor, ${\gamma }_b^R=I_b^R+i{Q}_b^R$. Subsequently, a moving mean filter is applied to extract the higher frequency component of the reflected signal and remove the low-frequencies contribution of the direct signal. The complex reflected filtered signal is given by ${\gamma }_a^R=I_a^R+i{Q}_a^R$. The phasor ${\gamma }_a^R$ is passed through a retracking module based on [28] to correct the time variations of the transmitter and receiver trajectory, the elevation angle, and the reflecting surface height.

The signal retracking starts by modeling the phase difference $\psi$ from the path difference model ${\mathsf{\Delta }}_p$. With $\lambda$ as the wavelength of the GPS L1 signal (~0.1904 m), the phase difference is obtained by:

$${\psi }_p=\mathrm{mod}\left(\frac{2\pi {\mathsf{\Delta }}_p}{\lambda },2\pi \right)$$

(10)

The phasor representation of the phase difference is built from:

$${\gamma }_p=\exp \left(-i{\psi }_p\right)$$

(11)

Finally, the retracked reflected signal denoted as the residual phasor ${\gamma }^R$ reads:

$${\gamma }^R={\gamma }_a^R {\gamma }_p^{*}$$

(12)

$$\Rightarrow {\gamma }^R=I^R+i{Q}^R$$

where * denotes the complex conjugate of the modeled phasor of the residual phase.

3.1.3. Spectral Retrievals

In addition to the signal as amplitude in the time domain, the PSD is used to reveal the power distribution of the residual retracked signal in the Doppler domain. We proceed as described in [34], and the PSD is defined by the Fourier amplitude as a function of the Doppler shift $f$ of the reflected signal, $\mathsf{\Gamma }\left(f\right):=\left|\mathfrak{F}\left\{ {\gamma }^R\right\}\left(f\right)\right|$, where $\mathfrak{F}\left\{\right\}$denotes the Fourier transformation.

The power spectral density is computed every 10 s for each satellite. Figure 6 shows the PSD of the reflected signal after the filtering $\left({\gamma }_a^R\right)$ and after retracking step $({\gamma }^R)$ for the GPS satellites PRN 30 $\left(E: ~9°\right)$ and PRN 8 $\left(E:~ 65° \right)$ on two different days, 17 July 2019 and 19 July 2019, the lowest and highest sea state, respectively.

In satellites with low elevations, a sharp spectrum with a noticeable peak is depicted. In contrast, high-elevation satellites show a spread spectrum with resulting peaks in the noise level. The effect of the sea state can be seen in PRN 30 on 19 July. In the presence of a higher sea state, the spectrum gets wider, losing its defined peak even when the elevation of this satellite is lower than on 17 July.

From each spectral retrieval integrated over 10 s, the five highest peaks of the spectrum ${\hat{\mathsf{\Gamma }}}_j\left(f_j\right)$ are selected, and the respective Doppler frequencies $f_j$ are retrieved from the x-axis location of the peaks in the spectrum. The Doppler spread ${\sigma }_f$ is estimated by computing the standard deviation of the five frequencies obtained in the step before. Figure 7 shows the power and the residual Doppler frequencies of the retracked signal on a logarithmic scale, including the elevation of the 8 satellites for each day. Low-elevation satellites concentrate in the center of the spectrum toward the relative frequency much more prominently on days of the lower sea state. The opposite is the case with high-elevation measurements, which become more distant. In addition, there is a reduction of peak power on days when the sea state is higher.

3.1.4. Residual Phase Retrieval and Tropospheric Residual Model

Based on the phasor ${\gamma }^R$ in Equation (12), a residual phase path after retracking is obtained. The observed residual phase ${\mathsf{\Delta }}_{\varphi }$of the retracked reflected signal is computed as follows:

$${\mathsf{\Delta }}_{\varphi }=\frac{\lambda }{2\pi } \left(\mathfrak{U}\left\{\arg \left\{{\gamma }^R\right\}\right\}+N\right)$$

(13)

where $\arg \left\{\right\}$ denotes the phase argument (defined between $-\pi$ and $+\pi$) and $\mathfrak{U}\left\{\right\}$ is the unwrapping operator to reconnect the phase argument by adding $\pm 2\pi$ when jumps greater than $\pi$ occur. The ambiguity $N$ remains unknown, and the factor $\frac{\lambda }{2\pi }$is applied to obtain the residual phase in path difference representation. The phase residual has been corrected for geometric components of the differential path (reflected signal relative to a direct one). It can be expected that troposphere residuals remain in the phase after retracking. Especially, the rather long propagation paths in the troposphere at grazing elevation (<15°). The modeled tropospheric residual ${\mathsf{\Delta }}_{tro}$ is obtained by the difference of the optical difference path, ${\mathsf{\Delta }}_{opt}$, using ray tracing as described in [23], assuming here the refractivity of a standard atmosphere [35] and the geometrical path difference ${\mathsf{\Delta }}_p$.

$${\mathsf{\Delta }}_{tro}={\mathsf{\Delta }}_{opt}-{\mathsf{\Delta }}_p$$

(14)

3.1.5. Zenith Total Delay Inversion

Similarly, as in GNSS processing for positioning and navigation applications, tropospheric propagation delay is an error source in GNSS-R. As presented in [36], the difference path delay between the path length of the reflected and the direct signal includes the zenith tropospheric delay that, multiplied by a mapping function, gives the differential tropospheric delay. According to [24], the tropospheric delay at position $P$ can be expressed as:

$${\rho }_{tro}^P=m_{hz}\mathsf{·}ZHD+m_{wz}\mathsf{·}ZWD\simeq m_{hz}\mathsf{·}ZTD$$

(15)

where $m_{hz}$ and $m_{wz}$ are hydrostatic, and the wet mapping functions, $ZHD$ and $ZWD$, are the zenith hydrostatic and wet delays. Therefore, ${\rho }_{tro}^P$ can be approximated by the product of the hydrostatic mapping factor and the zenith total delay, $ZTD,$ as the hydrostatic part accounts for 80–90% of the total delay. Assuming spherical symmetry in the atmosphere, the differential tropospheric delay $\widehat{{\rho }_{tro}}$ is twice the delay experienced between the specular point and the receiver [24]. An estimation of the delay at the specular point level is needed since there is no receiver in that location. Therefore, it is assumed that receiver height variations imply variations in the $ZTD$, so a height-dependent factor with exponential vertical decay, $h_f$, is included in the modeling. Finally, the approximation of the differential tropospheric delay is given by:

$$\begin{gathered} \widehat{{\rho }_{tro}}=2m_{hz}{h}_{f}ZTD-N \\ =2m_{hz}\left(1-e^{-\frac{H}{h_{scale}}}\right)ZTD-N \end{gathered}$$

(16)

where $H$ is the height of the receiver, and $h_{scale}$ is the scale height of the troposphere, assumed to be $h_{scale}=7160$ m [24]. The hydrostatic mapping factor $m_{hz}$ has been calculated from the Vienna mapping function (VFM3) [37,38], a discrete mapping function that employs the zenith distance $\left(90°-E\right)$, coordinates of the receiver, and day of the year.

As stated above, after the retracking step, troposphere residuals are expected in the residual phase from coherent observations. Thus, phase observations are used to inverse the $ZTD$ by using the linear regression approach in the form $y=ax+b$. The dependent variable $y$corresponds to phase observations ${\mathsf{\Delta }}_{\varphi }$. From Equation (16), the term $2m_{\backslash hz}{h}_f$ represents the independent variable $x$; the $ZTD$ is represented by $a$, which is the slope of the fitting line after the regression; and $b$ is the intercept, which is linked to the unknown phase ambiguity $N$ (see Equation (13)).

4. Results

4.1. Results on Residual Doppler Spread

The distribution of the Doppler spread on the four respective days is shown in Figure 8. For WS and SWH higher than 2.92 m/s and 0.26 m, respectively, the mean values of the Doppler spread are between 1.64 and 2.53 Hz with a dispersion that increases as the sea state increases. On the other hand, on 17 July, the mean Doppler spread value is 0.5 Hz, and the estimates are located close to this value, that is, a relatively low Doppler dispersion according to a calm sea for that date.

As discussed in Figure 7, there is a dependence between the Doppler shift and sea state depending on satellite elevation. Similarly, as presented by [39], to reduce the satellite elevation effect and represent sea-state impact on the residual Doppler shift, $f$ is multiplied by a mapping factor of $1/\sin \left(E\right)$ to obtain the mapped Doppler shift $f_v$. Figure 9 shows the mapped Doppler spread ${\sigma }_{f_v}$ for each day of measurements in the study area. The results show agreement between ${\sigma }_{f_v}$ and the sea state for each of the days. On 17 July, a blue response is shown for most of the satellite’s tracks, which corresponds to relatively low ${\sigma }_{f_v}$ levels $(0 Hz<{\sigma }_{fc}\le 2 Hz$). The mapped Doppler spread for the other days is higher, reaching values up to 10 Hz for some satellites.

The mapped Doppler spread is the estimate to correlate with the ancillary data from the ERA5 model for determining the relationship between the sea state and the Doppler spreading. To compute the Pearson correlation coefficient between ${\sigma }_{f_v}$ and the ancillary data, the mean value of ${\sigma }_{f_v}$ over 120 s is considered. The corresponding values of the ERA5 model parameters WS and SWH are interpolated at the specular point location from the grid files downloaded from Copernicus Climate Change Service (C3S) [40] with a spatial resolution of 25 and 50 km, respectively.

The Pearson correlations are presented in

Table 2. Pearson correlations between ERA5 parameters and mapped Doppler spread.

Parameter	Low	Mid	High
Wind Speed	0.88	0.66	0.58
SWH	0.75	0.58	0.56

. The events are classified depending on the elevation. Low events are satellites with $E\le 10°$, mid events are satellites with $10<E\le 30°$, and high events are satellites with $E>30°$.

Low-elevation satellites show a high correlation between ${\sigma }_{fc}$ and sea-state parameters. Although for mid- and high-elevation events, the degree of correlation remains moderate, it is significantly reduced with increasing elevation. Two aspects of our experiment may reduce the sensitivity at mid and high elevations: 1. The antenna gain decreases towards higher-elevation events that arrive at the aircraft closer to the nadir far off the antenna’s upward tilted boresight. 2. The high-elevation events lie close to or even on the beach where the sea-state effect is small.

4.2. Results on Carrier Phase Retrieval

In addition to the sea-state impact on Doppler spread, also a coherent carrier phase retrieval can be affected by sea state. As shown in [41], coherent GNSS-R observations tend to have continuous carrier phase measurements. The first analysis consists of examining the observed residual phase angle of the retracked signal computed as${\mathsf{\Delta }}_{\varphi a}=atan2\left(Q^R,I^R\right)$. Figure 10 depicts the ${\mathsf{\Delta }}_{\varphi a}$ of the satellites PRN 7, 30, 11, and 8 on two different days, including in the title the elevation angle range for each satellite. The first row shows the observations on 17 July, and the second row on 15 July. From initial inspection, PRN 7 and 30 show a continuous residual phase retrieval on 17 July. Continuous phase observations are lost when the elevation increases. At higher sea state (15 July), visual inspection suggests that no coherent observations are present except for the segment in the minute of the day from 760 to 763 for PRN 11, which corresponds to elevation angles between 5.9° and 7.5°.

The relation between Doppler spread and residual phase retrievals is presented in Figure 11 and Figure 12, on 17 and 15 July, respectively. As above, the same satellites with low, mid, and high elevations are analyzed on both days. The first row contains the Doppler spread over time ${\sigma }_f$ in blue for PRN 7, 30, 11, and 8. The red line is an empirical threshold established as 0.5 Hz based on the ${\sigma }_f$ mean value of the day with the lowest sea state. The second row contains the modeled and observed path differences. ${\mathsf{\Delta }}_{tro}$ denotes the tropospheric residual model calculated in meters, and ${\mathsf{\Delta }}_{\varphi }$ is the observed residual phase converted into a path as illustrated in Equations (14) and (13), respectively. The offset between the modeled and the observed path difference is due to the unknown ambiguity of the retrieved phase.

With a calm sea on 17th July, the low-elevation events $(E<10°)$ have Doppler spread below 0.5 Hz (PRN 7 and 30). These events present a smooth ${\mathsf{\Delta }}_{\varphi }$ that shows agreement with the tropospheric residual ${\mathsf{\Delta }}_{tro}$ in the path model, indicating coherent reflections. For satellites with mid and high elevations (PRN 11 and 8), this behavior does not remain and ${\mathsf{\Delta }}_{\varphi }$ is noisier, indicating a lack of coherent reflection.

On the other hand, with a rougher sea state on 15 July, the Doppler spread is significantly higher even for events below $10°$ elevation. Coherent phase retrieval fails except for a short period at the beginning of the event PRN 11 when the smooth phase appears at the lowest elevations $\left(~6°\right)$.

4.3. Results on Zenith Total Delay Inversion

From the coherent residual phase results, an agreement with the tropospheric residual model is found in the observations on 17 July 2019 at events with grazing elevation angles. Based on the observations and model correspondence, the $ZTD$ is estimated from linear regression as explained above. Figure 13 shows the residual phase observation versus the independent variable (scale factor) computed by using VFM3 and the height variation factor on 17 July 2019 and on 15 July 2019. The figure includes the linear fitting equation after the regression and the model uncertainty, which is the standard deviation of the observed minus predicted values.

The typical zenith total delay is about 2.3 m for a receiver located at mean sea level [42]. With smooth sea surface conditions (on 17 July 2019) at low-elevation angles, the estimated $ZTD$ for PRN 7 and 30 is 2.437 and 2.853 m, respectively, with an uncertainty of 0.02 m in both cases. Those events correspond to phase-coherent observations from which the estimated $ZTD$ only exceeds the expected value at 14 and 55 cm, respectively. Once the elevation angle increases and the phase coherence is lost, the $ZTD$ is far from the regular value. On the other hand, for events on the rougher sea surface (on 15 July 2019), there is no signature of phase coherence for a complete satellite track, so the slope of the fitting line after regression differs completely from the typical $ZTD$ value even for low-elevation observations (PRN 11).

5. Conclusions

The results show that loss of coherence in phase observations is accompanied by a Doppler spread of more than 0.5 Hz. The results also indicate a major influence of sea state in this respect depending on the elevation angle. As surface roughness and elevation angles increase, diffuse reflections dominate, and coherence in observation is lost. These findings apply to the here studied conditions of airborne observations over coastal waters. Similarly, as presented in [43] for reflectometry satellite measurements, it can be seen in this study that coherent observations respond to very calm waters and low-elevation angles. A Doppler spread ${\sigma }_f$ beyond 0.5 Hz causes a loss of coherent observations. In the given setup, the 0.5 Hz threshold is reached for SWH of ~0.3 m, and wind speeds of about 3 m/s when elevation angles go beyond $10°$.

However, even under coastal conditions, the coherent observations from airborne platform results are limited. For the 4-day measurement campaign, a total of 2646 Doppler spread estimations (every 120 s) are obtained. Out of these estimates, only 15% correspond to coherent observations below the ${\sigma }_f$ threshold. These retrievals are distributed as 10% and 5% for events at low and mid elevations, respectively, as shown in

Table 3. Total of coherent reflections below the Doppler spread threshold at low, mid, and high elevation.

Threshold	Low	Mid	High
${\sigma }_f\le 0.5 Hz$	10%	5%	0%

. None of the observations at high-elevation angles present coherence.

Nevertheless, it appears that even beyond the threshold, the Doppler spread carries information on sea state. By analyzing the Doppler spread of PRN 8 on 17 July, it is important to note the dependence on surface roughness even for high-elevation angles. Figure 14a shows ${\sigma }_f$ plotted with a vertical time axis aligned to the specular point track on the map. The highlighted section in yellow, which corresponds to the minute of day 806 until 809, represents a part of the track that crosses over very calm water near the port of Boulogne-sur-Mer (b) created by the breakwater built at the port location (c). Even though the elevation angle is $~70°$ in this track section, the smoother surface produces a significant reduction in the Doppler spread, nearly reaching the threshold of 0.5 Hz for coherence, allowing us to identify sensitive changes in sea state.

We found that the Doppler spread is an indicator for coherent observations, and it can furthermore carry sea-state information. However, it is not only the sea state that can affect the Doppler spread but also residuals that remain after correction with the geometric path model, in particular, variations of the receiver trajectory, which could have an effect. As discussed in the results sections, phase retrieval is possible under coherent reflection conditions below the 0.5 Hz threshold. Coherent phase observations over the sea surface are suitable for altimetric inversion. As described in [9,44], phase-coherent observations in GNSS-R allow altimetric retrievals from the carrier phase altimetry technique. Those studies have shown promising results in open ocean waters with high precision. This study may lead to potentially obtaining altimetric retrievals from airborne GNSS-R in coastal zones. However, when phase measurements are used, cycle slips caused by a low signal-to-noise ratio, ionospheric conditions, or aircraft dynamics need to be considered. Methods for cycle-slip correction are discussed in [24,39] by using airborne and satellite reflectometry data, respectively.

The standard GNSS antennas and receivers are capable of receiving RHCP signals only. At grazing angles, the reflected signal maintains its RHCP polarization. However, the reflected signal will change its polarization to left-hand circular polarization (LHCP) when it is coming from close to the nadir. A different antenna(s) setup, such as uplooking RHCP and downlooking LHCP antennas, as configurated in [32], may help capture the direct and reflected signals as they have a better field of view that may improve the final results.

GNSS reflectometry coherent phase observations also have the potential for atmospheric and ionospheric sounding, for example, to retrieve total electron content [45]. In this study, the comparison of coherent phase residuals and excess path model (tropospheric contribution) shows agreement. However, an offset remains due to phase ambiguity. The $ZTD$ estimations from linear regression and the expected zenith total delay value present differences at a centimetric level. Therefore, at low-elevation angles and smooth sea surface conditions, it is possible to retrieve zenith tropospheric delay from reflectometry airborne measurements. The literature suggests that the typical zenith total delay is approximately 2.3 m at mean-sea-level locations. On coastal calm waters (SWH: ~0.3 m, WS: 3 m/s), for elevation angles from 4.5° to 8.8°, the $ZTD$ estimation exceeds only 5% of the expected value. Once the elevation increases, from 8.3° to 10°, the difference is up to 24%. Future studies may use the sensitivity of phase observations to tropospheric contribution to retrieve parameters, for example, atmospheric water vapor, based on coherent reflectometry observations.