Ongoing Development of the Bass Strait GNSS/INS Buoy System for Altimetry Validation in Preparation for SWOT

Abstract

GNSS equipped buoys remain an important tool in altimetry validation. Progressive advances in altimetry missions require associated development in such validation tools. In this paper, we enhanced an existing buoy approach and gained further understanding of the buoy dynamics based on in situ observations. First, we implemented the capability to separate the ambiguity fixing strategy for different constellations in the processing software TRACK. A comparison between GPS and GNSS solutions suggested up to 3 cm reduction in the root mean square of the buoy minus co-located mooring SSH residuals over the selected sidereal periods. Then, comparison between double differencing and precise point positioning solutions suggested a possible common mode error external to GNSS processing. To assess buoy performance in different ocean conditions and sea states, GNSS and INS observations were used during periods where external forcings (waves, current and wind) were not interacting substantially. For the deployments investigated, no significant relationship was found, noting the maximum significant wave height and current velocity was ~2.3 m and ~0.3 m/s, respectively. In the lead up to the validation required for the SWOT mission, these results place important bounds on the performance of the buoy design under real operating conditions.

1. Introduction

Ocean altimetry has evolved considerably from pulse limited radar altimetry (TOPEX/Poseidon and Jason series missions from 1992) to Synthetic Aperture Radar (SAR) (Sentinel series missions since 2016), and now, the first swath-based SAR interferometric approach (SAR-In) with the anticipated Surface Water and Ocean Topography (SWOT) mission scheduled for launch in late 2022. Observations of ocean topography from altimetry missions have revolutionised the understanding of large-scale ocean circulation, eddy dynamics, and ocean tides [1,2,3,4,5,6]. Ongoing calibration and validation of altimetry has been an integral component of each mission design [7,8,9]. The iterative cycle of improved understanding of altimeter observations achieved through the validation process has assisted the improved geophysical interpretation that has followed. Over the progress of Jason and Sentinel series missions, the mission performance requirement for sea surface height (SSH) has advanced up to 3 cm [7,8,10]. This requirement, though defined differently for SWOT, has an ambitious goal of 1.2 cm for SSH at its minimal wavelength resolution [11].

Validation approaches have evolved to keep pace with developments in altimetry with the common aim of ensuring observations from the on-going altimetry missions remain in good quality. Among them, in situ approaches involving tide gauges (TG), moored oceanographic sensors and surface platforms (e.g., buoys, wave gliders) equipped with Global Navigation Satellite System (GNSS) receivers are widely used in validation activities [12,13,14,15,16,17]. The in situ observing platforms, such as moorings, unmanned surface vehicles, and buoys, can offer measurements at offshore comparison points (CP) and allow comparison between in situ SSH and altimeter-derived SSH to achieve a direct geometric validation, while off-site approaches observe SSH at a coastal location which is transformed to the offshore CPs for altimetry validation [18].

At the Bass Strait validation facility, a nearby TG, moored oceanographic sensors and GNSS equipped buoys are used to assess the fundamental accuracy in altimetry missions. Commencing with TG observations [19], the validation activity in Bass Strait has evolved to use a combination of in situ moorings and buoys [13,20]. In recent years, progress was made by Zhou et al. [17] to refine the buoy-derived SSH solutions. An Inertial Measurement Unit (IMU) was integrated into the buoy system addressing the previously ignored dynamic platform orientation. An initial attempt to address the bias induced by the tether of the buoy was also presented. They also examined the precision of the buoy derived SSH series by comparing solutions with the sea level record derived from the moored oceanographic instruments. While the study offered some insight into the precision of Global Positioning System (GPS) equipped buoys, the integration of non-GPS constellations was not considered, and the basic dynamics of the buoy were not investigated.

This paper advances the work by Zhou et al. [17] to further improve the understanding of GNSS and Inertial Navigation System (INS) equipped buoy measurements of SSH for altimeter validation. Using datasets from four buoy deployments of ~3–7 day duration, we aim to quantify the improvement gained by the enhanced treatment of GPS and GLObalnaya NAvigazionnaya Sputnikovaya Sistema (GLONASS) processing in the TRACK analysis suite [21]. We compare Double Differencing (DD) and Precise Point Positioning (PPP) solutions prior to characterising the dynamics of the buoy motion and investigating residual uncertainty with respect to the moored oceanographic sensors. Despite the difficulty of isolating external forcings from each other, we attempt to evaluate the performance of the buoy in different ocean conditions/sea states during each deployment and search for correlations between buoy minus mooring residuals and forcings such as wind, waves and current. Together, this work seeks to advance GNSS/INS-based determination of SSH that is required in order to keep pace with developments in satellite altimetry and make a meaningful contribution to validation of future missions.

2. Data

2.1. GNSS/INS Data

