1 Introduction and motivation

Fitting gravimetric geoid or quasigeoid models to GPS-levelling data is a technique now used in many countries. In Australia, this is necessary because of distortions in the Australian Height Datum (AHD) (e.g. Roelse et al. 1975; Morgan 1992; Featherstone 2004, 2006; Featherstone and Filmer 2008; Filmer and Featherstone 2009, 2012). This paper describes the fitting of the Australian Gravimetric Quasigeoid 2017 model (AGQG 2017; Featherstone et al. 2018a) to a nationwide GPS-levelling dataset (Featherstone et al. 2018b) to provide a model of the separation between the GDA2020 ellipsoid and the AHD, thus enabling a direct transformation between ellipsoidal and AHD heights (cf. Featherstone 2000). We term this a combined gravimetric–geometric model, though others may term it a “hybrid quasigeoid model”.

Combined gravimetric–geometric models are affected by errors in one, more or all of the GPS ellipsoidal heights, approximations in the computation of heights, errors in the levelling data and systematic errors in the geoid (e.g. Kotsakis and Sideris 1999). We briefly describe the computation of the geometric component and its accompanying grid of uncertainty values using least squares prediction (LSP; cf. Mysen 2014). We have developed a new AUSGeoid model, AUSGeoid2020, a combined gravimetric–geometric model accompanied by a 1′ by 1′ grid of uncertainty values propagated from all the input data into the LSP-gridded surface, providing users with location-specific uncertainty estimates.

There has been a change to the Australian geodetic datum, from GDA94 to GDA2020. GDA2020 is based on the International Terrestrial Reference Frame 2014 (ITRF2014; Altamimi et al. 2016) extrapolated to epoch 2020.0 using Australian GPS station velocities (ICSM 2018), under the assumption that the motion is linear. This datum change causes horizontal geodetic coordinates to move ~ 1.8 m in a north-easterly direction, primarily due to plate tectonic motion between 1994 and 2020, and ellipsoidal heights to decrease by ~ 0.090 m, due to improvement and refinement of the International Terrestrial Reference System between ITRF92 and ITRF2014. These changes have made it necessary to recompute AUSGeoid because its predecessor was fitted to AHD and GDA94 using LSP (Brown et al. 2011).

The location-specific uncertainty grid allows for the estimation of uncertainty from GPS-derived AHD heights when used in an absolute sense (i.e. using PPP or AUSPOS), but when GPS baselines are used to determine AHD height differences, correlated errors cancel (Kearsley 1986). This is because nearby AUSGeoid2020 model values are derived from similar gravity/GPS-levelling data. By subtracting them from one another to obtain height differences, common errors cancel. The closer two points are (horizontally), the more common data they are based on, and the more errors cancel (Fig. 1). We have therefore devised a decorrelation function (Sect. 4) that can be used to compute the relative uncertainty between GPS-AUSGeoid2020-derived AHD heights, similar to allowable misclosures used in differential levelling.

Fig. 1
figure 1

AUSGeoid2020 values (dots at centre of circles) are computed using data within a given radius (circles). The closer the dots are, the greater the amount of common data in the overlapping circles and the therefore the greater the correlated errors

Full size image

2 Fitting AGQG2017 to the AHD

Least squares prediction (LSP: Moritz 1980) was used to model the geometric component of AUSGeoid2020 from 7624 co-located GPS-levelling points across Australia; data were described and mapped in Featherstone et al. (2018b). After performing an initial cross-validation test (Sect. 3.1), 40 points were found to misfit AUSGeoid2020 by more than five standard deviations from the mean, so were removed and the geometric component recomputed.

A Gaussian empirical covariance function for the geometric component was computed using GEOCOL’s (Tscherning et al. 1992) empcov.f, with the covariance determined from GPS-levelling data aggregated over 5′ bin sizes (Fig. 2). The values determined by least squares for the Gaussian analytical covariance function are: variance 0.0092 m2, correlation length 75 km and constant 0.78. Beforehand, a tilted plane was removed from the GPS-levelling data using the remove linear function in x and y command in GEOCOL’s geogrid.f to remove the dominant north–south tilt in the AHD (cf. Featherstone and Filmer 2012). This tilt was added back after LSP to deliver the whole geometric component (cf. Fig. 3).

