Retrieval of Sea Surface Height from CYGNSS Data with Tropospheric Delay

The Global Positioning System (GPS) which was first built for positioning, navigation and timing has been used to measure seismic tectonic motions, Earth orientation and polar motion, gravimetry, neutral atmospheric temperature and water vapor profiles, ionospheric electron density profiles and global monitoring (Beutler et al. ^[1]). Hall and Cordey^[2] first proposed that Global Navigation Satellite System-Reflectometry (GNSS-R) was used to infer a number of geophysical parameters, such as sea surface wind speed and roughness (Garrison et al. ^[3]; Katzberg et al.^[4]; Gleason et al. ^[5]; Clarizia et al.^[6]; Foti et al.^[7]) and sea surface height (SSH) which is used to study sea level rise, a critical factor in understanding Earth's dynamic climate, and used to study hurricane intensity, tsunami dynamics, El Ni o Southern Oscillation, eddy dynamics, ocean boundary currents, coastal and shallow water tides, as well as weather and climate forecasting (Lowe et al.^[8]; Hajj and Zuffada^[9]; Rius et al.^[10]; Cardellach et al.^[11]). Since the GNSS-R in Low Earth Orbit (LEO) with an antenna pointed toward the Earth's surface can track about 10 GPS reflections simultaneously, it provides a coverage that is of order denser than nadir viewing altimeters. Such dense coverage can be translated into a higher temporal and spatial resolution.

CYGNSS which was launched in Dec. 2016, uses a constellation of eight microsatellites to act as observatories, each carrying a Delay Doppler Mapping Instrument (DDMI). We attempt to retrieve the SSH from the DDMI data; however, the retrieval accuracy is affected by the following factors (Mashburn et al.): (1) orbit error and reflecting surface models; (2) ionospheric delay and tropospheric delay; (3) calibration of satellite timing and measurement biases. The ionospheric delay estimation methods are commonly based on the thin-shell ionospheric model or single-layer model (SLM) and mapping function, which assume the ionosphere as horizontally homogeneous and the ionospheric height as a fixed value. However, it was demonstrated in a number of studies that the height of the ionospheric shell is very important for accurate ionospheric delay estimation. Komjathy and Langley^[14] studied the variation in the ionospheric shell height, based on the International Reference Ionosphere 1990 (IRI90) and IGS data, and indicated that the ionospheric shell height ranges from 300 to 500km. By considering the temporally and spatially varying ionospheric shell height as opposed to a commonly adopted fixed shell height, they found that the differential delays were up to 1 Total electron content unit (TECU) in a solar minimum and could be larger during high solar activity conditions. Birch et al. ^[15] suggested that the thin shell height between 600 and 1200km is preferred, which is much higher than the commonly adopted value of 350km. Li et al.^[16] studied the ionospheric shell height in China and suggested it ranges between 450 and 550 km. That is to say, the optimal ionospheric shell height depends on temporal variation, geographic latitudes, and solar activities. Therefore, the ionospheric delay estimation is complicated, due to the limitation of the length of the article, the ionospheric delay estimation method will be published in another article. The information about the orbits and accuracy of the CYGNSS is better determined, since the DDMs are acquired with a slightly higher antenna gain (14.5 dBi) and lower noise figure than TechDemoSat 1 (TDS-1). Among these, the tropospheric delay is one of the main error sources in Global Position System-reflections (GPS-R) measurement. At present, the most effective method to study tropospheric delay correction is to divide tropospheric delay into two parts: a hydrostatic part, which constitutes more than 90% of the total troposphere delay, and a wet part, which is usually less than 10% of the delay (El-Mowafy^[17]). Scholars have proposed a variety of models for the zenith hydrostatic delay (ZHD), such as the Saastamoinen^[18-19], Hopfield ^[20-21], Baby ^[22], Davis^[23], and Askne and Nordius ^[24]. The zenith wet delay (ZWD) models proposed include the Saastamoinen ^[18-19], Hopfield ^[20-21], Baby ^[22], Askne and Nordius ^[24], Ifadis ^[25], and Berman ^[26]. The mapping functions which are used to map the zenith delay to the elevation-angle-dependent slant delay (Langley^[27]), include the Davis^[23], Ifadis^[25], Chao^[28-29], Herring^[30], Niell^[31-32], Global and Vienna mapping functions (Boehm, et al. ^[33]). Overall, Tuka and El-Mowafy^[34] thought the Davis, and Askne and Nordius models performed evenly with Saastamoinen model, as they gave the exact values for the ZHD.

