Abstract: The new satellite gravity missions will lead to a major improvement of the knowledge of the gravity field of the Earth. Primarily the long-wavelength part will be improved (CHAMP, GOCE) and we will also be able to take time variations into account (GRACE).

In order to take advantage of this improvement when computing the quasi-geoid a number of computational procedures will have to be improved. The satellite-only gravity models should be used in remove-restore procedures. Mean gravity values should be computed so that the varying altitude of the point data (on land) is taken into account. Gridding should be avoided. The varying data quality must be taken into account. Downward continuation should be avoided, and if needed only harmonic downward continuation should be used. It must be possible to mix in a rigorous way different gravity field quantities like gravity gradients, gravity anomalies, gravity disturbances and height anomalies observed on the ground, in aircrafts or by a satellite. We must also be able to solve for systematic errors. Also new precise test possibilities have to be established so that we can verify the quality of quasi-geoidal computations based on the new gravity field solutions.

The new satellite gravity missions, CHAMP, GRACE and GOCE, will lead to important improvements of our knowledge of the long wavelength part of Earths gravity field, and thereby of the long wavelengths of the quasi-geoid. A factor 10 improvement compared to EGM96 (Lemoine et al. 1997) is expected, ( Tscherning et al. 2000). GRACE, furthermore, will determine time variations of the gravity field.

Similar improvements in the determination of the shorter wavelength part of the quasi-geoid will be possible. The same improvements of the geoid will be achieved where we here define the geoid as the surface which is constructed from the quasi-geoid using a convention for the determination of the geoid – quasi-geoid separation (see Torge, 1991). (The definition here differs from the one used by many geodesists who still thinks that one may determine the gravity potential inside the Earth accurately). In the following section 2 we will first describe the role of the satellite-only global gravity field models, which will be one of the results of the satellite missions. We then discuss possibilities for the improvement of existing local and regional quasi-geoid determination procedures (section 3), of the methods for global quasi-geoid determination (section 4) and subsequently the challenges we will encounter when we have to convince our customers about the high quality of the new quasi-geoids (section 5).

The determination of the quasi-geoid in an absolute sense is a problem, which only can be solved having a global coverage of for example gravity anomalies. However, for many applications the knowledge of quasi-geoid or height anomaly differences is sufficient. This led to the development of methods like astro-gravimetric geoid determination, where gravity data represented the local information and deflections of the vertical were hoped to represent the global information.

After the advent of the first satellite gravity models, it became clear that they could be used to represent the global, distant information. In the remove-restore technique (Forsberg & Tscherning, 1981) the contribution of the gravity model is first subtracted from the data, a so-called residual geoid is computed, and the contribution of the gravity model is added to the solution.

It was easily proved that the higher were the degree and order of the (error-free !) spherical harmonic coefficients of the gravity model, the smaller an area was needed with local data. This lead to the development of global combined satellite and ground gravity solutions such as for example GPM2, GPM98A,B,C (Wenzel, 1985, 1999), EGM96 (Lemoine et al. 1996) and GRIM5-C1 (Gruber et al., 2000).

Basically two types of methods were used for the determination of the quasi-geoid: Spherical harmonic analysis by integration and least-squares approximation. The integration method requires that the gravity data be downward continued to zero level, while this is not required by the approximation technique. Hence, least-squares approximation solutions have in the experience of the author been better than the other solutions if one takes into account the data base used in the calculations.

But there are other problems associated with the combined solutions. The integration methods, for example, could not use non-uniform weights on the data. The approximation methods have errors due to incorrect assignment of errors to the mean gravity anomalies used as weights in the least-squares adjustment. Typically the error estimate reflects the data coverage and not the fact that the error also depends on the magnitude of the data variance (and covariance). Furthermore mean values are computed by non-optimal procedures, which do not take into account that the data are distributed in varying altitudes, see 3.1.

In many cases did these early combination solutions not use in a satisfactory manner the high quality data available regionally. Methods for updating a global spherical harmonic model using only data from a local area were therefore developed and used in Europe (Basic et al. , 1989).

On the other hand, global combination solutions are still needed, and possibilities for improving existing procedures will be briefly discussed in section 4.

