Improvement of downward continuation values of airborne gravity data in Taiwan

An airborne gravity survey was carried out to fill gaps in the gravity data for the mountainous areas of Taiwan. However, the downward continuation error of airborne gravity data is a major issue, especially in regions with complex terrain, such as Taiwan. The root mean square (RMS) of the difference between the downward continuation values and land gravity was approximately 20 mGal. To improve the results of downward continuation we investigated the inverse Poisson’s integral, the semi-parametric method combined with regularization (SPR) and the least-squares collocation (LSC) in this paper. The numerically simulated experiments are conducted in the Tibetan Plateau, which is also a mountainous area. The results show that as a valuable supplement to the inverse Poisson’s integral, the SPR is a useful approach to estimate systematic errors and to suppress random errors. While the LSC approach generates the best results in the Tibetan Plateau in terms of the RMS of the downward continuation errors. Thus, the LSC approach with a terrain correction (TC) is applied to the downward continuation of real airborne gravity data in Taiwan. The statistical results show that the RMS of the differences between the downward continuation values and land gravity data reduced to 11.7 mGal, which shows that an improvement of 40% is obtained.


Introduction
There are many regions in the world where airborne gravity surveys have been carried out [1][2][3][4][5][6][7][8][9][10][11][12][13][14]. The Taiwan airborne gravimetry project was conducted over many hills and high mountains. In order to generate the best results, many issues were identified in data processing, such as kinematic positioning error, error models of measurement system, downward continuation error, etc., by Hwang et al. [3]. Among those issues, the downward continuation error is the biggest factor affecting the geodetic and geophysical applications of airborne gravimetry. In airborne gravimetry, gravity anomalies are collected above the mean sea level and commonly downward continuation to the earth surface before exploiting it in geodesy and geophysics is needed [15][16][17][18][19][20][21][22][23][24]. Without a proper method, downward continued gravity data will contain false gravity signals that would lead to erroneous geophysical interpretations [1].
Several downward continuation methods were studied in the literature [15][16][17]. On the one hand, the inverse Poisson's integral is a classical way of downward continuation [17]. The grid spacing and downward continuation height are critical factors in applying such an approach [1]. A direct approach, based on band limitation, was proposed by Novak and Heck [18] because band limitation has a significant positive impact on the stability of numerical evaluation of downward continuation values. Numerically, downward continuation methods can be classified into space-domain and frequency-domain modes. Almost all of these methods are derived from the Poisson's integral. In addition, the LSC method has been applied to the downward continuation of airborne gravimetry by Forsberg et al. [2]. The comparison of Alberts and Klees [15] using synthetic data showed that the LSC method performs slightly better than the Poisson's integral-based approaches in terms of RMS errors. The terrain effect always should be considered as the main error impacting the downward continuation in mountainous areas [4,5]. The residual terrain modeling (RTM) method is a classical approach for TC in the downward continuation [4,5]. In total, all the mentioned problems compose the contemporary issue of downward continuation of airborne gravity. The RMS of the differences between the downward continuation values and land gravity in Taiwan is about 20 mGal [3]. It should be noted that the main research objective in Hwang et al. [3] is not the downward continuation. The RMS given by Hsiao and Hwang [23] is 18.4 mGal by the Fourier transform method, with removal of the topographic gravity effect. Therefore, the downward continuation of airborne gravity in Taiwan deserves further study.
In this paper, we aim to improve the downward continuation of airborne gravity data in Taiwan. We investigated the inverse Poisson's integral, SPR and LSC methods based on synthetic data in the Tibetan Plateau. The Tibetan Plateau experiments reduce the synthesized gravity anomalies at 5500 m to gravity anomalies at sea level, which is similar to the downward continuation height of 5156 m in Taiwan. The outline of this paper is as follows. In Section 2, the 3 approaches are formulated and the function model is defined. Section 3 shows and analyzes the numerical results from the Tibetan Plateau and Taiwan. Section 4 discusses the findings of this research and Section 5 contains the main conclusions.

