Topography-Assisted downward Continuation of Airborne Gravity : An Application for Geoid Determination in Taiwan

We investigated the Fourier transform-based downward continuation (DWC) of airborne gravity anomalies around Taiwan assisted by topographic information. The topographic data are from the latest collections of elevations and ocean depths. The DWC employs a remove-compute-restore (RCR) procedure in which the topography is removed prior to computation and then restored to achieve stable solutions. The topographic gravity effect is evaluated point-wise using the Gaussian quadrature. A Gaussian filter with an optimal smoothing parameter reduces the noise-amplifying effect of DWC. Use of topography in DWC leads to improvements of 3 to 6 mgal of gravity on land. Surface and downward continued gravity anomalies are used to determine geoidal heights by least squares collocation (LSC) in a similar RCR procedure through the same topographic data. The accuracy of the geoidal heights at Taiwan’s first-order leveling benchmarks is improved by 1 to 2 cm due to inclusion of downward continued airborne gravity.


InTroDuCTIon
With the advent of airborne gravimetry in Taiwan (Hwang et al. 2007), an optimal combination of airborne gravity and surface gravity (the latter is defined as any gravity measurements on the ground and at the sea surface by a static or a moving platform) for geophysical and geodetic applications is increasingly important.Airborne gravity data were collected above sea level and they require downward continuation (DWC) before merging with surface gravity.DWC is an ill-posed problem that amplifies data noises.Without a proper method, downward continued gravity data will contain false gravity signatures that lead to erroneous geophysical interpretations.Several DWC methods are found in the literature, e.g., those given by Wang (1997), Sun and Vaníček (1998), Li (2000), Sjöberg (2001), Novák and Heck (2002), Trompat et al. (2003), Cooper (2004), and Tziavos et al. (2004).These methods were largely derived from the Poisson integral (Heiskanen and Moritz 1967), and were presented in some suitable forms for numerical evaluations.It is also possible to derive formulae of DWC based on the boundary value problem of physical geodesy, as in Novák and Heck (2002).Numerically, these DWC methods can be classified into a space-domain approach and frequency-domain approach.Because of limited spatial coverage of airborne gravity, a good result for DWC can be achieved using the so-called remove-compute-restore procedure (RCR), where a geopotential model and topography information are used to model low-and high-frequency gravity variations in DWC, respectively (see section 6).Use of topography in DWC is particularly helpful around Taiwan because the gravity values are fast-varying, especially over the Central Range and the ocean trenches east of Taiwan.Methods to compute the topographic gravity effect in DWC can be found in, e.g., Forsberg (1985), Omang and Forsberg (2000), Sjöberg (2000), Nahavandchi and Sjöberg (2001), and Flury (2006).In fact, use of topography in DWC of gravity in Taiwan has not been investigated before, and it is not clear how topography will affect the DWC result here.Furthermore, any method of DWC will be eventually assessed in its rigor and efficiency by comparison between downward continued gravity and surface gravity.
An important geodetic application of a combined airborne and surface gravity is geoid modeling, which can be challenging over Taiwan due to its rough terrain and bathymetry that leads to a high-frequency variation of the geoid.Because the gravity data around Taiwan were collected by different gravimeter types and different platforms, they have varying spatial resolutions and accuracies (section 2), posing a substantial difficulty in optimally combining them for geoid modeling.Like DWC, either a space-domain or a frequency-domain based method can be used to determine a local geoid model [A good summary can be found in Sanso and Rummel (1997)].An early work on geoid modeling in Taiwan using Stokes' integral implemented by fast Fourier transform (FFT) is Tsuei (1995).Hwang (1997) and Hwang et al. (2006) used a least-squares collocation to determine geoid models of Taiwan.These results are based only on surface gravity data, and show a geoid accuracy of several cm over coastal plains, and 1 to 2 cm over high mountains.Geoid determination with airborne gravity over Taiwan was first carried out by Hwang et al. (2007), but the airborne gravity has been downward continued without using the topographic information.The objective of this paper is to investigate the role of topography in DWC around Taiwan, and to develop an improved geoid model by combining downward continued gravity and surface gravity data.Accuracy assessments of downward continued gravity and the modeled geoid will be presented.

