ON THE COMPUTATION OF A PRECISE GEOID – TO – QUASIGEOID SEPARATION

In geodesy, orthometric and normal heights are considered as basic height systems on the earth. The reference surfaces for these heights are the geoid and quasigeoid respectively. Taking advantage of GNSS measurements, one can achieve a precise solution for the geoid and for the quasigeoid. Two methods, called direct and indirect, are worked out in this research for the computation of separation between geoid and quasigeoid in a mountainous region in the USA. The area selected for this purpose is mountainous and rough enough in order to be able to show the effect of roughness of topography in the sought quantity. The results of the two methods and testing them against GNSS-Levelling on 445 known points indicates an accuracy of 1.3 cm in RMS scale with the direct method, where there is 7 cm as an average difference between the observed geoid and quasigeoid separation and the same quantity derived from the direct method. Using Chi-squared goodness of fit test showed that the distribution of the residual quantities are normally distributed in the test area.


INTRODUCTION
Two concepts of Stokes-Helmert (Helmert, 1890) and Molodensky (Molodensky et al. 1960) approaches are used to solve geodetic boundary-value problem for geoid and quasigeoid heights both referred to the normal ellipsoid.The geoid and quasigeoid are in turn reference surfaces to measure orthometric and normal heights on the earth respectively.The separation between geoid and quasigeoid is usually needed, in practice, for transformation of one system of height to another.There is, however, an approximate formula in literature (Heiskanen and Moritz, 1967) to compute Geoid-to-Quasigeoid Separation (GQS) using Bouguer gravity anomaly; and it is also based on the known orthometric height at a point of interest.Precise estimation of Bouguer gravity anomaly requires more detailed topographical height information around the point.So far, many studies have been conducted for the determination of GQS.Sjoberg (1995) uses a model of orthometric height of higher order precision in the approximate formula (Heiskanen and Moritz, 1967) in order to reduce the terrain effect uncertainties.Rapp (1997) used precise height anomaly into the approximate formula.Nahavandchi (2002) used the Rapp technique as an indirect method for numerical evaluation on Iran's region GNSS-Levelling data.Sjoberg (2006) offered a precise formula, including terrain correction term and a term with regard to lateral variation of topographical densities for conversion from normal height to orthometric height.Tenzer et al. (2006) estimated the GQS by a correction model based on the mean gravity disturbance along the vertical within the topographic masses.Flury and Rummel (2009) derived the magnitude of the GQS equal to 24 and 48 cm in two test areas of the Alps.Sjoberg (2010) presented a strict formula for the GQS with two terms, terrain and gravimetric corrections.The result of his study improved the accuracy of Flury and Rummel formula up to 1 cm.Sjoberg (2012) used an arbitrary gravity correction model to improve the rough topographic effects accuracy on the GQS.Sjoberg and Bagherbandi (2012) evaluated the magnitude of this separation using the EGM08 and DTM2006 models complete to harmonic degree and order 2160 in the Tibet plateau and Indian Ocean to 5.47 and 0.11 m, respectively.Bagherbandi and Tenzer (2013) estimated the value of the GQS by GOCO02S model complete to degree 250 on the Tibet plateau and the Himalayan Mountains ranging from 0.15 to 3.62 m.The comparison of their results with EGM08 model shows ±20 cm differences.Sjoberg (2015) presented a rigorous scheme to estimate GQS using gravity disturbance expansion to Taylor series along plumbline.In his study, he used three different correction methods, applying the Bouguer gravity disturbance, an arbitrary compensation model and analytical continuation technique.
In this study two methods called direct and indirect are worked out to accurate estimation of GQS.Topographic masses above the geoid are a major obstacle in the gravity field determination.Hence, in two mentioned methods, different models for topographic effects are used.Finally, to show the precision superiority, in each one of two methods, the magnitude of geoidal height extracted from methods are compared with geometric geoidal height obtained GNSS-Levelling stations in test area.