Methods
In Section 1, we showed with the literature that the inverse Poisson's integral, the LSC, and the SPR methods all perform well for downward continuation. Therefore, we chose the three methods for our study. The basic principles of these methods will be described in the following.

The Inverse Poisson's Integral Method
Downward continuation can be viewed as the inverse operation of the Poisson's integral. In the spherical approximation, the upward continuation Poisson's integral yields [17,19]: where ∆g air (r, ϕ, λ) is the gravity anomaly at a field point, r, ϕ, λ are the spherical coordinates of the computation point, r = R + h, h is the height of upward continuation, and R is the integration radius. ∆g g (R, ϕ , λ ) is the gravity anomaly at the variable integration point on the geoid. l(r, ϕ, λ; R, ϕ , λ ) is the Euclidean distance between the fixed field point and the variable integration point. The distance l between the field point and the source point is calculated by l = r 2 + R 2 − 2rR cos ψ, where ψ is the spherical distance between the two points on the surface of a unit sphere, represented by (ϕ,λ), and (ϕ ,λ ), so that cos ψ= sinϕ sin ϕ + cos ϕ cos ϕ cos(λ − λ ) Equation (1) is called the Poisson's integral formula. Normally, the discrete Poisson's integration is used for the downward continuation. We use the matrix-vector notation for a grid representation of the Poisson's integral model which can be written as follows: where ∆g air is a vector of dimension M of the gravity anomalies measured in the air. ∆g g is a vector of length N of the gravity anomalies, and the M × N matrix B consists of elements b ij explained below. The diagonal entries of B are given by: where ∆σ j is the surface element centered at the block σ j . N C is the number of data within the spherical cap of radius 1 • , ψ 0 is the radius of the inner zone. The off-diagonal elements of B are given as follows:

The SPR Method
Based on the inverse Poisson's integral, a regularization parameter is added to suppress high-frequency errors. Additionally, the SPR has been developed to estimate the systematic error during the downward continuation process [1].
The nonparametric relation becomes: where L i is the observation including the error. b i is the ith row elements of matrix B. X is the estimated value. s(t i ) is a systematic error which is the nonparametric function. ∆ i is a random error vector. We consider as s(t), and as an approximation to the original systematic error function, the cubic smoothing spline, which for a given α S ≥ 0, minimizes the following relationship.
s (t) is the second order derivative of s(t). We call this the semi-parametric regression model or the "penalized likelihood" estimation. The first term of Equation (7) penalizes the lack of goodness of fit of the function to the data, and the second term penalizes the lack of smoothness of the approximating function. The solution to Equation (7) is unique at every data point t i . More importantly, the semi-parametric model can estimate the systematic error without the need for external data (such as upward continuation of surface gravity data [3]), because the "penalized likelihood" estimation allows the data "to speak for itself" without any a priori knowledge. The key point is to use α S to balance the two terms relative to each other. By adjusting the parameter α S , the smoothness of s(t) is varied. At the extremes, when α S goes to infinity, the left-hand side of Equation (8) is forced to be linear over the whole range of t values and becomes the best least squares line through the data. When α S goes to infinity, y tends to be an interpolating function for the data, fitting every data point exactly.
The generalized cross validation (GCV) method of the semi-parametric model is used to calculate the smoothing parameter α S . Based on Equations (6) and (7), the formulas are derived as follows: where M = (P + α S R s ) −1 P. R s is a semi-parametric matrix as in [1]. The Tikhonov regularization is a classic regularization method. Minimizing the Tikhonov cost function for the downward continuation yields: where α R is the regularization parameter and I is the unit matrix. The regularization parameter in this study is determined by GCV.
For the purposes of dealing with the systematic and random errors, we combine the two methods for the downward continuation: The semi-parametric model is used to estimate the systematic errors; then, we apply the regularization method to suppress the random errors.
The final equation is as follows. See Zhao [1] for more details.
where R ≤ r.
The major processing steps are as follows: (1) R s is generated by the natural cubic splines; (2) R s and the initial value of α S are used to calculate α S ; (3) R s and α S are used to estimate the systematic errorŝ; (4) The airborne gravity anomalies without the systematic errors can be obtained by subtractingŝ from the airborne gravity anomalies δg air ; (5) The airborne gravity anomalies without the systematic errors are used to calculate α R ; (6) Finally, the ground gravity anomalies are obtained.