GrAvITy AnD TopoGrAphIC DATA
In the remove-compute-restore procedure (section 6), the Bouguer anomaly was obtained from raw gravity minus the topographic gravity effect.Since most high-frequency variations of gravity are caused by a topographic effect, the Bouguer anomaly will be smoother than raw gravity leading to improved downward continuation.Virtually all existing gravity data gathered around Taiwan from the 1980s to early 2000s were used in this study.Figure 1 shows the locations of free-air gravity anomalies from these sources.Land gravity data (Fig. 1a) were collected by Academia Sinica (Yen et al. 1990;Yen et al. 1995), Base Survey Battalion of Taiwan, National Chiao Tung University and Ministry of the Interior (MOI) of Taiwan (Hwang et al. 2002).Part of the shipborne gravity data (Fig. 1b) were collected by National Central University, Taiwan, using a gravimeter R/Vl' Atalante KSS30 in 1996 (Hsu et al. 1998), and the rest were provided by the National Geophysical Data Center (NGDC), National Oceanic and Atmospheric Administration (NOAA), USA.In order to fill the data gaps in the existing ground gravity coverage (mainly in inaccessible areas), the Ministry of the Interior (MOI) of Taiwan sponsored an airborne gravity survey over the period from May 2004 to May 2005.
The project, including the field work and software development and was carried out by National Chiao Tung University, Taiwan, and National Survey and Cadastre (KMS), Denmark.The results are published in Hwang et al. (2007).The airborne survey area covered all of Taiwan and its surrounding seas (Fig. 1c).The airborne survey lines consist of 64 north-south, 22 east-west, 10 northeast-southwest, and 6 northwest-southeast trending lines with spacings of 4.5, 20, 5, and 30 km, respectively.A LaCoste and Romberg (LCR) Air-Sea Gravity System II (serial number: S-133) was used to collect the airborne gravity data at a 1-Hz sampling rate at an average altitude of 5156 m.The standard deviation of the crossover differences of gravity anomalies after a bias and drift correction is 2.88 mgal.The airborne gravity anomalies were compared with surface gravity data upward continued to the flight level, showing small differences (few mgals) at flat areas and large differences (over 10 mgals) at mountainous areas.As shown in Fig. 1b, there are still data gaps over the oceanic areas, and such gaps were filled by altimeter-derived gravity anomalies (Fig. 1d).In this paper, we chose to use the KMS02 gravity anomalies, which were given on a 2 2 # l l grid.The KMS02 gravity anomalies were determined from altimeter data from the repeat and geodetic missions of GEOSAT and ERS (Andersen et al. 2003).The DNSC08 gravity grid was not used here because it is not derived purely from altimeter data (Andersen et al. 2008).
The topographic data will be used to aid the DWC and geoid determination.The elevations on land were made available on a 3 3 # m m grid (with a horizontal resolution of approximately 80 m), which was constructed using the photogrammetric method by the Aerial Survey Office of the Forest Bureau, Taiwan (Hwang et al. 2003a).The accuracy of the 3 3 # m m elevation grid is approximately 4 m as determined by comparing with the orthometric heights at hundreds of first-order leveling benchmarks.The ocean depths around Taiwan were from ETOPO1, which is a 1 1 # l l bathymetry database generated by National Geographic data Center of USA (NGDC).Elevations on land and at sea from these different sources were decimated into three coarser grids sampled at 9 9 # m m, 90 90 # m m, and 6 6 # l l intervals, which will be used to compute the topographic gravity effects for DWC and geoid modeling.

MeThoD for DownwArD ConTInuATIon
Several methods can be used for DWC of airborne gravity values.There are also techniques for stabilizing DWC and for reducing the noise-amplifying effect.The choice of an optimal DWC method and an optimal stabilizing technique is not the subject of this paper.Rather, we will select one method of DWC to experiment with the use of topography in DWC.We expect that the conclusions on the topography effect on DWC drawn from this study will apply to parallel methods.
The basic relationship between the gravity value at an altitude of h, g (x, y, h), and the gravity value at sea level, g (x, y, 0), is (Heiskanen and Moritz 1967;Buttkus 2000) which can be expressed as the convolution g * , , , , , g x y h w x y is the impulse response function for upward continuation (UPC) from the z = 0 plane to the z = h plane.Its Fourier transform is = + is the radial (isotropic) frequency.Thus, in the frequency domain, Eq. ( 1) becomes where are the Fourier transforms of gravity at h and at sea level.Inverting Eq. ( 5) and applying a filter S G (f), we obtain the expression for DWC as Without a filter DWC is an operation that amplifies highfrequency components, leading to an unstable result.In this paper, we choose to use the isotropic Gaussian filter for DWC, defined as where k G is a positive number governing the smoothness of the filter.The unit of k G is to make the result k G f r dimensionless.Also, the degree of smoothing increases with the value of k G .