DIRECT METHOD IN DETERMINATION OF GEOID TO QUASIGEOID SEPARATION
The orthometric and normal heights H and   of a point on the earth are derived to be, (Heiskanen and Moritz, 1967), (1) where, C is the geopotential number, ̅ and ̅ are the mean actual and mean normal gravities along the actual plumb line from the geoid up to the point of interest and normal plumb line from the normal reference ellipsoid up to the telluroid respectively.The separation between the geoid and the quasigeoid (GQS) as the difference between normal and orthometric heights is given by, (ibid), where, Zetta, is the height anomaly and N is the geoidal height at the computation point.The evaluation of ̅ alone in practice is problematic since the information required about density distribution inside the earth around the point of interest is less reliable.The suggested solution is to approximate the differential ̅ − ̅ term with the Bouguer gravity anomaly Δ   at a computation point P. Therefore, the estimation of GQS is given by the formula (ibid).

𝑁 − 𝜁 ≈
Δ    ̅  (4) The formula above does not provide the required accuracy in mountainous areas due to the uncertainty in terrain correction (Helmert, 1890;Niethammer, 1932).Using Bruns formula, the disturbing potential is converted to geoid height and/or to height anomaly through normal gravity (Bruns, 1878;Molodensky et al., 1960).
where subscripts  and  denote locations on the Earth and geoid surfaces respectively,  0 and   represent the normal gravities on the reference ellipsoid with geocentric radius  0 and on the telluroid surface with   radius, respectively.Radius   equals the sum of  0 radius and normal height   .In Eq. ( 5) we can represent the disturbing potential T in form 1 Analytical continuation bias of spherical harmonic series as a (, Ω) point according to the following expansion that is complete to M harmonic degree.
where   and   are fully normalized spherical harmonics and spherical harmonics coefficients of disturbing potential of degree  and order , respectively,  is the maximum degree of expansion, Ω = (, ) the spatial angle at point of spherical coordinates co-latitude  and longitude .Considering Eq. ( 5), the expression of a strict model for GQS with the terrain correction term as following formula is given (Sjoberg, 2006).
where the ACB 1 is the topographic bias.This bias is evaluated applying external series expansion of earth gravity potential inside the topography.The following formula estimates the magnitude of the ACB to fourth powers of elevation according to spherical approximation and a constant density for topographic masses (Agren, 2004;Sjoberg, 2007).
where  = ,  is gravitational constant,  being the density of crust and for topographic height powers with its coefficients the following expansions exists = ∑    ,   (10) In Eq. ( 7), by using downward continuation technique and Taylor series, the following conversion between the telluroid and reference ellipsoid surfaces is expressed for normal gravity (Molodensky et al., 1960).
where by applying following spherical approximation and in Eq. ( 11) regardless of higher terms  = 1 , we have where  equal to In above Eq. is correction that is used for transferring normal gravity between telluroid and reference ellipsoid surfaces.By applying Eq. ( 13) in Eq. ( 7), final formula for the GQS computation is presented in the following form In Eq. ( 12) by assuming that   ~5000 m and R as the mean Earth radius approximating the mean sea level, correction quantity  for term of  = 1 equals to 0.9 and for the term of  = 2 this quantity equals to 1. Ergo, regardless of terms higher than  = 1, the amount of the GQS in Eq. ( 15) will not be affected substantially.On the other hand, in Eq. ( 12) the maximum error created by applying spherical approximation will be less than 1 mm (Sjoberg, 2004).Moreover, created error in Bruns formula is close to the 1.5 mm.This is due to the omission of the second and higher terms in the expansion of Taylor series.