The LSC Model
The classical methods of downward continuation also include the LSC model [2,4]. Provided that the auto-covariance function of gravity anomalies is known, the classical collocation formula iŝ whereŝ is the signal to be determined and x is the vector of observations. In this case, the signal is the disturbing potential at ground level, and the measurements are the gravity anomalies at flight level. C XX and C SX are the signal auto-covariances and signal cross-covariances, respectively, which are taken from an analytical covariance model. The matrix D describes the error covariances, which can be computed by using the prior noise information [15].
The LSC downward continuation model is expressed as where D nn is the (diagonal) noise matrix. The C hh is the covariance matrix following Kaula's rule. ∆g h is the air gravity anomaly. ∆g 0 is the ground gravity anomaly.
The covariance model used in this study is the attenuated planar logarithmic covariance model of Forsberg [35], which was designed for airborne gravimetry. The auto-covariance function for gravity anomalies at points P and Q at flight level h is given as where C h is the variance of gravity anomalies at flight level h. The parameters a i , z i and r i are given as where D i is the high-frequency attenuation parameter and T is the low-frequency attenuation parameter. See [35] for more details.

Simulation Data in Tibetan Plateau
To test the performance of three methods above, a simulation study on the Tibetan Plateau was carried out using gravity data simulated from the global geo-potential model (GGM) EGM2008 [36]. The test area is located from 27 • N to 31 • N latitude and from 92 • E to 93 • E longitude. For the study area, where surface gravity data are not available, the airborne gravity anomalies and the ground gravity anomalies were obtained using the EGM2008 from degree/order 360 to 2160. This was enough to roughly approximate the true field, which, nevertheless, was consistent both with the topography and at the flight height. EGM2008 is, of course, not a perfect approximation of the true field for the test region. Gravity anomalies were simulated at 5.5 km and 0 km (sea level) altitudes. Figure 1 shows the distribution of airborne gravimetry lines in the Tibetan Plateau. The simulated error data consist of random and systematic errors. The standard deviation (STD) of the former is 2 mGal, and the latter contains a 3 mGal bias, drift and periodic errors, the rate of drift is 0.002 mGal/s. The frequency of the periodic error is 0.003 Hz, and the amplitude is 6 mGal. The sine function is used to represent the periodic error. The errors in airborne gravimetric surveys can adequately be described as a linear combination of the random error, the bias, the linear drift, and the periodic error. These parameters are obtained based on an analysis of the relationship between the gravity measurement system and systematic errors.

Real Airborne Gravity Data in Taiwan
The Taiwan airborne gravity data was acquired by the National Chiao Tung University [3]. The Taiwan airborne gravimetry was conducted over Taiwan using a LaCoste and Romberg (LCR) System II air-sea gravimeter equipped with a with global positioning system (GPS) data sampled at 1 Hz [3]. This gravimeter has a nominal resolution of 0.01 mGal and an accuracy of better than 1 mGal for the gravity observations. There are four different survey lines with respect to four directions (e.g., north-south, east-west, northeast-southwest, and northwest-southeast), which are shown in Figure 2 (provided by Cheinway Hwang through private contact), the limits of which are from 21.7 • N to 25.5 • N in latitude and from 119.5 • E to 122 • E in longitude. The numbers of survey lines in the four directions are 64, 22, 10, and 6, respectively. The cross-line spacing is 4.5 km for all survey lines, except the west-east lines, which were spaced at 20 km. The west-east lines were used mainly for crossover analyses, which reduces the mis-fitting of the observed gravity signal that could be mistakenly considered as colored noise [24]. The cross-line spacing of 4.5 km equals to the theoretical resolvable wavelength of the gravity anomaly at the flight altitude, which was approximately 5156 m. The RMS of crossover differences was 2.88 mGal, corresponding to a 2 mGal STD error in terms of the gravity anomaly. Figure 3 shows the airborne gravity anomalies and land gravity anomalies in Taiwan. The land gravity anomalies were used to validate the airborne downward continuation values. The background of Figure 3 is the topography of Taiwan. According to Figure 3, the maximum values of gravity anomalies are located over the high mountain area. The minimum values are located in the east trench of Taiwan. Visually, the spatial distribution in Figure 3a shows that there is a good fit between the airborne gravity anomalies and the topography of Taiwan. The details of the statistics of airborne gravity anomalies and land gravity anomalies in Taiwan are given in Table 1. According to Table 1, the STD of airborne gravity anomalies and land gravity anomalies are 64.3 mGal and 67.1 mGal, respectively, which means that airborne gravity anomalies are smoother than land gravity anomalies.

