Computers & Geosciences 32 (2006) 1573–1584
Data reduction in scalar airborne gravimetry: Theory, software
and case study in Taiwan
$
Cheinway Hwang
, Yu-Shen Hsiao, Hsuan-Chang Shih
Department of Civil Engineering, National Chiao Tung University, 1001 Ta Hsueh Road, Hsinchu 300, TaiwanReceived 3 February 2005; received in revised form 17 February 2006; accepted 17 February 2006
Abstract
The objective of this paper is to present theories of scalar airborne gravimetry and to publish the related software to interested readers for free access. Numerical procedures and techniques are developed to compute velocities and accelerations from GPS-determined positions. A method based on the correlation analysis of raw gravity reading and vertical acceleration of aircraft is used to correct for gravimeter times. A method and a computer program for crossover adjustment of gravity values along survey lines are developed. This method allows flexible selection of a fixed survey line to overcome the rank defect problem. A computer program is developed to compute gravity anomalies while applying corrections of gravimeter position and filtering. Upward and downward continuations of gravity anomalies are performed using Fast Fourier Transform. Using observed airborne gravity data along a survey line and simulated data in Taiwan, all the programs have been validated and produced reliable results.
r2006 Elsevier Ltd. All rights reserved.
Keywords: Airborne gravimetry; Crossover adjustment; Filtering; Taiwan; Upward continuation
1. Introduction
Airborne gravimetry is a tool for mapping local gravity fields using a combination of airborne sensors, aircraft and positioning systems. It is suitable for gravity survey over difficult terrains and areas mixed with land and ocean. An example of such places is Taiwan, which is surrounded by the Pacific Ocean to the east and the Taiwan Strait to the west. About 75% of Taiwan is covered by hills
and high mountains, which make ground gravity survey rather difficult and expensive.
The principle of airborne gravimetry has been described in the geophysical literature. Examples of scalar gravimeters with a damping system are LaCoste and Romberg (LCR) Air–Sea Gravity System II (LCR, 2003) and Sea–Air–Gravimeter System Kss30/31. Other such gravimeters are discussed in Torge (1989). A number of papers, reports and lecture notes have discussed critical issues in airborne gravimetry.Torge (1989, Chapter 7) and Schwarz and Li (1997) summarize these issues. Torge (1989) also points out the accuracy requirements for aircraft positions, latitude, velo-cities and accelerations in order to achieve a sub-mgal accuracy in airborne gravimetry. Case studies
www.elsevier.com/locate/cageo
0098-3004/$ - see front matter r 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.cageo.2006.02.015
$
Code available from server athttp://www.iamg.org/CGEditor/ index.htm
Corresponding author. Tel.: +886 3 5724739; fax: +886 3 5716257.
E-mail address:[email protected] (C. Hwang).
of airborne gravimetry such as Schwarz and Li (1996), Bell et al. (1999), Childers et al. (2001),
Forsberg et al. (2003), Verdun et al. (2003) also mention important aspects in airborne gravity data reduction and/or gravity field modeling. In sum-mary, the most critical issues are precise positions, velocities and accelerations from Global Positioning System (GPS), an optimal filter for a best noise reduction and a best spatial resolution, a method for crossover analysis, and an optimal method for downward continuation of gravity.
This paper focuses on the above-mentioned issues and presents theories and computer programs for data reduction in airborne gravimetry. Simulated and observed airborne gravity data in Taiwan will be used to assess the performance of the developed computer programs.
2. Principle of scalar airborne gravimetry
In this paper, only scalar airborne gravimetry will be considered. For a complete review of vector gravimetry and scalar gravimetry, one can consult, e.g.,Schwarz and Li (1997). The principle of scalar airborne gravimetry can be expressed as
gz¼fzfb€u þ g0þ 2wecos f þ ve ðRNþzÞ ve þ v 2 n ðRMþzÞ , ð1Þ
where z is the flight altitude above sea level, gz the gravity value at z, fzthe gravimeter reading at z, fb the gravimeter reading at the airport, called base reading, €u the vertical acceleration of aircraft (positive when pointing to zenith), wethe rotational velocity of the earth (7.292115 105rad s1) f the latitude, RN,RM are radii of curvatures along the prime vertical and the meridian (Torge, 1989, p. 36), practically can be replaced by the Earth’s mean radius (about 6371 km), ne,nn are east and north components of velocity and g0is the gravity value at the airport.
The last two terms in (1) form the Eo¨tvo¨s effect. From (1), it is clear that the gravimeter readings (fzand fb) form only one part of airborne gravimetry. The accuracies of the aircraft position, velocity and vertical acceleration will also govern the accuracy of the gravity value. With the advances of GPS, it is possible to determine aircraft positions with a cm-level accuracy (Goad and Yang, 1997). With an adequate numerical technique, highly accurate
velocities and accelerations can be obtained from the positions. However, in a turbulent flight, sudden changes of aircraft positions will occur and the vertical accelerations of the aircraft might well exceed gravity values themselves (E9.8 m s2). Thus, only in a ‘‘smooth’’ flight, useful gravity values can be obtained. Thus, it is clear that gravimeter, GPS and aircraft form the three essential parts of airborne gravimetry.
3. Data reduction: theory
3.1. Numerical differentiation of position for determining velocity and acceleration
In a typical airborne gravity survey, GPS kine-matic positioning is used to determine aircraft position, see, e.g., Goad and Yang (1997) and
Verdun et al. (2003). Commercial software like GPSurvey or scientific software like Bernese Version 5.0 (Beutler et al., 2004) can be used for GPS kinematic positioning. With a GPS data sampling rate of 1 or 2 Hz and a careful data processing, the precision of GPS-determined aircraft positions can be at cm level (Goad and Yang, 1997). However, due to problems mainly in integer ambiguity resolution, GPS kinematic positioning may produce long wavelength errors up to several decimeters (Goad and Yang, 1997). If the coordinates of the aircraft are expressed in the rectangular system (x, y, z, see Fig. 1), the velocity and acceleration vectors in the local topographic system (n, e, u) can
h north (n) east (e) up (u) gravity vector equator geocenter z aircraft x, Greenwich Meridian y
Fig. 1. Geocentric coordinate system (x, y, z) and local topo-graphic coordinate system (n, e, u).
be computed as _rt¼ _n _e _u 2 6 4 3 7 5 ¼ vn ve vu 2 6 4 3 7 5 ¼ R _ x _y _z 2 6 4 3 7 5 ¼ R_rg, (2) €rt¼ €n €e €u 2 6 4 3 7 5 ¼ R € x €y €z 2 6 4 3 7 5 ¼ R€rg, (3)
where R is a rotation matrix defined as (Seeber, 1993, p. 19)
R ¼
sin f cos l sin f sin l cos f
sin l cos l 0
cos f cos l cos f sin l sin f 2 6 4 3 7 5, (4)
with l being longitude. In this case, the vectors _rg
and €rg are obtained by numerical differentiations of
rg, and rotations to _rt and €rt using R. The vertical
component in €rn is the value €u needed in (1).
Alternatively, the velocity and acceleration of the aircraft can be computed as
_rt¼ ðRMþhÞ _f ðRNþhÞ cos f_l _ h 2 6 4 3 7 5, (5) €rt¼ €n €e €u 2 6 4 3 7 5 ¼ _vn _ve _vu 2 6 4 3 7 5, (6)
where h is ellipsoidal height, and _f, _l and _h are velocities of latitude, longitude and ellipsoidal height, respectively. In this case, _f, _l and _h are first obtained by numerical differentiations of _f, _l and _h, and then converted to _rtusing (5). The acceleration
vector €rt is simply obtained by numerical
differ-entiation of _rt.
We use DERIV of the International Mathema-tical and StatisMathema-tical Library (IMSL) to perform numerical differentiations. DERIV first computes spline interpolants to the input functions (i.e., x, y and z, or f, l and h) and then differentiates the spline interpolants to obtain derivatives. For numerical differentiations, a user-supplied function for interpolations is needed. We use the lth degree Newton–Gregory forward polynomial (Gerald and Wheatley, 1994) to construct such a function. Given evenly spaced source data points (ti, fi), the value of
this function at a given epoch t is computed as
PnðtÞ ¼ f0þ Xl k¼1 s k Dkf0, (7) where s ¼ ðt t0Þ=d, (8) s k ¼sðs 1Þ ðs k þ 1Þ k! , (9) Dkf0 ¼fkþX k j¼1 ð1Þj k j ! fkj. (10)
Here, d is the sampling interval, i.e., ðtiþ1tiÞ. A
large degree (l) will introduce instability and a small degree will degrade the interpolation accuracy. According to our tests, l ¼ 14 yields the best result in terms of stability and interpolation accuracy. The method of numerical differentiation described above is chosen based on accuracy analysis in satellite gravimetry. By numerical integration of the equations of motion of a satellite orbit, we can compute satellite’s positions and velocities, see, e.g.,
Hwang and Lin (1998). The positions can be used to compute satellite velocities by numerical differentia-tion. Such velocities are then compared with the ‘‘true’’ velocities and the differences are of the order of 106m s1.Bruton et al. (1999)present a detailed analysis of numerical differentiation for velocity and acceleration, and subsequent use of low-pass filtering in different occasions.
3.2. Filtering of along-line values: a iterative Gaussian filter
In this paper, an iterative Gaussian filter is used to filter along-survey-line data values. Given data points ðxi; yiÞ; i ¼ 1; . . . ; n, where xiis abscissa and yi
is ordinate, the filtered values are computed as
^yk¼ X i¼i_ max i¼i_ min wk;iyi ,i¼i_ maxX i¼i_ min wk;i, (11)
where i_min and i_max are the lower and upper indices corresponding to a given filter width, s is 1/6 of the filter width (Wessel and Smith, 1995) and wk;i
is a weight function defined as
wk;i¼exp
ðxkxiÞ2
s2
The iterative filtering is done by first computing the standard deviation of the differences between the raw and the initially filtered data. If a difference exceeds three times of the standard deviation, the corresponding raw value is flagged as an outlier and its weight is reduced exponentially. Further filtering is carried out with newly assigned data weights, and the filtering stops when no more outliers are found. To see the performance of the Gaussian filter and the effect of filter width on the filtering results, we performed the following tests using a static, 60-km baseline between the National Chiao Tung Uni-versity (NCTU) GPS station and the Mt. Yangming GPS station. By holding fixed the coordinates of the Mt. Yangming station, the coordinates of the NCTU station were computed at a 1-s interval. The velocities and accelerations of the NCTU station were then computed using numerical differ-entiations (see Section 3.1). A total of 4300 sets of coordinates were determined. In theory, the velocity and acceleration of the NCTU station should be zero, and the standard deviations of the time series of velocity and acceleration are indicators of GPS positioning accuracy and filtering results. Table 1
shows the standard deviations at various filter widths. In Table 1, we assume that the aircraft velocity is 300 km h1, and the equivalent spatial resolution of airborne gravity data is just the multiplication of the aircraft velocity and the filter width. In general, with a filter width greater than 60 s, the standard deviations of velocities are at the 0.0001 m s1 level and the standard deviations of vertical accelerations are at the mgal level. A filter width larger than 200 s yields a sub-mgal standard deviation. This result is based on a static baseline. In an actual flight, the aircraft may experience turbulence and GPS data may contain cycle slips,
so the filter width to be used will vary from one case to another.
As discussed in Section 3.1, a long wavelength error may exist in the GPS-determined position along a gravity survey line. Such a long wavelength error poses little damage on velocity and accelera-tion, because differentiation will reduce or eliminate such a difference. This is also shown inTable 1. An example is as follows. Let the velocity be computed by the approximate differentiation vhDpDt, where
Dt ¼ 1 s (for an 1-Hz sampling interval), and Dp is the difference in along-line horizontal distance between two consecutive data points. Thus, in theory, an error of 1 cm in horizontal position will translate to an error of about 0.01 m s1 in horizontal velocity. However, the standard devia-tions of horizontal velocities in Table 1 are much smaller than 0.01 m s1(see the case with zero filter width). This implies that differentiation of position does reduce the effect of long wavelength error of position on velocity and acceleration. If one uses the phase approach of Jekeli and Garcia (1997) and
Kennedy et al. (2002) to compute velocity and acceleration, the effect of integer ambiguity is automatically removed because of differentiation of phase observable. Therefore, as far as the effect of integer ambiguity on velocity and acceleration is concerned, differentiation of position is equivalent to differentiation of phase observable.
3.3. Correction of gravimeter time
The gravimeter time associated with gravity read-ings comes from the clock of the computer attached to the gravimeter, and it will not be as accurate as the GPS time. Therefore, it is necessary to correct the gravimeter time. To find the correction of gravimeter
Table 1
Standard deviations of filtered north and east velocity components and vertical acceleration at a static point
Filter width (s) Equivalent spatial resolution (m) North velocity (m s1) East velocity (m s1) Vertical acceleration (m s2)
0 83 0.004257 0.004004 0.019871 60 5000 0.000175 0.000146 0.000036 80 6667 0.000143 0.000116 0.000023 100 8333 0.000120 0.000097 0.000016 200 16667 0.000068 0.000057 0.000005 300 25000 0.000051 0.000040 0.000003 400 33333 0.000042 0.000031 0.000002 500 41667 0.000037 0.000025 0.000001 600 50000 0.000033 0.000021 0.000001
time with respect to the GPS time, time series of raw gravity readings and vertical aircraft acceleration can be used (Olesen, 2003). The vertical aircraft accel-erations as determined from GPS are time-tagged to the GPS time. The lag between these two time series is the correction. For this purpose, a correlation analysis approach is used. First, the correlation function between two signals s1and s2is
z tð Þ ¼ Z 1
1
s1ð Þst 2ðt þ tÞdt, (13)
where t is the lag. The computation of (13) is carried out in the frequency domain as
zðtÞ ¼ F1ðS
1S2Þ, (14)
where S1and S2are the Fourier transforms of s1and s2, ‘‘*’’ is the conjugate operator, and F
1 is the inverse Fourier transform. The correction is defined to be the t value that produces the largest correlation value z(t). In practice, the z(t) values are determined at an even time interval. This time interval is the resolution of the correction and is the maximum of the sampling intervals of the GPS time and the gravimeter time.
3.4. Upward and downward continuation of gravity value
To evaluate the result of an airborne gravity survey, it is common practice to compare upward continued terrestrial gravity values with airborne gravity values at the flight altitude. For various geophysical applications, it is necessary to down-ward continue gravity values at the flight altitude to gravity values at sea level. In the space domain, the relationship between gravity value at an altitude of h above sea level and gravity value at sea level is
gðx; y; hÞ ¼ h 2p Z 1 1 Z1 1 gðx; Z; 0Þ h2þ ðx x2Þ þ ðy ZÞ2 3=2dx dZ; (15) where x and y are planar coordinates, g(x,y,0) is gravity value at sea level and g(x,y,h) is gravity value at h. In the frequency domain, this relation-ship becomes Gðkx; kyÞjz¼h¼e2ph ffiffiffiffiffiffiffiffiffiffi k2xþk 2 y p Gðkx; kyÞjz¼0, (16)
where kx and ky are spatial frequencies, Gðkx; kyÞjz¼0, and Gðkx; kyÞjz¼h are the Fourier
transforms of gravity values at sea level and at h, respectively. The relationship in (16) can be conveniently used for upward continuation (from
Gðkx; kyÞjz¼0 to Gðkx; kyÞjz¼h) and for downward
continuation (from Gðkx; kyÞjz¼h to Gðkx; kyÞjz¼0).
Since downward continuation is an ill-posed pro-blem, filtering must be applied:
Gðkx; kyÞjz¼0¼e 2ph ffiffiffiffiffiffiffiffiffiffik2 xþk2y p Gðkx; kyÞjz¼hSðkx; kyÞ, (17) where S(kx,ky) is a filter. Filtering of downward continued gravity values can also be carried out in the space domain. In this case, gravity values are first downward continued to sea level using Eq. (16), which are then filtered by a spatial filter such as the Gaussian filter.
Since gravity data cannot be given over an infinite domain, direct application of (15) is implemented over a finite domain and errors will occur. Errors in upward and downward continuations can be reduced by using the remove–restore procedure. In this procedure, reference gravity values at the needed level (sea level or an altitude of h) computed from a long wavelength gravity model such as EGM96 (Lemoine et al., 1998) are subtracted from the raw gravity values. The residual gravity values are then downward/upward continued to values at a new level. The values at the new level are added to the reference gravity values at this new level to yield the final gravity values.
In this paper, the program that does the down-ward/upward continuations is designed based on the frequency domain approach and the remove–r-estore procedure. Filtering of downward continued gravity values is done in the space domain, rather than in the frequency domain. For example, several filters in the GMT package (Wessel and Smith, 1995), such as the Gaussian, median, boxcar or mean filters, can be used for filtering. It is noted, there are many methods of upward/downward continuation in the literature, e.g., the methods of least-squares collocation (Moritz, 1980) and normal free-air gradient (Bayoud and Sideris, 2003). Cur-rently, our software implements only the Fourier transform method.
3.5. Crossover adjustment of gravity values
Crossover differences of gravity values from intersecting survey lines can be used to assess the quality of gravity values in an airborne gravity survey. Crossover adjustment of observed gravity values can reduce the effect of bias and drift due to both GPS and the gravimeter. For marine gravity
survey, a method of crossover adjustment can be found in Wessel (1989). Let ¯gq
r be the observed
gravity value at point r along survey line q; ¯gq r is
corrupted by a bias and a drift, plus random error, so that
¯gqr ¼gqrþaqþbqtqr þeqr, (18) where gq
r is the true gravity value, eqr is the random
error, aq and bqare the bias and the drift belong-ing to survey line q, and tq
r is the time at point r
relative to the beginning time of line q. At all intersecting points, observation equations can be formulated as
vklp þxklp ¼akþbktkpalbltlp,
k ¼ 1; i; l ¼ 1; m; p ¼ 1; n, ð19Þ
where xkl
p is the differenced gravity value at
cross-over point p pertaining to lines k and l, xkl p is the
residual, i+m and n are the numbers of survey lines and crossover points, respectively. Using a matrix representation, (19) can be expressed as
V þ L ¼ AX, (20)
where V, L and X are vectors containing residuals, observations and parameters (bias and drift), and A is the design matrix. The row vector formed at point p and lines k and l contain only four non-zero elements and it is easy to see that matrix A has a rank defect of order 2. To overcome this rank defect, at least one survey line must be held fixed. This line can be chosen as follows. For each candidate line, we first check the base reading, the
Fig. 2. (a) Given biases and drifts along north–south survey lines (beginning with NS) and west–east survey lines (beginning with EW), and (b) differences between given and estimated biases and drifts.
turbulence condition of the flight and the GPS positioning result. The gravity anomalies along this line are compared with upward continued ground data to see if any apparent bias and drift exist. The line with the least bias and drift is chosen as the fixed line.
In order to be flexible in selecting the fixed survey lines, we solve for X using
^
X ¼ ðATPA þ PxÞ1ðATPL þ PxLxÞ, (21)
where ^X contains the estimated biases and drifts, P is the weight matrix of differenced gravity values, Px is a diagonal matrix containing large weights (normally 1010time of the diagonal values of ATPA) pertaining to the fixed survey lines and zero diagonal elements elsewhere, and Lx contains the given bias and drift of the fixed survey lines and zeroes elsewhere.
Since no sufficient observed airborne gravity data were available, we tested the crossover adjustment program using simulated data. The existing ground gravity anomalies in Taiwan (see also Section 5.2.) were first upward continued to an altitude of 3000 m over the area 119.61–122.41E, 21.61–25.41N. In this area, gravity anomalies along 34 north–south survey lines and 21 west–east survey lines were interpolated from the gridded, upward continued gravity values at a 1-Hz sampling rate. The speed of flight is assumed to be 300 km h1. The cross-track intervals of the north–south and the west–east survey lines are 10 and 20 km, respectively. A total of 167 crossover points were formed in this case. The gravity values along each survey line were corrupted by a given bias and drift, plus random errors originating from a 1-mgal standard error.
Fig. 2 shows the given biases and drifts, and the differences between the given and the estimated biases and drifts. Table 2 shows the statistics of crossover differences before and after the crossover adjustment. From Fig. 2 and Table 2, clearly the crossover adjustment has successfully recovered the
biases and drifts and significantly reduced the crossover differences.
3.6. Correction of gravimeter position
The position determined by GPS refers to the GPS antenna phase center. In reality, the position of the gravimeter center (the position corresponding to the gravimeter reading) should be used in data reduction. First, a vehicle-fixed coordinate system is defined as follows (Fig. 3). The origin is at the antenna phase center, x is parallel to the fuselage and is positive to the flight direction, z is normal to x and positive to the zenith. Finally, y is normal to both x and z, and x, y, z form a right-hand system. Let the coordinates of the gravimeter center in the vehicle-fixed system be Dx, Dy and Dz. Ignoring the effect of pitch, yaw and roll of the aircraft, the geodetic coordinates of the gravimeter center ðfg; lg; hgÞ can be computed from the coordinates
of the antenna phase center ðfa; la; haÞ as
(Fig. 3) fg¼faþS cosðaabÞ Re , (22) fg¼faþ S sinðaabÞ Recos fa , (23) hg ¼haþDz, (24)
where Re is the Earth’s mean radius, S ¼pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiDx2þDy2, b ¼ tan1ðDy=DxÞ and aa is
theazimuth of the flight direction and is computed Table 2
Statistics of crossover differences of simulated airborne gravity values (in mgal) before and after crossover adjustment
Maximum Minimum Mean Standard deviation RMS Before 19.489 19.625 0.963 7.370 7.427 After 0.785 0.650 0.001 0.093 0.093 N GPS antenna a g gravimeter E ∆x ∆y S fuselage
Fig. 3. Geometric relationship between phase center of GPS antenna and gravimeter center.
by aa¼tan1 ve vn . (25)
4. Data reduction: software
Based on the theories presented above, we have developed computer programs in FORTRAN for data reduction in scalar airborne gravimetry. Program agp contains modules for interpolations and computations of velocities and accelerations from GPS-determined positions. Program agx
determines the locations of crossover points and performs crossover adjustment to estimate the bias and drift of each survey line. Program agr reads the output of agp and raw gravimeter readings, corrects the gravimeter time, corrects the gravimeter posi-tions, and finally generates filtered gravity values and gravity anomalies. Program agc reads gridded gravity anomalies and performs upward or down-ward continuation. In the case of downdown-ward continuation, agc does not filter the downward continued values, so filtering should be applied to the output of agc in a separate computation using a filtering program such as ‘‘grdfilter’’ of GMT
120° 120° 121° 121° 122° 122° 22° 22° 23° 23° 24° 24° 25° 25° NS14 CCK airport
GPS reference station at Shui-Nan
(Wessel and Smith, 1995). Appendix A shows sample jobs of data reduction using these programs and simulated and observed data in Taiwan.
5. A case study in Taiwan
5.1. Data
From May 2004 to May 2005, an airborne gravity survey is carried out to map the local gravity field around Taiwan. In this paper, the gravity and GPS data from a survey line (number NS14, collected on August 17, 2004) in this survey are used to assess the performance of the data reduction programs de-scribed above. Fig. 4 shows the locations of the survey line, the Ching-Chung Kang (CCK) airport (the airport for this survey) and the fixed GPS station for the differential positioning of the air-craft. An L&R Air–Sea Gravity System II gravi-meter (Serial number: S-130) onboard a King-Air Beechcraft 200 aircraft was used to collect the data. This gravimeter is owned by Ministry of the Interior, Taiwan. The flight altitude is approxi-mately 5150 m above sea level and the speed of flight is 300 km h1. The sampling rates of gravity reading and GPS reading are 1 and 2 Hz, respectively. The GPS receivers on the ground and onboard the aircraft are both Trimble5700. The coordinate components for position correction (Fig. 3) are Dx ¼ 2.0 m, Dx ¼ 0.0 m and Dx ¼ 1.5 m.
5.2. Results
Using the software developed in this paper, gravity values and gravity anomalies along line NS14 were computed. Fig. 5shows a time series of vertical aircraft accelerations and a time series of raw gravimeter readings. The two time series look similar, but a phase lag exists. Using the correlation method, we find that the gravimeter time is 30 s ahead of the GPS time. In order to compare the airborne gravity anomalies with terrestrial gravity anomalies, the surface gravity anomalies computed byHwang and Wang (2002)were upward continued to an altitude of 5150 m (orthometric height).Fig. 6
shows the distribution of surface gravity data and gravity anomalies at sea level and at 5150 m. Upward continuation is equivalent to low-pass filtering, so gravity anomalies at 5150 m are smooth compared to those at sea level. For the gravity anomaly computations, approximate orthometric heights of the aircraft are needed and they are
obtained by subtracting geoidal undulations implied by a Taiwan geoid model from the GPS-derived ellipsoidal heights of the aircraft.
Table 3 shows the statistics of the differences between the upward continued terrestrial gravity anomalies and the filtered airborne gravity anoma-lies at different filter widths of filtering. The RMS value of upward continued gravity anomalies along line NS14 is 11.5 mgal. Table 3 shows that the difference between the two sets of gravity anomaly decreases with increasing filter width. Both errors in the terrestrial and the airborne gravity anomalies contribute to the differences. As seen inFig. 6, the coverage of terrestrial and marine data is rather sparse along line NS14 and this incomplete coverage will introduce errors in the upward continuation.
-25000 -20000 -15000 -10000 -5000 0 5000 10000 15000 20000 25000 0 500 1000 1500 2000 2500 time (second)
vertical acceleration (mgal)
-15000 -10000 -5000 0 5000 10000 15000 20000 25000 30000 35000 0 500 1000 1500 2000 2500 time (second)
gravimeter reading (mgal)
-15000 -10000 -5000 0 5000 10000 15000 20000 25000 30000 35000 0 500 1000 1500 2000 2500 time (second)
gravimeter reading (mgal)
(a)
(b)
(c)
Fig. 5. (a) Vertical accelerations of aircraft, (b) raw gravity readings and (c) time-corrected gravity readings.
120° 120° 121° 121° 122° 122° 22° 22° 23° 23° 24° 24° 25° 25° 0 0 0 0 0 100 100 100 200 200 200 120° 120° 121° 121° 122° 122° 22° 22° 23° 23° 24° 24° 25° 25° -100 0 0 0 100 100 200 120° 120° 121° 121° 122° 122° 22° 22° 23° 23° 24° 24° 25° 25° (a) (b) (c)
Fig. 6. (a) Distribution of terrestrial and marine gravity data, (b) contours of gravity anomaly at sea level and (c) contours of gravity anomaly at 5150 m above sea level. Unit of contours is mgal.
The noises in the gravimeter readings and vertical accelerations (Fig. 5) may have corrupted the gravity signal to the extent that only gravity signal components beyond a certain wavelength can be recovered via filtering. Also, upward continuation is a smoothing process that will only retain gravity signal components at long wavelengths. Therefore, the outcome in Table 3 is what we have expected: the agreement between terrestrial and airborne gravity anomalies is better at long wavelengths than at short wavelengths (wavelength increases with filter width).
6. Conclusions
The goal of this paper is to develop methods for scalar airborne gravity data reduction. The numer-ical differentiation is highly accurate, delivering a commission error of the order of 106m s1for the velocity. The iterative Gaussian filter can effectively remove outliers and reduce the noise to the mgal level with a filter width large than 60 s. A cross-correlation technique is developed to correct for gravimeter time shift, and the accuracy of the computed time shift is the maximum of the sampling intervals of the GPS time and the gravimeter time. Our crossover adjustment method is flexible in that any line can be chosen as a fixed line, and the method can fully recover the biases and tilts in the simulated data set, and the accuracies of recovered parameters are subject to data noises. The Fourier transform method is used to upward continue ground gravity anomalies, which differ from filtered airborne gravity anomalies at the mgal level. Our software is freely available and can be improved by adopting new methods of downward/ upward continuation, filtering and other functions. Readers interested in the programs developed in this
paper can send requests to the authors at hwang@ geodesy.cv.nctu.edu.tw.
Acknowledgments
This research is funded by National Science Council of Taiwan, under project ‘‘Development of data reduction techniques in airborne gravime-try’’ (NSC-92-2211-E-009-047), and by Ministry of the Interior, Taiwan, under project ‘‘Airborne gravity survey in Taiwan’’. The authors are grateful to the pilots of King Air BE200 for their excellent flights.
Appendix A. Sample batch jobs of data reduction
(1) Computation of velocities and accelerations from determined positions: ‘‘agp’’ reads GPS-determined geodetic coordinates of the aircraft along survey line NS14 (ns14-p.txt) and compute velocities and accelerations (ns14-pva.txt). The starting and ending times are 11400 s and 13724 s of day 2 of the given GPS week.
agp Ins14p.txt –Ons14pva.txt S11400 E13724 -W2
(2) Reduction of gravity readings to gravity values and gravity anomalies: ‘‘agp’’ reads gravimeter readings (ns14-meter.txt), GPS data (ns14-pva.txt) and a Taiwan geoid model (twgeoid.grd3) and uses a 300-s window to filter data and to produce gravity values and gravity anomalies (ns14-gra.txt).
agr -Ins14-pva.txt –Gns14-meter.txt-Ntwgeoid.grd3-Ons14-gra.txt-W300-M1-T1-F2/0/-1.5-H1.9
(3) Crossover adjustment of along-line gravity values: ‘‘agx_simu’’ (a variant of agx, for simula-tion) reads simulated airborne gravity values (simulated_ag.txt) and fixed lines (fix-lines.txt), estimates biases and drifts (est_par.txt) and creates corrected gravity values (ag_corr.txt) and locations of crossover points (xover_loc.txt).
agx_simu simulated_ag.txt -Ffix-lines.txt -Gag_ corr.txt -Pest_par.txt-Xxover_loc.txt4agx.out
(4) Upward continuation: ‘‘agc’’ reads gravity anomalies (res_gravity.grd3) and upward continue them to values at 5156 m. ‘‘grd3toz’’ and ‘‘xyz2grd’’ convert the anomaly grid (res_gravity_5156.grd3) to a GMT .grd grid (res_gravity_5156.grd).
agc res_gravity.grd3-G res_gravity_5156.grd3-H5156 grd3toz res_gravity _5156.grd3|xyz2grd-Z-R119/123/21/26-I1 m-G res_gravity_5156.grd Table 3
Statistics of differences between upward continued and filtered airborne gravity anomalies at different filter widths
Filter width of filtering (s) Maximum (mgal) Minimum (mgal) Mean (mgal) Std. dev. (mgal) 200 10.3 20.5 1.8 3.7 300 6.3 6.1 1.4 2.1 400 5.2 4.1 1.7 1.4 500 4.2 2.6 1.7 1.2
References
Bayoud, F.A., Sideris, M.G., 2003. Two different methodologies for geoid determination from ground and airborne gravity data. Geophysical Journal International 155, 914–922. Bell, R.E., Childers, V.A., Arko, R.A., Blankeship, D.D.,
Brozena, J.M., 1999. Airborne gravity and precise positioning for geologic applications. Journal of Geophysical Research 104, 15281–15292.
Beutler, G., Bock, H., Brockmann, E., Dach, R., Fridez, P., Gurtner, W., Habrich, H., Hugentobler, U., Ineichen, D., Meindl, M., Mervart, L., Rothacher, M., Schaer, S., Springer, T., Urschl, C., Weber, R., 2004. Bernese GPS Software Version 5. Astronomical Institute, University of Bern, Switzerland 349pp.
Bruton, A.M., Glennie, C.L., Schwarz, K.P., 1999. Differentia-tion for high-precision GPS velocity and acceleraDifferentia-tion determination. GPS Solutions 2, 7–21.
Childers, V., McAdoo, A.D.C., Brozena, J.M., Laxon, S.W., 2001. New gravity data in the Arctic Ocean: comparison of airborne and ERS gravity. Journal of Geophysical Research. 106, 8871–8886.
Forsberg, R., Olesen, A.V., Keller, K., Moller, M., 2003. Airborne gravity survey of sea areas around Greenland and Svalbard 1999–2001. Technical Report 18, Kort & Matrikel-styrelsen. National Survey and Cadastre, Demark 55pp. Gerald, C.F., Wheatley, P.O., 1994. Applied Numerical Analysis,
fifth ed. Addison Wesley Pub. Co., New York 748pp. Goad, C.C., Yang, M., 1997. A new approach to precision
airborne GPS positioning for photogrammetry. Photogram-metry Engineering and Remote Sensing 63, 1067–1077. Hwang, C., Lin, M.J., 1998. Fast integration of low orbiter’s
trajectory perturbed by earth’s nonsphericity. Journal of Geodesy 72, 578–585.
Hwang, C., Wang, C.G., 2002. New gravity anomaly grid of Taiwan. Journal of Surveying Engineering 44, 1–22 (in Chinese). Jekeli, C., Garcia, R., 1997. GPS phase accelerations for
moving-base vector gravimetry. Journal of Geodesy 71, 630–639. Kennedy, S., Burton, L., Schwarz, K.P., 2002. Improving DGPS
accelerations for airborne gravimetry: GPS carrier phase
accelerations revisited. In: Adam, J., Schwarz, K.P. (Eds.), Vistas for Geodesy in the New Millennium, International Association of Geodesy Symposia, vol. 125, pp. 211–216. Springer, Berlin, 622pp.
LCR, 2003. L&R Model ‘‘S’’ Air-Sea Dynamic Gravity Meter System II, LaCoste & Romberg.http://www.lacosteromberg. com.
Lemoine, F.G., Kenyon, S.C., Factor, J.K., Trimmer, R.G., Pavlis, N.K., Chinn, D.S., Cox, C.M., Klosko, S.M., Luthcke, S.B., Torrence, M.H., Wang, Y.M., Williamson, R.G., Pavlis, E.C., Rapp, R.H., Olson, T.R., 1998. The development of Joint NASA GSFC and the National Imagery and Mapping Agency (NIMA) geopotential model EGM96. Report NASA/TP-1998-206861. National Aeronautics and Space Administration Greenbelt, MD 575pp.
Moritz, H., 1980. Advanced Physical Geodesy. Abacus Press, Tunbridge Wells, Kent, UK 500pp.
Olesen, A.V., 2003. Improved airborne scalar gravimetry for regional gravity field mapping and geoid determination. Technical Report 24, Kort & Matrikelstyrelsen. National Survey and Cadastre, Demark 55pp.
Schwarz, K.P., Li, Y.C., 1996. What can airborne gravimetry contribute to geoid determination? Journal of Geophysical Research 101, 17873–17881.
Schwarz, K.P., Li, Z., 1997. Introduction to airborne gravimetry and its boundary value problems. In: Sanso, F., Rummel, R. (Eds.), Geodetic Boundary Value Problems in View of the One Centimeter Geoid. Lecture Notes in Earth Sciences, vol. 65. Springer, Berlin 588pp.
Seeber, G., 1993. Satellite Geodesy: Foundations, Methods, and Applications. Walter de Gruyter, Berlin 531pp.
Torge, W., 1989. Gravimetry. Walter de Gruyter, Berlin 465pp. Verdun, J., Klingele, E.E., Bayer, R., Cocard, M., Geiger, A.,
Kahle, H.G., 2003. The alpine Swiss-French airborne gravity survey. Geophysical Journal International 152, 8–19. Wessel, P., 1989. XOVER: a cross-over error detector for track
data. Computers & Geosciences 15 (3), 333–346.
Wessel, P., Smith, W.H.F., 1995. New version of the Generic Mapping Tool released. EOS Transactions. American Geo-physical Union 76, 329.