MeThoD for CoMpuTInG TopoGrAphIC GrAvITy effeCT
The role of topography in DWC is to provide the highfrequency gravity component in gravity.In this paper, the computation of the topographic gravity effect is based on the Gaussian quadrature (Hwang et al. 2003b).Using a planar approximation, the topographic gravity effect at a point (x p , y p , h p ) can be computed as , , ( , , ) where G is Newtonian constant and t is density, X 1 , X 2 , Y 1 , and Y 2 are the geographic borders of the effective area of the topographic mass, and

R x x y y h h
In Eq. ( 9), the kernel function f (x, y, h) is a function of the coordinates and the elevation of the contributing topographic mass.Using the Gaussian quadratrure, the integral in Eq. ( 8) can be numerically evaluated as where f ij is the value of the kernel function at the nodes along longitude (for index i) and latitude (for index j), and w i , v j are the weights of the quadrature.These nodes are determined according to the sampling interval of elevation.A detailed description of the Gaussian quadrature in a one-dimensional case is offered by Press et al. (1989) and is extended to a two-dimensional formula found in Eq. ( 11), where the contribution along each latitude belt is first evaluated, followed by summing such contributions from the southernmost latitude belts to the northernmost most belts.
At the remove stage of the RCR procedure (section 6), h p is the flight altitude; at the restore stage, h p is the elevation at the land surface whose geodetic latitude and longitude are identical to those of the airborne gravity measurement.On land, the density in Eq. ( 8) was set to 2.57 × 10 3 kg m -3 , which is the average density of rock in Taiwan.The density in Eq. ( 8) is the density contrast between rock and sea water (about 1.03 × 10 3 kg m -3 ).For a computational point at the flight altitude [i.e., p in Eq. ( 8)], the contribution of a column of land mass, counting from a mean sea level (or the geoid) to the land surface along the plumb line, is positive, while the contribution of a column of oceanic mass is negative.If the computational point is at the land surface (with elevation = h p ), a column of land mass makes a positive or negative gravity contribution depending on its elevation is larger or smaller than 2h p .

MeThoD for GeoID MoDelInG: leAsT-squAres ColloCATIon
The downward continued gravity anomalies, along with surface gravity data, are to be used for geoid determination.In this study, least-squares collocation (LSC) is used as follows.LSC can accept hybrid data with different noises and spatial resolutions to estimate any functional of the earth's gravity field.The basic formula of LSC is (Moritz 1980) where s and l are vectors of signals and observations, respectively, C ll is the covariance matrix of l and C sl the covariance matrix between s and l, and D is the covariance matrix of the contained in the vector l, which functions as a filter.In the context of this paper, s, l, and D in Eq. ( 12) represent geoid undulation, gravity anomaly and error variance of gravity anomaly, respectively.C sl and C ll denote the covariance matrices for geoid-gravity anomaly and gravity anomaly-gravity anomaly.In addition, C sl and C ll were both determined by using the combination of EGM08 geopotential model up to degree 360 and a degree variance model of Tscherning and Rapp (1974).
For geoid modeling using LSC, various covariance functions are needed.The "origin" of all covariance functions is the covariance function of the earth's anomalous potential, which can be represented in a Legendre series as: where spherical angle between the radial vectors to P and Q.The needed covariance functions are constructed using the law of covariance propagation (Moritz 1980).In this paper, the following covariance functions are needed: (1) The covariance function between a gravity anomaly at P and gravity anomaly at Q: (2) The covariance function between the geoidal height at P and gravity anomaly at Q: In Eqs. ( 14) and ( 15), P c is the normal gravity at P. These covariance functions were tabulated at a 0.01° interval.The actual covariance values were determined by linear interpolations from the tables.