Fig. 2
figure 2

Gaussian analytical model fitted to the empirical covariance function computed for 5′ bins

Full size image
Fig. 3
figure 3

Geometric component of AUSGeoid2020. Units in metres

Full size image

AUSGeoid2020 (Fig. 4a) is the sum of the geometric (Fig. 3) and AGQG2017 1′ × 1′ grids (Fig. 12a in Featherstone et al. 2018a). The accompanying 1′ × 1′ uncertainty grid (Fig. 4b) was computed by SLOPOV-propagating that of the geometric component and AGQG2017 (Fig. 12a in Featherstone et al. 2018a), assuming independence. The uncertainty in AUSGeoid2020 is dominated by the uncertainty in the AHD heights, which results in a larger error estimate in AUSGeoid2020 than AGQG2017. We acknowledge there are neglected correlations amongst the data, e.g., because AHD heights have been used in the reduction of gravity observations in AGQG2017, and there is a covariance between the GPS-levelling error values because the levelling data are interconnected. This will be discussed further in Sect. 3.1.

Fig. 4
figure 4

a The AUSGeoid2020 combined gravimetric–geometric model and b its associated location-specific error model. Units in metres

Full size image

The public-domain AUSGeoid2020 model (http://www.ga.gov.au/scientific-topics/positioning-navigation/geodesy/ahdgm/ausgeoid2020) has been trimmed (values set to NaN) to a distance of ~ 50 km beyond the Australian coastline and its territorial islands. It was not trimmed exactly to the coastlines because this would prevent values being bi-cubically interpolated from the grid in near-coastal regions. This is administrative rather than scientific, but done deliberately so as to remind users that AUSGeoid2020 it is a model of the AHD (an onshore-only vertical datum) rather than the classical geoid or quasigeoid. For offshore users, the public-domain AGQG2017 is also available at the above website.

As described in “Appendix B”, post-publication inspection of AGQG2017 around the coastline revealed a spurious feature in Port Phillip Bay (near Melbourne in SE Australia), which introduced an unexplained localised bias between AGQG2017 and GPS-AHD data on land (Fig. 7a). We have patched this region using alternative altimetry data (“Appendix B”).

In regions where no GPS-levelling data are within 75 km of each other, the uncertainty of the geometric component (Fig. 3) is at its largest, 96 mm, which is the standard deviation (one sigma) of the GPS-levelling data. Given the highly correlated nature of the GPS-levelling data at short wavelengths beyond approximately 25 km (Fig. 2), AUSGeoid2020 thus transitions to AGQG2017 with a trend (cf. Fig. 3).

3 Testing AUSGeoid2020

3.1 Cross-validation

Cross-validation (aka “leave-one out”) in LSP was performed using GEOCOL’s geogrid.f program (Tscherning et al. 1992) as per our previous studies (Featherstone and Sproule 2006; Brown et al. 2011). The method is to remove one point from each LSP run and then use LSP to predict the GPS-levelling value at that removed point. The LSP value is then subtracted from actual GPS-levelling value to produce a residual. This process was repeated for all 7624 GPS-levelling points (after removal of the 40 outliers mentioned earlier), and the set of all residuals was used to assess the uncertainty of the combined model by looking at their distribution.

However, this cross-validation in LSP is not fully independent because the geometric component computed with each point removed is still based on AHD heights from connected levelling loops and ellipsoidal heights which may have been processed together with nearby points and also have correlated error. This can make the result of “leave-one out” cross-validation optimistic. The descriptive statics are max 228 mm, min 226 mm, mean − 2 mm and standard deviation (STD) ± 38 mm.

3.2 Pseudo-independent data

An additional 8409 GPS-AHD values were supplied by three states (Western Australia, New South Wales and Victoria). These data were not used in the determination of the geometric component (Sect. 2) because the GPS observation time was less than a recommended 6 h (ICSM 2017). However, the GPS heights still have an expected standard deviation of less than ± 20 mm (Jia et al. 2014), making them adequate for testing of the AUSGeoid2020 model given the formally propagated expected uncertainty in the model (Fig. 4b).

AUSGeoid2020 was bi-cubically interpolated to the locations of these pseudo-independent data to evaluate the differences. The descriptive statics are: max 249 mm, min 287 mm, mean − 7 mm and STD ± 27 mm. This lower STD, which is less than the formally propagated error in AUSGeoid2020, should, however, be treated sceptically. This is partly because these data are not fully independent in the strictest sense. The AHD heights (determined by levelling) are referenced to the Australian National Levelling Network (ANLN; Roelse et al. 1975), which means the same localised errors in the ANLN are common to both the data used to determine the geometric component and the pseudo-independent data. Another potential cause is the use of single covariance function to compute the geometric component across the continent. As a result, lower-quality data affect the variance of the data set and therefore falsely increase the uncertainty in regions with higher-quality data. Ideally, location-specific covariance functions could be developed (cf. Knudsen 2005), but this remains a topic of research for future models.

4 Decorrelation function to compute the relative uncertainty

When computing AHD height differences with GPS/GNSS baselines and AUSGeoid2020, correlated errors cancel (cf. Kearsley 1986). This fundamental principle behind relative GPS baseline processing is also applicable to AUSGeoid2020 because it is determined using similar gravity/GPS-AHD data (cf. Fig. 1). As such, the relative uncertainty of AHD heights derived from AUSGeoid2020 over GPS/GNSS baselines should be better than the absolute uncertainty, particularly for points closer together (cf. Fig. 1).

To assess the relative uncertainty of AUSGeoid2020 in this scenario commonly used by GPS/GNSS surveyors, height differences between the levelling data were compared to height differences derived using all the pseudo-independent GPS data (Sect. 3.2), minus bi-cubically interpolated AUSGeoid2020 values. The height difference residuals (i.e. difference between AHD and GPS-AUSGeoid2020) were grouped by the horizontal distance (calculated using Vincenty’s (1975) formula for geodesics) between points in bins of 1 km.

The STD of the height difference residuals increases as the geodesic distance between the points increases (Fig. 5). This demonstrates that relative heights determined using AUSGeoid2020 over shorter baselines are more accurate than over longer baselines, confirming that more of the correlated errors cancel. The approach used here was first presented in Smith and Roman (2001), Fig. 9). Our results show improvement, largely due to improvements in geoid modelling since circa 2001.

