Airborne gravity disturbances in sequential multipole ...

1 downloads 0 Views 717KB Size Report
gravimetric measurements onshore the Azores Islands provided by the ..... 0.11. 0.17 NSMA1(int)-NGPS. -0.79. 0.53. 0.23. 0.16. NFFT-NGPS. S. Jorge. 53. 0.00.
Airborne gravity disturbances in sequential multipole analysis for geoid determination and its test over the Azores D. Marchenko(1), U. Meyer(1), L. Bastos(2) (1)

GeoForschungsZentrum Potsdam, Department. 1: Geodesy and Remote Sensing, Telegrafenberg A17, 14473 Potsdam, Germany (2) Faculdade de Ciências, Universidade do Porto Observatório Astronómico, Monte da Virgem, Portugal

Summary In this paper we introduce a new method to use gravity disturbances in sequential multipole analysis to determine a regional geoid. The basics of the sequential multipole analysis are briefly described and the formulae to use gravity disturbances for input data are derived. It is shown that this method produces better results than using free-air anomalies and that airborne data can significantly improve regional geoid computations. For the test region, the newly derived solutions are compared to an existent FFT derived geoid model. Keywords: regional geoid determination, sequential multipole analysis, gravity disturbances, airborne gravity.

Introduction A new method to directly approximate gravity disturbances in geoid determination is developed and tested on a regional set of gravity data combining satellite altimetry derived free-air gravity anomalies from Kort- og Matrikelstyrelsen in Copenhagen (KMS), gravimetric measurements onshore the Azores Islands provided by the Observatório Astronómico in Porto, plus marine free-air data from RV Haakon Mosby and airborne gravity disturbance data observed during the AGMASCO project. The limits of the test area are defined by 43o N to 34o S in latitude and by 34o W to 22o W in longitude. The method is based on the remove-restore procedure in sequential multipole analysis (using EGM96 (360,360)). The developments are based on a modification of the sequential multipole analysis kernel and the connected recursive formulas for application on gravity disturbances. Details and the derivation of the formulae are given in this paper. The advantage of this method is improvement in accuracy of the final product by avoiding reductions (in this case: free-air reduction) on the airborne gravity disturbances. Thus creating cleaner processing. A preliminary model consisting of about 1400 multipoles was derived using the KMS data set, revision 1999, in combination with land and marine based observations. The surface data was reduced to Faye anomalies (Heiskanen & Moritz, 1967) using the ETOPO'2 terrain model. For the final model the observed gravity disturbances are used for re-adjustment. The direct approximation of airborne gravity disturbances by means of the radial multipole potentials needed new theoretical developments on the required recursive formulas and essential parameters. The applied method for computation allows geoid solutions with accuracies of about 10 cm in the test area. It is shown that the usage of gravity disturbances leads to better results than the usage of free-air anomalies in case of airborne gravity data. The accuracy estimation is based on independent comparisons between the obtained sequential multipole analysis geoid solution and geoid heights (observed in parallel to the land gravity stations) derived from GPS/leveling data. For comparisons to the FFT-model, bicubic spline interpolation is used.

1

1 Basic Theory We suppose that the location of a given point P is described by the coordinates x p , y p , z p within a rectangular Earth-fixed coordinate system, which is defined “in the usual way: the origin is at the center of mass (the geocenter); the z-axis coincides with the mean axis of rotation, the x-axis lies in the mean Greenwich meridian, the y-axis is normal to xz-plane and directed so that xyz system is right-handed; the xy-plane is thus the (mean) equatorial plane" (Moritz, 1980). According to Heiskanen and Moritz (1967, p. 84-88) we will use the following definition for the gravity disturbance, (1.1) g  g ( P)   ( P)  g P   P , taking the following expressions for g and , U W U , , (1.2)   g  n n n where W is the actual gravity potential and U is the normal gravity potential of the ellipsoid. n and n' symbolize the normals at the geoid and ellipsoid respectively. By inserting (1.2) into (1.1) and using T = W – U to describe the anomalous or disturbing potential we get, T  W U  g    , (1.3)   n  n n  P the basic formula for the gravity disturbance since, according to Heiskanen and Moritz (1967, p. 85), “the directions of the normals n and n almost coincide”. Because the normal n to reference ellipsoid is close to the direction of the radius vector r, we derive the approximation T . (1.4) g   r By regarding r as the radius vector of the spherical Earth, the expression (1.4) represents the first term of the well-known expression for a gravity anomaly in the spherical approximation. For further computations of the gravity disturbance functionals in the frame of the removerestore technique we will derive all necessary formulae for the spherical approximation only because the difference between the geodetic and geocentric latitude has a maximal value of about 12 (Marchenko, 1998, p.15). Note also, that if g is known we may solve the corresponding inverse problem to find T and the geoid undulation N above the reference ellipsoid.