Clarizia et al.^[35] used the GPS-R data from the TDS-1 satellite to obtain the SSH, the results show the good agreement with the global DTU10 mean sea surface height. Kucwaj et al. ^[36] applied the circular regression technique to estimate the SSH based on the GNSS-R phase delay observations, but the performances are not assessed with the real data. Mashburn ^[37] developed a tool-kit to produce state-of-the-art altimetric retrievals from the observed cross-correlation waveform of aircraft and spacecraft, the specular reflection timing was derived from the delay of the 70% peak correlation power (HALF method), the waveform leading edge peak first derivative (DER method), or the delay associated with a best fit function approximating the nominal waveform shape (PARA3 method). The tropospheric delays on the reflected signal path are accounted for using the UNB3m model (Leandro et al. ^[38]). It was found the HALF method produced the most precise measurements for a 5 second integration time with a standard deviation of σ = 0.6 meters. In the paper in order to find a more accurate method to get the SSH from the CYGNSS data, we attempt to exploit the leading edge derivation method (LED) while considering the Saastamoinen and numerical tropospheric delays.

The data set used for this study is described in Section 2. Section 3 determines the specular reflection point. Section 4 shows the tropospheric delay, including the numerical tropospheric delay method (NM) and Saastamoinen method (SM). The SSH estimation algorithm, LED, is given in Section 5. Section 6 shows the results and comparison of the SSH obtained from CYGNSS with that by Jason-2. Section 7 draws the conclusions and discusses some possible future work.

The data set of the geo-located Delay Doppler Maps (DDMs) are available to the public on a website (https://podaac-opendap.jpl.nasa.gov/opendap/allData/cygnss/L1/v2.1). These Level 1A DDMs are provided at a spatial resolution of 17 delay × 11 Doppler bins, corresponding to a surface area of about 50 km², at a time resolution of 1 Hz. Full details about CYGNSS can be found in Clarizia and Ruf^[35], for example, the Specular of DDM collected over sea from the CYGNSS on May 10, 2017 is shown in Fig. 1.

Figure 1.  The specular of DDM collected over sea from the CYGNSS on May 10, 2017.

The DDMs show the expected"horseshoe"-shaped characteristic (Fig. 2) of spaceborne GNSS-R correlated power over the ocean with the DDM peak corresponding to the area around the specular point on the ocean surface.

Figure 2.  An example DDM measured by the CYGNSS, showing the spatial distribution of the ocean surface scattering. Scattering cross section is plotted as a function of Doppler shift (x-axis) and relative propagation time of flight (y-axis), which is measured in units of Coarse Acquisition GPS Code, or"Chips".

The SSH from the OSTM /Jason-2 shall be used as the reference truth with a globally averaged RMS accuracy of 3.4 cm (OSTM / Jason-2 Products Handbook^[39]), assuming 1 second averages. The Jason-2 products generated and distributed by EUMETSAT in near real-time are also archived in the multi-mission EUMETSAT Archive. Any user can access Jason-2 data from the EUMETSAT archive upon registration. Or they can be directly downloaded from ftp://ftp.nodc.noaa.gov/ pub/data.nodc/jason2/gdr/gdr/.

The point on the Earth between the transmitter and receiver where the signal undergoes a specular reflection will satisfy several conditions. This point is commonly called the specular point and can be characterized by the following properties (Wu and Young ^[40]):

(1) The total path between the transmitter, specular point and the receiver will be the minimum of all possible travel paths.

(2) The specular point must lie on the surface of the Earth. For reflections from the ocean the specular reflection point can be reasonably assumed to lie on the WGS84 Earth geoid.

(3) The specular reflection must satisfy Snell's Law, or the angle between the incoming wave and reflected waves with respect to the surface normal must be equal.

In order to find the point on the Earth's surface that satisfies the above conditions, we first need to represent the signal path magnitude as a function of the unknown specular point location,

where $\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over T} $= the transmitter (i. e., GNSS satellite) location in the WGS84 reference frame, $\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over R} $= the receiver location in the WGS84 reference frame, and $\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over S} $= the specular point location in the WGS84 reference frame.

This expression for the path traveled can be expanded into three dimensions as follows,

Normally, the receiver location is known from the standard navigation output from the GPS receiver, and the transmitter location is either calculated during this same navigation solution or afterwards using data from the International GPS Service. It will be assumed that the transmitter and receiver locations are known, making the specular point location the only variable. As this equation is non-linear, an iterative method based on an initial guess will be used. In order to minimize this path we first take the partial derivatives of the specular point $\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over S} $ with respect to x, y and z. The partial derivative with respect to S_x is shown below with the results being identical with respect to S_y and S_z.

It can be noted that the denominators above are the incoming and reflected vector magnitudes, respectively. Simplifying the above equation and expanding it to include three dimensions follows as,