Fig. 5
figure 5

Comparison of the uncertainty in computing the difference in height from levelling and GPS using AUSGeoid2020

Full size image

From Fig. 5, over distances greater than approximately 3 km, the uncertainty of relative height differences computed from GPS observations and AUSGeoid2020 is less than the allowable misclose of traditional third-order levelling in Australia, which is equivalent to 12 mm of misclose per square root km. At distances less than ~ 3 km or when compared with more accurate levelling techniques over longer distances (e.g. allowable misclose of 2 mm or 4 mm per square root km of levelling), height transfer by levelling is superior.

The variance associated with relative AHD heights, determined by taking the difference pairs of AUSGeoid2020 and GPS ellipsoidal height observations \( \left( {h_{1} - \zeta_{1} ) - (h_{2} - \zeta_{2} } \right) \), is given by

$$ \begin{aligned} \sigma^{2} \left( {h_{1} - \zeta_{1} - h_{2} + \zeta_{2} } \right)&= \sigma^{2} \left( {h_{1} } \right) + \sigma^{2} \left( {\zeta_{1} } \right) + \sigma^{2} \left( {h_{2} } \right) + \sigma^{2} \left( {\zeta_{2} } \right) \\ &\quad\quad- 2\text{cov} \left( {\zeta_{1} ,\zeta_{2} } \right) - 2\text{cov} \left( {h_{1} ,h_{2} } \right) \\ & \quad\quad - 2\text{cov} \left( {h_{1} ,\zeta_{1} } \right) + 2\text{cov} \left( {h_{1} ,\zeta_{2} } \right) \\ &\quad\quad+ 2\text{cov} \left( {h_{2} ,\zeta_{1} } \right) - 2\text{cov} (h_{2} ,\zeta_{2} ) \\ \end{aligned} $$
(1)

The cross-covariance terms in Eq. (1) (i.e. terms involving a GPS \( h \) and AUSGeoid2020 \( \zeta \) values) are zero since they are independent; this leads to a simplified expression