Why do we then need the satellite solutions ? First because there are many unsurveyed areas of the Earth. On the oceans we have today gravity anomalies determined from satellite altimetry, but the effects of the ocean circulation pollute this data since we today only have estimates of the (quasi-stationary) sea-surface topography at very long wavelengths. And where we have gravity data, we do not always know the reference system used. Sometimes Bouguer-anomalies and free-air anomalies are mixed up. So what is important is that the satellite only solutions are not polluted by surface data like the existing high resolution spherical harmonic models based on data combination.

Is everything then so very nice ? Naturally there are problems hidden under the carpet. For GOCE, which is expected to give high resolution (80 km) results, the polar gaps are not covered. Consequently it is quasi-geoid differences, which are well determined (at latitudes between -82 degrees and 82 degees), (Rummel et al., 1995). We therefore have to be careful and take this fact into account when we determine the local/regional geoids using remove-restore techniques.

Similarly, when propagating the errors of the global models into the estimated local/regional geoid heights, we should consider that the errors are correlated. This is a difficult matter to deal with.

Quasi-geoid determination has been one of the favoured activities of geodesists for more than a century. So is it still necessary to improve computational procedures, if not theory ? Yes, the prospects of achieving higher order of precision makes it necessary to improve procedures. Furthermore the use of space methods (GPS/GLONASS/Galileo) for position determination with high precision makes it necessary that the quasi-geoid and thereby height-anomalies are determined with the same precision. Again, what is needed is not THE geoid, but height-anomalies at the surface of the Earth.

3.1. The data.

If data are correct and well distributed nearly any method will work and the results will agree. So the basis for high quality quasi-geoid determination is good, well distributed data. In mountains, especially, gravity has only been collected along roads in valleys. Removing and restoring topographic information, see 3.2, can only partly account this for. New airborne gravity surveys with regularly spaced tracks must be carried out like done in Greenland (Forsberg et al., 1999)

We know that some older gravity surveys are not converted correctly to IGSN71. Frequently the gravity value is precisely measured, but the height and the positions have errors. These can be systematic (datum errors, especially height datum errors) or random. If there is any doubt about the data, then this doubt must be expressed through parameters, which e.g. express a height datum bias, and the parameters must be estimated. The new satellite gravity models should form a very good basis for such a parameter determination if the extend of the survey is larger than the minimum wavelength represented by the satellite model. For sea-gravity a comparison between ship-borne gravimetry and altimeter-derived gravity may be a possibility. This requires that the area being considered only have week currents and (water) density variations or that improved estimates of the sea-surface topography (as obtained from the satellite only gravity field models and satellite radar altimetry) have been used to re-compute altimeter-derived gravity anomalies.

A much disturbing error source is the mixing of anomalies of different types. Many data-bases contain Bouger-anomalies, where in geoid determination free-air anomalies are the basic quantity. The data should be screened, and the correlation with the topographic heights must be determined. If the correlation is zero, something is wrong !

Random errors are also frequent. BGI and NIMA have both devoted a large effort into screening data, but errors will remain, or what is more usual: the users do not know about or use the screened data. Therefore, methods for geoid-determination must be used which filter the data, and which permit that the residuals can be inspected. (The integration methods (see 3.3) may have included a filtering, but the residuals are generally not calculated).

In some methods mean values are used, typically in remote areas. But many mean values have been calculated using simple averaging without taking the varying distribution and quality into account. Especially the distribution in height must be accounted for. Preferably the mean values should be calculated based on upward continued values at an altitude equal to the maximal altitude of the heights in the area in order to avoid any hypothesis about the density inside the topography.

It is also important to be able to take into account all types of data, (see also 3.6). Airborne gravimetry today produces gravity disturbances, which might have to be mixed with gravity anomalies. Gravity disturbances are furthermore obtained when the height is determined using GPS. It is very important to be able to use this data type, because it is not polluted by an unknown height datum zero value. (Naturally these disturbances are like gravity anomalies affected by the error in the zero degree and order coefficient of the reference potential used).