INDIRECT METHOD IN DETERMINATION OF GEOID TO QUASIGEOID SEPARATION
The other model for the estimation of GQS is the Eq. ( 4) with a correction topographic term of the second power of orthometric height added, (Sjoberg, 1995).
where the vertical gradient of free air gravity anomaly equals, (Heiskanen and Moritz, 1967), In above Eq., ∆  is free air gravity anomaly,  is unit sphere and  is the spatial distance between the computation point P and integration point.Rapp (1997) used the Taylor series expansion formula for the height anomaly on the Earth's surface as for the evaluation of GQS, he suggested the following equation in 1997, where and  0 is the value of the height anomaly at the point with ellipsoidal radius   and it is equal to, (Heiskanen and Moritz, 1967), where  is semi major axis of the ellipsoid and GM is geocentric gravitational constant.By using above expansion and applying spherical approximation  ℎ ⁄ = −2  3 ⁄ , the  1 correction term will be calculated through two following expansions.
In the first section of  2 correction term, by assuming the constant crust density and applying the Bouguer reduction δ  = −2 to free air gravity anomaly ∆  , the Bouguer gravity anomaly ∆  from and by applying spherical approximation, the mean normal gravity  from is computed, (Heiskanen and Moritz, 1967).
Based on existing definition, in the second section of  2 correction term, the vertical gradient of the free air gravity anomaly equals to Eq. ( 17).By applying the planer approximation, we can express Eq. ( 17) in the following form (Heiskanen and Moritz, 1967;Bian, 1997).
where ∆ 0  illustrates the value of free air gravity anomaly at the computation point and [, ] are horizontal parameters for planer coordinates of running point.Furthermore, in above Eq.2∆   ⁄ term is disregarded in the light of its being little amount that demonstrates the indirect effect of free air gravity anomaly.
In this situation, to solve the integration of Eq. ( 27), the Newton-Cotes formula is used.This formula in innermost area of computation points would have a solution in the following form (Sadiq et al., 2010) where for  1 , … ,  5 we have In above Eq.,  is the spacing of the data on a grid with −2 ≤  ≤ 2 and −2 ≤  ≤ 2 for integration.This integration is performed for  = 4 interpolated nodal points.

NUMERICAL EVALUATION
The results of this section is obtained from a numerical calculation in a mountainous region of the USA located in an area with latitudes 37 °N to 41 °N and longitudes 104 °W to 109 °W.This region includes the Rocky Mountains with roughness topographic to 3776 m high.The harmonic coefficients of the EGM08 model are completed to 2190 degree/order.A Digital Terrain Model (DTM) arranged in a regular grid data of 6 arc-min spacing is obtained from an average data of SRTM model with 3 second interval.Global height harmonic models of  2 ,  3 and  4 are used.The information of 445 GNSS-Levelling stations is available in the test area.The GRS80 normal gravity field, (Moritz, 1980), is used to generate the normal gravities.In direct method, the GQS is calculated from Eq. ( 15).In this equation, the disturbing potential is obtained on the Earth's surface   and geoid   considering height ℎ referred to the ellipsoid.The disturbing potential on the surfaces is computed by the expansion of spherical harmonic series in Eq. ( 6) and by using the EGM08 model.The Somigliana formula was used in order to calculate the value of normal gravity on the reference ellipsoid  0 (Somigliana, 1929). Quantity which describes the correction of the normal gravity between the reference ellipsoid and telluroid, was determined by Eq. ( 14).By accepting  = 2670 kg m −3 as the value of constant density for topographic masses, terrain effect as the analytical continuation bias was calculated from the expansion of Eq. ( 8) to maximum degree 2190.In this expansion, the second to fourth powers of topographic height were used.In Eq. ( 10) the multiple topographic powers are computable from height coefficients.The global height harmonic models make these coefficients accessible.The value of terrain effect on the GQS, obtained from various powers of topography has been illustrated in Table 1.
Table1.The statistical parameters of terrain effect on the GQS with topographic powers dissociation, unit: cm.The results show that the contribution of second power of topography is high and about 1.640 m.The contribution of the third power is in opposite sign.The fourth power effect is still considerable compared to the third power effect.It is shown that the even power effects compared to their earlier odd powers effects play important roles in modelling topographic masses potential.Figure (1-c) shows the total topographic effect on the GQS maximum to 2190 harmonic degrees.By determining right side of quantity in Eq. ( 15), separation between the geoid and quasigeoid is obtained.In Figure (1-d) the amount of this separation is illustrated.Table (2) represents the numerical results of direct method from separation determination between the geoid and quasigeoid.