$$\begin{aligned} &\sigma^{2} \left( {h_{1} - \zeta_{1} - h_{2} + \zeta_{2} } \right)\\ &\quad= \sigma^{2} \left( {h_{1} } \right) + \sigma^{2} \left( {\zeta_{1} } \right) + \sigma^{2} \left( {h_{2} } \right) + \sigma^{2} \left( {\zeta_{2} } \right)\\&\qquad - 2\text{cov} \left( {\zeta_{1} ,\zeta_{2} } \right) - 2\text{cov} \left( {h_{1} ,h_{2} } \right)\end{aligned} $$
(2)

We have estimated the shape of the covariance term \( 2cov\left( {\zeta_{1} ,\zeta_{2} } \right) \) by considering the variances of the 1-km-binned relative height difference residuals. The empirical binned variance values can be modelled with a semi-variance exponential function (Eq. 9; Fig. 6) of the form

$$ S = \sigma^{2} \times \left( {1 - \rho e^{{ - \frac{3l}{\alpha }}} } \right) $$
(3)
Fig. 6
figure 6

Variance of the 1-km-binned height difference residuals (blue circles) and a fitted exponential semi-variance function (red line)

Full size image

In Eq. (3), \( \alpha \) is the distance parameter that represents the point at which the correlation in the errors is no longer statistically significant and \( \sigma^{2} \) is the far-field variance, often referred to as the “sill”. The value \( \rho < 1 \) prevents the function from being too small for very small values of \( L \). We fitted these parameters to the binned relative height difference residual variances by least squares. The fitted values are \( \alpha = 63.151 \) km and \( \rho = 0.68 \). We thus propose that the covariance term, \( 2cov\left( {\zeta_{1} ,\zeta_{2} } \right) \), can be approximated using this fitted function such that

$$ 2cov\left( {\zeta_{1} ,\zeta_{2} } \right) = \left( {\sigma^{2} \left( {\zeta_{1} } \right) + \sigma^{2} \left( {\zeta_{2} } \right)} \right) \rho e^{{ - \frac{3l}{a}}} . $$
(4)

Therefore, the standard deviation of the uncertainty associated with AHD height differences determined using relative GPS and AUSGeoid2020 can be approximated by

$$ \sigma \left( {h_{1} - \zeta_{1} - h_{2} + \zeta_{2} } \right) = \sqrt {\sigma^{2} \left( {h_{1} } \right) + \sigma^{2} \left( {h_{2} } \right) - 2\text{cov} \left( {h_{1} ,h_{2} } \right) + \left( {\sigma^{2} \left( {\zeta_{1} } \right) + \sigma^{2} \left( {\zeta_{2} } \right)} \right)(1 - \rho e^{{ - \frac{3l}{a}}} )} $$
(5)

In Eq. (5), \( \sigma^{2} \left( {\zeta_{1} } \right) \) and \( \sigma^{2} \left( {\zeta_{2} } \right) \) can be determined by interpolating and squaring values from the AUSGeoid2020 error grid (http://www.ga.gov.au/scientific-topics/positioning-navigation/geodesy/ahdgm/ausgeoid2020) and \( \sigma^{2} \left( {h_{1} } \right) \), \( \sigma^{2} \left( {h_{2} } \right) \) and \( cov\left( {h_{1} ,h_{2} } \right) \) should be extracted from GPS processing software. “Appendix A” provides a worked example of the application of Eq. (5) to the outputs of some GPS processing software.

5 Conclusion

AUSGeoid2020 provides the method to convert between GPS-derived ellipsoidal heights and AHD heights (and vice versa) and to compute AHD height differences. AUSGeoid2020 is also accompanied by a location-specific combined uncertainty of the GPS, AHD and quasigeoid components used in its construction. Pseudo-independent testing of AUSGeoid2020 has shown that it is capable of assisting users of GPS to convert from ellipsoidal heights to AHD heights with an absolute uncertainty of ± 27 mm. We also presented a method to compute the relative uncertainty between AUSGeoid2020-derived heights based on the principle of correlated errors cancelling between two points. Our analysis implies AUSGeoid2020 is equivalent to, or better than, traditional third-order levelling in Australia at distances beyond ~ 3 km.