1

Height Determination Using GPS Data, Local Geoid

and Global Geopotential Models

by

Khairul Anuar Abdullah

Panel for Industrial Measurement and Hydrography

Faculty of Surveying and Real Estate

Universiti Teknologi Malaysia

e-mail: Khairnl@fu.utm.my

Abstract

Orthometric heights are normally derived using the spirit levelling. This requires the spirit
level equipment to be set up from point 10 point along a levelling line which is a time
consuming and tedious task. GPS offers a new altematiV)O in orthometric height
determination very accurately over a comparatively short period. The ellipsoidal height
derived from GPS technique can be transformed into orthometric height if we know the
geoidal height normally derived from a gravimetric geoid of the area. Unfortunately, we
have yet to compute an accurate gravimetric geoid for such purpose, largely due to non­
existence of gravity data for a larger part of the Country. As an alternative, a study was
undertaken to look into the feasibility of using a simple local geoid solution and a global
geopotential geoid model solution The data used in this study consist of GPS data and
known orthometric height ofseveral GPS points. This paper present some of the results
obtained so far in estimating orthometric height from GPS data in local en,ironment.

1.0 L"lTROD(JCTION

The position of points derived from GPS measurements are lIsually computed in a three­

\ dimensional Cartesian coordinate system, and are then transformed into the more recognisable