The satellites, especially GOCE, will produce gravity gradient data, which may be used locally with a higher resolution that the one represented in a global gravity field model. This might especially be true in mountainous areas where the signal to noise ratio is large. We should also be prepared to include such data in our computations.

Height anomaly differences obtained from GPS/levelling must also be used. This means that if orthometric heights are used, they must be converted to normal heights using the same procedure which was used to convert the potential differences obtained from the levelling into orthometric heights.

Unfortunately height anomalies are often considered as some kind of absolute quasi-geoid values, but they are only valid as differences due to the levelling datum zero point problem. However, they represent long-wavelength information, which may aid in reducing the error of not having (or using) data from remote areas. Furthermore they have to be used so that calculated height anomalies may be converted to the local datum (see Tscherning et al.,2001).

Some methods require gridded data at zero altitude. But in the gridding process information is lost (Tscherning & Forsberg, 1992), and the data noise becomes correlated in a complicated manner. Another complication is that the grid generally is needed at zero altitude, i.e. in the continents (expect parts of Holland) inside the masses. If the highest precision has to be achieved both gridding and downward continuation must be avoided. Furthermore if one insists on using gridded data one must carefully select the gridding method. Fortunately gridding is generally done using LSC, which (like Kriging) permit parameters and data noise to be taken into account. Also, if downward continuation is unavoidable, harmonic downward continuation should be used following the ideas of Bjerhammar.

At present it is impossible to avoid using integration methods when dealing with areas of continental size (Denker et al., 1997, Forsberg et al., 1996, Smith & Milbert, 1999). However an alternative to the gridding used when applying FFT to evaluate the Stokes integral is ring integration (see Kearsley, 1985, 1987). Here mean values of integration elements are used. Possibly this method can be tuned, using a flexible strategy for selecting the integration elements, so that the information loss is less than when using gridding.

3.2. Remove-restore procedures.

As recalled above, quasi-geoid determination requires global data. But a global gravity field model may represent data far outside the area of interest. The best would be (as mentioned above) to use one of the unpolluted new satellite only gravity fields.

So if we use the remove-restore method we will to a large extend (but not at all completely !) have removed our need for gravity field data outside our area of interest. This is now a standard procedure used with nearly all sensible quasi-geoid determination methods. What we should remember here is that we have also removed the effect of the topography and its possible compensation corresponding to the wavelengths represented in the global model. Hence, when dealing with topographic effects, we must take this into account.

There is another effect of removing (and restoring) a global gravity field model. The removal of the model will make the local mean value much smaller, close to zero. (If this is not the case, the data probably has a systematic error, Bouger anomalies are used instead of free-air anomalies or a too small area is used for the quasi-geoid determination. Furthermore, if the global model is based on data combination it is also possibly that erroneous data was used in the data combination).

If the short wavelength part of the topography is removed the residual quantities are much smoother. It thereby becomes much easier to interpolate (predict) in between data points. In general the error goes down with a factor equal to the ratio between the standard deviation before and after the removal of the global field and the topography. Also the errors of other linear processes (such as computing the height anomalies !) are reduced with the same factor.

In areas with a strongly varying topography (on land or at the ocean bottom) the removal and subsequent restoration of the potential of the topography is very important. On land the ground gravity is most often only available along valleys, and the remove-procedure makes these values more representative. If the data are converted to Faye-anomalies like used for the computation of GEOID99 (Smith and Milbert, 1999) the most important part of the topographic information is lost.

What also has to be remembered in this step is that when removing the long-wavelength part using a spherical harmonic expansion, a large part of the topography (and its isostatic compensation) has already been removed, see e.g. (Smith & Milbert, 1999, Appendix 2). Long-wavelength errors may arise if this is not properly taken into account.

3.3 Methods for quasi-geoid determination.

Basically two methods are used today for local/regional geoid determination: Stokes integration and least-squares approximation (including least-squares collocation, LSC). While Stokes integration can be performed very fast using FFT on gridded data, least squares approximation requires a larger computational effort if the base functions are not selected associated with a regular grid (of point masses for example), see 3.8.

Stokes (or Hotine) integration generally only use one data type with a uniform noise. Methods have however been found to mix two sets of gridded data, see Sideris (1995). However, generally only gravity anomaly data is available in a quantity, which make gridding sensible.