Numerical Simulation Test in Tibetan Plateau
In the experiments, we downward continue the gravity anomalies from approximately 5.5 km to 0 km (sea level) based on the three methods described in Section 2. Before LSC is used, the parameters of the covariance function must be estimated. The parameters of LSC in the Tibetan Plateau are determined by GPFIT of GRAVSOFT [27]. For the SPR method, the smoothing parameter is 1.007, and the regularization parameter is 0.504. These parameters are generated by GCV method. The differences between downward continuation values and gravity anomalies are shown in Table 2.
According to Table 2, the inverse Poisson's integral method yields downward continued gravity anomalies containing larger random and systematic errors because the corresponding STD is the largest one and the mean value is 4.4 mGal. The SPR method is more efficient than the inverse Poisson's integral. The mean value 4.8 mGal illustrates that the results from LSC are impacted by systematic errors. However, the LSC generates the best result in terms of RMS.

Downward Continuation of Airborne Gravity Data in Taiwan
The LSC approach can generate the best results in the Tibetan Plateau. The topography of the Tibetan Plateau is as rough as it is in Taiwan. In addition, the inverse Poisson's integral and the SPR methods are usually applied in grid mode, not in a point-to-point mode [3]. The LSC can downward continue values in point-to-point mode. Therefore, we chose LSC for the downward continuation of airborne gravity in Taiwan. Before the LSC was used, the parameters of the covariance function had to be estimated. The parameters in Taiwan were set according to the parameters fixed by GPFIT from GRAVSOFT [37]. The Satellite only model GO_CONS_GCF_2_DIR_R5 (DIR_R5) [38] was used to remove-restore the long wavelength information in airborne gravity and surface gravity.
To reduce the downward continuation error, we also remove the TC effects in gravity anomalies. According to Forsberg [5], only the short wavelength variations of TC effects is taken into account, which is done by choosing a smooth mean elevation surface, and computationally remove masses above this surface and fill up valleys. This is called the RTM terrain reduction method [5], which is used to compute the short wavelength part of the TC effects in this paper by using the prism integration, and assuming a constant density of 2.67 g/cm 3 for all topographic masses. Considering the task of the computations, RTM terrain effects were performed in 0.5 • blocks and expanded with 0.3 • overlaps. The 30" × 30" Shuttle Radar Topography Mission (SRTM) data is used here, and the data processing is same to the one used in [4,6]. After this step, the downward continuation values of airborne gravity in Taiwan are estimated by GPCOL1 in GRAVSOFT [35,37].
The differences between downward continued values and land gravity data are shown in Table 3. The RMS of the differences corresponding to the LSC without TC is 27.881 mGal, and is reduced to 11.653 mGal when TC was applied. Hwang et al. [3] studied the downward continuation using Fourier transform and LSC, and reported that the Fourier transform method delivered the optimal results, which are also shown in the Table 3. The RMS is 18.4 mGal when applying Fourier transform method, which removed the topographic gravity effect [23]. According to Table 3, the improvement of the downward continuation results from LSC is mainly related to the TC effects. Therefore, we analyze the relationship between the heights of land gravity points and the downward continuation errors from the LSC approach with and without the TC applied (see Figures 4 and 5).