The reMove-CoMpuTe-resTore proCe-Dure
As a standard practice in many fields of physical geodesy, we use the RCR procedure for DWC and geoid modeling.As shown in Fig. 2, there are two cases of DWC: one that removes the long-wavelength part of gravity only (Case A), and the other removing both the long-wavelength part of gravity and the topographic gravity effect (Case B).The long-wavelength part is attributed to a geopotential model.In Fig. 2, g  level than that for the remove stage) were then added back to the downward continued residual values.
In the RCR procedure for geoid modeling, the longwavelength component and short-wavelength components are first removed from the raw gravity anomalies to obtain residual gravity anomalies as The long-wavelength component, gl D , is computed from a geopotential model, and the short-wavelength component, gs D , is evaluated using the mass between the actual land surface and a mean surface, the later called the "residual terrain model (RTM)" (Forsberg 1985).The residual gravity anomalies were then used in Eq. ( 12) to compute the residual geoidal heights.Finally, the long-and short-wavelength components of the geoidal height based on the same geopotential model and the same RTM were restored to obtain the final geoidal heights.As in DWC, the geopotential model EGM08 truncated at degree and order 360 is used to compute the long-wavelength effects for both DWC and geoid modeling (see also section 5).The error degree variances of EGM08 enter Eq. ( 13).In addition, the RTM gravity and geoid effects were estimated by FFT (Schwarz et al. 1990).Two elevation grids with different resolutions -9 9 # m m, 90 90 # m m, were used to compute the RTM gravity effects for the inner zone and the outer zone.The residual elevation is with respect to a mean elevation surface modeled by a 6 6 # l l elevation grid.

ACCurACy AssessMenT of DwC GrAvITy AnoMAlIes
Downward-continued gravity anomalies were computed over the area 119.5 -122.1°E and 21.8 -25.5°N on a 2 2 # l l grid for Cases A and B (Fig. 2).We experimented with several choices of k G in the Gaussian filter [Eq.( 7)], and it was found that the best accuracy of downward continued gravity anomalies, based on comparisons with surface gravity data, was achieved when using k G = 10 km. Figure 3a shows the airborne gravity anomalies at 5156 m, which is highly correlated with elevation and bathymetry (Fig. 3b).Figures 3c and d show the topographic gravity effects at 5156 m and at sea level; the statistics of these effects are given in Table 1.In general, the topographic gravity effect decreases with elevation.The range (maximum value -minimum value) of the topographic gravity effect at sea level is about 682 mgal, compared to the maximum free-air gravity anomaly of about 350 mgal in Taiwan.The airborne gravity anomalies over flat regions are moderate; however, their magnitudes become large over the Central Range and certain areas east of Taiwan, reaching extreme values of approximately 250 and -150 mgal, respectively.Figure 3e shows the downward-continued gravity anoma-lies from Case A. Compared to Fig. 3a, the signatures of the downward-continued gravity anomalies are more prominent.After DWC, the extreme gravity values increase to 300 and -200 mgal.The differences in gravity anomaly between Cases A and B are given in Fig. 3f.The largest differences occur over high mountains, reaching approximately 50 mgal.In general, the difference increases with the absolute magnitude of gravity anomaly, and both the largest difference and the largest magnitude of gravity anomaly occur over the Central Range.
The downward continued gravity anomalies were then compared with the surface gravity data given in Figs.1a and b over the same area where DWC was carried out.We employed the GMT package (Wessel and Smith 1995) to interpolate the land and shipborne gravity values on a 2 2 # l l grid, and subsequently filter the grid at a wavelength of 3 km by the Gaussian filter.Figure 4 shows the differences between the downward continued and surface (land and shipborne) gravity anomalies on 2 2 # l l grids.Each black dot in Fig. 4 represents a grid.In general, on land (elevations > 0 ) both the absolute differences in Cases A and B increase with elevation.In particular, the differences are highly correlated with elevation over areas with 0 < elevation < 1000 m.However, there is virtually no correlation between the gravity differences and the depths at sea (elevations < 0 ).Table 2 shows the statistics of the differences between the surface and downward continued gravity anomalies.At sea, the standard deviations of the differences in Cases A and B are almost identical.Thus, the depth data in fact do not improve the accuracy of DWC.However, the advantage of using topography becomes evident on land: compared to Case A, the standard deviations and the means of the gravity differences are significantly smaller in Case B (with the topographic data).The best improvement of DWC due to topography occurred over areas with 1 km < elevation < 2 km.In summary, use of topography in DWC has led to improvements of 3 to 6 mgal on land, but no improvement at sea.