Approximation methods can use mixed data types and take varying noise into account. Since the methods will produce a gravity field model, and not just the quasi-geoid, gravity can also be computed from the model. This can be used to compute the residuals between observed and model (predicted) values, which can be used in outlier detection. The methods also permit the determination of parameters like biases due to datum definition (see Tscherning et al. 2001).

Stokes integration does function very efficiently using some kind of gridded data because the integral may be evaluated using FFT. In fact the key problem here is to limit the loss of information, which takes place when interpolating in order to obtain the grid values as discussed above. Smoothing the observations can do this. Typically the residual topography is removed and later restored – before or after the Stokes integration is performed.

If removed before the Stokes integration is performed, the potential is to a good approximation harmonic down to zero level. Harmonic downward continuation can then be used to construct values at zero level. The advantage of this procedure is that we take advantage of the very smooth field in this step.

3.4 Downward continuation to zero-level.

When using approximation methods (with functions harmonic outside the surface of the Earth as base functions) no downward continuation is required. However in order to stabilize the approximations when infinite dimensional base spaces are used, spaces which contain functions harmonic down to a surface inside the Earth (e.g. a Bjerhammar-sphere) are used. If the topography is removed as completely as possible using a remove-restore method then no negative consequences of this have been found (Tscherning & Forsberg, 1992). In some way the approximation methods have a build-in downward continuation procedure. The downward continued values are simply obtained by evaluating the approximation at zero height.

When using integration methods the same philosophy can be used. Using FFT a preliminary solution is computed, which in the spectral domain is converted to vertical gravity gradients. These gradients are then used to harmonically down-ward continue the gravity values, and the process is repeated until no changes are found. The only problem here is in case the gravity field is strongly varying and that this variation is not represented in the topographic model. The problem has been significant in deep fjords with not precisely known depths (Dahl & Forsberg, 1998).

If the effect of the (residual) terrain is not removed, we have to make some kind of downward continuation through the masses like using the strongly varying and not representative Faye anomalies (Smith & Milbert, 1999). Faye anomalies do not give sufficient smoothing, since only the terrain correction, and not the Bouguer plate is removed. But in general non-harmonic downward continuation should be avoided. See however the results in (Omang & Forsberg, 2000).

3.4 Spherical or Ellipsoidal reference surface ?

The original Stokes and Hotines integration formulae were developed for a spherical earth. We know the errors committed are of the order of the flattening, i.e. 1/300. If, after the removal of a long wavelength field as well as the residual topography, the residual geoid has a magnitude of 0.3 m, then the error is 1 mm. Hence, we are in the interval where we might consider using an ellipsoidal reference surface. Fortunately procedures have been found to make an ellipsoidal correction so that the remaining error is reduced to a quantity of the order of the flattening squared (Fei & Sideris, 2000).

For approximation methods, which take place in a finite dimensional space of functions harmonic outside the Earth, the reference surface plays no role. But if a Bjerhammar sphere is used (as in present implementations of LSC) this has the statistical implication that the gravity variance is smaller at Equator than on the Poles.

3.5. Improving the statistically based methods.

The statistically based approximation method, least-squares collocation (LSC), should in theory give the best solution, where "best" is to be understood as the minimum mean square variance of all data-configurations which may be obtained by moving and rotating the point configuration all over the unit sphere (see Sanso’, 1986). However, this is in the case where a so-called isotropic covariance function or norm is used. It is possible to use an an-isotropic covariance function, but the problem is how to find such a function, see Tscherning (1999).

Another problem is also how to represent the error estimates of the coefficients of the global gravity model, which is used in the remove-restore process. The problem goes partly back to the manner in which the errors of the global models were determined. The errors represent more the data distribution than the gravity field variation. (The error should be large in mountains and small in Denmark). In practical implementations the errors of the gravity field models are scaled so that they are multiplied with a factor, which represents how well the global model, fits the local data. But this is a procedure, which is not well founded in theory.

6. Mixing data types.