TC Analysis
To illustrate why the TC helps downward continuation, the statistical results of both the airborne gravity anomalies and land gravity anomalies after removing the DIR_R5 model up to d/o 220 degree [38] and the TC effects are reported in Table 4 and shown Figure 6. According to Table 4, the Max/Min/Mean/STD values before and after removing the effects of the GGM and TC are of the same order of magnitude but lower in value. Compared with the values in Table 1, after removing GGM, the STD of the free-air airborne and land gravity anomalies reduced from 64.3 mGal and 67.1 mGal to 59.6 mGal and 44.2 mGal, respectively. After removing GGM and TC, the STD of the free-air airborne and land gravity anomalies reduce to 39.0 mGal and 53.7 mGal, respectively. Because the high frequency contribution from the residual terrain can be removed in the remove-compute-restore procedure by RTM method, the STD of airborne gravity anomalies reduces from 44.2 to 39.0 mGal after removing TC [2,4,5]. According to Figure 6, the spatial pattern of free-air anomalies after removing GGM and TC is smoother than that in Figure 3. Figure 7 shows the terrain and bathymetry around Taiwan.

Discussion
According to Table 3, the RMS of the LSC method is reduced by more than one half after removing the TC effects computed by RTM method. Compared with Figure 5, the distribution of the points in Figure 4 is more dispersed and shows a tilted feature, especially for the points with heights higher than 500 m. The most significant improvement is that the maximum error reduces from 191.1 mGal in Figure 4 to 67.5 mGal in Figure 5. The big improvement shows that removing TC effects in the complex terrain region is very effective for the LSC method on the one hand, and about 3/4 of downward continuation errors are caused by residual terrain effects on the other hand. Moreover, the LSC method with TC improves the accuracy of downward continuation by approximately 40% compared with the optimal method in [3,23]. Additionally, the improvement indicated in Table 3 and Figure 5 also shows that using the RTM method to compute terrain effects is very effective [2,[4][5][6].
Generally, TC effects in the mountainous areas are larger than those in the low elevation areas. Therefore, after removing the TC effects, the spatial pattern of gravity anomalies will become smoother than before. According to Table 4, the airborne gravity anomalies are nicely smoothed by removing the GGM and the TC. Over 75% of Taiwan's topography is hills and high mountains. The variation of topography over Taiwan and its adjacent ocean area is large, ranging from approximately −6000 m in the trench east of Taiwan to approximately 4000 m in the Central Range of Taiwan. Because the complex topography contributes greatly to the high-frequency part of the gravity anomalies [4][5][6], after removing the TC, the spatial pattern of gravity anomalies in Taiwan in Figure 6 is smoother than the one in Figure 3. The smoother gravity anomalies are more suitable for executing reliable gravity functional computations [26] in general, and for performing a better downward continuation of gravity anomaly signals based on the LSC method [2,[4][5][6].

Conclusions
To obtain the optimal results of downward continuation in Taiwan, we analyzed the theoretical equations of the three downward continuation methods explained in Section 2. In the experiment step, we conducted comparison tests in the Tibetan Plateau, which showed that the least-squares collocation method generates the best results compared to the inverse Poisson's integral and the semi-parametric method combined with the regularization method.
We applied the least-squares collocation method with a terrain correction to the downward continuation of real airborne gravity data in Taiwan. The terrain correction was accounted for by the residual terrain modeling approach. The downward continuation results show that, with the help of terrain correction, the RMS of the differences between downward continuation values and land gravity values is 11.7 mGal in this paper. An improvement of approximately 40% was gained compared with the results given by Hwang et al. [3] and Hsiao and Hwang [23]. For the airborne gravity data in the region with complex topography, the terrain correction must be considered in the downward continuation. After removing the terrain correction effects, the magnitude of gravity anomalies decreases and its spatial pattern becomes smoother, which contribute greatly to reducing error in downward continuation. In the future, for further improvement in the downward continuation of airborne gravity, the influences of underwater topography, and more accurate and high-resolution topography data should be considered in data processing. Moreover, some new methods could be tested, such as radial basis functions, machine learning, and etc.