AUSGeoid2020 combined gravimetric–geometric model: location-specific uncertainties and baseline-length-dependent error decorrelation

AUSGeoid2020 is a combined gravimetric–geometric model (sometimes called a “hybrid quasigeoid model”) that provides the separation between the Geocentric Datum of Australia 2020 (GDA2020) ellipsoid and Australia’s national vertical datum, the Australian Height Datum (AHD). This model is also provided with a location-specific uncertainty propagated from a combination of the levelling, GPS ellipsoidal height and gravimetric quasigeoid data errors via least squares prediction. We present a method for computing the relative uncertainty (i.e. uncertainty of the height between any two points) between AUSGeoid2020-derived AHD heights based on the principle of correlated errors cancelling when used over baselines. Results demonstrate AUSGeoid2020 is more accurate than traditional third-order levelling in Australia at distances beyond 3 km, which is 12 mm of allowable misclosure per square root km of levelling. As part of the above work, we identified an error in the gravimetric quasigeoid in Port Phillip Bay (near Melbourne in SE Australia) coming from altimeter-derived gravity anomalies. This error was patched using alternative altimetry data.


Introduction and motivation
Fitting gravimetric geoid or quasigeoid models to GPSlevelling data is a technique now used in many countries. In Australia, this is necessary because of distortions in the B N. J. Brown Nicholas.Brown@ga.gov.au 1 Australian Height Datum (AHD) (e.g. Roelse et al. 1975;Morgan 1992;Featherstone 2004Featherstone , 2006Featherstone and Filmer 2008;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 prop- Fig. 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 agated 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.

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 m 2 , 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). 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. The public-domain AUSGeoid2020 model (http://www.g a.gov.au/scientific-topics/positioning-navigation/geodesy/a hdgm/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).

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 Fig. 4 a The AUSGeoid2020 combined gravimetric-geometric model and b its associated location-specific error model. Units in metres 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.

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 AUS-Geoid2020, 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.

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 AUS-Geoid2020 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 AUS-Geoid2020 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.
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 (h 1 − ζ 1 ) − (h 2 − ζ 2 ), is given by − 2cov(h 1 , ζ 1 ) + 2cov(h 1 , ζ 2 ) + 2cov(h 2 , ζ 1 ) − 2cov(h 2 , ζ 2 ) (1) Fig. 6 Variance of the 1-km-binned height difference residuals (blue circles) and a fitted exponential semi-variance function (red line) The cross-covariance terms in Eq. (1) (i.e. terms involving a GPS h and AUSGeoid2020 ζ values) are zero since they are independent; this leads to a simplified expression We have estimated the shape of the covariance term 2cov(ζ 1 , ζ 2 ) 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 In Eq. (3), α is the distance parameter that represents the point at which the correlation in the errors is no longer statistically significant and σ 2 is the far-field variance, often referred to as the "sill". The value ρ < 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 α 63.151 km and ρ 0.68. We thus propose that the covariance term, 2cov(ζ 1 , ζ 2 ), can be approximated using this fitted function such that Therefore, the standard deviation of the uncertainty associated with AHD height differences determined using relative GPS and AUSGeoid2020 can be approximated by In Eq. (5), σ 2 (ζ 1 ) and σ 2 (ζ 2 ) can be determined by interpolating and squaring values from the AUSGeoid2020 error grid (http://www.ga.gov.au/scientific-topics/positioning-nav igation/geodesy/ahdgm/ausgeoid2020) and σ 2 (h 1 ), σ 2 (h 2 ) and cov(h 1 , h 2 ) 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.

Conclusion
AUSGeoid2020 provides the method to convert between GPS-derived ellipsoidal heights and AHD heights (and vice versa) and to compute AHD height differences. AUS-Geoid2020 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. in the AHD. However, the data used to produce Fig. 5 are corrected for correlated AHD errors since they are residual height differences. This results in a smaller σ value in Eq. (3) for the residuals, but retains the general shape of the spatial dependence of the covariance for the function fitting.

Appendix B: Patch correction to AGQG2017
Post-publication inspection of AGQG2017 around the coastline revealed a spurious feature in Port Phillip Bay (near Melbourne in SE Australia), which introduced a localised bias between AGQG2017 and GPS-AHD data on land (Fig. 7a). Upon inspection of the data used to compute AGQG2017, the feature is attributable to a large positive The patch to the AGQG2017 model has reduced the residuals between GPS-levelling and AUSGeoid2020 in the region close to the patch. a Before patch was applied, b after patch was applied gravity anomaly in the Sandwell et al. (2014) satellitealtimetry-derived gravity anomaly grid (Fig. 7b). The same gravity anomaly was not evident in the DNSC08, DTU10 and DTU15 Knudsen 2009, 2016;Andersen et al. 2010 satellite altimetry gravity anomaly grids. We believe this is not a true gravity anomaly, but an error in the Sandwell model, likely due to the known problems computing gravity anomalies from satellite altimetry data in coastal regions. We computed a new quasigeoid model using DTU15 gravity data, in place of the Sandwell et al. (2014) model, over a 5°× 5°patch centred over the region containing the anomalous feature in AGQG2017. This quasigeoid computation used the same integration parameters as AGQG2017 (i.e. modification degree of 40 and integration cap radius of 0.5°) and EGM2008. We then blended this 5°× 5°degree patch with the AGQG2017 model by linearly combining the two grids using cosine tapered weights to ensure there were no discontinuities. DTU15 does not have an accompanying error grid, so the AGQG2017 error model was left unchanged in this region. Figure 7c shows the patched quasigeoid, and Fig. 7d shows the gravity anomaly data where DTU15 is used offshore.
As shown in Fig. 8, the patch resolved the localised bias between AUSGeoid2020 and GPS-AHD data. The residuals beyond the limits of the patch have not changed, while the residuals near the patched region have reduced. At the point closest to the patch at approximately − 38°00 , 144°30 , the residual has been reduced from 50 to 0 mm.