2 Representation of gravity disturbances by potentials of radial multipoles The representation of the disturbing potential T by potentials of non-central radial multipoles was developed by A. Marchenko (1998). Each multipole represents a special point object located at the point i inside the Bjerhammar sphere (Marchenko et al. 2001) (see fig. 1). The potential of a non-central radial multipole is characterized by the degree ni and by the geocentric spherical coordinates di , i , i for geocentric distance, latitude and longitude, respectively.

2

Multipole location

P

ri

i

di  i

r

O

Figure 1: Location of radial multipole (i) inside the Bjerhammar sphere

After Marchenko (1998, pp. 125) the disturbing potential T can be represented at any external point P with the geocentric spherical coordinates r, , , by a convergent series of non-orthogonal harmonic functions n

 

where  in

GM  n  a  i (2.1) T ( P)  T (r ,  ,  )    i   v n (r ,  ,  ) , r i 1  r  are the dimensionless coefficients of the expression (2.1) or the dimensionless

multipole moments. vni (r,  ,  ) is the dimensionless potential of a multipole at the point P with each potential v ni having the appropriate degree n  ni . G and M represent the gravity constant and the mass of the Earth. The harmonic function v ni of degree n is now defined by the following expression, 1   ,  qi  where s i is the relative geocentric distance between the origin and the multipole, 1 n v  n! sin i n

di , r qi is the relative distance between the multipole and the point P, r qi  i  1  si2  2si  cos i , r and  i is the geocentric spherical distance between the multipole and point P, cos i  sin  sin  i  cos  cos  i cos(  i ) . si 

(2.2)

(2.3)

(2.4)

(2.5)

The functions v ni can be computed by means of the following recursion formula (Marchenko, 1998, pp. 126): 1  v0i  ,  qi  cos i  si  i (2.6) v1  ,  3 qi  nqi2 v ni  (2n  1)(cos i  si )v ni 1  (n  1)v ni  2   based on the recursion formula for Legendre polynomials (Heiskanen and Moritz, 1967, pp. 23).

3

To derive the gravity disturbance via the potentials (2.2) of radial multipoles we apply the expressions (1.4) and (2.1), assuming a uniform convergence of the series (2.1). Straightforward differentiation of the series (2.1) with respect to r gives the final result,

g ( P )  g ( r ,  ,  )  

 T    ing~ni (r ,  ,  ) , r i 1

(2.7)

or - interchanging summation and differentiation – gives, n 1   T   n n1 g    GM  a   i   v ni (r ,  ,  )   r r  r i 1 

   1   GM  a     i 1 r   r  n



n 1

n i

 v ni (r ,  ,  )  .  

(2.8)

The coefficients of the expansion (2.7) are,

 in   in ,

(2.9)

and the functions g~ni (r ,  ,  ) are the result of the differentiation in (2.8),

g~ni (r ,  ,  )  

GM r2

a   r

n

 v ni  GM  r  (n  1)v ni (r ,  ,  )    2 r  r 

n

a i   g n (r ,  ,  ) . (2.10) r

The relationship (2.10) can be simplified applying the following basic relation, v ni  (n  1)vni . si To solve (2.10) the following substitution will be applied, v ni s v i  i n . r r si Therefore, with v ni s   i (n  1)v ni 1 , r r

(2.11)

(2.12)

(2.13)

our expressions (2.10) and (2.7) become: n



