Preparation of a geodetic calculator for earth gravity field computation

دسته ژئوفیزیک
گروه سازمان زمین شناسی و اکتشافات معدنی کشور
مکان برگزاری بیست و ششمین گردهمایی علوم زمین
نویسنده سید روح الله عمادی
تاريخ برگزاری ۰۲ اسفند ۱۳۸۵



An accurate global model of the Earth’s gravity field is required for many applications in geodesy, geophysics and engineering, A practical application of such accurate geopotential models for surveying is to convert GPS determined ellipsoidal heights to heights referring to the geoid or quasigeoid.(GPS/LEVELLING).High resolution geopotential models have been developed to degree ۱۸۰ for the first time by Rapp (۱۹۷۸) and Lerch et al. (۱۹۷۸). Numerous efforts in refinements of the used data and algorithms have been performed to increase the resolution to degree ۳۶۰. A collaboration of NASA, US National Imagery and  Mapping Agency (NIMA) and Ohio State University has succeeded in computing the most accurate high resolution geopotential model EGM۹۶.(G.Wenzel,۱۹۹۸)



The Earth’s gravity potential field contains infinitely many level surfaces, each parallel to the others. The geoid is one of the equipotential surfaces of the earth’s gravity field. As such, it is defined mathematically by the following prescription:

           W (

) =                                                                                      (۱)               

Where, W stands for the earth’s gravity potential. The value of the  is selected so that this equipotential surface represents the mean sea level (MSL) surface as well as possible. The idea is due to J.K.F.Gauss (۱۷۷۷-۱۸۵۵) who was the first to recognize the geoid as a mathematical figure of the earth.The geoid is usually described by the separation between itself and a mathematical surface. This separation is called geoidal height or geoidal undulation. The mathematical surface is a biaxial ellipsoid with the dimension chosen so that it describes the Earth’s shape as closely as possible and it is the same ellipsoid that is used as a reference surface for measuring latitude and longitude.The geoid is also the reference surface for most height networks since leveling gives the height above the geoid. In geodesy these heights are called orthometric heights.





With introduction of space techniques in geodetic measurements the geoid more important than before. The reason for this is that the Global Positioning System (GPS) and other space techniques are measuring in coordinate system related to the reference ellipsoid. The heights one will get are therefore not the orthometric heights H but the ellipsoidal heights h. There is a need of a precise geoid model to convert ellipsoidal heights to orthometric heights.


       H=h-N                                                                                                         (۲)         


Deflection of vertical

There are three kinds of deflection of vertical in geodesy. The first is geoidal deflection of vertical that defined as the spatial angel between the normal gravity vector on the reference ellipsoid and the actual gravity vector on the geoid (θ).

It can also be interpreted as the maximum slope of the geoid with respect to the reference ellipsoid at the point of interest.


If the angel is reckoned between the normal gravity vector on the reference ellipsoid and the actual gravity vector on the earth’s surface, we define it the surface deflection and denote it by θ’. These two deflections differ by an amount dictated by the curvature of the actual plumb line.


We can compute the surface deflection of vertical with following equations:





( ): astronomical coordinate at interest point