Anomalies, disturbances, GPS/levelling and gravity gradients (observed in aircrafts or in GOCE) may be mixed in a rigorous manner in approximation methods. In integration methods the mixing is done by fitting the gravimetric quasi-geoid to height anomaly (differences) obtained from GPS/levelling. (This is typically done without taking into account the varying quality of the quasi-geoid, neither the irregular distribution of the GPS/levelling.)

If this fit is parameterised as a translation it is equivalent to a datum transformation, which may or may not be justified (see Smith & Milbert, 1999). It is much better to include the GPS/levelling with appropriate weights in the set of observations used to determine the least-squares approximation, (see Tscherning et al, 2001). It is, especially in the future, important to be able to mix gravity anomalies (from ground data at points with known orthometric heights) and gravity disturbances (at points with known ellipsoidal height, contingently observed in an aircraft).

3.7 Using non-homogeneous noise.

Using non-homogeneous noise is important if we have data from old and new surveys. Old data may be important in areas with no other data. The use of non-homogeneous noise is unproblematic in approximation methods, but it still requires development before it can be taken into account in the integration methods.

8. Solving big systems of equations.

One of the drawbacks of the approximation methods are the fact that as many unknowns as parameters have to be solved for. This is especially a problem for LSC because the number of equations (parameters) is equal to the number of observations. However for gridded data, the so-called fast collocation method can be used (Bottoni & Barzaghi, 1993), which takes advantage of the repetitive structure of the grid.

If data are of the same kind and is associated with points which altitude does not differ more than a few thousand meters, then the base functions (covariances) may be substituted with functions which are zero outside a given distance (Moreaux et al. 1999, Moreaux, 2001). This makes it possible to use sparse matrix methods combined with the conjugant gradient technique in order to compute a solution. So this problem of the approximation methods is currently disappearing due to the improvements in theory and naturally also due to the immense improvement of computer performance in the recent years.

It is obvious that there is a possibility for improving the satellite-only models using ground data. Many of the considerations in section 3 regarding regional and global models are also valid for the mixed global models.

First of all data at the polar gaps can be used. Procedures exist now for an efficient combination of satellite data and gridded ground gravity data (Sanso’ and Tscherning, 2001), but general approximation methods will also be able to handle such a data-combination such as for example the method used to create EGM96.

The most important improvement should come from improvements of the ground data. As mentioned above mean gravity anomaly computation procedures should be improved . Altimeter derived anomalies may also be improved taking advantage of new SST estimates. Or procedures which use directly the altimeter heights (minus SST) should be used.

At present the data used in the mixed global models are exclusively mean gravity anomalies (and satellite-only spherical harmonic coefficients). Maybe the use of mean-values can be avoided ? In this case error-correlations like those present in air-borne gravity could be taken into account.

The use of the topography in a remove-restore procedure like used for regional modelling is not possible. The global topography contains very large errors which may effect the solution. However it should be used (as for EGM96) in the "densification" of the information in mountains.

The increase in computer performance will also make it possible to use non-uniform weights. It will be possible to solve very big systems of equations.

Other types of base-functions than spherical harmonic (or ellipsoidal harmonic) functions should be investigated. But the base-functions should be harmonic (point-masses, covariance functions/reproducing kernels, (Tscherning, 2001)).. Otherwise equations representing the harmonicity in space must be introduced. (The author does not think wavelets are worthwile trying to use).

Several approximation methods (and especially LSC) gives (formal) error estimates and error-covariances reflecting the data distribution and varying quality, Stokes integration does only have associated an overall error-estimate arising from the mean data error. At least using LSC one can see where there were no data !

The general method for the evaluation of a geoid has therefore been a comparison with external data. Geoid height differences are compared with differences obtained from GPS/levelling. Here a most important set of points is those of the European (Torge) calibration line (Torge et al., 1989) and the North American set of data (Smith and Milbert, 1999).

Here we have agreements at the 10 cm level. But today we should hope for better agreements, and some of the disagreement is surely due to the GPS and the levelling, to datum errors and not just the geoid error. We should therefore re-observe the Torge line, and possibly extend it as originally planned south and east-west. The best would naturally be to have precise GPS positions made at all the points of the European Unified Levelling Network.