Iterating on $\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over S} $ using (4) will then result in a convergence to the minimum path. However, if you think about the true minimum location between $\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over R} $ and $\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over T} $, you will realize it will lie at the midpoint of the line connecting these two points and a great distance away from the Earth's surface, where the actual reflection occurred. This is the reason for restraining the correction to the Earth's surface as stated in the second condition above. This can be done using a simple scaling procedure at each new estimate of $\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over S} $. The radius of the Earth according to the WGS84 model can be calculated as a simple function of the specular point estimate z coordinate (i.e. latitude) as follows,

where a_WGS84 = 6378137 meters and e_WGS84 = 0.08181919084262 are the semi-major axis and the eccentricity of the WGS84 Earth geoid, respectively. The point on the Earth's surface that satisfies the three conditions listed above is then solved for iterative use of equations (6) and (7) below. A correction gain K has been added to quicken the convergence, considering that the initial guess for the location of $\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over S} $ (for example, the sub-$\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over R} $ Earth surface point) will be a great distance from the final solution.

where $\hat s = \frac{{{\rm{d}}\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over S} }}{{|{\rm{d}}\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over S} |}}$ is the directional unit vector for the correction.

This intermediate value is then converted to a unit vector and scaled by the Earth radius, giving us the new estimate for $\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over S} $ to be used during the next iteration.

The specular point can be considered found when the difference between the old and new values of $\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over S} $ falls below a specified tolerance after several iterations. Finally, as a last sanity check that the value of $\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over S} $ is correct, we can test the third condition listed above, which specifies that the Snell's law is satisfied with respect to the incoming and reflected wave directions.

Figure 3.  Geometry to determine the specular reflection point.

When the ZTD is considered, the numerical method (NM) is exploited as well as the Saastamoinen method (SM) with the best performance.

The ZTD in meters can be expressed as

where z is the height above the surface in meters, and N is the scaled-up refractivity (Smith and Weintraub, 1953^[41]); the refractivity of air is scaled by a factor of 10⁶, which can be expressed as

where p_v is partial pressure of water vapor, T is temperature, and p is the total pressure. The constant values a and b are 77.6 K hPa^-1 and 3.73 3×10⁵ K² hPa^-1, respectively. It is assumed that N is constant between the model levels.

Equation (10) can be used to calculate N_i for each model level i, using

where q is the specific humidity, T_i is the mean temperature in layer i centered on the i^th model level, and ε is the ratio of molecular mass of water and dry air. The ZTD starts at the highest model level and iterates downward, adding the delay calculated for each layer to the column total. Assuming the top of the model is high enough that the water vapor contribution to the delay is negligible, we can drop the second term in Eq. (10) at the top model level. Assuming hydrostatic equilibrium at the top of the model, it exists

where ρ_d is the density for dry air, p_top is the pressure at the model top, and ρ_∞ can be considered to be zero. Applying the ideal gas law for dry air therefore yields