GM a g~ni (r ,  ,  )   2 (n  1)  vni (r ,  ,  )  si vni 1 (r ,  ,  ) r r 

GM r2 n



n

a i   g n (r ,  ,  ) , r



(2.14)



GM  a g ( P)  g (r ,  ,  )   2   in (n  1)  vni (r ,  ,  )  si vni 1 (r ,  ,  ) , r i 1 r

(2.15)

where the dimensionless function g~ni (r ,  ,  ) is given as

g ni (r,  ,  )  (n  1)vni (r,  ,  )  si vni 1 (r,  ,  ).

4

(2.16)

3 Multipole potential parameters determination based on gravity disturbance data Let gravity disturbance data be given at some discrete set of the points Pk  and largest absolute value of the data be located at the point A  Pi such that,

 i  max  k .

(3.1)

k

Thus, we can postulate this external point as the epicenter of a multipole vni and put

i   A  (3.2) . i   A  After Marchenko (1998) a global maximum or the first essential parameter of the functions vni and g ni exists at this epicenter: n 1

 1   , ν ( i  0)    1  si  n  2si  1 g ni ( i  0)  . (1  si ) n  2 i n

(3.3) (3.4)

The global maximum of the dimensionless function g ni (P) will be at the epicenter of this multipole as well. Considering the expression (2.16) for  i  0 , inserting the relationship (3.3) into (2.16) for the required degrees n, (n+1) gives the first essential parameter of the function g ni ( i ) : n2

 1   . (3.5) g ( i  0)  (n  1) 1  s i   The geocentric distance d i , the degree n and the moment  in of one multipole (located at the point i) may be obtained from the empirical isotropic function (i.e. independent of the azimuth; in further reading abbreviated EIF), which is introduced as any discrete function f ( i ) of the spherical distance  i computed by averaging initial gravity disturbance data over the azimuth. The functions g ni ( i ) , f ( i ) have their global extremes at the point A ( f (0)  0 by definition). Under the assumption that the gravity disturbance EIF can be approximated in the vicinity of an epicenter by the function (2.14) based on one multipole potential of degree n, we can determine a preliminary value of the moment  in by applying the relationships (2.14) and (2.16) as, n

r2  r  f (0)   , (3.6)   i GM  a  g n ( i  0) if the values of gravity disturbance data are used for the construction of EIF. Now, with a preliminary known value of the moment  in , we approximate EIF with the n i

analytical isotropic function (AIF) given by one of the functions g ni ( i ) . We determine the relative geocentric distance of the corresponding multipole by least squares in several iterations, if some preliminary value of si is given. A preliminary value of si can be found by applying the second essential parameter or the so-called decreasing length (Marchenko, 1998; Marchenko et al., 2000). If the decreasing length    emp is derived numerically from the normalized EIF (for example by means of the inverse linear interpolation) we obtain the standard non-linear equation, 5

1 i (3.7) f n ( i  0) , 2 regarding the relative distance si. This equation is solved numerically by iterations for the function g ni ( i ) (n0). f ni ( i   emp ) 

4 Data description Within the frame of the AGMASCO project (Airborne Geoid Mapping System for Coastal Oceanography, EU-MAST III, MAS-CT95-0014, Timmen et al., 2000) an aerogravimetric survey took place over the Azores islands in 1997. This data set supplied the airborne gravity disturbance data to test the above described, newly derived technique for geoid determination. In addition to the airborne data marine gravity data assembled in the same project was used in conjunction with KMS’99 gravity anomalies inverted from satellite altimetry (Andersen & Knudsen, 1998). A total of 391 land based observation points distributed over seven of the Azores islands were available with gravity and GPS/leveling measurements parallel to each other. The onshore gravity data was used in the model but the connected GPS/leveling was left out for later independent comparison. The complete data distribution is shown in the fig. 2.

Figure 2: Gravity data distribution; (solid green lines – aerogravity profiles, solid blue lines – marine gravity profiles, red points – land based gravity stations)