The best quasi-geoid results are lately obtained after the gravimetrically determined geoid has been corrected for a bias and tilted. In fact, it is probably the heights, which should have been changed. A tilt is very dangerous, since this hides the long-wavelength errors in the quasi-geoid, as well as systematic errors in the levelling. The users must not be cheated by quoting results obtained after tilting. The (in the opinion of the author: bad) excuse has generally been that the tilted surface could be used to reproduce the levelling results.

If approximation methods are used, it is possible to evaluate any quantity, which can be expressed as a linear, functional applied on the result. If we do not use all data, an evaluation can be made comparing with these data. This again favours the use of approximation methods over integration methods.

The advent of pure satellite gravity models will greatly improve our possibilities for computing precise quasi-geoidal differences, globally, regionally or locally. We must revise the procedures we use both for calculating and for evaluating the quasi-geoid and its errors.

Systematic errors must be removed and data with random errors flagged. Error estimates must be taken into account. Mean values must be re-computed in a correct way. The loss of information due to gridding must be avoided. Data should be used in the points where it is observed.

Methods, which imply explicit non-harmonic downward continuation inside masses, should be avoided. Masses must as far as possible be completely removed before any downward continuation is made.

Error estimates must be available both in terms of formal error estimates produced by error-propagation and in the form of comparisons with external data, not used when computing the quasi-geoid. The Torge calibration line should be re-observed to reach the sub-cm level.

If fitting a quasi-geoid to GPS/levelling is done in order to construct a surface, which permits the users to reconstruct heights in a specific levelling datum one have to be very careful not to cheat the user by delivering a false product. There are many – random and systematic – errors in levelling data, or for that sake also in GPS data (especially missing precise excentricity information). The levelling data should be converted to normal heights.