where R_d is the gas constant for dry air and g is the acceleration due to gravity. Using the mean acceleration due to gravity at the surface, g₀ (9.80665 m s^-2), and the delay is added to the total delay from the other model layers. The mapping function uses the Rocken et al. ^[42] direct mapping technology based on the NWP. The numerical model data come from the NCEP reanalysis data (https://rda.ucar.edu).

In the Zenith angle direction, the delay is the smallest. In general, the relation between zenith delay and slant delay can be written as follows:

where θ is the elevation of the satellite, m_h (θ) is a hydrostatic mapping function, and m_w (θ) is a wet mapping function. The expression of the mapping function is as follows:

where θ is the satellite elevation angle at the specular point. For the hydrostatic mapping function m_h (θ), each coefficient a, b, c is obtained from the mean value a_avg and the amplitude a_amp,

where ϕ is the geographic latitude at the specular point, t is the day-of-year of the date, T₀= 28, and the coefficient is linearly interpolated from the data given in Table 1 and Table 2 according to the geographical latitude of the specular point.

The ZHD $\Delta {L_{{\rm{zh}}}} = \frac{{2.2767{p_s}}}{{f(\phi )}}$, here p_s is the surface pressure. f (ϕ) = 1 - 0.00266cos2ϕ indicates the the variation of the gravity acceleration with the geographical latitude ϕ on the sea surface. And ZWD $\Delta {L_{{\rm{zw}}}} = 0.002277(\frac{{1255}}{{{T_s}}} + 0.05){e_s}$, here T_s is the sea surface atmospheric temperature, and e_s is sea surface water vapor pressure.

Our SSH estimation algorithm is based on the leading edge derivation (LED) approach described in Hajj and Zuffada^[9]. The algorithm is applied to the delay waveforms. The delay difference Δτ is converted to SSH with knowledge of the measurement geometry, as given by

where h is the estimated SSH, c is the speed of light, R_t and R_r are, respectively, the transmitter and receiver ranges from the predicted specular point (in section 3) on the Earth's WGS84 ellipsoid, and H_t and H_r are, respectively, the transmitter and receiver altitude with respect to the plane tangent to the WGS84 Ellipsoid at the specular point.

The procedure of our work is as follows,

(1) Obtaining the DDM data, SSH from the Jason-2 and NCEP reanalysis data. A landmask is applied, which filters out the samples over land;

(2) Determining the specular reflection point as Section 2;

(3) Estimating the tropospheric delay by the SM and NM in Section 3, respectively. Since the reflected signal passes through the troposphere twice below the receiver altitude, and the direct signal does not pass through troposphere, the tropospheric delay is double;

(4) Calculating the delay difference by subtracting the tropospheric delay from the measured delay on the DDM;

(5) Retrieving the SSH by the LED method;

(6) Comparing the retrieved SSH from the CYGNSS with that from Jason-2 within 100 kilometers with the time window being 3 hours.

First, we conduct the SSH retrieval using the data from 06:29:53 UTC to 07:26:06 UTC, 28 November, 2018 by following above-mentioned steps. The resulted SSH bias is shown in Fig. 4. Here the NO scheme indicates that the tropospheric delay is not corrected.

Figure 4.  The SSH bias from 06:29:53 to 07:26:06, 28 November, 2018.

From Fig. 4 we can see that the bias from the NO scheme is the largest. Mean bias is about -2.3 m, while the result from the NM scheme is comparable with that from the SM scheme; however, SM scheme has the abnormal result in some special time and is inferior to the NM scheme as a whole.

The parameters in the SM scheme with annual mean and amplitude for temperature, pressure, and water vapour pressure varying with respect to latitude and height are computed for a particular latitude and day of year using a cosine function for the annual variation and a linear interpolation for latitude under the condition of the hydrostatic balance. If the atmosphere is nonhydrostatical, the SM scheme has greater error. However, the NM scheme performance has been investigated using reanalysis data, which is applicable to the hydrostatic or nonhydrostatical condition. Since the limitation of the SM scheme, the performance of the NM scheme is superior to that of the SM scheme.

16348 DDMs collected have been employed to retrieve the SSH from Mar. 2017 to Dec. 2018. We calculated the minimum, maximum, mean, median, mode, standard deviation and range of SSH bias (units: m) which are listed in Table 3.

Overall if the tropospheric delay is corrected, the SSH bias can decrease. The resulted SSH bias from the Jason-2 SSH by the LED retrieval method is of order decimeter. The SSH bias more obviously decreases with the NM scheme than that with the SM scheme. Although the SM scheme has the better performance, the NM scheme is superior to the SM scheme.

Since the transmitter and receiver clocks drift can exist, and the positions of the transmitter and the receiver have errors, the SSH bias is of order decimeter. If we track 8GPS reflected signals simultaneously, a receiver in LEO will observe 0.2million 4-s measurements of ocean height in one day. These are separated by about 25km in the direction of the reflection point motion and an average of 100 km between tracks. By dividing the ocean into small areas, mean sea height estimates, obtained from different GPS signal reflection within the same area, can be averaged to reduce the random errors by square root of the number of measurements.

The qualitative assessment of the results is promising; however, the bias is quite large but not surprising given that the characteristics of the DDMI on board CYGNSS are not optimized for SSH, and several large sources of error can be identified. The accuracy of the delay difference between direct and reflected signal provided in the DDM metadata is typical within one to two DDM pixel. Also, the error due to the limited receiver bandwidth and the error due to the noise affecting the waveforms are estimated to be between 3.5m and 7m depending on wind speed and geometry. Another significant source of error comes from the uncertainty in orbit knowledge. The receiver position calculated onboard differs from the one provided in the metadata by 5m or more. In addition, our errors also include the sea state bias, solid and ocean tides. In particular, the errors caused by the ionospheric delay also are nonnegligible. Finally, all the specular points are not colocated in the nadir point of Jason-2 at the same time, and this can cause the SSH bias.

The ability of CYGNSS to provide SSH information has been demonstrated using DDM data from the CYGNSS satellites. After considering the tropospheric delay with SM or NM method, SSH is derived from DDMs using the LED algorithms. The comparison between the SSH from the CYGNSS and Jason-2 is carried out. The overall SSH bias is comparable with precision estimates available in the literature. The NM scheme is superior to the SM scheme. Since the space-time coverage of the CYGNSS become quite dense over ±38° latitude, they will provide more useful SSHs for meteorologists and oceanologists. The work can be applied to the Beidou system in the future.