I geodetic latitudes (<1», longitudes (A.), ellipsoidal heights (h) or in tenn of geodetic coordinate
differences A<j>, AA and Mi. GPS ellipsoidal heights are very useful for deformation and
I subsidence studies and other applications where the emphasis is Ilot so much ill 10catiJJg a precise
point in space as in the relative change of height from one time epoch to another. It is, however,
the case that the ellipsoidal heights delivered by GPS are not the same as those historically
obtained with spirit levelling (providing' orlhometric heights). Conventionally, topographic maps.
engineering design ODd construction project plans, usually depict relief by means ()f orthometric
height. Thus. the application of GPS will be f\!flher extended if acew:ate transformations between
GPS ellipsoidal height differences ODd the orthometric height differences can be realised. This can
be accomplished on the condition that we know the geoidal heighl, or rather, the geoidal height
difference which relates the orthometric height difference to the GPS ellipsoidal height difference.
Hence today, a great deal ofintcrest is being shown in the development ofthe geoid models which
are important to provide the necessary geoid height to transform GPS ellipsoidal heights to
orthometric heights.

2.0 GPS HEIGHTlNG

The basic results of the precise differential GPS survey of a baseline are the Cartesian coordinate
differences AX, AY ODd 62. Baselines connecting the observed GPS points are then put through a
network adjustment such as UD-HEIGHT (Khairul, 1994). The resnIting X, Y, and Z coordinates
of the GPS points are then transformed, using a reference ellipsoid, into geodetic coordinates in
terms of latitude (.p), longitude (A) and ellipsoidal height (h). The orthometric height (II) is related
to the ellipsoidal height (h) by the following relation:

H=h-N or, [I]

H= bnps - NMODEL [2]

where,
bnps is the GPS derived ellipsoidal height, and,
NMODEL is the geoidal height derived from a chosen geoid model,

or by the relative approach, the orthometric height difference between two GPS points may be
deduced from: I
Llli = tihaps • ANMODEL [3]
t
From the expressions above, the errors in H in eqn. [2] will depend upon the accuracy of the
parameters used in its evaluation. It is generally known that, the differences in h between two
points measured simultaneously by GPS are much more precise than h at either ofthe points. This
is because of the presence of systematic errors which, being significantly the same at the two
points, cancels in the difference. Similarly, ANMODEL, is much more precise than the geoid
height at either points. This means that for detennination requiring highest precision, the approach
implied in eqn. [3] is preferred to that in eqn. [2]. .

3.0 DERIVAnON OF GEOIDAL HEIGHT

The most precise method of obtaining accurate geoidal beight is using gravimetric observations.
Numerical integration of gra'imetric observations using Stokes integral equation provides us with
a local gravimetric geoid solution. Unfortunately, we have yet to compute our very own precise
gravimetric geoid. The reasons behind it are many but the most significant one is related to the
difficulties encountered in observing gravity values in large parts of the country due to the
topography. However, plans are underway to utilised modem approach in gravimetric
measurement such as using an airborn gravimeter. So when no local gravimetric geoid solution is
available, a local geoid surface fittiug solution may be employed as an alternative for small area of
gentle undulation, but for large area, a global geoid solution should be used instead.

3.1 Local Geoid Surface Fitting Solutions

This solution involves the use of a local geoid sutface model using a SUIface fitting procedure. The
fundamental theory of this solution is based on the following three assumptions:

(i) the adjusted GPS observations are of very high quality and considered to be exact,
(ii) the orthometric heights of at least three GPS stations in the network are known and also

considered to be exact,

(iii) the area involved is small and that the geoid features does not vary rapidly,

If we have three or more GPS points with known orthometric heights (referred here as Height
Control Point, (Hep), then a local geoid surface solution using a surface fitting model can be
employed. This is accomplished by taking Ni as a function of the position of each HCP in the
network. The sutface fitting model would take the following mathematical form:

z
 {4]

where;
ao, ai, a2, a3 are the unknown model coefficients,
xi = (4)i - 4>o)pmll!
Yi = ('-1 - Ao)vm cos <Pm/I!
4>m = (4)i + 4>0)/2
and;
<Pi, Ai are the latitude and longitude of station i,
4>0' Ao are the latitude and longitude of a point chosen as the origin,
Pm is the prime vertical radius of cwvature at mid-latitude,
vm is the meridian radius of cwvatU)"e'at mid-latitude,
4>m is the mid-latitude between 4>i and 4>0' and,
I! = 206264.806"
When using 3 HCPs, eqn. [4] will represent a simple plane SUl"face passing through all the three
points on the surface of the geoid. The east-west and the north-south tilt of the plane surface are
designated by the model coefficients ao,
a I and a2. Once these coefficients are computed using
the least squares method, the unknown orth6metric height of the other GPS points can then be
derived. If more than three HCPs are available, a more complex (cwve) SUl"face can be modelled
in place ofthe plane surface.

3.2 Global Geoid Solutions

Global geoid solutions are obtained from global geopotential models which are given as a set of
coefficients consisting of a series of spherical harmonic functions. The coefficients of the various
terms in the series are determined using a combination of satellite orbit analyses (for the long
wavelength geoid features), terrestrial gravity (medium to short wavelength features) and geoidaJ
heights measured by satellite altimetry over the ocean (medium to short wavelength features).

Some global geopotential models are derived from satellite observations only and are known as
sateUite-only solutions. These models such as GEM9(Lerch et al., 1979), GEM-L2(Lerch et al.,
1982} and GEM-T2(Marsh et al., 1989) involve only low-order spherical harmonics and thus
contain relatively few coefficients. Other global geopotential models which are known as
combined solutions are obtained by adding surface gravimetry and altimetry data to the satellite­
only solutions and they usually contain more coefficients. For example, the development of the
GEMIOB(Lerch et a1.,l981) model are based on the GEM-9 satellite-only model. OSU86F(Rapp
and Cruz, 1986) is based on the GEM-L2 model, while OSU89B(Rapp and Pavlis, 1990) and
OSU9IA(Rapp et al., 1991) are based on the GEM-T2 model.

The geoidaJ height NGM from a global geoid solution is computed from a set of normalised
geopotential coefficients using the following equation:

N GM GM
= - -
ry
'~aJ'
-
n=2
~
L
r m""O
[-.
C~ - j- (
cosm2+ S"m sinm2 P,rn sin (J ) [5]

where,
nMAX is the maximum degree at which the coefficients are known.

3
 r

-* -
Cnm are the Cnm less the ronal coeffidents of the the normal potential of the

selected reference ellipsoid.

G is the gravitational constant

M is the mass of the earth, including the atmosphere.

a is the earth's equatorial radius.

r is the distance from the earth's centre ofmass.

<1>, A are the geocentric latitude and longitude.

P.... (sin ¢) is the nonnalised associated Lagendre ftmction.

g is the nOlmal gravity.

n, m are the degree and order respectively.

Generally, the more coefficients there are in a model, the more detailed the model usually is since
it contains shorter wavelength infonnation of the earth's gravity field. This means that in general,
the best solution to use is one that has determined up to the maximum degree and order of360. In
this study, the global geopotential model adopted is the OSU9lA which was developed using 30'
by 30' mean gravity anomalies derived from terrestrial and altimetric data (Rapp et aI., 1991).
From previous slndy conducted by Ahmad and Kearsley(I994), it was found that this model is
best suited for Malaysia One of the reason for this may be due to the utilization of gravity data
observed over various parts of Malaysia in the OSU91A solution. Once the geoidal heights NGM
have been derived, the ortbometric height of GPS stations in the network can simply be computed
using eqn. [I].

4.0 THE EXPERIMENT

4.1 The GPSlHeight Network
In ordcr to test and evaluate the proposed method of height determination, a network of points
with known heights is clearly needed. As such a network of 10 points with known heights was
established within UTM campus for this purpose. The distribution of the points is shown in
Figure[I.Oj. The heights of the 10 points are derived using the conventional levelling method from
two bench marks (i.e., NI84 and N185) previously established in the campus and this is shown in
Table[2.0]. The GPS observations were made using three Ashtechn ' and one Topcon™ receivers.
A total of 24 baselines were processed using the GPPSTM post-processing softwares.

Table l.() Height of Known Points

STATION HElGlIT (m)

BM02 18.478
BM03 25.897
BM04 32.139
BMOS 39.412
BM06 33.958
BM07 13.813
BM08 09.315
BM09 28.893
BM 10 25.688
BM II 18.%9
Be 10 145.190

4
,.....

4.2 The Tests

Two approaches described previously in (3.1) and (3.2) were tested using the GPS network

Using geoid surface fitting teclmigue

With reference to eqn.[4], a number of surface fitting models can he formed subject to the number
of known HCP available. For this study three models were used. The models takes the following
form:

Model PI : Ni = ao + alxi + aZYi [6]

Model PZ : Ni ~ ao + a I Xi + aZYj + a3XjYi [7]

Model P3: Ni ~ ao + aIxi + aZYi + a3Mlj [8]

In the model P3, Mli= hi - Ito, hi being the ellipsoidal height of the GPS point and hO is the
ellipsoidal height at the point selected as the origin. The rational of adding the height term Ml in
the model is based on the assumption that the geoid undulation generally follow approximately
the topography ofthe area

The first model (PI) requires at least 3 HCP to form the surface model. This model is basically a
plane surface that passes throught all the three points on the surface of the geoid. It is generally
expected that this kind of model is most suitable for small area where the terrain is considered flat
in nature. Tllble [2.0] shows the estimated orthometric height derived using this model. The HCPs
used consist of BCIO, BM08 and BMI I selected according to their location and forming a
triangle enclosing the test area. The r.m.S calculated from the difference between the true and
estimated height for six points is found to be 10.8 em.

Table 2.0 Solution using model PI and three known heights

GEODETIC COORDINATE TRUE ESTIMATED DIFF.

STATION (WGS84) ORTH. HEIGHT(m) (m)
HEIGHT
cb )" It (m)
BM03 013334.1931 103 3801.1370 32.0257 25.897 25.906 0.009
BM04 01 3354.6876 103 3743.2444 38.1366 32139 31.947 0.192
BM05 01 33 54.0404 1033801.839l 45.4659 39.412 39.327 0.085
BM06 01 3329.9165 1033813.1441 40.4659 33.958 34.037 0.079
BM07 01 33 22.0442 103 3826.5890 19.9838 13.813 13945 0.132
BM 10 013346.2812 103 3821.2902 31.8193 25.688 25.740 0.052
RMS- 10.8 cm

The second model PZ was tested to see if it fits more closely with the actual geoid surface of the
area. This is made apparent from the use of the fourth term in the model which will represent any
curvature of the geoid surface. In this model at least four known heights are needed to fit the geoid
surface. The four points chosen are BCIO, BM02, BM09 and BMU. The computed orthometric
heights are shown in Table[3.0]. The r.m.s computed from this model is 3.5 em which is an
improvement over model PI.

5
 Table 3.0 Solution using model P2 and four known heights

GEODETIC COORDINATE TRUE ESTIMATED DJFF.

STATION (WGS84) ORTH. HEIGHT(m) (m)
HEIGHT
d> A h (m)
BM03 013334.1931 103 38 01.l370 32.0257 25.897 25.908 0.011
BM04 01 33 54.6876 103 3743.2444 38.1366 32.139 32.187 0.048
BM05 01 33 54.0404 103380l.8391 45.4659 39.412 39.459 0.047
BM06 013329.9165 1033813.1441 40.4659 33.958 33961 0.003
BM07 01 33 22.0442 103 38 26.5890 19.9838 13.813 13.773 0.040
BM 10 01 3346.2812 103 38 21.2902 31.8193 25.688 25.723 0.035
RMS= 3.5cm

The third model P3, which is a slight variation from model P2 was also used in this experiment.
This model contain the term Mi which refers to the difference between the ellipsoidal height of
each GPS point and the average ellipsoidal height within the test network.This is to reflect the
assumption suggesting the geoid generally follows the terrain. Again similar computation steps as
above were taken. Table [4.0] shows the result computed using 4 HCPs consisting of BCIO,
BM02, BM09 and BMII. The r.m.s computed is found to be 6.6 em and it is apparent that the use
of a height term in the model does not contribute any significant improvement to tbe result
compared to those given by model P2.

Table 4.0 Solution using model P3 and four known heights

GEODETIC COORDINATE TRUE' ESTIMATED DIFF.

STATION (WGS84) ORTH. HEIGHT(m) (m)
HEIGHT
cI> A h (m)
BM03 013334.1931 1033801.1370 32.0257 25.897 25.924 0027
BM04 01 33 54.6876 1033743.2444 38.1366 32.139 32.200 0.061
BM05 01 33 54.0404 1033801.8391 45.4659 39.412 39.499 0.087
BM06 01 3329.9165 10338 13.1441 40.4659 33.958 33.960 0.002
BM07 01 33 22.0442 1033826.5890 19.9838 13813 13.793 0.020
BMlO 013346.2812 103 3821.2902 31.8193 25.688 25.806 0.118
RMS- 6.6cm

Using Global Geoootential Model

Global geopotential model as discussed previously can be used to derive thc geoidaJ height to
correct the ellipsoidal height to give us the orthometric height. The geoidal height is computed
using eqn.[5] using a set of coeflicients. As discussed previously there are quite a number of
global geopotential model available that can be utilised but for this study OSU9IA(Rapp et al.)
coeflicients were used. Table[5.0] shows the orthometric heights derived using the OSU9lA
coefficients. The r.m.s computed from this model is about 82 em and this accuracy is not suitable
for most engineering application.

6
 Table S.O Solution using global geopotential model (OSU9tA) only

GEODETIC COORDINATE TRUE ESTIMATED DIFF.

STATION (WGS84) ORTH. HEtGHT(m) (m)
HEIGHT
cI> A h (m)
BMOJ OJ 3334.1931 103 3801.1370 32.0257 25.897 25.109 0.788
BM04 01 33 54.6876 103 37 43.2444 38.1366 32.139 31.242 0.897
BM05 013354.0404 103 3801.839\ 45.4659 39.412 38.550 0.862
BM06 013329.9165 1033813.1441 40.4659 33.958 33.188 0.770
BM07 0\ 33 22.0442 1033826.5890 19.9838 13.813 13.036 0.777
BM \0 013346.2812 1033821.2902 31.8\93 25.688 24.880 08~
RMS ­ 81.8 em

A strategy to improve the orthometric height estimation using the global geopotential model was
attempted. This is because the. geoidal height computed using the solution described above may
contain biases due to several factors, such as the problem arising from the differences in the GPS
and geoid model datums. These biases can be reduced or absorbed by implementing some kind of
transformation procedure such as that used by Forsberg et al. (1990). The geoid change (N' ­
NMODEL) due to these biases can be expressed in geodetic coordinates in the form of a
regression formula (ibid.):

By using at least four known geoidal heights, N', in the above equation, the four coefficients in the
regression model can be computed. These coefficients are then used·in computing the 'correction'
that will be applied to NMODEL in deriving the geoid height and the orthometric heights at the
other points. Table[6.0] shows the orthometric heights of GPS points computed in this manner
using 4 HCPs(BClO, BM02, BM09 and BMII). The resulting r.m.s of 10.4 em signifies a
significant improvement in the height estimation.

Table 6.0 Solution using OSU9tA and four known heighIII

GEODETIC COORDINATE TRUE ESTIMATED DIFF.

STATION (WGS84) ORTH. HEIGHT(m) (m)
HEIGHT
cI> A h (m)
BM03 013334.\931 1033801.1370 32.0257 25.897 25.895 0.002
BM04 013354.6876 103 3743.2444 38.1366 32.139 32.048 0.091
BM05 0\ 33 54.0404 103 3801.8391 45.4659 39.412 39237 0.175
BM06 01 33 29.9165 10338 13.1441 40.4659 33.958 33.949 0.009
BM07 01 33 22.0442 10338265890 19.9838 13.813 13.831 0.018
BMIO 0133462812 1033821.2902 31.8193 25.688 25.525 0.163
RMS= 10.4 cm

7
 Be 10

BM02

Figure 1.0 The Test GPS Network

5.0 CONCLUSIONS
Based on the tests described above, the following conclusions can be made:
• For small area such as that used in this study, the use of a simple polynomial model (such as
model P2) to fit the local geoid surface is adequate in providing estimated orthometric height.
The expected accuracy in height determination using this approach is about 10 em or better.
This level of accuracy is more than adequate for most engineering applications.
• For small area, using geoidal height from a global geopotential model may give an accuracy
of about 80 ,:m to I metre to the height determination. This level of accuracy is below the
requirement of many engineering applications.
• Using a linear regression to overcome biases that may arise from using a global geopotentail
model does contribute a Significant improvement on the height estimation, but the accuracy
level gained is almost on par with that using a simple plane model depicted by model PI.
• For some engineering applications, the use of GPS data in conjunction with additional known·
height of several points has the potential of replacing the conventional spirit levelling for
height determination.

8
r

I
I
I Acknowledl:ements

I Assistance received from a OlJlllber of people in implementing this study is greatly appreciated Unit
Penyelidikan dan Pembangunan. l!niversiti Tekllologi Malaysia Cor funding this project, Tajul Aritrin Musa
(or his assistance in the data processing; Zuber Che Tak. Hassan Ithnin and Bakri Dahalan for their help in
performing the field observation. and Dr R H Rapp tor providing the OSU91 A gcopotential coefficients

Reference

Ahmad, Z. and A.H.W. Kearsley (1994). Test on High Order Geopotential Models lor Geoid
Evaluation Over the Malaysia Pemnsular Region. The Surveyor. Vol. 29 '\0. :1, The Institution of
Surveyors Malaysia.

Khairul A. Abdullah (1994). UD-HEIGHT : A Computing Tool for GI'S Heighting. Presented
atlh~ hr.\1 IlIrklsh Srll1[1olillll1 (in /)eti>r.lI1ot/(ln<,
Istanhul, Turkey.

Lerch, F.J., S. \1. Klosko, R. E. I.aubscher, and C. A. Wagner (1979). Gravity 'dodel
Improvement using GEOS:1 (GE\19 and 10). J (i('(}[lhrs. f(~I, 84(B8), :1897-:1916.

Marsh, S. M., el al. (1989). The GEM-T2 Gravitational Model. NASA lech M~II1(/. f()()7./Ii.

RaPII, R. H., and J. Y. Cruz (1986). Sphelical Hannonic Expansion of the Earth's Gravitational
Potential to degree :160 using :1()"" mean anomalies. I<e[lori No. ,i:'Ii, Dept. of Geodetic Science, O.
S. U
Rallll, R. U., and N.K. Pavlis (1990). The Development and Analysis of GeopOlenliai
Coefficients Model to Spherical lIannonie Degree 360. J (;etJ[lIH'.'lwl lleasearch, 95(B U).
21885-21911.

RaplI, R.H., Y.VI. Wanl:, and :-.I.K. Pavlis (1991). The Ohio Stale 1991 Geupulenlial and Sea
Surface Topugraphy Hannonic Coetlieienl Models. 11<'[101'/ NII./ltl ot the {Jc[li. lit Ci~(/detic
Science, The Ohio State Uni\ ersity, Columbus, Ohio. 91 pp.

<)