5 Geoid processing and results The geoid determination was performed in two steps. A preliminary geoid model was constructed using sequential multipole analysis (SMA). The SMA approach is based on the remove-restore technique using EGM’96 (360,360) as the global gravity model. The resulting SMA model consisted of 1400 multipole moments and positions with an estimated accuracy of 1.5 mGal and based on shipborne, land based and KMS’99 gravity anomalies. For the terrain correction the ETOPO2 (NOAA, 2001) digital elevation model was used. In the second step the derived model was re-adjusted with gravity disturbance data using the direct approximation method (use of gravity disturbances at the measurement altitude) described here. As shown in the subsequent tables, the application of the method results in an improved accuracy of the final model by avoiding reductions (herein assumptions or not exactly known reduction levels, in this case the free-air gradient and the geoid level) applied on the gravity disturbances and creating a cleaner processing. The resulting geoid model (fig. 3) was compared with GPS/leveling data, the global gravity model and also with the existing geoid solution for Azores area computed using FFT (Fernandes et al, 2000). For the comparison of 6

FFT results with GPS/leveling we used a bicubic spline interpolation to restore geoid heights from the FFT model to the points of the GPS/leveling. From tests we derived an accuracy of about 1.5 cm for the bicubic spline interpolation. To check how large the errors of this interpolation are, the same comparison was done with a SMA model, marked as NSMA(int), instead of direct comparison from multipole model. The statistics of this comparison are given in the tab. 2 to 4. Tab. 5 shows the statistic of the comparison of four different SMA models to the GPS/leveling data. Model SMA1 is based on KMS’99 gravity anomalies, marine gravity, land based data and airborne gravity disturbances (fig. 3). Model SMA2 is based on KMS’99 gravity anomalies, marine gravity, land based data and airborne free-air gravity values. Model SMA3 was computed without using any land based gravity data at all. A fourth model, SMA4 is based on KMS’99 gravity anomalies, marine gravity and land based data, leaving out all airborne data. An overview is given in tab. 1. SMA1 (AGD): SMA2 (FAA): SMA3 (no land): SMA4 (no air): FFT (Fernandes et al.)

KMS'99 free-air anomalies KMS'99 free-air anomalies KMS'99 free-air anomalies KMS'99 free-air anomalies KMS'99 free-air anomalies plus AZOMSS99 mean sea surface data

Marine free-air anomalies Marine free-air anomalies Marine free-air anomalies Marine free-air anomalies Marine free-air anomalies (RV Haakon Mosby plus extended marine surveys)

Land based Faye anomalies Land based Faye anomalies Land based Faye anomalies Land based anomalies

Airborne gravity disturbances Airborne free-air anomalies Airborne gravity disturbances Airborne free-air anomalies

Table 1: Sequential multipole analysis model overview

Islands

Points

Graciosa

26

S. Miguel

79

St. Maria

24

Terceira

76

Pico

72

Faial

61

S. Jorge

53

Min 0.00 0.09 -0.03 0.22 -0.35 -0.40 -0.33 -0.71 -0.14 -0.11 -0.17 -0.10 -0.31 -0.53 -0.32 -0.07 -0.35 -0.90 -0.36 -0.18 -0.80 -0.94 -0.83 -0.79 0.00 0.10 0.00 0.26 -0.80

Max 0.28 0.25 0.28 0.47 0.23 0.30 0.24 0.23 0.14 0.34 0.14 0.13 0.45 0.47 0.45 0.74 0.45 0.48 0.48 0.93 0.41 0.44 0.43 0.53 0.27 0.17 0.27 0.53 0.45

Average (m) 0.07 0.07 0.09 0.36 0.00 0.08 0.00 -0.12 0.03 0.00 0.02 0.07 0.03 0.02 0.03 0.31 0.06 0.02 0.06 0.28 0.07 0.10 0.11 0.23 0.07 0.05 0.08 0.40 0.05 7

 (m) 0.05 0.06 0.06 0.06 0.10 0.13 0.11 0.27 0.05 0.05 0.05 0.05 0.13 0.14 0.14 0.14 0.12 0.15 0.14 0.20 0.15 0.16 0.17 0.16 0.05 0.05 0.06 0.05 0.11