The first cm-geoid was computed for an area around Hannover in 1987 (Denker and Wenzel, 1987) using least-squares collocation. The method has earlier and was subsequently used in many countries with success see for example (Andreu & Simo, 1992, Barzaghi et al., 1992, Benciolini et al., 1984, Sevilla, 1997, Tscherning, 1982, Tscherning et al., 2001). Approximation using point masses has also been used successfully ((Vermeer, 1994, Kukums, 1993). Meanwhile integration has been used for the recent continental size geoids and a number of regions, see for example (Boziane, 1996, Brunje et al. 1997, Denker et al. 1997, Forsberg et al., 1996, Forsberg, 1998, Kuroshi, 1995, Smith & Milbert, 1999). Now that methods have been found for reducing the large numerical effort required when using LSC then this method should in the future be used for local, regional, if not global, quasi-geoid computation so that the best possible use can be made of the new satellite gravity data and models. A description of software, which fulfils many of the requirements summarized above, can be found in Tscherning (1994).

Acknowledgement. This paper war prepared at the occasion of the 70 birthday of Wolfgang Torge, June 4, 2001.

Thanks to H.Denker and R.H.Rapp for important suggestions for the improvement of the paper.

Andreu, M.A. and C.Simo: Determination del geoide UB91 a Catalunya. Institut Cartografic de Catalunya, Monografies tecniques, Num. 1, 1992.

Basic,T., H.Denker, P.Knudsen, D.Solheim and W.Torge: A new geopotential model tailored to gravity data in Europe. Presented IAG Gen. Meeting, Symp. 103, Edinburgh, Aug., 1989.

Barzaghi,R., M.Brovelli, F.Sanso' and C.C.Tscherning: Geoid Computations in the Mediterranean Area. Mare Nostrum - GEOMED Rep. 1, pp. 7-36, Politecnico di Milano, 1992.

Barzaghi R., Brovelli M. & Sguerso D.: The geoid in the central Mediterranean Area: determination and comparisons. IGeS Bulletin n. 2, pp. 7-24, 1993.

Benciolini, B., L.Mussio, F.Sanso, P.Gasperini and S.Zerbini: Geoid Computation in the Italian Area. Boll. di Geodesia e Sc. Affini, Vol. XLIII, n. 3, pp. 213-244, 1984.

Bottoni, G.P. & R.Barzaghi: Fast Collocation. Bulletin Geodesique, Vol. 67, pp. 119-126, 1993.

Bouziane, M.: The new Algerian quasi-geoid 96 computed using spherical FFT. Laboratory of Geodesy, CNTS, Algeria, 1996.

Brunje, A.J.T, R.Haagmans, E.J. de Min: A preliminary North sea geoid model GEONZ97. MD-rap., MDGAP-9735, Delft, 1997.

Denker, H.: Hochaufloesende regionale Schwerefeldbestimmung mit gravimetrischen und topograpischen Daten. Wiss. Arb. Fachrichtung Vermessungswesen der Univ. Hannover, nr. 156, 1988.

Denker, H. and W.Torge: The European Gravimetric Geoid 1997 - An IAG supported enterprise. IAG Scientific Assembly proceedings, Rio de Janeiro, 1997.

Denker, H. and H.G.Wenzel: Local geoid determination and comparison with GPS results. Bulletin Geodesique, Vol. 61, pp. 349-366, 1987.

Denker, H., D.Behrend & W.Torge: The European Gravimetric Quasigeoid EGG96. IAG Symp., Proc. GraGeoMar96, Springer Verlag, 1997.

Fei, Z.L. and M.G.Sideris: A new method for computing the ellipsoidal correction for Stokes's formula. J. of Geodesy, Vol. 74, pp. 223-231, 2000.

Forsberg, R., D. Solheim and J. Kaminskis: Geoid of the Nordic and Baltic area from gravimetry and satellite altimetry. Proc. Int. Symposium on Gravity, Geoid and Marine Geodesy, Tokyo, pp. 540-548, Sept. 1996.

Forsberg R.: Geoid tayloring to GPS -with example of a 1 cm geoid of Danmark. Proc. 2nd Continental Workshop on the Geoid in Europe, Budapest, Hungary. Reports of the Finnish Geodetic Institute, 98:4, pp. 191-198, March 1998.

Forsberg, R. and C.C.Tscherning: The use of Height Data in Gravity Field Approximation by Collocation. J.Geophys.Res., Vol. 86, No. B9, pp. 7843-7854, 1981.

Forsberg, R., A. V. Olesen, K. Keller: Airborne gravity survey of the North Greenland continental shelf. IAG proceedings volume of Gravity, Geoid and Geodynamics conference, Banff, 2000 (in print).

Gruber, T., A. Bode, C.Reigber, P.Schwintzer, G.Balmino, R.Biancale and J.-L. Lemoine: GRIM5-C1: Combination solution of the global gravity field to degree and order 120. GRL, Vol. 27, pp. 4005-4008, 2000.

Kearsley,A.H.W.: Towards the optimum evaluation of the inner zone contribution to geoidal heights. Aust. J. Geod. Phot. Surv., No. 42, pp. 75-98, 1985.

Kearsley, A. H. W.: The computation of Geoid height differences using ring integration. Boll. di Geodesia e Sc. Aff., 1987.

Kenyeres, A.: GPS-Gravimetric geoid determination based on combination of GPS/Levelling and Gravity Data. Pres. 1. Cont. Workshop on the Geoid in Europe, Prague, May 11-14, 1992.

Kukums, K.: Latvian geoid determination with mass point frequency domain inversion. Rep. Finnish Geodetic Institute, 93:3, Helsinki, 1993.

Kuroishi, Yuki: Precise Gravimetric Determination of Geoid in the Vicinity of Japan. Bulletin of the Geographical Survey Institute, Vol. 41, 93 pp., 1995.

Lemoine, F.G., D.Smith, R.Smith, L.Kunz, E.Pavlis, N.Pavlis, S.Klosko, D.Chinn, M.Torrence, R.Williamson, C.Cox, K.Rachlin, Y.Wang, S.Kenyon, R.Salman, R.Trimmer, R.Rapp and S.Nerem: The development of the NASA GSFC and DMA joint geopotential model. Proc.Symp. on Gravity, Geoid and Marine Geodesy, Sept. 30 - Oct. 5, 1996. The University of Tokyo, Tokyo, 1996.

Moreaux, G., C.C.Tscherning & F.Sanso': Approximation of Harmonic Covariance Functions by non Harmonic Locally Supported Ones. Journal of Geodesy, Vol. 73, pp. 555 - 567, 1999.

Moreaux, G.: Some preconditioners of harmonic spherical spline problems.Inverse Problems, Vol. 17, pp. 157-177, 2001.

Omang, O.C.D., and R.Forsberg: How to handle topography in practical geoid determination: three examples. J. of Geodesy, Vol. 74, pp. 458-466, 2000.

Sanso, F.: Statistical methods in physical geodesy. In: Suenkel, H.:Mathematical and Numerical Techniques in Physical Geodesy. LectureNotes in Earth Sciences, Vol. 7, pp. 49-155, Springer-Verlag, 1986.

Sanso', F. and C.C.Tscherning: Fast spherical collocation. Paper prepared for IAG2001, Budapest, Sept. 2001.

Sevilla, M.: A high-resolution gravimetric geoid in the Strait of Gibraltar. J. of Geodesy, Vol. 71, pp. 402-410, 1997.

Sideris, M.G.: On the use of heterogeneous noisy data in spectral field modelling methods. J. of Geodesy, Vol. 70, pp. 2 - 12, 1995.

Smith, D.A. and D.G.Milbert: The GEOID96 high resolution geoid height model for the United States. J. of Geodesy, Vol. 73, pp. 219 - 236, 1999.

Torge, W.: Geodesy. 2. ed., de Gruyter, Berlin, 1991.

Torge, W., G.Weber and H.-G. Wenzel: Computation of a high resolution European gravimetric geoid (EGG1). Proc. 2nd. Int. Symp. on the Geoid in Europe and the Mediterranean Area, Rome, Sept. 1982. pp. 437-460, Istituto Geografico Militare Italiano, Firenze, 1983.

Torge, W., T.Basic, H.Denker, J.Doliff and H.-G. Wenzel: Longe Range Geoid Control through the European GPS Traverse. DGK Reihe B, Heft Nr. 290, Muenchen 1989.

Tscherning, C.C.: Determination of a (quasi) geoid for the Nordic Countries from heterogeneous data using collocation. Proceedings of the 2nd International Symposium on the Geoid in Europe and Mediterranean Area, Rome 13-17 Sept. 1982, pp. 388-412, Istitu to Geografico Militare Italiano, Firenze, 1983.

Tscherning, C.C.: Geoid determination by least-squares collocation using GRAVSOFT. Lecture Notes "Int. School of the Determination and Use of the Geoid", Milano, Oct., 1994, pp. 135 - 164, published by International Geoid Service, 1994.

Tscherning, C.C.: Construction of an-isotropic covariance-functions using Riez-representers. Journal of Geodesy, Vol. 73, pp. 332-336, 1999.

Tscherning, C.C.: Computation of spherical harmonic coefficients and their error estimates using Least Squares Collocation. Accepted Journal of Geodesy, 2001.

C.C.Tscherning, Anwar Radwan, A.A.Tealeb, S.M.Mahmoud, Abd El-monum Mohamed, Ramdan Hassan, El-Syaed Issawy and K. Saker: Local geoid determination combining gravity disturbances and GPS/levelling: A case study in the Lake Naser area, Aswan, Egypt. Accepted Journal of Geodesy, March 2001.

Tscherning, C.C., D.Arabelos and G.Strykowski: The 1-cm geoid after GOCE. Prepared for GGG2000, Banff, August 2000.

Tscherning,C.C. and R.Forsberg: Harmonic continuation and gridding effects on geoid height prediction. Bulletin Geodesique, Vol. 66, pp. 41-53, 1992.

Vermeer M.: A fast delivery GPS-gravimetric geoid of Estonia. Report 94:1 Finnish Geodetic Institute, Helsinki. 1994.

Wenzel,H.-G.: Hochaufloesende Kugelfunktionsmodelle fuer das Gravitationspotential der Erde. Wiss. Arb. Fachrichtung Vermessungswesen der Universitaet Hannover, no. 137, 1985.

Wenzel, H.G.: Ultra hochaufloesende Kugelfunktionsmodelle GMP98A und GMP98B des Erdschwerefeldes. Proceedings Geodaetische Woche, Kaiserslautern, 1998.