ACCurACy AssessMenT of GeoID MoDel
In this section, two cases of geoid modeling were carried out.Case 1 uses land, shipborne and altimeter-derived gravity anomalies (Fig. 1).Case 2 uses additional the downward continued gravity anomalies from Case B (section 7).For geoid modeling by LSC, we assigned data noises of 0.1, 1.0, 5.0, and 3.0 mgal to land, shipborne, altimeter-derived, and airborne gravity anomalies, respectively.Figure 5a shows the geoid model from Case 2. Figure 5b shows the difference between the geoid models of Cases 1 and 2. On land, large differences occurred over high mountains and offshore areas, reaching a maximum of about 30 cm.At sea, large differences occurred off the southeast coast.These large differences are mainly due to the airborne grav-   ity anomalies that were excluded in Case 1 and included in Case 2.Over high mountains (e.g., the Central Range) and areas with sparse surface gravity data (e.g., central Taiwan Strait), the differences were prominent.This is because over these areas the airborne gravity anomalies help to pick up high-frequency gravity signatures which in turn enhance the spatial resolution and the accuracy of the modeled geoid.
The "observed" geoidal heights at selected GPS/leveling benchmarks (Fig. 5a) were compared with the modeled geoidal heights.At theses benchmarks, 3-h GPS observations were collected to compute ellipsoidal heights accurate to 1 ~ 2 cm, and their orthometric heights were determined by precision leveling accurate to a few mm.An observed geoidal height is the difference between the ellipsoidal height and the orthometric height.Figures 6a and b show the differences between the observed and the modeled geoidal heights from Cases 1 and 2, respectively, and the statistics of the differences are given in Table 3.For both Cases, the geoidal differences increase with elevations.At elevations < 2 km, the mean and standard deviation of the difference are around -0.250 and 0.085 m in Case 1, and -0.120 and 0.080 m in Case 2. However, they reach -0.443 and 0.121 m in Case 1, and -0.298 and 0.115 m in Case 2 at elevations > 2 km.In other words, the dependence upon differences in height becomes evident over areas with elevation > 2 km. Figure 7 shows the improvement percentages of geoid at the benchmarks when the airborne gravity anomalies were used (Case 2) or not used (Case 1).An improvement percentage at a benchmark is computed as   < 500 m, the improvement percentages fluctuate between negative and positively values.Over areas with elevation > 500 m, the airborne gravity anomalies do enhance the geoid accuracy, despite some few negative improvement percentages.The individual averages of improvement percentages in Fig. 7 are 30.5, 29.8, 36.5 at heights of 0 -1000, 1000 -2000, and > 2000 m, respectively.In conclusion, compared to Case 1, the geoid accuracy in Case 2 is improved by 1 to 2 cm over both high mountains and flat regions.

ConClusIon
The primary contribution of this paper is to show the value of topography in improving the accuracy of DWC,   and also to present a Taiwan geoid model and its accuracy assessment.The improvement of DWC due to topography is more evident over high mountains than flat regions.At sea, no improvement was found.The downward continued gravity anomalies improve the geoid accuracy by 1 to 2 cm, and the improvement increases with elevation.The geoid model (Case 2, section 7) obtained in this study is expected to aid in GPS leveling, ocean circulation determination east of Taiwan, and linking Taiwan's height datum based at Keelung Harbor to the height data of the offshore islands around Taiwan.The gravity data used in this study are from a variety of sources.These gravity data have varying spatial resolutions and noise levels.It is possible to derive a more rigorous method for combining such data than the current method in this paper.For example, an optimal combination can be based on frequency-dependent weighting of data.The method for geoid modeling can also be improved using, e.g., band-limited covariance functions in LSC or band-limited kernel functions in Stokes' integral (Novák and Heck 2002).With more gravity data in Taiwan to be collected in future campaigns, such improvements will be highly desirable and are important subjects for future investigations.
Fig.2.The remove-compute-restore procedures of DWC without (Case A) and with (Case B) topographic gravity effect.The long-wavelength gravity contribution is considered in both cases.

Fig. 3 .
Fig. 3. (a) airborne free-air gravity anomalies at 5156 m, (b) terrain and bathymetry around Taiwan, (c) topographic gravity effect at 5156 m, (d) topographic effect at sea level, (e) downward-continued free-air anomaly from case A, and (f) difference of downward-continued free-air gravity anomalies between cases A and B. The horizontal color scale is for gravity anomaly and the vertical one for elevation.

Fig. 6 .
Fig. 6.Differences between the observed and the modeled geoidal heights from (a) Case 1 and (b) Case 2.

Fig. 7 .
Fig. 7. Improvement percentages of geoid model due to use of downward continued gravity.

Table 1 .
Statistics of topographic gravity effects (in mgal) at 5156 m and sea level over 119.5 -122.1°E and 21.8 -25.5°N.

Table 2 .
Differences (in mgal)between downward continued and surface gravity anomalies.

Table 3 .
Differences (in meter) between the observed and modeled geoidal heights at first-order leveling benchmarks.