Method NSMA1-NGPS NSMA2-NGPS NSMA1(int)-NGPS NFFT-NGPS NSMA1-NGPS NSMA2-NGPS NSMA1(int)-NGPS NFFT-NGPS NSMA1-NGPS NSMA2-NGPS NSMA1(int)-NGPS NFFT-NGPS NSMA1-NGPS NSMA2-NGPS NSMA1(int)-NGPS NFFT-NGPS NSMA1-NGPS NSMA2-NGPS NSMA1(int)-NGPS NFFT-NGPS NSMA1-NGPS NSMA2-NGPS NSMA1(int)-NGPS NFFT-NGPS NSMA1-NGPS NSMA2-NGPS NSMA1(int)-NGPS NFFT-NGPS NSMA1-NGPS

Total

391

-0.94 -0.83 -0.79

0.48 0.48 0.92

0.07 0.05 0.20

0.16 0.13 0.25

NSMA2-NGPS NSMA1(int)-NGPS NFFT-NGPS

Table 2: Comparison of different solutions of geoid models with GPS/Leveling data.

Min Max Average (m)  (m)

NSMA1 – NEGM’96 -0.70 1.38 0.00 0.12

NSMA1 – NFFT -0.52 1.23 0.60 0.16

NFFT – NEGM’96 -0.61 2.21 0.34 0.19

Table 3: Comparison of SMA model with EGM’96 global model and FFT model.

Total points

Points 391

Min -0.76

Max 1.50

 (m) 0.46

Average (m) 0.19

Table 4: Comparison of EGM’96 model with GPS/Leveling data

SMA 1 (AGD) SMA 2 (FAA) SMA 3 (no land data) SMA 4 (no air data)

Min -0.80 -0.94 -0.34 -0.97

Max 0.45 0.48 1.25 0.50

Average (m) 0.05 0.07 0.60 0.05

 (m) 0.11 0.16 0.29 0.17

Table 5: Comparison of SMA models with GPS/Leveling data.

Figure 3: SMA geoid model for Azores area (0.5 meter contours).

8

Figure 4: Differences between SMA1 and FFT (Fernandes contours).

et al, 2000) geoid solutions (0.1 meter

Conclusions The method of direct approximation of airborne gravity disturbances embedded in the sequential multipole analysis was successfully tested. By this, a cleaner data processing gives best results for regional geoid determination (see tab. 5). Whereas the FFT accuracies (from GPS/leveling comparison over the different islands) are within the decimeter range, SMA solutions mostly are within the sub-decimeter range. Even leaving out the land based data set, SMA still leads to a stable solution. Differences between SMA1 and FFT solution given at the fig. 4. From this example we can safely deduce that airborne data of good quality can significantly improve the SMA solution. Direct approximation of gravity disturbance data might be applied to collocation methods and other methods for geoid determination as well, resulting in improved better resolution of regional or local geoid models.

References: Andersen, O. B. ; Knudsen, P., 1998: Global marine gravity field from the ERS-1 and Geosat geodetic mission altimetry, J. Geophys. Res. Vol. 103 , No. C4 , p. 8129 (97JC02198) Fernandes, M. J., Bastos, L., Forsberg, R., Olesen, A., Leite, F., 2000, International Association of Geodesy Symposia, Vol. 121, Schwarz (ed.), Geodesy Beyond 2000 – The Challenges of the First Decade, Springer Verlag,

9

Heiskanen, W. A., Moritz, H., 1967: Physical Geodesy, Freeman and Company, San Francisco and London Marchenko, Barthelmes, F., Meyer, U., Schwintzer, P., 2000: Regional Geoid Determination : An Application to Airborne Gravity Data in the Skagerrak. Scientific Technical Report ; 01/07, GeoForschungsZentrum Potsdam Marchenko, A. N., 1998: Parametrization of the Earth's Gravity Field: Point and Line Singularities, Lviv Astronomical and Geodetical Society, Lviv Timmen, L., Bastos, L., Forsberg, R., Gidskehaug, A., Meyer. U., 2002: Airborne Gravity Field Surveying for Oceanography, geology and Geodesy – Experiences from AGMASCO, IAG, Volume 121, Springer Verlag

10