- Email: [email protected]
levelling data with an example in the north part of Algeria
ARTICLE IN PRESS
Computers & Geosciences 32 (2006) 1733–1745 www.elsevier.com/locate/cageo
A procedure for modelling the differences between the gravimetric geoid model and GPS/levelling data with an example in the north part of Algeria$ S.A. Benahmed Dahoa,, S. Kahlouchea, J.D. Fairheadb a
National Center of Space Techniques, Geodetic Laboratory—BP 13 Arzew 31200, Algeria b Department of Earth Sciences, GETECH-University of Leeds, Leeds LS2 9JT, UK Received 8 June 2005; received in revised form 2 April 2006; accepted 2 April 2006
Abstract Local gravimetric geoid models from a combination of a global geopotential model and local gravity data generally contain errors of dm-level on the long wavelengths and sometimes they may be significantly higher in areas lacking accurate gravity information like Algeria where only few gravity data from Bureau Gravime´trique International (BGI) have been included in the development of the recent geopotential models. Consequently, these models do not have the required accuracy to transform the GPS ellipsoidal heights to orthometric heights. One of the main causes for this is the limited precision of the global and detailed DTM models. On the other hand, we can now measure by means of the space techniques, on land through a combination of GPS positioning and precise levelling and at sea through satellite altimetry, the geoid on some points on the earth’s surface with very high absolute accuracy. These points can be used to correct the systematic effects, the medium and longer wavelength errors in the gravimetric geoid. The main goal of this study is to propose a procedure, for combination of available gravimetric geoid and external data from GPS and levelling in an optimal way and for estimating the gravimetric geoid accuracy using the collocation approach. So, the question is to find what is the adequate functional representation of the correction that should be applied to the gravimetric geoid? Several functions have been tested and the most suitable will be selected in test area from a statistical testing procedure. For this purpose, the improved Algerian gravimetric geoid computed by the Geodetic Laboratory of the National Center of Space Techniques from the gravity data supplied by the Geophysical Exploration Technology Ltd. (GETECH), and the precise GPS data collected from the international TYRhenian GEOdynamical NETwork (TYRGEONET), ALGerian GEOdynamical NETwork (ALGEONET) projects with baseline length ranging from about 1 to 1000 km have been used. The comparisons based on different GPS campaigns provide after fitting a RMS of the differences 71.9 cm and prove that a good fit in experimental area between the gravimetric geoid and GPS/levelling data using the seven-parameter model transformation has been reached. Moreover, the analysis of statistics shows that the residuals in benchmarks are due principally to gravimetric geoid errors. The main outlines of the Algerian geoid computation, the available GPS/levelling data, the developed procedure and the obtained results will be presented. r 2006 Elsevier Ltd. All rights reserved. Keywords: Corrector surface; Collocation approach; TYRGEONET project
$
Code available from server at http://www.iamg.org/CGEditor/index.htm.
Corresponding author. Tel.: +213 4147 2217; fax: +213 4147 3665.
E-mail address: [email protected] (S.A. Benahmed Daho). 0098-3004/$ - see front matter r 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.cageo.2006.04.003
ARTICLE IN PRESS 1734
S.A. Benahmed Daho et al. / Computers & Geosciences 32 (2006) 1733–1745
1. Introduction The problem of height determination presents one of the most interesting research areas related to GPS. It is well known that ellipsoidal (GPS) height cannot be used in the same way as conventional heights derived from levelling. In order to be able to use GPS technique to obtain accurate conventional heights related to mean sea level (orthometric or normal heights), geoid heights or height anomalies must be known with sufficient accuracy. The crucial part of this method is the geoidal height, which normally is obtained with a lower accuracy than the ellipsoidal height and thus affecting the accuracy of the orthometric height. This geoid is still contaminated by errors in the lower frequencies of the selected global geopotential model. The magnitude of the error is in the order of a decimeter and sometimes it may be significantly higher in areas lacking accurate gravity information like Algeria. If the computed gravimetric geoid is compared to external and independently derived precise geoid information, we may see the error pattern of the long wavelengths. Besides gravimetric approach, precise geoid heights can also be determined geometrically by means of GPS/levelling combination on land area, and altimetric techniques at sea. If a relatively large number of well distributed points where ellipsoid and orthometric heights are available, the geometric geoid derived by GPS and levelling can be determined and the correction to be applied to the gravimetric geoid can be modelled and computed. In this way, the possibility exists to develop an empirical surface (corrector surface) which relates a given gravimetric geoid model to the reference system of GPS ellipsoidal heights, and to the vertical datum of one’s orthometric height system. But, the question is to find what is the adequate functional representation of the correction that should be applied to the gravimetric geoid? Several functions have been tested and the most suitable will be selected in test area from a statistical testing procedure. The main goal of this paper is to propose a procedure, for combination of available gravimetric geoid and external data from GPS and levelling in an optimal way and for estimating the gravimetric geoid accuracy using the collocation approach. For this purpose and in order to test the developed methodology, the improved Algerian gravimetric geoid computed by the Geodetic La-
boratory of the National Center of Space Techniques from the gravity data supplied by the Geophysical Exploration Technology Ltd. (GETECH) and the precise GPS data collected from the international TYRhenian GEOdynamical NETwork (TYRGEONET) and ALGerian GEOdynamical NETwork (ALGEONET) projects have been used.
2. Modelling procedure The procedure we followed can be divided into several steps:
compute detailed gravimetric geoid from gravity data, fit the empirical function to the residuals between the gravimetric and external geoid, apply a statistical test procedure for finding the best representation.
So, an external evaluation for the quality of a gravimetric geoid model can be performed by comparing its interpolated values (N) at a network of GPS benchmarks with the corresponding GPS/ levelling-derived geoid heights (NGPS). Such a comparison is traditionally based on the following model: NGPS Ni ¼ hi Hi Ni ¼ aTi x þ vi , i
(1)
where x is an n 1 vector of unknown parameters, ai is an n 1 vector of known coefficients, and vi denotes a residual random noise term. The parametric model aTi x is supposed to describe the systematic errors and datum inconsistencies inherent in the different height data sets. Its type varies in form and complexity depending on a number of factors. In practice, the various wavelength errors in the gravity solution may be approximated by different kinds of functions in order to fit the geoid to a set of GPS levelling points through an integrated least squares (LS) adjustment. Several models can be used ranging from a simple linear regression to a more complicated seven-parameter similarity transformation model (Personal communication, Kotsakis et al., 2001). The parametric models tested in this paper are: Polynomial model: This model may be adequate to absorb non-periodic large scale systematic errors,
ARTICLE IN PRESS S.A. Benahmed Daho et al. / Computers & Geosciences 32 (2006) 1733–1745
such as systematic levelling errors X aTi x ¼ x0 þ xp;q ðUi U0 Þp p;q
ðki k0 Þq cosq Ui ;
p; q ¼ 0; 1; 2; . . . ,
where U0 , k0 are the mean values of U, k, respectively, and xp,q are the polynomial coefficients. Three-parameter model: aTi x ¼ cos Ui cos ki x1 þ cos Ui sin ki x2 þ sin Ui x3 .
Four-parameter model: aTi x ¼ cos Ui cos ki x1 þ cos Ui sin ki x2 þ sin Ui x3 þ x4 .
1735
reliable geoid assessment (Kotsakis and Sideris, 1999). 3. Data description This section gives an overview on the data sets used in the Algerian geoid computation. However, it is well known that due to the window effect, the results of the FFT techniques in the windowed zone and its immediately adjacent zone have to be totally disregarded because they have unacceptable large errors occurring due to the periodic nature of the discrete Fourier transform. Therefore, the test area used in the computations is larger than the effective area for the computed geoid. The used test window is bounded by the limits 16pjp401 and 101plp141 while the effective window for the geoid results is bounded by the limits 181pjp371 and 81plp121.
Five-parameter model: aTi x ¼ cos Ui cos ki x1 þ cos Ui sin ki x2 þ sin Ui x3 þ sin2 Ui x4 þ x5 . Seven-parameter model: aTi x ¼ cos Ui cos ki x1 þ cos Ui sin ki x2 þ sin Ui x3 cos Ui sin Ui cos ki þ x4 Wi cos Ui sin Ui sin ki þ x5 Wi sin2 Ui þ x6 þ x7 , Wi where Ui and ki are the horizontal geodetic coordinates of the network or baseline points, e is the eccentricity of a common datum ellipsoid and pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Wi ¼ 1 e2 sin2 Ui . The vector of unknown parameters x is solved by minimising the quantity vTv. The adjusted values for the residuals vi give a realistic picture of the level of absolute agreement between the gravimetric geoid and the GPS/levelling data. The previous parameter model absorbs most of the datum inconsistencies among the available height data sets, as well as possible the long wavelength geoid errors. In such comparisons, it should always be kept in mind that the final residual values vi are not purely gravimetric geoid determination error, but contain levelling and GPS positioning errors as well and that it will be separated into its individual components for a more
3.1. Gravity data Between 1986 and 1988, GETECH undertook a comprehensive gravity compilation study of Africa, consisting of over 750,000 land and 1,500,000 marine gravity values constructed into a coherent data set. These data have been acquired in the framework of African Gravity Project (AGP) from Bureau Gravime´trique International (BGI), Defence Mapping Agency (DMA) and from oil companies and many national academic and nonacademic organisations in different countries. The data set is referred to as the IGSN71 gravity datum, processed using the GRS80 gravity formula and terrain corrected to 167 km. All data have been checked and duplicate points removed in a consistent manner. The remaining data have been verified for intersurvey compatibility and the presence of rogue points identified by viewing an isometric 3D triangulated network of the data values on an Intergraph workstation. Erroneous data have been removed from GETECH database (Fairhead et al., 1988). Fig. 1 gives a graphical representation of the gravity data coverage in the computation area. From this figure, it becomes clear that the coverage with gravity observations is not sufficient for some land areas particularly in the south of Algeria and new measurements are needed to accomplish a homogeneous coverage. These data are irregularly distributed and were interpolated onto a regular grid using the minimum curvature technique. The GETECH’s derived
ARTICLE IN PRESS 1736
S.A. Benahmed Daho et al. / Computers & Geosciences 32 (2006) 1733–1745
Fig. 1. Geographical distribution of GETECH gravity data in computation area.
gravity product is in the form of 50 50 grid of freeair, Bouguer and Isostatic gravity anomalies. This does not mean that the observed data set has 50 50 data resolution (Fairhead et al., 1988). For this work, we have used the pre-processed 50 50 grid of the free air anomalies covering the area bounded by the limits 161pjp401 and 101plp141. This grid contains 289 289 points which has been provided to us by GETECH through the agreement between the National Centre of Space Techniques/Geodetic Laboratory and University of Leeds/GETECH without any information on the accuracy of different values. However and since the history of GETECH gravity data processing is not clearly understood, it is very difficult to estimate the accuracy of the resulting geoid heights. However and considering that the 12,472 BGI validated gravity data on the Algerian territory have been included for the generation of the previous 50 grid of anomalies, we have proceeded to assess the prediction accuracy by compar-
Table 1 Statistics of differences between BGI free-air gravity anomalies and predicted ones (mGal) Anomalies
Minimum
Maximum
Mean
St. dev
Dg (Obs) Dg (Pred.) Differences
82.59 89.650 29.923
136.20 107.825 69.666
4.887 4.380 4.657
26.845 26.184 0.502
ing each of the validated BGI free-air gravity anomalies with the corresponding value predicted from the previous GETECH grid using the Spline interpolation. Table 1 summarises the statistics of the differences. The analysis of the results shows the large discrepancies between the BGI free-air gravity anomalies and the predicted ones and proves that gridded gravity set was not derived for precise geoid determination purpose, but for other geophysical applications such as regional geological interpretations.
ARTICLE IN PRESS S.A. Benahmed Daho et al. / Computers & Geosciences 32 (2006) 1733–1745
3.2. Digital terrain model In gravity reduction to correct for terrain effects, and in gravity interpolation/prediction to smooth the gravity field, Digital Terrain Model (topography and bathymetry data) of the region is required. However and besides gravity data, GETECH also provides a 50 50 grid of topography and bathymetry for the whole world called Global DTM5. In constructing this data set, various elevation data sets were used: ETOPO5, height data from gravity stations, topographic maps, national DTM, bathymetry, satellite derived heights and shoreline. Even though it covers land and sea globally, it does not mean that the height of grid points in the data sets refer to a unique vertical datum. It is suspected that the height component of the DTM is adopted directly from the original data source. Another weakness of this data set is the resolution of the data. With 50 50 or about 10 km 10 km grid, the high frequency information of topography cannot be recovered (GETECH, 1995). 3.3. Geopotential coefficient model One of the most important kinds of information for gravity field approximation today is a good reference field supplied by geopotential model. By subtracting the reference field information from the observations, the prediction based on the residual observation is of improved accuracy. This is particularly true for the fast Fourier Transformation method. From the numerical computation point of view, it may reduce the leakage error caused by the existence of the long wavelength parts, which cannot be resolved using only gravity data of limited extent. The problem of choosing a geopotential model which best fits the gravity field in Algeria has not solved definitely. Several studies have shown that the geopotential models tailored to regional or local gravity data are best suited for high precision geoid computations. For the determination of the long wavelength part of the Earth’s gravity field, various models are now available. However, in this work, a global model OSU91A, is completed to a degree and an order 360, was used (Rapp et al., 1991). Replacing the global model OSU91A by the most recent model EGM96 does not improve significantly the present solution, because no new gravity data from Algerian territory was incorporated in EGM96. An analysis based on the gravity data
1737
supplied by the BGI, was carried out to define the geopotential model which fits best the gravity field in Algeria. The following global geopotential solutions GPM2, OSU81, OSU91A and EGM96 are used in the comparison. The result obtained for the geopotential models comparison shows that the OSU91A model fits best the BGI gravity data in Algeria (Benahmed Daho and Kahlouche, 2000).
3.4. Geometrically derived geoid height There are several GPS/levelling points distributed over some regions of Algeria principally in the north part of the country. The distribution is fairly good but the total number of the GPS stations is too small in relation to north part of Algeria’s area. The GPS data used in this work are collected from the international TYRGEONET, ALGEONET projects and some local GPS/levelling surveys with baseline length ranging from about 1 to 1000 km. Started in 1990, the TYRGEONET project constituted by fifty points with some VLBI and SLR stations, was initiated by the Istituto Nazionale di Geofisica e Volcanologia (INGV) and the University of Bologna whose initial objective was the oceanographical study and the geodynamics control of the Italian peninsula. The extension of the TYRGEONET project contributes to set up the new ALGEONET project, composed of ten stations with two points TYRGEONET (Arzew and Algiers), and the essential objective is the geodynamical and seismic control of active zones of the north of Algeria. A first GPS observing campaign of the ALGEONET network was conducted during 6 days in June 1998, simultaneously with the TYRGEONET campaign. The data of this campaign was processed with the Bernese GPS software version 4.2 developed at the Astronomical Institute of the University of Bern (Beutler et al., 2001). Code pseudo ranges were used to estimate receiver a priori coordinates and clock synchronisation with respect to GPS time. Precise satellite ephemerides and satellite clock file information provided by the International GPS Service (IGS) were used in the computation. The final station coordinates are obtained by the free adjustment of all the session solutions for the campaign by the NETGPS software keeping fixed the same station (Algiers) with the same a priori coordinates. This software allows for a rigorous a posteriori statistical analysis because it takes into
ARTICLE IN PRESS 1738
S.A. Benahmed Daho et al. / Computers & Geosciences 32 (2006) 1733–1745
account the complete covariance matrix of the parameters (Crespi, 1996). The ALGEONET networks have been used later for the densification and establishment of the new local geodetic networks. Many GPS campaigns have been carried out in the past years in Algeria. The GPS observations were performed using ASHTECH Z-12 dual frequency receivers with a length between 3 and 12 h and were processed with the Bernese GPS software version 4.2 using the precise ephemerides supplied by IGS. The computed ellipsoidal heights were referred to WGS84 system and their standard deviations do not exceed 3 cm. All GPS stations have been connected to the national levelling network, which consists of orthometric heights. The connection of the GPS stations to the levelling network has been carried out by traditional levelling. The accuracy of the levelling heights may be estimated to about 6 cm depending on the type of connection measurements because some GPS points used in this work as benchmarks are located in mountainous regions in which the spirit levelling would be impractical. 4. Algerian geoid computation The present gravimetric quasigeoid solution is built up in the usual ‘‘remove–restore’’ technique by three terms z ¼ z1 þ z2 þ z3 ,
(2)
where z1 gives the contribution of the geopotential model, while z3 gives the contribution of residual gravity anomalies with the effect of the geopotential model and the terrain removed. z2 gives the indirect effect of the terrain reduction. Terrain effects have been removed in a consistent RTM terrain data reduction, taking into account the topographic irregularities relative to the mean height surface with a resolution of approximately 100 km. The residual gravity anomalies are converted into
residual geoid undulations by spherical FFT evaluation. Pre-processing of geoid computation involves the computation and removal of the geopotential model and terrain contributions from the free-air gravity anomaly. The terrain corrections were computed on a 50 50 grid using DEM data supplied by GETECH. The area considered for DTM is two degree larger than the computation area. In this study the Fortran program TC of the GRAVSOFT package was used to compute the topographic gravimetric correction and subsequent RTM reduction using two-dimensional FFT and strict numerical integration without any approximation. The residual gravity anomalies were calculated by simply subtracting the contribution of OSU91A model complete to degree and order 360 and the terrain effect from the original gravity anomalies. Table 2 shows the impact of the reference field and RTM terrain effect in the considered area. The residuals after reduction are significantly smoother than the original data. The gridded residual gravity anomalies have subsequently been converted to heights anomalies by using spherical FFT method. The FFT was carried out on a grid of 289 289 points excluding zeropadding. The Fortran program SPFOUR of the GRAVSOFT package written by Rene Forsberg was used (Tscherning et al., 1992). Fig. 2 represents the residual quasigeoid heights computed and their statistics with the final Algerian quasigeoid ones tabulated in Table 3. The major contributions to the final quasigeoid are coming from the OSU91A geopotential model. The standard deviations of the contributions from gravity data and DEM are 71.4 m and 70.007 m, respectively. However, it has been noted that in the Algerian territory, the indirect quasigeoid effect is significant only in mountainous areas where it reaches values from one to a few decimeters. On the remaining territory it is on the level of a few
Table 2 Statistics of grid of reduced gravity data in computation area (mGal) Anomalies
Min
Max
Mean
St. dev
DgObs DgOSU91A DgRTM DgObs DgOSU91A Dgr ¼ DgObs DgOUS91A DgRTM
172.956 122.150 295.074 122.900 120.117
218.840 126.422 85.505 192.028 277.075
2.432 4.723 0.729 2.291 1.561
28.807 25.782 5.719 14.936 14.340
ARTICLE IN PRESS S.A. Benahmed Daho et al. / Computers & Geosciences 32 (2006) 1733–1745
1739
Fig. 2. Residual quasigeoid (m).
Table 3 Geoid heights statistics Unit: meter
Minimum
Maximum
Mean
St. dev
zr (Residual effect from FFT) zGM (Global effect from GM) zRTM (RTM effect) Final quasigeoid Nz
7.054 15.847 4.956 17.695 1.446
5.689 56.004 0.854 60.648 0.054
0.000 35.187 0.004 35.183 0.060
1.415 9.184 0.007 8.831 0.132
millimeters. The final quasigeoid heights were obtained by adding three effects. Fig. 3 shows the final quasigeoid surface referred to the GRS80 ellipsoid (Benahmed Daho and Fairhead, 2004). The geoidal height N is derived from the height anomaly z using the expression N z¼
g¯ g¯ Dg H B H, g¯ g0
(3)
where g¯ and g¯ are mean actual and mean normal gravity, respectively, DgB are the Bouguer anomalies and H is the orthometric height. The statistics the difference between N and z are summarised in Table 3. The maximum absolute values of the separation in the computation area mentioned above and in the Algerian territory are 1.45 and 0.54 m, respectively and therefore highly significant.
ARTICLE IN PRESS 1740
S.A. Benahmed Daho et al. / Computers & Geosciences 32 (2006) 1733–1745
Fig. 3. Quasigeoid solution in Algeria (m).
5. Fortran 90 program for combined adjustment between the gravimetric geoid model and GPS/ levelling data A program, called ADJ_GLG, in FORTRAN 90 was developed at National Centre of Space Techniques for combined adjustment between the gravimetric geoid model and GPS/levelling data (Benahmed Daho, 2004). This program had two main objectives. The first aim was to adapt the gravimetric geoid to a set of levelled GPS reference points and to select among the parametric models (see Section 2), the appropriate functional representation of the correction which should be applied to the gravimetric geoid capable to describe all the datum/systematic distortions between the available height data sets. However, the second aim was to proceed to meticulous and reliable analysis gravimetric geoid accuracy according to the developed methodology presented by Kotsakis and Sediris using two modelling approaches called a pure
deterministic parametric model and a hybrid deterministic and stochastic model (Kotsakis and Sideris, 1999). Special attention has been paid to the organisation of data of gravity geoid grid and the GPS/ levelling points to establish the observation equation matrix, and to save memory and to speed up computation. The input data necessary to perform the ADJ_GLG program were: The grid of the geoid in standard and binary format, the levelled GPS benchmarks coordinates, and the variances–covariance matrices of the GPS/levelling and Geoid networks if they are available, otherwise the uniform accuracy of the ellipsoidal, orthometric and geoid heights will be used (Benahmed Daho, 2004). 6. Analysis of results of the numerical tests The accuracy of a geoid model can be assessed in two ways: first, an internal error budget can be
ARTICLE IN PRESS S.A. Benahmed Daho et al. / Computers & Geosciences 32 (2006) 1733–1745
1741
Fig. 4. Geographical distribution of GPS stations (D: Benchmark point, D: Control point).
evaluated by looking at the inherent accuracy of the gravity and topographic data and at the imperfections in the algorithm by which they are transformed to geoid height; secondly, there can be an external comparison with a completely independent method of geoid determination. For our study and in order to conduct comparative evaluation, the gravimetric geoid estimate was compared on 80 points with GPS derived undulations. All of these points are located in the north of Algeria and distributed as in Fig. 4 of which the most are close to the station of Arzew. In these 80 double points, both h (ellipsoidal height) and H (orthometric height) are known so that NGPS/ Lev ¼ hH can be computed. Thus, the NGPS/Lev values can be compared with the gravimetric estimate to assess its precision. However, among these 80 GPS levelled stations only 16 well distributed GPS points are used as benchmarks and the remaining stations considered as control points were excluded in order to estimate the real accuracy given by the comparison between the adjusted values and the known ones. The statistics of the absolute differences in benchmarks
Table 4 Statistics of absolute differences in benchmarks Unit: meter
Minimum
Maximum
Mean
St. dev
NGPS/Lev N gravimetric Differences
43.766 42.374 0.063
48.451 47.939 1.531
47.559 46.915 0.644
1.216 1.458 0.377
are summarised in Table 4 and there indicate that the overall agreement between the gravimetric and GPS/Levelling derived height anomalies is about 38 cm in terms of standard deviation and there exist systematic biases between these two kinds of geoid representations with a mean value of about 0.65 m. These systematic biases are due to the systematic difference between the gravimetric geoid and the orthometric datum as well as the long wavelength error in the gravimetric geoid. In order to minimise these deviations, several tests by varying the type of parametric model (see Section 2) were conducted to find the adequate functional representation of the correction that should be applied to the gravimetric geoid. In this
ARTICLE IN PRESS 1742
S.A. Benahmed Daho et al. / Computers & Geosciences 32 (2006) 1733–1745
Table 5 Comparison of gravimetric with GPS/levelling derived undulations at benchmarks in experimental area after fit Unit: meter
Minimum
Maximum
Mean
St. dev
Linear trend model Second-order polynomial model Third-order polynomial model Three-parameter model Four-parameter model Five-parameter model Seven-parameter model
0.679 0.214 0.027 0.688 0.461 0.255 0.030
0.254 0.094 0.018 0.243 0.454 0.086 0.052
0.000 0.000 0.000 0.000 0.000 0.000 0.000
0.224 0.095 0.019 0.226 0.208 0.097 0.023
Table 6 Comparison of gravimetric with GPS/levelling derived undulations at control points Unit: meter
Minimum
Maximum
Mean
St. dev
Linear trend model Second-order polynomial model Third-order polynomial model Three-parameter model Four-parameter model Five-parameter model Seven-parameter model
0.135 0.231 0.418 0.128 0.287 0.355 0.087
0.357 0.333 0.171 0.358 0.300 0.209 0.201
0.132 0.029 0.000 0.134 0.091 0.056 0.029
0.171 0.088 0.084 0.172 0.141 0.116 0.057
paper, we have used various models in order to determine the possible and improved transformation model that can describe more effectively the general trend of the discrepancies between the GPS/ levelling and the gravimetric geoid. The statistics of the differences after fitting out the systematic biases and tilts between the gravimetric geoid and the GPS/levelling data using the previous parameter models transformation are summarised in Table 5. These statistics in benchmarks show that a good fit between the Algerian geoid and GPS/levelling has been reached when the third-order polynomial was used as corrector surface. Also, we can see that the results from this model are very close to those given by the seven-parameter similarity transformation model in terms of the standard deviation. Levelled GPS points with global residual module larger than 4 times RMS were rejected. No outliers can be identified. This does not necessarily mean that outliers do not exist but more exhaustive statistical means must be used to identify the possible blunders on gravimetric geoid, levelling and GPS networks by using their individual components residuals. Since this is not the purpose of this work, it was assumed that all GPS/levelling data points were valid. In addition, and before any special tests are performed for checking particular features of the general adjustment model, some global statistical
test must be used first. For our case, the Chi2 (w2) test with significance level a ¼ 1% was employed. This test is positive and indicates that the model is not problematic. However, we do not believe that the previous RMS value is the real accuracy of the determined geoid, this provided proof that the combined adjustment can optimally fit the gravity geoid to the GPS levelling point in the least square sense. So and in order to estimate the real accuracy, the GPS/ levelling undulations at control points are compared at adjusted ones computed using the estimated parameters for each model transformation. The statistics of these differences for the all parameter models transformation are summarised in Table 6. Fig. 5 shows the differences between the third-order polynomial and the seven-parameter model transformation at control points. We can see that the combination of the Algerian geoid with GPS/levelling using the last model as the corrector surface gives the best results. The standard deviation of the discrepancies between these two representations at control points is about 76 cm after fitting. This fact confirms the good fit in the experimental area between them. In the next step and in order to estimate the gravimetric geoid accuracy, a common adjustment using the collocation approach of all the 3 data sets
ARTICLE IN PRESS S.A. Benahmed Daho et al. / Computers & Geosciences 32 (2006) 1733–1745
1743
Fig. 5. Differences between third-order polynomial and seven-parameter model transformation at control points.
has been attempted with the gravimetric geoid and the levelled GPS network described above. Moreover, and since the variance–covariance matrix of the GPS and levelling networks adjustments necessary for this kind of combined adjustment are not available, we have used the a priori uniform accuracy of the ellipsoidal, orthometric fixed according to networks accuracy. Furthermore and in absence of the prior geoid error model, we have used a unit weight matrix and have got an estimate for the a posteriori unit geoid variance. On the other hand and for the collocation approach, in which the signal part is considered as additional stochastic parameters, an analytical expression of the covariance function of reduced signals is more convenient. In this work, three different analytical covariance functions broadly used in physical geodesy were tested for the selection of the appropriate analytic covariance model of the reduced signals capable to describe their spatial behaviour: the simple exponential function, Hirvonen model and the second-order Markov. The free parameters in these models are variance of the reduced signals and correlation distance. These parameters are estimated by fitting these models to a number of empirical covariance values using the least square adjustment. The empirical covariance function of the reduced signals was computed using the following formula: C ss ðcÞ ¼
1 X si sj , M
(4)
where M is the number of combinations, c is the spherical distance between Qi and Qj and s is the reduced signal in benchmarks points. The summation was made for all the combinations of the data points Qi and Qj whose distance was comprised between ðc Dc=2Þ and ðc þ Dc=2Þ, and Dc denotes the sampling interval size. By comparison between the adjusted values of parameters of the different models and empirical ones, the best result, in the computation area, are globally obtained by
exponential model. However, it is necessary to note that these results remain tributaries of the density and the quality of data used. Finally, the analyses of the results from different tests using the previous models transformation as the corrector surface shows that residuals in benchmarks are due in the major part to gravimetric geoid errors. For example, Fig. 6 illustrates the individual components of GPS, levelling and geoid residuals relative to appropriate seven-parameter model transformation. 7. Conclusion In this paper a procedure for the choice of the corrector parametric model to be incorporated in Least Squares adjustment problems of linked GPS/ levelling/geoid vertical networks able to absorb the longer wavelength errors in the gravimetric geoid, the systematic errors and datum inconsistencies inherent in the different height data sets was presented. In this way, the possibility exists to develop an empirical surface (corrector surface) which relates a given gravimetric geoid model to the reference system of GPS ellipsoidal heights, and to the vertical datum of one’s orthometric height system. Several models have been tested and the most suitable will be selected in test area from a statistical testing procedure. This procedure has been successfully applied to determine—as a first attempt—a consistent geoid for the north part of Algeria with an accuracy acceptable for the low order levelling network densification. The analysis of the results come out that the seven-parameter model transformation used for modelling the deviations between the gravimetric geoid and the available GPS/levelling points gives the best results in the experimental area and has proven that the four-parameter which corresponds to a simplified version of the similarity (Helmert) transformation model usually employed in most studies dealing with the combined
ARTICLE IN PRESS 1744
S.A. Benahmed Daho et al. / Computers & Geosciences 32 (2006) 1733–1745
Fig. 6. Individual components of GPS, levelling and geoid residuals.
adjustment of GPS, levelling and geoid data is not always sufficient for removing the systematic errors introduced by the datum differences of the three height data sets. Furthermore, a common adjustment of all the 3 heights data sets using the collocation approach and the a prior uniform accuracy of the ellipsoidal, orthometric and geoid heights fixed according to networks accuracy instead of the variance–covariance matrices of different height data sets has been performed. For this first attempt, the results of the numerical tests for all parameter models transformation show that the estimated signals values are insignificant and the residuals in benchmarks are due mainly to gravimetric geoid errors in the experimental area. However, the difficult step in the application of general adjustment model using the collocation approach is the estimation of the covariance function of reduced signals and subsequently the selection of its corresponding analytic representation. In this context, special attention will be given for the optimal choice of a local covariance model of the reduced signals capable to describe their spatial behaviour. Moreover, the global statistical test Chi2 (w2) applied here gives an overall impression of the model validity but further testing is always required in order to validate the adjustment results. The next series of tests to be traditionally carried out are those for individual outliers in the observations, using data snooping and the individual variance
components for each network. Such tests are not discussed here, their usefulness however should be always exploited after having eliminated all possible modelling errors. In perspective and in order to increase the accuracy of the orthometric height values that are derived from GPS observations, the gravimetric geoid model for Algeria should be improved by including the maximum and new surrounding sea and land gravity data, topographic information for terrain modelling, and the new GPS/levelling data with a homogeneous coverage and enough density for a more reliable evaluation of the Algerian geoid accuracy. Acknowledgements The authors gratefully thank the many Institutions and persons who provided data and software, making this work possible: Professors R. Forsberg, C.C. Tscherning and P. Knudsen for the GRAVSOFT package, The Geophysical Exploration Technology Ltd., and the Bureau Gravime´trique International, for the gravity set on the Algerian territory. Appendix A. Supplementary materials Supplementary data associated with this article can be found in the online version at doi:10.1016/ j.cageo.2006.04.003.
ARTICLE IN PRESS S.A. Benahmed Daho et al. / Computers & Geosciences 32 (2006) 1733–1745
References Benahmed Daho, S.A., 2004. A computer program for an adjustment of combined GPS/levelling/geoid networks: case of study: North of Algeria. Newton’s Bulletin No. 2, A Joint Bulletin of the Bureau Gravime´trique International and of the International Geoid Service—ISSN 1810-8547, pp. 42–51. Benahmed Daho, S.A., Fairhead, J.D., 2004. A new quasigeoid computation from gravity and GPS data in Algeria. Newton’s Bulletin No. 2, A Joint Bulletin of the Bureau Gravime´trique International and of the International Geoid Service—ISSN 1810-8547, pp. 52–59. Benahmed Daho, S.A., Kahlouche, S., 2000. Geopotential models comparison in Algeria—Bulletin No. 10 of the International Geoid Service (IGeS)—ISSN 1128-3955, pp. 85–90 Beutler, G., Bock, H., Brockmann, E., Dach, R., Fridez, P., Gurtner, W., Hugentobler, U., Ineichen, D., Johnson, J., Meindl, M., Mervart, L., Rothacher, M., Schaer, S., Springer, T., Weber, R., 2001. Bernese GPS Software Version 4.2 Manual. Astronomical Institute, University of Bern, Bern, Switzerland, (418pp). Crespi, M., 1996. A software package for the adjustment and the analysis of GPS control. In: Ungendoli, M. (Ed.), Reports on Surveying and Geodesy. Nautilus, Bologna, Italy, pp. 237–264.
1745
Fairhead, J.D., Watts, A.B., Chevalier, P., El-Haddadeh, B., Green, C.M., Stuart, G.W., Whaler, K.A., Whibdle, I., 1988. African gravity project. Technical Report, University of Leeds Industrial Services Ltd., Leeds, United Kingdom, 14pp. GETECH, 1995. Global digital terrain model; Global DTM5. Geophysical Exploration Technology, Department of Earth Sciences, University of Leeds, Leeds, UK, CD-ROM. Kotsakis, C., Sideris, M.G., 1999. On the adjustment of combined GPS/levelling/geoid networks. Journal of Geodesy 73, 412–421. Kotsakis, C., Fotopoulos, G., Sideris, M.G., 2001. Optimal fitting of gravimetric geoid undulations to GPS/levelling data using an extended similarity transformation model. Proceedings of the Annual Scientific Meeting of the Canadian Geophysical Union, Ottawa, Canada, May 14–17, 2001. Rapp, R.H., Wang, Y.M., Pavlis, N.K., 1991. The Ohio State 1991 geopotential and sea surface topography harmonic coefficient models. Report 410, Department of Geodetic Science and Surveying, Ohio State University, Columbus, OH, 94pp. Tscherning, C.C., Forsberg, R., Knudsen, P., 1992. Description of the GRAVSOFT package for geoid determination. In: Proceedings of First Continental Workshop on the Geoid in Europe, Prague, pp. 327–334.
Computers & Geosciences 32 (2006) 1733–1745 www.elsevier.com/locate/cageo
A procedure for modelling the differences between the gravimetric geoid model and GPS/levelling data with an example in the north part of Algeria$ S.A. Benahmed Dahoa,, S. Kahlouchea, J.D. Fairheadb a
National Center of Space Techniques, Geodetic Laboratory—BP 13 Arzew 31200, Algeria b Department of Earth Sciences, GETECH-University of Leeds, Leeds LS2 9JT, UK Received 8 June 2005; received in revised form 2 April 2006; accepted 2 April 2006
Abstract Local gravimetric geoid models from a combination of a global geopotential model and local gravity data generally contain errors of dm-level on the long wavelengths and sometimes they may be significantly higher in areas lacking accurate gravity information like Algeria where only few gravity data from Bureau Gravime´trique International (BGI) have been included in the development of the recent geopotential models. Consequently, these models do not have the required accuracy to transform the GPS ellipsoidal heights to orthometric heights. One of the main causes for this is the limited precision of the global and detailed DTM models. On the other hand, we can now measure by means of the space techniques, on land through a combination of GPS positioning and precise levelling and at sea through satellite altimetry, the geoid on some points on the earth’s surface with very high absolute accuracy. These points can be used to correct the systematic effects, the medium and longer wavelength errors in the gravimetric geoid. The main goal of this study is to propose a procedure, for combination of available gravimetric geoid and external data from GPS and levelling in an optimal way and for estimating the gravimetric geoid accuracy using the collocation approach. So, the question is to find what is the adequate functional representation of the correction that should be applied to the gravimetric geoid? Several functions have been tested and the most suitable will be selected in test area from a statistical testing procedure. For this purpose, the improved Algerian gravimetric geoid computed by the Geodetic Laboratory of the National Center of Space Techniques from the gravity data supplied by the Geophysical Exploration Technology Ltd. (GETECH), and the precise GPS data collected from the international TYRhenian GEOdynamical NETwork (TYRGEONET), ALGerian GEOdynamical NETwork (ALGEONET) projects with baseline length ranging from about 1 to 1000 km have been used. The comparisons based on different GPS campaigns provide after fitting a RMS of the differences 71.9 cm and prove that a good fit in experimental area between the gravimetric geoid and GPS/levelling data using the seven-parameter model transformation has been reached. Moreover, the analysis of statistics shows that the residuals in benchmarks are due principally to gravimetric geoid errors. The main outlines of the Algerian geoid computation, the available GPS/levelling data, the developed procedure and the obtained results will be presented. r 2006 Elsevier Ltd. All rights reserved. Keywords: Corrector surface; Collocation approach; TYRGEONET project
$
Code available from server at http://www.iamg.org/CGEditor/index.htm.
Corresponding author. Tel.: +213 4147 2217; fax: +213 4147 3665.
E-mail address: [email protected] (S.A. Benahmed Daho). 0098-3004/$ - see front matter r 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.cageo.2006.04.003
ARTICLE IN PRESS 1734
S.A. Benahmed Daho et al. / Computers & Geosciences 32 (2006) 1733–1745
1. Introduction The problem of height determination presents one of the most interesting research areas related to GPS. It is well known that ellipsoidal (GPS) height cannot be used in the same way as conventional heights derived from levelling. In order to be able to use GPS technique to obtain accurate conventional heights related to mean sea level (orthometric or normal heights), geoid heights or height anomalies must be known with sufficient accuracy. The crucial part of this method is the geoidal height, which normally is obtained with a lower accuracy than the ellipsoidal height and thus affecting the accuracy of the orthometric height. This geoid is still contaminated by errors in the lower frequencies of the selected global geopotential model. The magnitude of the error is in the order of a decimeter and sometimes it may be significantly higher in areas lacking accurate gravity information like Algeria. If the computed gravimetric geoid is compared to external and independently derived precise geoid information, we may see the error pattern of the long wavelengths. Besides gravimetric approach, precise geoid heights can also be determined geometrically by means of GPS/levelling combination on land area, and altimetric techniques at sea. If a relatively large number of well distributed points where ellipsoid and orthometric heights are available, the geometric geoid derived by GPS and levelling can be determined and the correction to be applied to the gravimetric geoid can be modelled and computed. In this way, the possibility exists to develop an empirical surface (corrector surface) which relates a given gravimetric geoid model to the reference system of GPS ellipsoidal heights, and to the vertical datum of one’s orthometric height system. But, the question is to find what is the adequate functional representation of the correction that should be applied to the gravimetric geoid? Several functions have been tested and the most suitable will be selected in test area from a statistical testing procedure. The main goal of this paper is to propose a procedure, for combination of available gravimetric geoid and external data from GPS and levelling in an optimal way and for estimating the gravimetric geoid accuracy using the collocation approach. For this purpose and in order to test the developed methodology, the improved Algerian gravimetric geoid computed by the Geodetic La-
boratory of the National Center of Space Techniques from the gravity data supplied by the Geophysical Exploration Technology Ltd. (GETECH) and the precise GPS data collected from the international TYRhenian GEOdynamical NETwork (TYRGEONET) and ALGerian GEOdynamical NETwork (ALGEONET) projects have been used.
2. Modelling procedure The procedure we followed can be divided into several steps:
compute detailed gravimetric geoid from gravity data, fit the empirical function to the residuals between the gravimetric and external geoid, apply a statistical test procedure for finding the best representation.
So, an external evaluation for the quality of a gravimetric geoid model can be performed by comparing its interpolated values (N) at a network of GPS benchmarks with the corresponding GPS/ levelling-derived geoid heights (NGPS). Such a comparison is traditionally based on the following model: NGPS Ni ¼ hi Hi Ni ¼ aTi x þ vi , i
(1)
where x is an n 1 vector of unknown parameters, ai is an n 1 vector of known coefficients, and vi denotes a residual random noise term. The parametric model aTi x is supposed to describe the systematic errors and datum inconsistencies inherent in the different height data sets. Its type varies in form and complexity depending on a number of factors. In practice, the various wavelength errors in the gravity solution may be approximated by different kinds of functions in order to fit the geoid to a set of GPS levelling points through an integrated least squares (LS) adjustment. Several models can be used ranging from a simple linear regression to a more complicated seven-parameter similarity transformation model (Personal communication, Kotsakis et al., 2001). The parametric models tested in this paper are: Polynomial model: This model may be adequate to absorb non-periodic large scale systematic errors,
ARTICLE IN PRESS S.A. Benahmed Daho et al. / Computers & Geosciences 32 (2006) 1733–1745
such as systematic levelling errors X aTi x ¼ x0 þ xp;q ðUi U0 Þp p;q
ðki k0 Þq cosq Ui ;
p; q ¼ 0; 1; 2; . . . ,
where U0 , k0 are the mean values of U, k, respectively, and xp,q are the polynomial coefficients. Three-parameter model: aTi x ¼ cos Ui cos ki x1 þ cos Ui sin ki x2 þ sin Ui x3 .
Four-parameter model: aTi x ¼ cos Ui cos ki x1 þ cos Ui sin ki x2 þ sin Ui x3 þ x4 .
1735
reliable geoid assessment (Kotsakis and Sideris, 1999). 3. Data description This section gives an overview on the data sets used in the Algerian geoid computation. However, it is well known that due to the window effect, the results of the FFT techniques in the windowed zone and its immediately adjacent zone have to be totally disregarded because they have unacceptable large errors occurring due to the periodic nature of the discrete Fourier transform. Therefore, the test area used in the computations is larger than the effective area for the computed geoid. The used test window is bounded by the limits 16pjp401 and 101plp141 while the effective window for the geoid results is bounded by the limits 181pjp371 and 81plp121.
Five-parameter model: aTi x ¼ cos Ui cos ki x1 þ cos Ui sin ki x2 þ sin Ui x3 þ sin2 Ui x4 þ x5 . Seven-parameter model: aTi x ¼ cos Ui cos ki x1 þ cos Ui sin ki x2 þ sin Ui x3 cos Ui sin Ui cos ki þ x4 Wi cos Ui sin Ui sin ki þ x5 Wi sin2 Ui þ x6 þ x7 , Wi where Ui and ki are the horizontal geodetic coordinates of the network or baseline points, e is the eccentricity of a common datum ellipsoid and pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Wi ¼ 1 e2 sin2 Ui . The vector of unknown parameters x is solved by minimising the quantity vTv. The adjusted values for the residuals vi give a realistic picture of the level of absolute agreement between the gravimetric geoid and the GPS/levelling data. The previous parameter model absorbs most of the datum inconsistencies among the available height data sets, as well as possible the long wavelength geoid errors. In such comparisons, it should always be kept in mind that the final residual values vi are not purely gravimetric geoid determination error, but contain levelling and GPS positioning errors as well and that it will be separated into its individual components for a more
3.1. Gravity data Between 1986 and 1988, GETECH undertook a comprehensive gravity compilation study of Africa, consisting of over 750,000 land and 1,500,000 marine gravity values constructed into a coherent data set. These data have been acquired in the framework of African Gravity Project (AGP) from Bureau Gravime´trique International (BGI), Defence Mapping Agency (DMA) and from oil companies and many national academic and nonacademic organisations in different countries. The data set is referred to as the IGSN71 gravity datum, processed using the GRS80 gravity formula and terrain corrected to 167 km. All data have been checked and duplicate points removed in a consistent manner. The remaining data have been verified for intersurvey compatibility and the presence of rogue points identified by viewing an isometric 3D triangulated network of the data values on an Intergraph workstation. Erroneous data have been removed from GETECH database (Fairhead et al., 1988). Fig. 1 gives a graphical representation of the gravity data coverage in the computation area. From this figure, it becomes clear that the coverage with gravity observations is not sufficient for some land areas particularly in the south of Algeria and new measurements are needed to accomplish a homogeneous coverage. These data are irregularly distributed and were interpolated onto a regular grid using the minimum curvature technique. The GETECH’s derived
ARTICLE IN PRESS 1736
S.A. Benahmed Daho et al. / Computers & Geosciences 32 (2006) 1733–1745
Fig. 1. Geographical distribution of GETECH gravity data in computation area.
gravity product is in the form of 50 50 grid of freeair, Bouguer and Isostatic gravity anomalies. This does not mean that the observed data set has 50 50 data resolution (Fairhead et al., 1988). For this work, we have used the pre-processed 50 50 grid of the free air anomalies covering the area bounded by the limits 161pjp401 and 101plp141. This grid contains 289 289 points which has been provided to us by GETECH through the agreement between the National Centre of Space Techniques/Geodetic Laboratory and University of Leeds/GETECH without any information on the accuracy of different values. However and since the history of GETECH gravity data processing is not clearly understood, it is very difficult to estimate the accuracy of the resulting geoid heights. However and considering that the 12,472 BGI validated gravity data on the Algerian territory have been included for the generation of the previous 50 grid of anomalies, we have proceeded to assess the prediction accuracy by compar-
Table 1 Statistics of differences between BGI free-air gravity anomalies and predicted ones (mGal) Anomalies
Minimum
Maximum
Mean
St. dev
Dg (Obs) Dg (Pred.) Differences
82.59 89.650 29.923
136.20 107.825 69.666
4.887 4.380 4.657
26.845 26.184 0.502
ing each of the validated BGI free-air gravity anomalies with the corresponding value predicted from the previous GETECH grid using the Spline interpolation. Table 1 summarises the statistics of the differences. The analysis of the results shows the large discrepancies between the BGI free-air gravity anomalies and the predicted ones and proves that gridded gravity set was not derived for precise geoid determination purpose, but for other geophysical applications such as regional geological interpretations.
ARTICLE IN PRESS S.A. Benahmed Daho et al. / Computers & Geosciences 32 (2006) 1733–1745
3.2. Digital terrain model In gravity reduction to correct for terrain effects, and in gravity interpolation/prediction to smooth the gravity field, Digital Terrain Model (topography and bathymetry data) of the region is required. However and besides gravity data, GETECH also provides a 50 50 grid of topography and bathymetry for the whole world called Global DTM5. In constructing this data set, various elevation data sets were used: ETOPO5, height data from gravity stations, topographic maps, national DTM, bathymetry, satellite derived heights and shoreline. Even though it covers land and sea globally, it does not mean that the height of grid points in the data sets refer to a unique vertical datum. It is suspected that the height component of the DTM is adopted directly from the original data source. Another weakness of this data set is the resolution of the data. With 50 50 or about 10 km 10 km grid, the high frequency information of topography cannot be recovered (GETECH, 1995). 3.3. Geopotential coefficient model One of the most important kinds of information for gravity field approximation today is a good reference field supplied by geopotential model. By subtracting the reference field information from the observations, the prediction based on the residual observation is of improved accuracy. This is particularly true for the fast Fourier Transformation method. From the numerical computation point of view, it may reduce the leakage error caused by the existence of the long wavelength parts, which cannot be resolved using only gravity data of limited extent. The problem of choosing a geopotential model which best fits the gravity field in Algeria has not solved definitely. Several studies have shown that the geopotential models tailored to regional or local gravity data are best suited for high precision geoid computations. For the determination of the long wavelength part of the Earth’s gravity field, various models are now available. However, in this work, a global model OSU91A, is completed to a degree and an order 360, was used (Rapp et al., 1991). Replacing the global model OSU91A by the most recent model EGM96 does not improve significantly the present solution, because no new gravity data from Algerian territory was incorporated in EGM96. An analysis based on the gravity data
1737
supplied by the BGI, was carried out to define the geopotential model which fits best the gravity field in Algeria. The following global geopotential solutions GPM2, OSU81, OSU91A and EGM96 are used in the comparison. The result obtained for the geopotential models comparison shows that the OSU91A model fits best the BGI gravity data in Algeria (Benahmed Daho and Kahlouche, 2000).
3.4. Geometrically derived geoid height There are several GPS/levelling points distributed over some regions of Algeria principally in the north part of the country. The distribution is fairly good but the total number of the GPS stations is too small in relation to north part of Algeria’s area. The GPS data used in this work are collected from the international TYRGEONET, ALGEONET projects and some local GPS/levelling surveys with baseline length ranging from about 1 to 1000 km. Started in 1990, the TYRGEONET project constituted by fifty points with some VLBI and SLR stations, was initiated by the Istituto Nazionale di Geofisica e Volcanologia (INGV) and the University of Bologna whose initial objective was the oceanographical study and the geodynamics control of the Italian peninsula. The extension of the TYRGEONET project contributes to set up the new ALGEONET project, composed of ten stations with two points TYRGEONET (Arzew and Algiers), and the essential objective is the geodynamical and seismic control of active zones of the north of Algeria. A first GPS observing campaign of the ALGEONET network was conducted during 6 days in June 1998, simultaneously with the TYRGEONET campaign. The data of this campaign was processed with the Bernese GPS software version 4.2 developed at the Astronomical Institute of the University of Bern (Beutler et al., 2001). Code pseudo ranges were used to estimate receiver a priori coordinates and clock synchronisation with respect to GPS time. Precise satellite ephemerides and satellite clock file information provided by the International GPS Service (IGS) were used in the computation. The final station coordinates are obtained by the free adjustment of all the session solutions for the campaign by the NETGPS software keeping fixed the same station (Algiers) with the same a priori coordinates. This software allows for a rigorous a posteriori statistical analysis because it takes into
ARTICLE IN PRESS 1738
S.A. Benahmed Daho et al. / Computers & Geosciences 32 (2006) 1733–1745
account the complete covariance matrix of the parameters (Crespi, 1996). The ALGEONET networks have been used later for the densification and establishment of the new local geodetic networks. Many GPS campaigns have been carried out in the past years in Algeria. The GPS observations were performed using ASHTECH Z-12 dual frequency receivers with a length between 3 and 12 h and were processed with the Bernese GPS software version 4.2 using the precise ephemerides supplied by IGS. The computed ellipsoidal heights were referred to WGS84 system and their standard deviations do not exceed 3 cm. All GPS stations have been connected to the national levelling network, which consists of orthometric heights. The connection of the GPS stations to the levelling network has been carried out by traditional levelling. The accuracy of the levelling heights may be estimated to about 6 cm depending on the type of connection measurements because some GPS points used in this work as benchmarks are located in mountainous regions in which the spirit levelling would be impractical. 4. Algerian geoid computation The present gravimetric quasigeoid solution is built up in the usual ‘‘remove–restore’’ technique by three terms z ¼ z1 þ z2 þ z3 ,
(2)
where z1 gives the contribution of the geopotential model, while z3 gives the contribution of residual gravity anomalies with the effect of the geopotential model and the terrain removed. z2 gives the indirect effect of the terrain reduction. Terrain effects have been removed in a consistent RTM terrain data reduction, taking into account the topographic irregularities relative to the mean height surface with a resolution of approximately 100 km. The residual gravity anomalies are converted into
residual geoid undulations by spherical FFT evaluation. Pre-processing of geoid computation involves the computation and removal of the geopotential model and terrain contributions from the free-air gravity anomaly. The terrain corrections were computed on a 50 50 grid using DEM data supplied by GETECH. The area considered for DTM is two degree larger than the computation area. In this study the Fortran program TC of the GRAVSOFT package was used to compute the topographic gravimetric correction and subsequent RTM reduction using two-dimensional FFT and strict numerical integration without any approximation. The residual gravity anomalies were calculated by simply subtracting the contribution of OSU91A model complete to degree and order 360 and the terrain effect from the original gravity anomalies. Table 2 shows the impact of the reference field and RTM terrain effect in the considered area. The residuals after reduction are significantly smoother than the original data. The gridded residual gravity anomalies have subsequently been converted to heights anomalies by using spherical FFT method. The FFT was carried out on a grid of 289 289 points excluding zeropadding. The Fortran program SPFOUR of the GRAVSOFT package written by Rene Forsberg was used (Tscherning et al., 1992). Fig. 2 represents the residual quasigeoid heights computed and their statistics with the final Algerian quasigeoid ones tabulated in Table 3. The major contributions to the final quasigeoid are coming from the OSU91A geopotential model. The standard deviations of the contributions from gravity data and DEM are 71.4 m and 70.007 m, respectively. However, it has been noted that in the Algerian territory, the indirect quasigeoid effect is significant only in mountainous areas where it reaches values from one to a few decimeters. On the remaining territory it is on the level of a few
Table 2 Statistics of grid of reduced gravity data in computation area (mGal) Anomalies
Min
Max
Mean
St. dev
DgObs DgOSU91A DgRTM DgObs DgOSU91A Dgr ¼ DgObs DgOUS91A DgRTM
172.956 122.150 295.074 122.900 120.117
218.840 126.422 85.505 192.028 277.075
2.432 4.723 0.729 2.291 1.561
28.807 25.782 5.719 14.936 14.340
ARTICLE IN PRESS S.A. Benahmed Daho et al. / Computers & Geosciences 32 (2006) 1733–1745
1739
Fig. 2. Residual quasigeoid (m).
Table 3 Geoid heights statistics Unit: meter
Minimum
Maximum
Mean
St. dev
zr (Residual effect from FFT) zGM (Global effect from GM) zRTM (RTM effect) Final quasigeoid Nz
7.054 15.847 4.956 17.695 1.446
5.689 56.004 0.854 60.648 0.054
0.000 35.187 0.004 35.183 0.060
1.415 9.184 0.007 8.831 0.132
millimeters. The final quasigeoid heights were obtained by adding three effects. Fig. 3 shows the final quasigeoid surface referred to the GRS80 ellipsoid (Benahmed Daho and Fairhead, 2004). The geoidal height N is derived from the height anomaly z using the expression N z¼
g¯ g¯ Dg H B H, g¯ g0
(3)
where g¯ and g¯ are mean actual and mean normal gravity, respectively, DgB are the Bouguer anomalies and H is the orthometric height. The statistics the difference between N and z are summarised in Table 3. The maximum absolute values of the separation in the computation area mentioned above and in the Algerian territory are 1.45 and 0.54 m, respectively and therefore highly significant.
ARTICLE IN PRESS 1740
S.A. Benahmed Daho et al. / Computers & Geosciences 32 (2006) 1733–1745
Fig. 3. Quasigeoid solution in Algeria (m).
5. Fortran 90 program for combined adjustment between the gravimetric geoid model and GPS/ levelling data A program, called ADJ_GLG, in FORTRAN 90 was developed at National Centre of Space Techniques for combined adjustment between the gravimetric geoid model and GPS/levelling data (Benahmed Daho, 2004). This program had two main objectives. The first aim was to adapt the gravimetric geoid to a set of levelled GPS reference points and to select among the parametric models (see Section 2), the appropriate functional representation of the correction which should be applied to the gravimetric geoid capable to describe all the datum/systematic distortions between the available height data sets. However, the second aim was to proceed to meticulous and reliable analysis gravimetric geoid accuracy according to the developed methodology presented by Kotsakis and Sediris using two modelling approaches called a pure
deterministic parametric model and a hybrid deterministic and stochastic model (Kotsakis and Sideris, 1999). Special attention has been paid to the organisation of data of gravity geoid grid and the GPS/ levelling points to establish the observation equation matrix, and to save memory and to speed up computation. The input data necessary to perform the ADJ_GLG program were: The grid of the geoid in standard and binary format, the levelled GPS benchmarks coordinates, and the variances–covariance matrices of the GPS/levelling and Geoid networks if they are available, otherwise the uniform accuracy of the ellipsoidal, orthometric and geoid heights will be used (Benahmed Daho, 2004). 6. Analysis of results of the numerical tests The accuracy of a geoid model can be assessed in two ways: first, an internal error budget can be
ARTICLE IN PRESS S.A. Benahmed Daho et al. / Computers & Geosciences 32 (2006) 1733–1745
1741
Fig. 4. Geographical distribution of GPS stations (D: Benchmark point, D: Control point).
evaluated by looking at the inherent accuracy of the gravity and topographic data and at the imperfections in the algorithm by which they are transformed to geoid height; secondly, there can be an external comparison with a completely independent method of geoid determination. For our study and in order to conduct comparative evaluation, the gravimetric geoid estimate was compared on 80 points with GPS derived undulations. All of these points are located in the north of Algeria and distributed as in Fig. 4 of which the most are close to the station of Arzew. In these 80 double points, both h (ellipsoidal height) and H (orthometric height) are known so that NGPS/ Lev ¼ hH can be computed. Thus, the NGPS/Lev values can be compared with the gravimetric estimate to assess its precision. However, among these 80 GPS levelled stations only 16 well distributed GPS points are used as benchmarks and the remaining stations considered as control points were excluded in order to estimate the real accuracy given by the comparison between the adjusted values and the known ones. The statistics of the absolute differences in benchmarks
Table 4 Statistics of absolute differences in benchmarks Unit: meter
Minimum
Maximum
Mean
St. dev
NGPS/Lev N gravimetric Differences
43.766 42.374 0.063
48.451 47.939 1.531
47.559 46.915 0.644
1.216 1.458 0.377
are summarised in Table 4 and there indicate that the overall agreement between the gravimetric and GPS/Levelling derived height anomalies is about 38 cm in terms of standard deviation and there exist systematic biases between these two kinds of geoid representations with a mean value of about 0.65 m. These systematic biases are due to the systematic difference between the gravimetric geoid and the orthometric datum as well as the long wavelength error in the gravimetric geoid. In order to minimise these deviations, several tests by varying the type of parametric model (see Section 2) were conducted to find the adequate functional representation of the correction that should be applied to the gravimetric geoid. In this
ARTICLE IN PRESS 1742
S.A. Benahmed Daho et al. / Computers & Geosciences 32 (2006) 1733–1745
Table 5 Comparison of gravimetric with GPS/levelling derived undulations at benchmarks in experimental area after fit Unit: meter
Minimum
Maximum
Mean
St. dev
Linear trend model Second-order polynomial model Third-order polynomial model Three-parameter model Four-parameter model Five-parameter model Seven-parameter model
0.679 0.214 0.027 0.688 0.461 0.255 0.030
0.254 0.094 0.018 0.243 0.454 0.086 0.052
0.000 0.000 0.000 0.000 0.000 0.000 0.000
0.224 0.095 0.019 0.226 0.208 0.097 0.023
Table 6 Comparison of gravimetric with GPS/levelling derived undulations at control points Unit: meter
Minimum
Maximum
Mean
St. dev
Linear trend model Second-order polynomial model Third-order polynomial model Three-parameter model Four-parameter model Five-parameter model Seven-parameter model
0.135 0.231 0.418 0.128 0.287 0.355 0.087
0.357 0.333 0.171 0.358 0.300 0.209 0.201
0.132 0.029 0.000 0.134 0.091 0.056 0.029
0.171 0.088 0.084 0.172 0.141 0.116 0.057
paper, we have used various models in order to determine the possible and improved transformation model that can describe more effectively the general trend of the discrepancies between the GPS/ levelling and the gravimetric geoid. The statistics of the differences after fitting out the systematic biases and tilts between the gravimetric geoid and the GPS/levelling data using the previous parameter models transformation are summarised in Table 5. These statistics in benchmarks show that a good fit between the Algerian geoid and GPS/levelling has been reached when the third-order polynomial was used as corrector surface. Also, we can see that the results from this model are very close to those given by the seven-parameter similarity transformation model in terms of the standard deviation. Levelled GPS points with global residual module larger than 4 times RMS were rejected. No outliers can be identified. This does not necessarily mean that outliers do not exist but more exhaustive statistical means must be used to identify the possible blunders on gravimetric geoid, levelling and GPS networks by using their individual components residuals. Since this is not the purpose of this work, it was assumed that all GPS/levelling data points were valid. In addition, and before any special tests are performed for checking particular features of the general adjustment model, some global statistical
test must be used first. For our case, the Chi2 (w2) test with significance level a ¼ 1% was employed. This test is positive and indicates that the model is not problematic. However, we do not believe that the previous RMS value is the real accuracy of the determined geoid, this provided proof that the combined adjustment can optimally fit the gravity geoid to the GPS levelling point in the least square sense. So and in order to estimate the real accuracy, the GPS/ levelling undulations at control points are compared at adjusted ones computed using the estimated parameters for each model transformation. The statistics of these differences for the all parameter models transformation are summarised in Table 6. Fig. 5 shows the differences between the third-order polynomial and the seven-parameter model transformation at control points. We can see that the combination of the Algerian geoid with GPS/levelling using the last model as the corrector surface gives the best results. The standard deviation of the discrepancies between these two representations at control points is about 76 cm after fitting. This fact confirms the good fit in the experimental area between them. In the next step and in order to estimate the gravimetric geoid accuracy, a common adjustment using the collocation approach of all the 3 data sets
ARTICLE IN PRESS S.A. Benahmed Daho et al. / Computers & Geosciences 32 (2006) 1733–1745
1743
Fig. 5. Differences between third-order polynomial and seven-parameter model transformation at control points.
has been attempted with the gravimetric geoid and the levelled GPS network described above. Moreover, and since the variance–covariance matrix of the GPS and levelling networks adjustments necessary for this kind of combined adjustment are not available, we have used the a priori uniform accuracy of the ellipsoidal, orthometric fixed according to networks accuracy. Furthermore and in absence of the prior geoid error model, we have used a unit weight matrix and have got an estimate for the a posteriori unit geoid variance. On the other hand and for the collocation approach, in which the signal part is considered as additional stochastic parameters, an analytical expression of the covariance function of reduced signals is more convenient. In this work, three different analytical covariance functions broadly used in physical geodesy were tested for the selection of the appropriate analytic covariance model of the reduced signals capable to describe their spatial behaviour: the simple exponential function, Hirvonen model and the second-order Markov. The free parameters in these models are variance of the reduced signals and correlation distance. These parameters are estimated by fitting these models to a number of empirical covariance values using the least square adjustment. The empirical covariance function of the reduced signals was computed using the following formula: C ss ðcÞ ¼
1 X si sj , M
(4)
where M is the number of combinations, c is the spherical distance between Qi and Qj and s is the reduced signal in benchmarks points. The summation was made for all the combinations of the data points Qi and Qj whose distance was comprised between ðc Dc=2Þ and ðc þ Dc=2Þ, and Dc denotes the sampling interval size. By comparison between the adjusted values of parameters of the different models and empirical ones, the best result, in the computation area, are globally obtained by
exponential model. However, it is necessary to note that these results remain tributaries of the density and the quality of data used. Finally, the analyses of the results from different tests using the previous models transformation as the corrector surface shows that residuals in benchmarks are due in the major part to gravimetric geoid errors. For example, Fig. 6 illustrates the individual components of GPS, levelling and geoid residuals relative to appropriate seven-parameter model transformation. 7. Conclusion In this paper a procedure for the choice of the corrector parametric model to be incorporated in Least Squares adjustment problems of linked GPS/ levelling/geoid vertical networks able to absorb the longer wavelength errors in the gravimetric geoid, the systematic errors and datum inconsistencies inherent in the different height data sets was presented. In this way, the possibility exists to develop an empirical surface (corrector surface) which relates a given gravimetric geoid model to the reference system of GPS ellipsoidal heights, and to the vertical datum of one’s orthometric height system. Several models have been tested and the most suitable will be selected in test area from a statistical testing procedure. This procedure has been successfully applied to determine—as a first attempt—a consistent geoid for the north part of Algeria with an accuracy acceptable for the low order levelling network densification. The analysis of the results come out that the seven-parameter model transformation used for modelling the deviations between the gravimetric geoid and the available GPS/levelling points gives the best results in the experimental area and has proven that the four-parameter which corresponds to a simplified version of the similarity (Helmert) transformation model usually employed in most studies dealing with the combined
ARTICLE IN PRESS 1744
S.A. Benahmed Daho et al. / Computers & Geosciences 32 (2006) 1733–1745
Fig. 6. Individual components of GPS, levelling and geoid residuals.
adjustment of GPS, levelling and geoid data is not always sufficient for removing the systematic errors introduced by the datum differences of the three height data sets. Furthermore, a common adjustment of all the 3 heights data sets using the collocation approach and the a prior uniform accuracy of the ellipsoidal, orthometric and geoid heights fixed according to networks accuracy instead of the variance–covariance matrices of different height data sets has been performed. For this first attempt, the results of the numerical tests for all parameter models transformation show that the estimated signals values are insignificant and the residuals in benchmarks are due mainly to gravimetric geoid errors in the experimental area. However, the difficult step in the application of general adjustment model using the collocation approach is the estimation of the covariance function of reduced signals and subsequently the selection of its corresponding analytic representation. In this context, special attention will be given for the optimal choice of a local covariance model of the reduced signals capable to describe their spatial behaviour. Moreover, the global statistical test Chi2 (w2) applied here gives an overall impression of the model validity but further testing is always required in order to validate the adjustment results. The next series of tests to be traditionally carried out are those for individual outliers in the observations, using data snooping and the individual variance
components for each network. Such tests are not discussed here, their usefulness however should be always exploited after having eliminated all possible modelling errors. In perspective and in order to increase the accuracy of the orthometric height values that are derived from GPS observations, the gravimetric geoid model for Algeria should be improved by including the maximum and new surrounding sea and land gravity data, topographic information for terrain modelling, and the new GPS/levelling data with a homogeneous coverage and enough density for a more reliable evaluation of the Algerian geoid accuracy. Acknowledgements The authors gratefully thank the many Institutions and persons who provided data and software, making this work possible: Professors R. Forsberg, C.C. Tscherning and P. Knudsen for the GRAVSOFT package, The Geophysical Exploration Technology Ltd., and the Bureau Gravime´trique International, for the gravity set on the Algerian territory. Appendix A. Supplementary materials Supplementary data associated with this article can be found in the online version at doi:10.1016/ j.cageo.2006.04.003.
ARTICLE IN PRESS S.A. Benahmed Daho et al. / Computers & Geosciences 32 (2006) 1733–1745
References Benahmed Daho, S.A., 2004. A computer program for an adjustment of combined GPS/levelling/geoid networks: case of study: North of Algeria. Newton’s Bulletin No. 2, A Joint Bulletin of the Bureau Gravime´trique International and of the International Geoid Service—ISSN 1810-8547, pp. 42–51. Benahmed Daho, S.A., Fairhead, J.D., 2004. A new quasigeoid computation from gravity and GPS data in Algeria. Newton’s Bulletin No. 2, A Joint Bulletin of the Bureau Gravime´trique International and of the International Geoid Service—ISSN 1810-8547, pp. 52–59. Benahmed Daho, S.A., Kahlouche, S., 2000. Geopotential models comparison in Algeria—Bulletin No. 10 of the International Geoid Service (IGeS)—ISSN 1128-3955, pp. 85–90 Beutler, G., Bock, H., Brockmann, E., Dach, R., Fridez, P., Gurtner, W., Hugentobler, U., Ineichen, D., Johnson, J., Meindl, M., Mervart, L., Rothacher, M., Schaer, S., Springer, T., Weber, R., 2001. Bernese GPS Software Version 4.2 Manual. Astronomical Institute, University of Bern, Bern, Switzerland, (418pp). Crespi, M., 1996. A software package for the adjustment and the analysis of GPS control. In: Ungendoli, M. (Ed.), Reports on Surveying and Geodesy. Nautilus, Bologna, Italy, pp. 237–264.
1745
Fairhead, J.D., Watts, A.B., Chevalier, P., El-Haddadeh, B., Green, C.M., Stuart, G.W., Whaler, K.A., Whibdle, I., 1988. African gravity project. Technical Report, University of Leeds Industrial Services Ltd., Leeds, United Kingdom, 14pp. GETECH, 1995. Global digital terrain model; Global DTM5. Geophysical Exploration Technology, Department of Earth Sciences, University of Leeds, Leeds, UK, CD-ROM. Kotsakis, C., Sideris, M.G., 1999. On the adjustment of combined GPS/levelling/geoid networks. Journal of Geodesy 73, 412–421. Kotsakis, C., Fotopoulos, G., Sideris, M.G., 2001. Optimal fitting of gravimetric geoid undulations to GPS/levelling data using an extended similarity transformation model. Proceedings of the Annual Scientific Meeting of the Canadian Geophysical Union, Ottawa, Canada, May 14–17, 2001. Rapp, R.H., Wang, Y.M., Pavlis, N.K., 1991. The Ohio State 1991 geopotential and sea surface topography harmonic coefficient models. Report 410, Department of Geodetic Science and Surveying, Ohio State University, Columbus, OH, 94pp. Tscherning, C.C., Forsberg, R., Knudsen, P., 1992. Description of the GRAVSOFT package for geoid determination. In: Proceedings of First Continental Workshop on the Geoid in Europe, Prague, pp. 327–334.