𝐀𝐂𝐁 γ
In indirect method, the GQS values are computed from Eq. ( 16).The values of correction terms  1 and  2 are introduced to Eq. ( 19).The term  1 is the correction to transfer height anomaly  0 from the ellipsoid to the Earth surface (telluroid surface).The height anomaly  0 from Eq. ( 22) expansion and correction term  1 from Eq. ( 23) and Eq. ( 24) were computed using the EGM08 model.In term  1 by applying spherical approximation, this term  ℎ ⁄ was approximated to −2  3 ⁄ .Correction term  2 indicates the separation amount between the geoid and quasigeoid.To be more precise, this term is calculable through Eq. ( 21).In this Eq.The Bouguer gravity anomaly ∆  assuming a constant amount of crust density of the Earth and applying the Bouguer reduction δ  to free air gravity anomaly ∆  can be obtained.In Eq. ( 25) the free air gravity anomaly is expanded to spherical harmonics series and it is computed by using the EGM08 model.In Figure (2-a) the Bouguer gravity anomaly amount has been displayed.The mean normal gravity quantity ̅ is calculated by applying spherical approximation from Eq. ( 26).In the second section of correction term  2 , the vertical gradient of the free air gravity anomaly is computed from integral Eq. ( 17).This Eq. is converted to Eq. ( 27) by applying planer approximation.Needless to say, the free air anomaly indirect effect amount (2∆   ⁄ ) has been ignored on account of its low value.Eq. ( 27) was computed by using the Newton-Cotes integral solution for  = 4 interpolated points in horizontal location 2 ≤  ≤ 2 and−2 ≤  ≤ 2.Eq. ( 28) is a solution to the integral Eq. ( 27) for the vertical gradient of the free air gravity anomaly.In this Eq. 1 , … ,  5 quantities are achieved Eq. ( 29) by computing the free air gravity anomaly.In addition, the free air gravity anomalies in a regular grid of data with spacing of  = 6 ° are derived from the EGM08 model.shows the vertical gradient of the free air gravity anomaly in test area.By calculating right side of quantity in Eq. ( 16) the GQS was obtained in the indirect method.Figure (2-c) illustrates this separation.In Table (3) the numerical results of this method were shown.
The GNSS-Levelling stations as the test points of known GQS values are used to evaluate the accuracies of the direct and indirect methods.Hence, the geoidal heights were extracted from direct and indirect methods at the test points.Table (4) illustrates the results of geoidal heights.As the results show, the geoid modelled in direct and indirect methods has on average 16 cm difference.This difference demonstrates the models are comparable in centimeter accuracy level.In the performed comparison between the obtained geoid from two methods with the geometric geoid of GNSS-Levelling stations, the direct method with 1.3 cm precision superiority in respect to the indirect method at RMS scale is obvious.It has been considered as an efficient method.In addition, this is indicated as a precise approach in determining GQS.The stimulating result of direct method is due to the precise modelling of terrain effects.Comparison of the results of direct method in the test area with the GNSS-Leveling results in the area shows an average difference of 7 cm in GQS values.Also, the Chisquared goodness of fit test shows the normal distribution of the residuals in the test area.

CONCLUSIONS
In this study, two different methods (direct and indirect) were used to reach at a reliable method to determine the Geoid-to-Quasigeoid Separation (GQS).These two methods used the EGM08 geopotential model in a mountainous area.
The results of this study show that the GQS values in this two methods are different in average of approximately 3.3 cm.This is an indication of the fact that these two methods are comparable.It shows obviously that the direct method with the ability of modelling terrain effect up to high degrees of topographic height approximates more closely the GQS values in the test region.Comparison of the results of direct method in the test area with the GNSS-Leveling results in the area shows an average difference of 7 cm in GQS values.Also, the Chi-squared goodness of fit test shows the normal distribution of the residuals in the test area.

Figure 1 .
Figure 1.(a) the topographic of test area, unit: meter, (b) GNSS-Levelling stations location in test area, (c) the total topographic effect on the GQS, unit: meter and (d) separation between the geoid to quasigeoid as a result of the direct method, unit: meter.

Figure 2 .
Figure 2. (a) The Bouguer gravity anomaly with 6 arc-min interval, unit: mGal, (b) the vertical gradient of the free air gravity anomaly with 6 arc-min interval, unit: mGal and (c) the GQS obtained from indirect method, unit: m. .

Table 2 .
The statistical characteristics of GQS based on direct method, unit: meter.

Table 4 .
The statistical parameters of the comparison of the geoidal height in direct and indirect methods with geometric geoidal height of GNSS_Levelling stations, unit: meter.