( : geodetic coordinate at interest point

These equations assume that the geodetic datum (G) is properly aligned with the average terrestrial coordinate system (AT).The third species of deflection used in geodesy is defined as the angel between the actual gravity vector on the earth’s surface and the normal gravity vector on the telluroid. This angle is called Molodenskij’s deflection and denoted by .

We can separated deflection of the vertical angle in two components, the north-south component along the meridian is denoted by ξ, and the east-west component in the prime vertical plane is denoted by η .Since the deflection components at the point of interest can be considered as slope of the geoid with respect to the reference ellipsoid, in the meridian and prime vertical direction we have, [Heiskanen and Moritz, ۱۹۸۱]:




Where r is the radius of the earth at the computation point. The negative signs are due to convention, and for deflection components as normally used in geodesy [Bomford, ۱۹۷۱]


Gravity anomaly

The gravity anomaly is defined as the scalar whose value is equal to the difference between the magnitude of the actual gravity on the geoid, , and the normal gravity on the geocentric ellipsoid, ( .



       : observed gravity on the earth surface

      : gravity gradient                                                                                                   

When  is defined in this way, it is called the gravity anomaly on the geoid, or just the

geoidal gravity anomaly. Depending on how the actual gravity on the geoid is obtained from the observed gravity, i.e., which formula is used for vertical gradient of gravity, we get different kinds of gravity anomalies. For geoid determination we use from Free-air gravity gradient:  

     ( mGal/m                                                                          (۱۰)      


A counterpart to the geoidal anomaly is the surface gravity anomaly , which is defined as the difference between the magnitude of the observed gravity taken on the earth’s surface and the normal gravity taken on the telluroid.



This kind of gravity anomaly does not require the knowledge of the vertical gradient of the actual gravity within the earth. The exact value of the normal gravity on the telluroid needed for this anomaly can be obtain from following equation: [Torge.Wolfgang, ۱۹۸۰]       

         +                                                                 (۱۲)                                    


         m=                                                                                             (۱۴)    

 f: flattening of the reference ellipsoid

 m: geodetic constant of the reference ellipsoid

: angular velocity of the earth

 a: semi-major axis of the reference ellipsoid

 b: semi-minor axis of the reference ellipsoid

normal gravity on the reference ellipsoid

and for   we can write: [Vanicek and Krakiwsky, ۱۹۸۶]

 Gal.                                                                                                                    (۱۵)                                     


Actual Gravity Potential                                     

The actual gravity potential of the earth W is the sum of the gravitational potential (W ) and the centrifugal potential (W ),

 W=W +W                                                                              (۱۶)      

The gravitational potential W  is a harmonic function outside the earth (excluding the mass of atmosphere) and can be expressed as an infinite harmonic function, formulated in a geocentric coordinate system as:


W (r, =                           (۱۷)   

W : gravitational potential of the Earth (on the earth or above it)

r, θ ,λ: geocentric radius, geocentric co-latitude, longitude of computation point

a: radius of the Brillouin’s sphere  

GM : gravity-mass constant of the Earth

: fully- normalized coefficient depending on the mass distribution of the earth

:  fully –normalized associated Legendre function of the first kind

The normalized associated Legendre functions of the first kind are defined by means of recurrence relation of type: [A.Ardalan and E.Grafarent, ۲۰۰۰]                       


For n  and m

With starting values                                                                                              (۱۹)  


Earth Geopotential model

The largest force a satellite orbiting the earth experiences is the earth gravity (within reasonable distance).This force due to the earth’s gravitational potential (i.e., the   geopotential) can be modeled using a geopotential model.Gravity field modeling is essentially based on the analysis of satellite tracking data, terrestrial gravity measurement, and satellite altimetry data.Gravity field models consist of spherical harmonic series for disturbing potential, extended  up to degree and order ۷۰(implying more than ۵۰۰۰ harmonic coefficients).The models provide harmonic coefficients(and mostly error estimates)and imply a set of coordinates for stations that were used to track the satellites. Gravity field model may be classified as:

۱-      Satellite-only models based only on satellite tracking data such as GRIM۵-S۱, EIGEN-۱s…

۲-      Combined models: using both, satellite tracking data and terrestrial gravity, such as GRIM۵-C۱, TEG۴,…

۳-      High resolution models: with a spherical harmonic representation, typically up to degree and order ۳۶۰(with more than ۱۳۰۰۰۰ coefficients).They combines the above gravity model with high resolution gravity data derived by satellite altimetry and terrestrial gravity data. Such as, GFZ۹۷, EGM۹۶, …


The EGM۹۶ geopotential model is a high resolution geopotential model of the Earth consisting of spherical harmonic coefficients complete to degree and order ۳۶۰. It is a composite solution, consisting of: (۱) a combination solution to degree and order ۷۰; (۲) a block diagonal solution from degree ۷۱ to ۳۵۹; and (۳) the quadrature solution at degree ۳۶۰. This model, just completed, is the result of collaboration between the National Imagery and Mapping Agency, the NASA۱:placename> Goddard۱:placename> Space۱:placename> Flight۱:placename> Center۱:placetype>۱:place>, and the Ohio۱:placename> State۱:placetype> University۱:placetype>۱:place>.

The resolution of the EGM۹۶-۳۶۰ model on the earth surface is about ۵۵ km, and its Omission Error is about ۰.۲۲۸ meter in height anomalies, and ۲۵.۲۷ mGal in gravity anomalies, and about ۳.۷۶ arc-secs in the vertical deflection [Wenzel, ۱۹۹۸]


If we use from a geopotential model for determination of the gravitational potential of the earth we must make some changes in equation (۱.۱۰) as following [Dru A.Smith, ۱۹۹۸]:

 W (r, =                 (۲۰)     

W : gravitational potential to degree N from geopotential model (i.e.EGM۹۶)

r, θ, λ: geocentric radius, geocentric co-latitude, longitude of computation point

a: equatorial scale factor of the geopotential model  

GM : gravity-mass constant of the geopotential model

: fully- normalized coefficient of geopotential model

: fully –normalized associated Legendre function of the first kind


۲.۱.۶ Normal Gravity Potential

Since the variation of the actual gravity potential (W) of the earth is very small compared to its magnitude, the bulk of the potential is expressed mathematically in such a way that remaining part, called the disturbing or anomalous potential T is small and thus easily handled using even linear approximations. The closest model potential to the actual potential, yet easy handle mathematically, is the normal potential U generated by the normal body (the reference ellipsoid) of the earth with the mass equal(ideally) to the mass of the earth, co-axial with the earth, and spinning with the same velocity as the earth  .

Normal potential same as actual potential is sum of   gravitational potential U  and

Centrifugal potential U ,

U=U +U                                                                                                     (۲۱)        U is a harmonic function and be expanded, in the geocentric coordinate system, in a

Infinite series of harmonic function [Heiskanen and Moritz, ۱۹۸۱]


U (r,θ)=                                                 (۲۲)    

U : normal gravitational potential

r, θ: geocentric radius, geocentric co-latitude

GM : gravity-mass constant of normal field

a : semi-major axis of the reference ellipsoid

: fully- normalized coefficient of normal field

For computation of these coefficients we have:  [Heiskanen and Moritz, ۱۹۸۱]     :



                                             (۲۴)             e: first eccentricity of reference ellipsoid    

: second degree coefficient (Dynamic form factor)

This coefficient can be determined from the four given constant [Bomford, ۱۹۷۱]

The coefficient

سید روح الله عمادی  ،دانشجوی دکتری ژئودزی دانشگاه خواجه نصیر طوسی،عضو هیات علمی دانشگاه آزاد اسلامی واحد تهران-جنوب

فرهاد صادقی،کارشناس ارشد ژئودزی،سازمان نقشه برداری کشور

مرتضی قنبری،کارشناسی ارشد ژئوفیزیک،سازمان زمین شناسی کشور


 اگر چه امروزه تعیین پارامتر های میدان ثقل زمین با استفاده از مدل های جهانی پتانسیل زمین برای کاربردهای مختلف امری رایج برای ژئودزین ها و ژئو فیزیکدان ها شده است ولی در بسیاری از موارد بدلیل تفاوت در ساختاربرنامه ها و تصحیحات لازم اختلافات قابل ملاحظه ای در جوابهای بدست آمده ملاحظه میشود.همچنین اکثر این برنامه ها به نحوی تهیه شده اند که برای عموم استفاده کنند گان قابل استفاده نبوده ،و میبایستی کاربر اطلاعاتی از زبان برنامه نویسی استفاده شده و روابط ریاضی مورد نظر داشته باشد.از این رو نویسندگان این مقاله با توجه به این نیاز و با در نظر گرفتن کلیه موارد علمی مورد نظر سعی نمودند ماشین حسابی را تهیه نموده تا کلیه پارامترهای میدان ثقل زمین را برای نقطه مورد نظر تعیین نماید.

Abstract :

Although, nowadays determination of the earth gravity field parameters for different applications is a usual matter for geodesy and geophysics scientists, but in many cases because of the variety in structure of programs and correction, considerable discrepancy in results were observed. Also most of them are not user-friendly and applicant need to familiar with program structure and related formula. So base on this demand we attempted and in consideration of all scientific corrections, have produced a geodetic calculator for computation  of the earth gravity field components at interesting point.

کلید واژه ها: ژئودزین ژئوفیزیک میدانثقل سایر موارد