To assess the performance of the new GNSS/INS approach on the precision of SSH estimates, four deployments in total were undertaken at the historic Jason comparison point (JAS CP, latitude: −40°38.43′S, longitude: 145°35.41′E). Buoys were horizontally tethered to a surface float that was anchored in ~52 m of water (see Watson et al. [20] for details). This configuration constrains the buoy to a circle at the surface with radius of ~100 m (Figure S1). Deployment durations ranged from ~3.1 to ~7.4 days (Figure 1a), yielding a total of ~20.2 days (~485 h) of 2 Hz GNSS data (i.e., GPS + GLONASS) and 100 Hz INS data (note buoys are powered by internal batteries—power draw per buoy for the GNSS/INS ©s 2.3 Watts at 12 V”. Th’ four deployments provide different combinations of weather and ocean conditions useful to assess the performance of the buoy under different dynamic regimes. In Figure 1a, ocean conditions are quantified via significant wave height (SWH), calculated over hourly segments as four times the standard deviation of the SSH solutions.

2.2. Mooring Data

The sub-surface mooring located at the JAS CP is equipped with a SBE26+ bottom mounted pressure gauge (at ~52-m depth) along with a temperature and salinity sensor (SBE37) at depth of 25 m through the column. A Nortek Aquadopp current-meter was also installed at ~17-m depth to measure the eastward ($u$) and northward ($v$) current velocity in the upper water column, sampling with a ~20 min interval.

In this paper, we assumed that the bottom pressure derived SSH was representative of the SSH in the area at the multiple 100-m scale required for satellite validation. To achieve this, the mooring processing involves some temporal smoothing necessary to alleviate any higher frequency biases that occur over shorter scales. It is, however, not in the datum of the satellite altimetry and the GNSS buoy allows for the datum connection to be made.

Following Watson et al. [20], the bottom pressure is converted to water height by correcting for atmospheric pressure and the dynamic height of the ocean column, yielding the resultant time series with 5 min interval. The water height is then transferred onto the datum determined over multiple previous buoy deployments, thus now representing a geodetic height in the same reference frame as of the buoy solution, hereafter referred to as mooring SSH.

The estimated uncertainty of the mooring solution is at 12 mm including contributions from pressure sensor (5 mm), dynamic height correction (5 mm) and overhead atmospheric pressure correction (10 mm) [20]. In this study, we use the mooring SSH as the ground truth in the assessment of the SSH solutions from the buoy. It should be noted that the systematic error within the mooring is not the focus of this paper, and we return to this issue in the discussion.

2.3. ACCESS-G Wind Stress

Despite limited temporal and spatial resolution (hourly, ~12 km, respectively), wind stress from the Australian Community Climate and Earth-System Simulator Global (ACCESS-G) model is used as another measure to quantify external forcings on the tethered buoy during each deployment. The wind stress variable considers the roughness of the sea surface which has been seen as a useful proxy of the impact of wind on the platform for our investigation of the buoy dynamics [17].

3. Methods

SSH solutions from the buoy are derived using both TRACK [17,22] software developed by MIT and GipsyX [23] software developed by NASA/JPL. As we describe below, we modified TRACK to include INS observations and enhance GLONASS processing capability. Within our TRACK solutions, GPS/INS and GNSS/INS solutions were processed at 2-Hz/100-Hz, while for GipsyX, 1-Hz GPS-only solutions were generated. A 25-min moving mean smoother is then applied to obtain the smoothed SSH timeseries—selection of such filtering window length is determined via investigation of the optimal window length for the buoy solutions (Figure S2). Figure 2 provides an overview of GNSS/INS data processing procedure via TRACK. Further details regarding the general settings, correction models applied during processing can be found in the supplement material (Table S1). In the following sections related to GNSS processing, we focus on our latest enhancements made to the TRACK suite.

3.1. Processing of GNSS/INS Observations

3.1.1. GNSS Ambiguity Resolution

The four deployments considered in this paper provide the opportunity to assess the inclusion of GLONASS data, noting the work presented by Zhou et al. [17] used a GPS-only approach. This adds 24 GLONASS satellites to the 32 from the GPS constellation, creating a more robust satellite geometry [24,25], especially during certain windows when the number of GPS satellites above the cut-off elevation is low or poorly distributed. In turn, parameter estimation and their uncertainty are expected to improve.

Since TRACK release version 1.31, GNSS processing capability is implemented by mapping functioning frequencies from different constellations uniformly onto the GPS L1 (1575.42 MHz) and L2 (1227.60 MHz) band. While this approach enables TRACK to include more satellites to form more robust geometry at observing epochs, it ignores inter-system and inter-frequency biases across constellations and receivers. This will cause failures in the ambiguity resolution (AR) process, which DD relies heavily on to provide quality solutions.

Here, focusing on the AR processing, we enhance the treatment of multi-GNSS observations in TRACK. We extended the dimension of the variables storing the GNSS Differential Code Bias (DCB) table in TRACK to address the code bias issue in observations across different frequencies and constellations. On the other hand, the receiver-dependant bias in the phase observations is considered to largely cancel intrinsically in the DD processing when all reference stations and buoys share the same type of receiver. For our buoy deployments, this criterion is satisfied given we solely use Septentrio^TM AsteRx-2 receivers for all land-based stations and buoys. Finally, we implemented constellation specific constraints in AR to achieve optimal GNSS solutions. These changes and set-ups remove most of the limitations within TRACK for GPS/GLONASS processing.

As mentioned above, the approach taken in TRACK to include GLONASS data is to map the GLONASS observations across different frequencies onto the GPS L1 and L2 bands. During the process, inter-system and inter-frequency biases are two major concerns in high precision GNSS processing [26,27,28]. As a result, apart from an updated DCB table and the use of same type of receiver across all stations, we also modified the ambiguity resolution from a unified setting for all constellations into a constellation-specific strategy. This allows the users to manually control how TRACK deals with ambiguities for satellites from different constellation.

Given the complexity associated with the phase bias issue in GLONASS being hardware and even station specific [29,30], we adopted an epoch-wise approach on all GLONASS ambiguity estimation, effectively forgoing any attempt to fix GLONASS ambiguities in the AR process. Meanwhile, the GPS ambiguities are set to be fixed when possible, using appropriate parameter settings as recommended within the TRACK suite. Averaged ambiguity resolution rates for GPS were ~80% across all 7-h processing sessions and varied insignificantly before and after addition of GLONASS.

3.1.2. Outlier Detection and Removal

From the raw output of TRACK solutions, we performed a two-step outlier detection/removal procedure. First, a differenced series between the SSH and the reconstructed tide series based on a decade long mooring SSH series at JAS CP was calculated. Any values outside three standard deviations from the median of the residual series were removed (hereinafter referred as “3-sigma approach”). Following this, over each 2.5 h segment of the differenced series, any values out of the 1.5 interquartile range (IQR) from the mean of the residual series was further removed (hereinafter referred as “1.5-IQR approach”). The resultant series is considered as clean SSH solutions.

When comparing buoy solutions from DD and PPP approach or against the mooring SSH, we performed another round of 3-sigma and 1.5-IQR editing on the differenced residual to further remove outliers. For the buoy-against-mooring analysis, an additional absolute cut-off of 0.06 m for gross errors was applied to exclude any potential GNSS processing quirks caused by failure of ambiguity resolution. The detailed gap percentage of buoy-mooring residual for each deployment is provided in Table S2.

3.2. Characterising Buoy Dynamics

An idealised GNSS buoy platform precisely follows the instantaneous sea surface and provides uninterrupted sky-view to the overhead GNSS constellation. In reality, the hydrodynamics of a sizable buoy platform will influence the precision and accuracy of buoy derived SSH. When the buoy and water is at rest, the offset between the water level and the GNSS antenna reference point (ARP) is known precisely [20]. Under dynamic conditions, the instantaneous position of the antenna with respect to water surface will vary. Such variations are driven by the constant disruption of the hydrostatic equilibrium given movement of the water around the tethered buoy platform. This translates to the GNSS buoy derived SSH being possibly affected by the waves/current (and their interaction) in a complex way. The existing assumption is that such variation occurs at high frequency (largely driven by the waves) and has negligible contribution to lower frequency sea level signals. This may not be the case if for example, a sustained wave or current state systematically elevates (or lowers) the antenna position with respect to the water surface. We seek to explore this assumption by investigating the observed dynamics of the buoy alongside the residual bias between buoy and mooring SSH time series.

The variation of the reference water level can be further dissected into rotational and translational motions. The rotational motion is caused by a combination of tether forcing and an imbalanced momentum due to the asymmetric design of the buoy platform (Figure 1c), while the translational motion occurs when an external force in the form of a wave induces a momentary acceleration to the buoy platform. Both the rotational and the translational motion of the buoy can be effectively measured by the IMU onboard. Meanwhile, the translation motion can also be described the GNSS derived quantities. In the following sections, we have determined three quantities to assist us to better understand its dynamics, namely the vertical acceleration, the wave direction from both GNSS coordinate timeseries and INS observations, and the wave magnitude derived from spectrum analysis of GNSS coordinate timeseries.

3.2.1. Vertical Acceleration of the Buoy

Acceleration is the fundamental indicator of a non-equilibrium system, where forcings within remain imbalanced. Such forcing will lead to non-zero velocity which, when integrated over a certain period, results in displacement. When assessing potential bias in buoy derived SSH, it is inevitable that the vertical acceleration must be considered.

The accelerometer in the IMU directly senses the vertical acceleration including the local gravity in the body frame. To transform these accelerations into a topocentric East-North-Up (ENU) frame and exclude gravity from the raw measurements to derive the “free” acceleration, we need to apply (1) to the 3-axis accelerometer readings:

$$a_U=\sqrt{a_{bx}^2+a_{by}^2}·\sin \left({\theta }_{tilt}\right)+a_{bz}·\cos \left({\theta }_{tilt}\right)-g_{local}$$

(1)

where $a$ generally denotes acceleration, subscript $U$ indicates the vertical component in the ENU frame, subscript $bx$, $by$, $bz$ denotes $x,y,z$ readings from accelerometer in body frame, ${\theta }_{tilt}$ indicates the tilt of the buoy platform, where tilt is defined as the angle between z-axis of the platform in body frame and the upright direction in the ENU frame, and $g_{local}$ represents the local gravity constant. Schematic derivation for Equation (1) can be found in Figure S4 in the Supplementary Material.

Alternatively, to describe the instantaneous acceleration at GNSS observational epochs without introducing any interpolation errors, we derive acceleration using three sequential instantaneous SSHs via Equation (2):

$$a_n=2·\frac{\left(S_{n+1}-S_{n-1}\right)·t_n-\left(S_{n+1}-S_n\right)·t_{n-1}-\left(S_n-S_{n-1}\right)·t_{n+1}}{\left(t_{n+1}-t_{n-1}\right)·\left(t_{n+1}-t_n\right)·\left(t_n-t_{n-1}\right)}$$

(2)

where subscript $\left(n-1\right), n, \left(n+1\right)$ denotes three consecutive epochs, $a$ denotes acceleration, $S$ denotes instantaneous SSH at corresponding epochs and $t$ denotes the time for the epochs. Schematic derivation for Equation (2) can be found in Figure S5 in the Supplementary Material.

3.2.2. Wave Direction Information as Sensed by the Tethered Buoy

Given the potential for interaction with ocean currents, an important sea state parameter that can be derived from the buoy-based IMU is the wave direction. The direction of the wave not only provides key information when examining the platform dynamics but also serves as a valuable derived quantity in oceanographic studies, with possible use in validating regional wave models [31,32]. However, it is worth mentioning that given the buoy is tethered rather than free drifting, wave directions derived from either GNSS or INS are not independent of the current. We will address implications of this later in the discussion.

Via Equation (3) [33], the wave direction from INS can be calculated:

$${\vartheta }_w=\mathrm{acos}\left(\sqrt{co{s}^2\left(\alpha \right)-si{n}^2\left(\beta \right)}\right)+\phi$$

(3)

where $\alpha ,\beta , \phi$ are the pitch, roll and yaw of the platform, ${\vartheta }_w$ is the wave direction. Both yaw and the wave direction are referenced to the geographic north. For comparison, using the GNSS coordinate time series, the direction of the smoothed horizontal displacement of the buoy over a given time interval (e.g., 5 min) in the ENU frame is derived and used as the indicator of the wave direction from GNSS.

3.2.3. Wave Magnitude Based on SSH Spectrum Analysis

For the magnitude of the waves, we rely on the power spectrum analysis of the GNSS/INS SSH solutions. The waves are categorized into two frequency bands to reflect the wave magnitude for swell and wind waves. The swell dominates the frequency band with period between 12 and 20 s, while the wind waves are considered between a period of 3 to 6 s. The wave magnitude is important for assessing the response from the buoy to the impact of the waves from both rotational and translational perspectives.

The procedure to derive wave magnitude can be described as follows: first, a spectral density estimation is performed using Welch’s method [34] on segments of 5 min detrended SSH solutions from the buoy. During the estimation, a Hamming window of 1 min is used to smooth the spectrum. Then, the average signal power density from the corresponding frequency bands is determined, yielding a magnitude time series on a 5 min interval for swell and wind waves, respectively.

4. Results

4.1. GPS/INS versus GNSS/INS Buoy Solution—AR Update

We first assess the enhancements made to GNSS processing. Figure 3 (upper panels) provides a visual comparison between segments of SSH over the same sidereal period of four consecutive days in May 2021, showing 2 Hz GPS/INS and GNSS/INS solutions. To highlight the differences between the two sets, solutions shown are smoothed using a 25-min moving average low-pass filter.

Several apparent spurious SSH signals can be seen in the four solution segments selected from representative deployment 66. At least centimetre-level biases between solutions are evident. Due to the intrinsic repeat cycle of the GPS constellation being close to 23 h 56 min 4 s (i.e., one sidereal day), it is observed that these relatively large biases occurred repeatedly over consecutive days of a week-long deployment, as is evident in Figure 3. Analysis of the GPS constellation observed at these times (Figure 3 lower panels) shows only five GPS satellites are in view above the elevation cut-off threshold (10°), one of which is of low elevation. With GLONASS introduced, up to 15 total satellites were tracked, forming a better geometry for estimation leading to solutions with noticeably reduced deviations from a smooth trajectory (Figure 3 upper panels). The different repeat period for the GPS and GLONASS constellations results in varied geometries at sidereal timescales—note the GPS satellites (red) shown on the sky plots in Figure 3 remain in almost identical locations.

Comparison with the mooring shows that the GNSS/INS solutions are improved over the GPS/INS solutions over these times. Summary statistics are shown in Figure 3 (upper panels), revealing a 17–45% reduction in root mean square (RMS) for these periods when adding GLONASS to the solutions.

Considering a substantially longer period spanning deployments 63–66, shows that the addition of GLONASS tends to improve the solutions, but not always. When compared against the in situ mooring SSH, the GNSS/INS solutions typically outperform the GPS/INS solutions. For deployments 63 and 64, the standard deviation of the buoy-mooring residual is reduced by 0.5 cm and 0.3 cm, respectively. For deployments 65 and 66, the two solutions are more comparable, with GNSS/INS showing a standard deviation 0.1 cm larger than GPS/INS. Nevertheless, it is noticeable in the time series (Figure 4a) a section of GNSS/INS solution is excluded as outliers for deployment 65—we will return to this in the discussion. Overall, the range of the buoy-mooring standard deviation using GNSS/INS is from ~1.7 cm to ~2.6 cm, showing some dependence on the average significant wave height over the respective deployment (Figure 4a). The histogram of Buoy-minus-Mooring (B-M) residuals for all deployments shows a marginal 0.1 cm improvement from GPS/INS to GNSS/INS solutions (Figure 4b). Importantly, the GPS/INS residuals are slightly negatively skewed, while the GNSS/INS histogram is more symmetric about zero.

In Figure 4c, the residual spectra show slightly reduced long-period noise in GNSS/INS solutions although with some increased high frequency noise. There is some evidence of a ~12 h periodical signal in both solutions; however, the energy at this band is reduced in the GNSS/INS solution. In the band up to 6 h, both residual series express coloured noise [35] with a slope between −2 and −1 but flattening at longer periods, similar to results from land-based stations in another study [36]. The GPS/INS solutions have a slightly higher noise level between periods of 1 to 2 h. We note this signal includes contribution from the uncertainty in the in situ mooring SSH, and we return to this in the discussion.

4.2. DD versus PPP Buoy Solution

In this section, SSH solutions are further assessed from a perspective of processing methods. At the Bass Strait facility, DD processing using TRACK is possible and potentially preferrable to PPP given the relatively close proximity (20–50 km) of land-based reference stations (Figure 1). The alternative PPP approach (e.g., via GipsyX) offers the advantage of positioning without the need for land-based sites [37]. Figure 5 shows a comparison between GPS-only solutions from both PPP and DD for deployment 66. For this deployment, PPP shows a marginally lower RMS of the residual computed against the in situ mooring. The residual timeseries also show similar temporal variability indicating some common mode error. This is potentially associated with the in situ mooring. When investigated in daily segments of the residuals (Figure 5a), as RMS of the segments decreases, the correlation coefficient between DD and PPP solutions increases reaching a maximum of 0.41.

Taking the RMS of the residual of each solution from the mooring values, it is shown that both solutions are at a precision of ~1.5 cm, which is consistent with previous studies in Bass Strait [17,20]. Figure 5c shows the spectra, indicating that the PPP solution has more energy around the 12-h band, while the DD solution includes a signal near 6-h period. For the high frequency bands, the PPP solution tends to have higher energy. The slope in the high frequency bands implies both solutions express coloured noise up to at least 6 h which we return to in the discussion.

PPP solutions for the other three deployments were more challenging yielding several periods where the solution filter had diverged resulting in periods of output data being flagged as outliers (amounting to 17% for deployment 63 and 9% for deployment 64, see Figure S6—while for DD solutions, 11% for deployment 63 and 6% for deployment 64, see Table S2). This is, in our experience, unusual, and data quality issues associated with a drifting receiver clock may potentially be the source of this issue (resolved in the DD approach through differencing)—this remains under investigation. When excluding these outliers, the buoy-minus-mooring SSH residuals are comparable between DD and PPP approaches, suggesting that solutions from both have similar quality, whereas DD solutions appears more robust in general for these four deployments when using GPS only solutions.

4.3. Measured Indicators of the Buoy Dynamics

We now present the observed buoy dynamics with the aim to provide further insight into the buoy performance and potentially further understand the characteristics of the B-M residual. In Figure 6, metrics representing possible drivers for the buoy motion are shown for deployment 66 (see supplementary material for the equivalent plots for other deployments, Figures S7–S9), including the magnitudes of wind waves, swell and current along with magnitude of accelerations. These metrics are shown without detailed quantification since the time correlation among these are the focus rather than their absolute values.

For INS acceleration, 100-Hz accelerometer readings are smoothed via a 2-Hz moving mean window to reduce the systematic noise of the unit, while GNSS accelerations is derived using Equation (2). Buoy tilt is in the unit of cosine of the tilting and is derived from adaptive ARP correction [17]. Magnitude of wind waves and swells are derived based on methods previously described and are normalized within [0, 1] based on their range separately, hence they are not suitable for quantitative inter-comparison. GNSS waves are derived as a combination of GNSS wave direction and wind stress amplitude from the model—interpolated to 2 h in Figure 6. Quivers of the wind stress indicates the direction and is scaled with its amplitude and presented at also 2-h resolution. Currents are observed by the in situ current-meter with a sampling rate of 20-min. The directions are further down sampled to 1 h in Figure 6. Quivers of the current showed tidal behaviour in both amplitude and direction.

In Figure 6, we see that accelerations from GNSS and from INS are consistent in time and in phase with the observed tilt of the buoy. Over the duration of the deployment, three distinct periods characterised by increased sea state were evident. To facilitate further investigation, these periods are designated: window #1 (SWH from 0.6 to 0.8 m), window #2 (SWH from 0.5 to 0.7 m) and window #3 (SWH from 0.7 to 1.3 m). Both accelerations and the buoy tilt correlate with these transitions in conditions.

Based on our derived wave information, the first two elevated periods of SWH were likely driven by local increases in wind speed given there was no associated increased energy over swell periods. In Window #3, it can be observed that a far-field energy propagated through the water as evidenced in the swell time series prior to Window #3. Within the window, a local frontal system was observed, and the amplitude of the wind waves started to increase accordingly. As a result, the buoy motion was altered in both acceleration and buoy tilt variables. The wind stress series also captures the progression of this event.

In the lower panel of Figure 6, B–M residuals are presented to provide a sense of temporal variability in response to changing conditions during the deployment. There is some indication within three windows that higher B–M residuals are associated with higher sea states. This again indicates the buoy dynamics are altered under the impact of external forcings given the underwater moored sensors are insensitive to surface varying conditions.

A close-up investigation of the buoy dynamics over a period of ~10 s was carried out to study the orientational/horizontal motion under impacts from external forcings. The study period was chosen void of any gaps and processing outliers to best present typical motion of the buoy observed by the GNSS/INS sensors during the deployment.

Figure 7a shows that the trend of the instantaneous unsmoothed SSH is in generally good agreement with the titling and the decomposed roll and pitch, despite the orientations being at a much higher sampling rate (i.e., 2 Hz SSH versus 100 Hz orientations). The yaw in Figure 7a describes a rather stable heading of the buoy—varying within 10 degrees, while the horizontal position from GNSS along with the black arrows showing epoch-wise movement in Figure 7b indicates a circular or orbital motion within a general south-eastward trend.

In Figure 7b, we also provided information of four horizontal directions related to the buoy motion: model wind stress, waves derived from each of GNSS and INS based on methods previously described, and ocean currents from the mooring. Waves derived from GNSS and INS instrument show good consistency with each other. The observed current on the in situ mooring provides information of the surface currents, which is the dominant driver of the buoy route in general over this period. Wind stress direction from ACCESS-G model, provides the impact from the wind on the buoy with the consideration of the sea surface roughness. Given the low resolution (hourly, 12 km) of the model, the agreement between the observation-based wave directions and the model is reasonable (GNSS/INS waves versus Model wind in Figure 7b).

4.4. Bias Identification by the Buoy Dynamics

Having inspected the metrics describing the buoy motion, there appears to be some evidence suggesting momentary correlation between the buoy dynamics and the external forcings (Figure 6 and Figure S7–S9) which warrants further investigation. In this study, the B–M residual is considered as an indicator of biases induced by the buoy dynamics under the impact from all the near-surface forcings, given the in situ mooring at the deployed depth is less sensitive to these forcings, especially at high frequencies. Nevertheless, this assumes the in situ mooring is error free which we return to in the discussion.

To determine if there is any significant correlation between the B–M residual and the dominant external forcings, two typical cases that isolate specific forcings were extracted from the four buoy deployments including: (1) low current with moderate wave state and (2) moderate current with low wave state. Over any given buoy deployment, such scenarios are rarely isolated from each other—the waves and currents are often simultaneously driven by passing oceanic/meteorological conditions. In our analysis, we define low current by a current velocity smaller than 0.05 m/s, while the low wave condition is determined by thresholding the tilting of the buoy platform below 10°.

Case #1—Low Current & Moderate Wave State

Under this scenario of external forcing setup, a quadratic fit between the buoy tilt and wave magnitude (calculated as described in Section 3.2.3) is shown in Figure 8. A window of 5 min is used to derive the standard deviation of all the segments of the buoy tilt series centred on time of each derived wave quantity for the full span of the buoy deployment. As expected, the scatter of the buoy tilt increases while the energy in the waves surrounding the buoy increases (Figure 8a). In other panels of Figure 8, three sets of quantities are fitted to the B–M residual to further investigate the relation between external forcings and the buoy dynamics.

Wind stress reflects some combination of wind and surface water forcings on the buoy which is tethered in place. In this scenario, a quasi-linear relationship is found between the wind stress projected to the tether direction using conventional 2-D vector projection and the B–M residual as shown in Figure 8b. For two different deployments with data fitting this scenario, the slope of the two fitted linear models appears the same, though the absolute values vary substantially (Figure 8b). Since both proxies for wind stress magnitude and direction from the model could potentially be substituted by quantities derived from observations, two further modified inputs are provided. First, we replaced the magnitude of the wind stress with GNSS-derived wave magnitude, maintaining the directions from the wind stress model. Second, we did likewise but for direction we substitute the observed wave direction from the INS data. Using these two models showed negligible relation (Figure 8c,d) with the B–M residual confirming that the original linear fit with modelled wind stress was driven by the small sample and insufficient spatial/temporal resolution of the model.

Case #2—Moderate Current & Low Wave State

The second case was used to determine if there is any specific relationship between the current and the B–M residual. To test such a case, the waves induced from either a far field or a near field source need to be absent, ideally leaving the current as the dominant external forcing. 180 min of data (B–M residual samples at 5 min intervals) from buoy deployment 64 was found to meet the input criteria.

As shown in Figure 9a, a linear relation between the current projected onto the tether and the B–M residual is observed. However, there is considerable scatter around this linear fit and some systematic errors are evident (e.g., the B–M residual remains consistent at around 5 cm when the current increases from ~0.2 m/s to its maximum at ~0.3 m/s in Figure 9a). A similar oceanic condition (lasting merely 40 min) during deployment 63, suggests further uncertainty since when the current increased from 0.08 to 0.27 m/s, the B–M residual remains within the range of 0 to 2 cm, as shown in Figure 9a.

In Figure 9b, when taking the “moderate current + high wave” condition into consideration as a contrasting case—samples over ~2 h from deployment 64 and ~1.5 h from deployment 66, the B–M residual stayed mostly positive. In Figure 9b, significant scatter of the B–M residual is observed when the current is measured within the range of 0.15 to 0.2 m/s suggesting that any relationship remains equivocal using data from these four deployments.

5. Discussion

5.1. Addition of GLONASS

GPS solutions from buoys have contributed to the altimetry validation activities in Bass Strait for almost two decades [13,17,20]. For our double differencing approach, past validation activities have shown limitations in solutions at times of poor satellite coverage. These relatively short periods of time (~40 min) repeat each sidereal day given the repeat period of the GPS constellation. The introduction of GLONASS, which includes 24 satellites and a repeat period of eight sidereal days, solves this issue by improving the solution geometry (Figure 3).

Figure 3 clearly demonstrates the improvements achieved when adding the GLONASS constellation during periods of poor GPS coverage. Without the addition of GLONASS, the detection of these periods of degraded solution quality is often problematic. If left undealt with, the buoy solution is clearly biased which propagates into the uncertainty of the altimetry validation. With GLONASS included, more satellites can be observed at each epoch, providing more confidence and improving geometry for the SSH solutions.

An indicator of improved quality of the GPS and GLONASS solutions is the significant drop in RMS when comparing against the in situ mooring. In Figure 3, a reduction of up to 3 cm in the RMS of the B–M residual over the selected windows is observed. On average, a ~2 cm improvement can be seen across the consecutive four-day period centred on the same sidereal period when using GNSS instead of GPS only. The improved solution geometry is clearly indicated in the sky plots shown in the lower panels of Figure 3. It is also obvious in the spectrum in Figure 4 where a noticeable flattened level of noise is observed in the high frequency domain below the 2-h band. In general, improved geometry is associated with reduced level of uncertainty when forming normal equations during filtering process, which will further help build a more robust estimation of the parameters.

The sidereal multipath on the carrier phase is likely to affect the base stations onshore more than the offshore buoys. The orientation of the buoy antenna with respect to the reflection surface (water) is dynamic and highly complex. The reflection surface for land-based sites is, however, static. Hence, in a DD approach, multipath from the reference site will influence the buoy parameter estimates. This impact could be one of many possible causes of the 6-h sidereal harmonic presence in the DD time series (Figure 5), as shown also in other studies [38,39]. When GLONASS is included in the processing, the impact of this is reduced given the overall geometry no longer repeats, and the contribution of the GPS multipath is diluted.

There remain areas for further improvement in a coupled GPS/GLONASS and INS solution. For some deployments, the standard deviation of the SSH solution has increased rather than dropped, e.g., deployment 65 from 2.4 cm to 2.5 cm. This is due to the fact that GNSS ambiguity resolution in TRACK is presently limited by the inter-frequency and inter-system bias in multi-constellation processing [27,30]. Moreover, the possible reason behind a large segment of outliers removed from GNSS/INS solutions for deployment 65 could be traced back to the ambiguity resolution process where the float GLONASS ambiguity with large uncertainty gets passed to the final adjustment phase, which further results in unacceptably large errors distributed to the positional state vector. Our inclusion of a separated bias-fixing strategy—integer for GPS whereas float for GLONASS, has so far worked well to avoid the ambiguities getting adjusted incorrectly.

5.2. Processing Strategy

Double differencing is the traditional way of positioning where a rover station is referenced against a fixed static station. Through inter-satellite and inter-station differencing, the method cancels common unwanted signal and noise. The DD approach becomes limited when the unwanted signal or error contribution to two stations becomes spatially decorrelated—typically beginning at baseline separations of 20–50 km. In the marine domain, this limits high precision applications to near coastal regions. On the other hand, PPP, removes the constraint of needing a reference station, noting it still relies on the global land-based network for computations of orbits and satellite clocks. The PPP approach typically requires an extended length of time to converge and, historically, the precision is not as high as the DD approach over short (<~20 km) baselines [40]. Moreover, because of its reliance on the functional model including and addressing all possible errors (rather than cancellation as occurs in DD), it relies heavily on the sustained quality of the observations. As such, periods of filter divergence were evident in PPP solutions for deployments 63, 64 and 65 (Figure S6). We note typical percentages of outliers from PPP was ~10% compared to 6% for DD (Table S2). We can potentially attribute these challenges in the PPP solution to issues associated with clock steering on our receivers, noting that DD is not similarly affected given differencing.

From the perspective of our buoy deployment, advantages of a traditional DD method include robustness and high precision, both of which depends on the cancellation of path delay errors between the buoy and the reference station onshore. However, the disadvantage of DD method is, as previously indicated, the constrained distance between the base station and the buoy and the addition of base station multipath. The former limits the potential deployment configuration to remain in relative proximity to land based reference sites.

The future interferometric mission, SWOT, requires two-dimensional validation with strict precision requirement for (sub-)mesoscale ocean processes. Some in situ deployment configurations for SWOT validation may stretch the traditional DD approach given its baseline distance constraint. This justifies our re-examination of both the DD and PPP method to derive solutions in the marine domain, since PPP has evolved quickly over the past decade [41] and its precision is now comparable [42]. From Figure 5, the two methods have approximately equivalent precision for the deployment configuration used here, achieving precision below the 2-cm level. In the frequency domain, the PPP may have a higher level of noise in higher frequency bands, but it is not impacted as much as the DD method by the sidereal harmonics in the lower frequency bands, likely driven by multipath effects at land-based reference sites. The parallel development of both approaches remains critical given the advantages of them becoming eminent in very different deployed areas—DD in coastal areas, while PPP more offshore adaptive.

Another key message from the comparison of the PPP and DD solutions against our in situ mooring (Figure 5) is that despite the intrinsic positioning mechanism differences, the two series show many in-phase signals (noting daily correlation values only reach 0.41). This suggests, however difficult to validate, that biases within the buoy derived SSH solutions as compared against the mooring are not entirely GNSS processing method dependant. Instead, error contributions are expected from the motion of a tethered buoy during the deployment, in addition to error contributions from the mooring SSH time series. The first issue reinforces the necessity to characterise and investigate the buoy motion using the on-board sensors, i.e., GNSS and INS.

5.3. Contribution of the Buoyancy Variation Due to Dynamics

To further understand the impact of buoyancy variation (i.e., change in antenna position with respect to water level) on the SSH solutions, we need to properly characterise the buoy dynamics. As a preliminary attempt, the IMU on-board is used in this study to quantify the vertical motion induced by external forcings that put the buoy platform in a non-equilibrium condition.

The assumption of a static reference waterline with respect to the antenna is not in reality true during buoy deployment. While the cylindrical design of the buoy platform and outer floats (Figure 1d) was conceived with the intent of minimising the generation of lift as a function of water velocity, the tether of the buoy inevitably introduced a rotational variation of the platform by inhibiting its translational motion when the tether is stretched to its full length. The buoy needs to work to maintain the equilibrium state, and this is when the nominal buoyancy level may shift the reference waterline (either vertically up or down depending on the hydrodynamics of the platform), which further induces a bias into the instantaneous SSH measurements. Chupin et al. [16] has attempted to deal with this issue by observing it via an acoustic sensor. Such a technique shows promise yet is increasingly problematic for sustained tethered deployments in higher sea states where the water in the vicinity of the platform is highly dynamic and often in the form of turbulent spray/wash.

The vertical translational and rotational motions are likely key variables in quantifying any biases induced by external forcings. In Figure 6, we confirm that both the GNSS and IMU can describe the vertical acceleration of the buoy well, judging from the near symmetric shape of the acceleration time series. The buoy tilt and derived orientations in both Figure 6 and Figure 7 further verify the ability of the IMU to characterise the rotational motion of the buoy. During large wind events (e.g., Window #3 in Figure 6), the buoy tilting and the vertical acceleration in response to the increased sea state needs proper treatment, otherwise the SSH derived over such a period is going to be biased and will have negative impact on altimetry validation.

The magnitude and direction of each external forcing is critical in order to understand any impact on SSH uncertainty. In Figure 7, from the representative segment shown, we see relatively good agreement among GNSS derived wave direction, INS derived wave direction and the modelled wind stress direction—noting the low temporal and spatial resolution of the model product. In particular, the GNSS and INS counterparts agree well despite not presenting the “true” wave direction given the impact of the tether. The current, as a low-frequency component, is constantly affecting the direction of the buoy motion, while the tether may inhibit the buoy motion causing waves to pass over or through the buoy rather than solely displacing it.

As for the modelled wind stress, due to its low temporal-spatial resolution, the use of it to decode any bias induced in the SSH at centimetre level is likely to be unreliable and will possibly introduce more uncertainty. On the other hand, while the waves are responsible for most buoy rotational and vertical motion, the observed current appears to be the dominant driver of buoy location about the central anchor (Figure S1). As a result, the current will often determine the heading of the buoy route in the water rather than its dynamics.

5.4. Relation between Dynamic Indicators and B–M Residuals

Our analysis of B–M residuals from four deployments offers insufficient evidence to support any bias in B–M series as a function of any of the observed or inferred quantities that reflect ocean conditions. In situ derived wave amplitudes seem to be widely scattered due to limited sample size, causing an uncertain estimation of linear models. The modelled wind stress seemingly shows a good linear fit; however, it is not a strong indicator given the apparent relationship is likely induced by the low temporal resolution and small sample size of the model input. The inconclusive relationship between B–M residual and wave quantities for the four deployments considered highlights that for the bias in B–M series to form spurious signals as in our previous study [17], the specific alignment of the external forcings on the tether is likely to be vital, since it will project its contribution onto the tether more consistently over time and in turn generate signal rather than noise.

Nevertheless, we do observe a non-linear relationship between buoy tilt and the magnitude of the waves derived from GNSS along the full-time span of the deployment (Figure 8a). From this, it can be learnt that buoy tilt series and the GNSS derived wave directions hold the potential to describe the dynamic performance of the buoy and its surrounding state of the sea. Links between the two are likely to be good indicators of the possible biases in the buoyancy position. However, (non)linear fitting is equivocal suggesting that further investigation is required with increased sample size. Meanwhile, it is worth noting that, without the assistance of indoor simulation with strictly controlled conditions (e.g., a towing tank experiment), isolation of external forcings from each other is extremely difficult to achieve in an actual buoy deployment. The impact of this is evident in the analysis of the second scenario (moderate current, low wave state—Figure 9a), where an unreliable linear model suggests a seemingly good fit to the observations, when in fact, large uncertainty can be observed towards higher current velocities. For example, in Figure 9a, B–M residual is observed in a range from 0 to 6 cm when the projected current velocity is at a moderate level of 0.2 m/s. This uncertainty is likely due to the coincidental alignment of the current with other forcings compounded with processing errors from multiple sources: GNSS buoy, Bottom Pressure (BP) gauge and the current meter.

So far, and within the limited scope of this study, we have considered the mooring SSH as the ground truth ignoring all potential error contribution from it towards the SSH residual. In reality, the mooring SSH derived from the BP gauge includes multiple error sources. For example, the BP sensor itself suffers from bias and drifting errors. A longer-term study is on-going to investigate these errors. Additionally, during the process of deriving mooring SSH, we used modelled atmospheric pressure to remove the atmospheric component from the BP record. Errors in the model could potentially be introduced into the solutions at sub-centimetre level (1 hPa = 1 cm under an inverse barometer assumption). As a result, we need to keep in mind that the uncertainty we observe when exploring the relation between buoy dynamics and the B–M residuals may also be generated from the mooring, hence the statistics we derive for the buoy system represent an upper bound for the errors involved.

5.5. Implication for High-Resolution Altimetry Validation

Advancing altimeters require improved in situ validation approaches. By way of an example, SWOT, the swath-based altimetry mission, requires: (a) higher precision in validation approaches; (b) extension of the traditional point-based approaches (e.g., buoy, mooring, and TG) towards a higher spatial scale; and (c) better understanding of the intra-swath variability of the ocean-atmosphere interactions. The progress of these three aspects is of equal importance, while in this study we focus solely on the precision aspect of the validation approach. Our ongoing development of GNSS/INS equipped buoys has delivered a number of enhancements, namely the integration of INS [17] and the addition and enhanced handling of GLONASS in GNSS processing. GNSS/INS buoy errors remain at the 1.5-cm level as assessed against in situ mooring. Notwithstanding likely error contributions from the mooring sensors, this requires further development.

By design, SWOT has a validation requirement at the level of 1.2 cm [11], aiming to resolve ocean processes up to a wavelength of 15 km. The performance of the buoy-mooring system needs to be better understood in order to compensate the remaining discrepancy between the precision of our buoy-mooring validation system and the mission requirement. First, while the buoy may not be able to validate the 7.5 km averaged observation from the SWOT mission, it could be useful for 15 km averaged samples or higher, despite some further understanding of the spatial scale of its errors in need. Second, a set of Current, Waves, Pressure Inverted Echo Sounder (CWPIES) is being prepared for future in situ comparison at Bass Strait validation facility. This could provide more insights into the error contribution from the BP gauge towards the overall 1.5-cm precision level. Finally, potential improvement can still be achieved for GNSS processing, such as adding a spatial constraint on the estimation of tropospheric delay in the buoy network, advanced ambiguity resolution method resolving inter-system/inter-frequency biases for multi-GNSS constellations (e.g., GLONASS, Galileo and the BeiDou Navigation System) and further improvement to the undifferenced multi-GNSS PPP approach.

6. Conclusions

We report on further developments made to GNSS/INS equipped buoys in response to increasingly stringent validation requirements for satellite altimetry. Using deployments at the Bass Strait validation facility, the benefit of introducing GLONASS into the processing workflow and a hybrid ambiguity fixing approach (i.e., integer GPS ambiguities + float GLONASS ones) are assessed via comparison between GPS/INS and GNSS/INS SSH solutions, and against our in situ mooring. Up to ~3 cm RMS reduction can be achieved for the B–M residual over the selected sidereal periods with the addition of GLONASS. Over the full deployment timeseries, a decrease in the standard deviation of the solutions ranging from 0.1 to 0.5 cm is commonly observed, highlighting a maximum ~20% reduction in the RMS, although in some cases the RMS increases by a small amount (up to 0.1 cm). The spectrum analysis performed indicates a drop in noise level across the frequency bands, indicating the solution benefits from the change to the constellation repeat period.

A comparison between DD and PPP processing methods was carried out using GPS-only solutions. Results show that both methods reach equivalent precision around ~1.5 cm when compared against our SSH series from an in situ mooring. The prevalence of outlier periods suggests that PPP is more dependent on the data quality in order to retain filter convergence, whereas the DD processing suffers from reference station multipath and inter-frequency biases when receivers from different manufacturers are used. Another important finding suggests that the B–M residual has some common signals independent of the processing methods, which are from two possible sources: one being the buoy dynamics as a result of external forcings, the other being the systematic errors in the mooring sensors. In this regard, the B–M statistical results in this study are therefore considered an upper bound to the buoy system error since they include a contribution from uncertainty in the mooring.

Evidence of temporal correlation is found among wave magnitudes, buoy tilts and accelerations as expected. Derived wave directions from GNSS and INS show good agreement and are also generally consistent with modelled wind stress from the ACCESS-G model (albeit at very different temporal resolutions). During selected periods when multiple forcings were isolated as much as possible from each other, no significant quantitative relation was found between them and B–M biases over four deployments, confirming the performance of the buoy design. However, more data is required to further understand the buoy dynamics, especially in higher sea states. Possible contribution from a hydrodynamic model of the buoy and experimentation in a towing tank could be beneficial since it allows possible isolation of external forcings from the environment.

The ongoing buoy development in Bass Strait seek to improve the SSH solutions continuously via addressing potential biases within the system. As part of the preparation for the soon to be launched SWOT mission, the buoy presented here will be scaled up to facilitate sustained operation in higher sea states. A buoy array is being developed and will enable better understanding of intra-swath variability of the ocean-atmosphere conditions in proper spatial scales for future validation activities.