Investigation of the Tikhonov Regularization Method in Regional Gravity Field Modeling by Poisson Wavelets Radial Basis Functions

0. INTRODUCTION

Regional gravity field modeling is one of the classical issues in physical geodesy, the precise determination of local gravity field is of considerable importance not only for surveying and mapping, but also for the research fields like investigating the structure of seismic activities as well as lithosphere in geophysics and geodynamics (Wu et al., 2017; Guo et al., 2015; Xu et al., 2015). With the launch of the state-of-the-art satellite gravity missions, especially the dedicated GRACE/GOCE ones, the gravity field at long-wavelength up to tens or hundreds kilometers has been significantly improved. The incorporation of more local high-quality heterogeneous gravity related data also contributes to further improvement of the regional gravity field, especially in the short-wavelength parts down to several kilometers. Over the years, radial basis functions (RBFs) has been extensively used in regional gravity field modeling by combining heterogeneous data sets (Wu and Luo, 2016; Luthcke et al., 2013; Simons and Dahlen, 2006; Chambodut et al., 2005; Albertella et al., 1999). However, as the heterogeneous data have different spatial coverage, various error characteristics as well as different spectral contents, the associated solutions based on RBFs are typically ill-conditioned, where regularization is mandatory for deriving the reliable results (Wittwer, 2010).

Typically, several approaches could be applied for ill-posed problems, e.g., Tikhonov regularization (Tikhonov, 1963), Truncated singular value decomposition (TSVD) (Xu, 1998; Hansen, 1993), Least squares collocation (LSC) (Rummel et al., 1979), biased estimation (Hoerl and Kennard, 1970) and ridge estimation (Xu and Rummel, 1994). Tikhonov regularization method does not need any apriori information, which is extensively used in geodesy. Two key elements affect the quality of Tikhonov regularization, i.e., regularization parameter and regularization matrix. L-curve (Hansen et al., 2007; Johnston and Gulrajani, 2000), generalized cross validation (GCV) (Xu, 2009; Girard, 1989) and variance component estimation (VCE) (Koch, 1987) are the three widely used methods for estimating the optimal regularization parameter. Kusche and Klees (2002) used L-curve and GCV methods for estimating the optimal regularization parameter in the global gravity field modeling based on simulated GOCE-derived gradients. Koch and Kusche (2002) proposed the VCE approach in dealing with the ill-conditioned problems, which could simultaneously estimate the variance components of various observation groups as well as the regularization parameter. Moreover, Kusche and Klees (2002) compared the performances of various regularization matrices based on first- and second-order derivative constraints, which showed there were no significant differences among the solutions derived from various constraints. However, according to the research by Xu (1998), the L-curve method tends to derive an oversmoothed solution, which may affect the quality of the solution. Besides, Xu et al. (2006) indicated the VCE approach proposed by Koch and Kusche (2002) may introduce significant biases and lead to the inappropriate estimation of the variance components. As a result, Xu et al. (2006) proposed a bias-corrected VCE approach for estimating the variance components as well as the regularization parameter, which may outperform the original VCE approach proposed by Koch and Kusche (2002) in many cases. As the L-curve and VCE approaches are the two typical methods that are most extensively used in the ill-conditioned problems in physic geodesy, we mainly investigate the performances of these two approaches in dealing with ill-conditioned problems for regional gravity field modeling.

Moreover, unlike using the diagonal regularization matrices in the global gravity field modeling by spherical harmonics, the regularization matrices derived from various constraints by RBFs in the regional gravity field are not totally diagonal any more. Thus, the choice of regularization matrices may affect the quality of the solutions. Together with the regularization matrices, the approaches for deriving the optimal regularization parameters also affect the quality of the solutions, which need further investigation. Given the fact that Poisson wavelets have the explicit analytical kernel, which could be fast computed (Holschneider and Iglewska-Nowak, 2007), in this research, the Poisson wavelets are chosen as the basis functions for regional gravity field recovery, and the fast synthesized method for computing various regularization matrices is introduced. As a case study, the locally-distributed terrestrial, shipboard as well as airborne gravity data are combined for regional gravity field modeling, and the Tikhonov regularization method is applied to solve the ill-conditioned problem.

1. METHODOLOGY

1.1 Functional Model and Parameters Estimation

Based on the remove-compute-restore (RCR) method, the residual disturbing potential T could be expressed as the linear combination of Poisson wavelets (Klees et al., 2008)

where z and y are the vectors that include 3-D coordinates of the observation and position of Poisson wavelets, respectively. K is the number of Poisson wavelets, β_i is the unknown coefficient of Poisson wavelets, which should be determined from the data. d is the order of Poisson wavelets, and W_y^{ext, d}(z) is the Poisson wavelets (Holschneider and Iglewska-Nowak, 2007), which is expressed as

with w=z–y is the vector difference between z and y, $\mathit{\boldsymbol{\hat z}}$ and $\mathit{\boldsymbol{\hat y}}$ are the unit vector of z and y, respectively. ρ is the radius of a sphere that is totally located inside the earth, which is typically chosen as Bjerhammar sphere. P_l is the fully normalized Legendre functions. The coefficient $\alpha _l^d$ is iteratively computed as

where

Terrestrial, shipboard gravity anomalies △g and airborne gravity disturbances δg are used for regional gravity field modeling, and the GPS/leveling data is applied to evaluate the quality of the solutions. After linearization and spherical approximation, the gravity anomalies and gravity disturbances could be linked to residual disturbing potential through

The geoidal height N can also be related to residual disturbing potential through Bruns formula

and γ is the normal gravity value.

For the particular type of gravity observation group p, we map the observations L_p to the residual disturbing potential as follows

where Δ_p are the corresponding observation errors, f_p is the functional that links the observations to the residual disturbing potential, J is the number of observation groups.

We rewrite Eq. (8) in terms of vector-matrix notation

And x is the K×1 vector of unknown coefficients of Poisson wavelets, A_p is the m_p×K design matrix of observation group p, m_p is the number of the gravity observations of group p, l_p is the m_p×1 corresponding observation vector, e_p is m_p×1 vector of corresponding stochastic errors.

The final observational equations could be expressed as

Then

The observations in different groups are assumed to be uncorrelated to each other, then the error variance-covariance matrix of the observations is

with $E\left({{\mathit{\boldsymbol{e}}_p}} \right) = 0$, and $D\left({{\mathit{\boldsymbol{e}}_p}} \right) = {\mathit{\boldsymbol{C}}_p} = \mathit{\boldsymbol{\sigma }}_p^2{\mathit{\boldsymbol{Q}}_p}$ is the error variance- covariance matrix of observations in group p, $\sigma _p^2$ is the variance factor and Q_p is the cofactor matrix.

Due to the data gaps in the gravity observations, various characteristics of error sources as well as the imperfectly designed network of Poisson wavelets, the least squares-derived equations are usually ill-conditioned, where the regularization is mandatory. Based on penalized least squares principle, for a given α (regularization parameter), we minimize the quadratic functional (Klees et al., 2008)

and

where R is the regularization matrix, then the minimum of the quadratic functional of $\Phi \left( \mathit{\boldsymbol{x}} \right)$ is obtained by

with

Variance component estimation (VCE) approach (Koch and Kusche, 2002) is used to properly determine the weight of the disjunctive observation group, where the variance factor of group p is computed as

and ${\mathit{\boldsymbol{v}}_p} = {\mathit{\boldsymbol{A}}_p}\mathit{\boldsymbol{\hat x}} - {\mathit{\boldsymbol{l}}_p}$ is the vector of residuals of observation group p, r_p is the corresponding redundancy number, which is expressed as

As the regularization is compulsory for deriving the reliable solutions, two key factors that affect the quality of the regularized solutions, i.e., the regularization matrix as well as regularization parameter, are studied in the following parts. It is also worthy to be mentioned that as the prior expectation of the residual regional gravity field is not equal to the zero vector, the estimation of unknown parameters through Eqs. (16) and (17) are supposed to be the biased ones (Xu et al., 2006; Xu and Rummel, 1994; Xu, 1992). However, according to the research by Xu (1992), the biases of the unknown parameters are proportional to the regularization parameter and the true values of the unknown parameters. We use the remove- compute-restore method to reduce the sizes of the unknown parameters such that they can be as close to zero prior mean values as possible (see the statistics of the residual gravity observations in Table 1). As a result of this extra effort, we can further reduce the effect of possible biases on our results.

Table  1.  Statistics of the residual gravity observations (units: mGal)

	Max	Min	Mean	Sd
Terrestrial	90.47	-165.52	-1.08	10.65
Shipboard	158.39	-78.93	-0.74	12.40
Airborne	49.22	-25.75	2.15	10.44

 | Show Table

DownLoad: CSV

1.2 The Choice of Regularization Matrix

For the selection of the optimal regularization matrix, the performances of zero- and first-order Tikhonov regularization matrices are investigated.

1.2.1   Zero-order Tikhonov regularization

Using the addition theorem of spherical harmonics, the Poisson wavelets is rewritten as

and ${\bar Y_{lm}}$ is the fully normalized spherical harmonics, l and m is its degree and order, respectively.

We assume the target function as the Poisson wavelets on the sphere, i.e., $s = W_\mathit{\boldsymbol{y}}^{ext, d}(\mathit{\boldsymbol{z}}).$ Then the element of zero-order Tikho- nov regularization matrix $\mathit{\boldsymbol{R}}_{P, Q}^0$ is expressed as the inner product of two different Poisson wavelets (Heiskanen and Moritz, 1967)

Based on the research by Holschneider and Iglewska- Nowak (2007), Eq. (21) is fast synthesized by

where $\mathit{\boldsymbol{y}}_P^* = \frac{{{\rho ^2}}}{{{{\left| {{\mathit{\boldsymbol{y}}_P}} \right|}^2}}}{\mathit{\boldsymbol{y}}_P}$ and $\mathit{\boldsymbol{y}}_Q^* = \frac{{{\rho ^2}}}{{{{\left| {{\mathit{\boldsymbol{y}}_Q}} \right|}^2}}}{\mathit{\boldsymbol{y}}_Q}$.

1.2.2   First-order Tikhonov regularization

The zero-order Tikhonov regularization derived above is suitable for functions that restrict to a sphere where σ_U belongs to L₂(σ_U). However, this function space is quite large, it means that quite unsmooth functions also belong to this space. Therefore, one may wish to introduce a scalar product, which is appropriate for smoother functions. Suppose the target function is the first-order derivative of Poisson wavelets, then the inner product of different target functions could be used to compute the elements of the first-order Tikhonov regularization matrix.

Assume the target function is the first-order radial derivative of Poisson wavelets, then elements of radial constrained regularization matrix $\mathit{\boldsymbol{R}}_{P, Q}^r$ is computed as (Heiskanen and Moritz, 1967)

As $\frac{\partial }{{\partial \left| \mathit{\boldsymbol{z}} \right|}}{\psi _l} = - \frac{{l + 1}}{{\left| \mathit{\boldsymbol{z}} \right|}}{\psi _l}$, see Eq. (20), then Eq. (23) is rewritten as

According to Holschneider and Iglewska-Nowak (2007), if the operator $\Lambda $defines through

then

Combining Eqs. (22), (24) and (26), $\mathit{\boldsymbol{R}}_{P, Q}^r$ is fast computed as

The target function could also be chosen as the horizontal derivative of Poisson wavelets, the corresponding elements of horizontal constrained regularization matrix $\mathit{\boldsymbol{R}}_{P, Q}^h$ is computed as

Similarly, Eq. (28) is computed through

1.3 Methods for Regularization Parameter Estimation

1.3.1   L-curve

L-curve is the plot of the residual norm $r\left(\alpha \right) = {\left\| {\mathit{\boldsymbol{A}}{{\mathit{\boldsymbol{\hat x}}}_\alpha } - \mathit{\boldsymbol{l}}} \right\|_P}$ versus the solution norm $s\left(\alpha \right) = {\left\| {{{\mathit{\boldsymbol{\hat x}}}_\alpha }} \right\|_K}$ for all the non-negative regularization parameters (Hansen et al., 2007). For discrete ill-conditioned problems, this plot displays as the 'L-shape' with the corner point, where a certain balance between the regularized solution and the fit to the data is achieved. Usually, the L-curve is drawn as log-log scale based on the two- dimensional matrix as follows (Hansen et al., 1993)

The corner point of L-curve is obtained based on the curvature of this curve, where the point derives the maximum curvature is treated as the corner point. The curvature of the L-curve is derived as (Hansen et al., 1993)

where μʹ, λʹ, μʺ and λʺ are the first- and second-order derivatives of μ(α) and λ(α), respectively.

It is noticeable that the L-curve method may lead to an oversmoothed solution, see Xu (1998), where a significant number of components were discarded when this approach was used to determine the proper cutting number of the TSVD approach. This result means the regularization parameters computed from the L-curve method may be larger than the proper ones, which may affect the quality of the solutions. Thus, we also introduce VCE approaches for cross validation and make necessary comparisons among the solutions derived from various methods.

1.3.2   Variance component estimation

VCE approach could also be used to properly determine the regularization parameter (Koch and Kusche, 2002). We rewrite Eq. (13) as (Klees et al., 2008)

and the "virtual" J+1th observation group (l_J+1=0) is added, which is treated as the prior information for the unknown parameters. The corresponding error equation and the error variance-covariance matrix of observations for this group are expressed as (Klees et al., 2008)

where I_K is the K × K identity matrix, σ²_J+1=1/α is the variance factor of observation group J+1, Q_J+1=R^-1 is the corresponding cofactor matrix. Similarly, we estimate the minimum value of the quadratic functional of $\Phi \left(\mathit{\boldsymbol{x}} \right)$ in Eq. (32), and the estimated unknown parameters could also be expressed as

where the normal equation matrix and right-hand side vector are

The Eqs. (34) and (35) have the same form as the Eqs. (16) and (17). And the reciprocal of regularization parameter could be estimated as the variance component of the J+1th observation group, which is expressed as (and Kusche)

And regularization parameter α is estimated as

However, as shown by Xu et al. (2006), the VCE approach derived by Koch and Kusche (2002) may introduce significant biases to the estimated variance components when regularization is applied to tackle the ill-conditioned least squares system. However, given the numerical efficiency associated with the Monte Carlo method embedded in the VCE approach proposed by Koch and Kusche (2002), especially for the regional gravity field modeling with millions of heterogeneous observations, we will limit ourselves to the original VCE approach by Koch and Kusche (2002) in this research, further comparisons with other methods may be left to a future work.

1.3.3   The minimum standard deviation approach

L-curve and VCE approaches estimate the regularization parameters only depend on the observations, and these two methods may have potential pitfalls as mentioned above. Thus, we also use an alternative approach for further validation, which may provide the users with more insight regarding the feasibiltiy of using the L-curve and VCE approaches. As the high- quality GPS/leveling data are available over the target area, the quality of the solutions based on various regularization parameters could be directly evaluated, which is the so-called minimum standard deviation approach (MSTD) (Wittwer, 2010). We assume there are m GPS/leveling points in total over the region, the difference between the geoidal height modeled from Poisson wavelets $N_{SRBF}^i$ and the GPS/leveling-derived one $N_{GPSL}^i$ at the i-th GPS/leveling point is expressed as

And the standard deviation of these differences is computed by (Wittwer, 2010)

and $\overline {{N_{res}}} $ is the corresponding mean value.

Thus, the regularization parameter that minimizes Eq. (39) is regarded as the optimal one, which could be used for cross validation when L-curve and VCE approaches are applied for regularization parameters estimation. It is noticeable that the principle of MSTD approach is similar to GCV method, and both of these two approaches estimate the unknown parameters based on the predicted and observed values. In our opinion, the MSTD method is still different from the GCV approach. GCV method uses the leave-one-out idea, where the proper regularization parameters are derived if the good estimation of the leaving out observation is obtained based on the other observations (Xu, 2009). To be more specific, GCV approach estimates the regularization parameters from residuals analysis only from observations. However, when applying the MSTD approach, the independent ground-based control data (e.g., GPS/leveling data) are used instead of the leaving out observation. In this way, the high-quality GPS/leveling data are treated as the true values, and the regularization parameter derives the solution that best fit to the control data would be regarded as the proper one. Thus, the prerequisite for using MSTD method is that the high-quality control data sets are available, while the GCV approach doesn't need this precondition. Moreover, the reliability of the proper parameters estimated by MSTD method largely depends on the quality of the control data sets, and this approach may lead to the inappropriate estimation of the unknown parameters if the quality of the control data is suspicious.

2. RESULTS

2.1 Data Sets

2.1.1   The global geopotential model and digital terrain model

Based on RCR methodology, the long- and short-wavelength of regional gravity field signals are represented by global geopotential model (GGM) and residual terrain model (RTM), respectively, and the residual part is recovered from the residual gravity data. The GGM is chosen as the combined GRACE/GOCE geopotential model, i.e., DGM1S complete to degree and order of 250 (Hashemi Farahami et al., 2013). A high-resolution and high- quality digital terrain model (DTM) is derived for RTM reduction, where three different data sources, i.e., EuroDEM, SRTM and GEBCO, are combined together, see Fig. 1a. To obtain the mean elevation surface (MES) used in RTM correction, a moving- average filter is applied to the combined DTM model (Hirt, 2013), which makes the spatial resolution of the derived MES consistent with the adopted GGM, i.e., DGM1S with the resolution of approximately 0.72 , see Fig. 1b. Given the difference of the mean crust density in land and at seas, the density parameter used in RTM is selected as ρ_land=2.67 g/cm³ and ρ_sea=1.64 g/cm³, and the RTM corrections were computed based on tesseroid integrals (Heck and Seitz, 2006).

Figure  1.  Digital terrain model (a) and mean elevation surface (b).

DownLoad: Full-Size Img PowerPoint

2.1.2   Heterogeneous gravity data sets

Terrestrial, shipboard and airborne gravity data are combined for regional gravity field modeling. Point-wise terrestrial data cover the whole Netherlands, Belgium, England as well as parts of Germany and France, which are derived from Bureau Gravimétri- que International (BGI), Bundesamt für Kartographie und Geod sie (BKG) and Nordic Geodetic Commission (NKG). The corresponding spatial resolution is approximately 5 km and the accuracy is roughly at 1 mGal level. Shipboard gravity data are provided by BGI, British Geological Service (BGS), Institut für Erdmessung (IFE), National Geophysical Data Center (NGDC) and NKG, the spatial resolution and precision of which is approximately 7 km and 2 mGal, respectively. Airborne gravity disturbances are obtained from BKG and NKG, the estimated accuracy is roughly at the magnitude of 2 mGal. Crossover adjustment is applied to remove the systematic errors in the shipboard and airborne data, and low-pass filter is used for reducing the effect of high-frequency noise in these two data sets. Moreover, DTU10-GRA model (Andersen, 2010) is applied for removing the outliers exist in the shipboard data through the threshold rules. Then, all the gravity data sets are referred to European Terrestrial Reference System 1989 (ETRS89) and European Vertical Reference Frame 2007 (EVRF2007) as the horizontal and vertical datum, respectively. Based on RCR methodology, the residual gravity anomalies/disturbances are computed by subtracting the short- and long-wavelength parts from the original data, which are shown in Fig. 2. The corresponding statistics are displayed in Table 1.

Figure  2.  Residual gravity observations. (a) Terrestrial gravity anomalies; (b) shipboard gravity anomalies; (c) airborne gravity disturbances.

DownLoad: Full-Size Img PowerPoint

2.2 Determination of Regularization Matrix

To choose the optimal regularization matrix in dealing with the ill-conditioned normal equation in regional gravity field modelling, the numerical experiment is designed as follows. The computational region is bounded by 50°N to 55°N in the latitude direction and 2°E to 8°E in the longitude direction. This region includes whole mainland of Netherlands, as well as parts of the North Sea, Belgium, UK, Germany and France, and the high- quality GPS/leveling data over Netherlands are served as evaluation data. Poisson wavelets are put on a grid that is located on a constant depth beneath the Bjerhammar sphere (Tenzer and Klees, 2008). The depth of Poisson wavelets is selected as 50 km, and the number of Poisson wavelets is chosen as 7 394 (the mean distance between Poisson wavelets is approximately 8.7 km). In this case, the regularization is mandatory as the derived least squares system is highly ill-conditioned. The performance of various Tikhonov regularization matrices is investigated, and the L-curve is applied to estimate the optimal regularization parameters.

Figure 2 shows the L-curve plots for various regularization matrices, where the optimal regularization parameter with zero-order regularization is between 10^-11 and 10^-10. While, for the first-order regularization (both for the radial and horizontal constrains), the optimal parameter located between 10^-12 and 10^-11. As shown in Fig. 4, although all the normalized regularization matrices are diagonal-dominated (the regularization matrices are shown in log scale), there are still differences between the zero- and first-order regularization matrix, which lead to the difference between the estimated regularization parameters. To further demonstrate the corner point of L-curve, we densify the regularization parameters around the intervals where the optimal parameters located in, see the red circles in Fig. 3, and the densified regularization parameters for various regularization matrices could be referred to Table 2 and Table 3. Based on Eq. (31), the estimated optimal regularization parameter for zero- and first-order regularization is α=10^-10.6 and α=10^-11.3, respectively, which coincide well with the estimated parameters derived from the MSTD method, see Table 2 and Table 3. In addition, there are no differences for the performances of two different first-order regularization matrices, i.e., the radial and horizontal constrained matrices, thus we don't distinguish these two matrices in the following research. However, when comparing the quality of the solutions derived from zero- and first-order regularization, we find the application of first-order regularization derives better results, and the accuracy of the corresponding geoid reaches 1.3 cm. For the solutions computed from zero-order regularization, the accuracy decreases to 1.5 cm. These results show the choice of the regularization matrices has nonnegligible effects on the solutions, and the first-order regularization matrix may be more preferable in the regional gravity field modeling based on Poisson wavelets.

Figure  3.  Optimal regularization parameters computed from L-curve methodology based on various regularization matrices. (a) Zero-order regularization, (b) first-order regularization (radial gradient constraint), (c) first-order regularization (horizontal gradient constraint).

DownLoad: Full-Size Img PowerPoint

Figure  4.  Various normalized regularization matrices. (a) Zero-order matrix, (b) first-order matrix (radial gradient constraint), (c) first-order matrix (horizontal gradient constraint).

DownLoad: Full-Size Img PowerPoint

Table  2.  Evaluation of the different geoids based on various regularization parameters using the zero-order Tikhonov regularization (Units: m)

Regularization parameters	10-14	10-13	10-12	10-11	10-10.8	10-10.7	10-10.6	10-10.4	10-10.2	10-10	10-9	10-8	10-7
	0.025	0.022	0.020	0.017	0.017	0.016	0.015	0.016	0.017	0.019	0.021	0.023	0.027

 | Show Table

DownLoad: CSV

Table  3.  Evaluation of the different geoids based on various regularization parameters using the first-order Tikhonov regularization (Units: m)

Regularization parameters	10-14	10-13	10-12	10-11.8	10-11.6	10-11.4	10-11.3	10-11.2	10-11	10-10	10-9	10-8	10-7
Radial constrained	0.021	0.018	0.018	0.015	0.015	0.014	0.013	0.014	0.014	0.015	0.016	0.018	0.020
Horizontal constrained	0.021	0.018	0.018	0.015	0.015	0.014	0.013	0.014	0.014	0.015	0.016	0.018	0.020

 | Show Table

DownLoad: CSV

2.3 Optimal Approach for Estimating Regularization Parameter

The L-curve, VCE and MSTD approaches are investigated for their performances in deriving the optimal regularization parameters at this part. The numerical experiment we study here is the same as the one designed above, and the first-order regularization is used as the case study. Figure 5 shows the estimated regularization parameter from VCE approach after five iterations, and the convergent parameter is equal to 0.500 1× 10^-11≈10^-11.3, which coincides well with the results derived from L-curve and MSTD approaches, see Fig. 3, Tables 2 and 3. Compared to the L-curve and MSTD approaches, which need to compute plenty of solutions with various regularization parameters, VCE methodology derives the convergent parameters in several iterations, typically no more than 7 iterations in the usual case, which is more efficient. Besides, the VCE method estimates the proper weight of disjunctive observation group as well as the regularization parameter simultaneously, which is more suitable in regional gravity field modeling from hetero- geneous data sets. Although the L-curve and VCE approaches work fine in this case study, users should always keep the potential pitfalls of these two approaches in mind (Xu et al., 2006; Xu, 1998), and the incorporation of various approaches for cross validation may be necessary.

Figure  5.  Optimal regularization parameter computed from VCE approach.

DownLoad: Full-Size Img PowerPoint

2.4 Regional Gravity Field Modeling

The design of Poisson wavelets' network has significant effects on the solutions, where the depth and number of Poisson wavelets are the two key factors (Wu et al., 2016; Klees et al., 2008). To our experience, the depth of Poisson wavelets' grid should not be chosen too shallow, which may lead to overfitting problem. Too deep Poisson wavelets may cause heavily ill- conditioned problem, where strong regularization is need, and the associated regularization errors may corrupt the solution. Together with the depth, the number of Poisson wavelets also affects the stability and the quality of the solution. The good choice of depth and number of Poisson wavelets should be the trade-off between the fit to the data and smoothness of the solution.

To find the optimal depth and number of Poisson wavelets, different depths and number of Poisson wavelets are combined together to form the various parameterizations of networks, based on which various geoids are computed. To select the best solution, the GPS/leveling points are served as evaluation data, where the combination of the depth and number of Poisson wavelets that obtains the best fit to the evaluation data is considered as the optimal parameters. The quality of this fit is determined by the standard deviation (STD) of the differences between the modeled geoidal heights and the ones derived from the GPS/leveling data. By trial and errors, we make spatial resolution of Poisson wavelets vary from 7.2 to 9.0 km with an increment of 0.3 km, and the depth of Poisson wavelets changes from 30 to 50 km with the interval of 5 km. To deal with the ill-conditioned problem, first-order regularization together VCE approach is applied. Figure 6 shows the accuracy of the different geoids based on various Poisson wavelets' network, where the smallest standard deviation (0.011 m) is obtained when the depth is 40 km and the number of Poisson wavelets is equal to 7 394 (the mean distance between Poisson wavelets is approximately 8.7 km), the detailed statistics could be referred to Table 4. It is also noticeable that there are still some unfavourable solutions even Tikhonov regularization is applied for tackling the ill-condition system, see the black grids in Fig. 6. The main reason is that the overlapping between the Poisson wavelets is significantly increased with the incorporation of more deep Poisson wavelets, which leads to highly ill-conditioned normal matrix. As a result, heavy regularization should be applied for deriving the stable solutions, and the quality of the solutions is corrupted with the associated over-regularization. Figure 7 shows the evaluation results of the gravimetric geoid computed from the optimal Poisson wavelets' network in the Netherlands.

Figure  6.  Determination of the optimal Poisson wavelets' network.

DownLoad: Full-Size Img PowerPoint

Figure  7.  Evaluation of the gravimetric geoid computed from the optimal Poisson wavelets' network.

DownLoad: Full-Size Img PowerPoint

Table  4.  Evaluation of different geoids computed from various Poisson wavelets' networks (units: m)

	6720	7394	8043	8694	9333	9970	10571
30	0.019	0.025	0.016	0.022	0.017	0.020	0.013
35	0.015	0.020	0.016	0.013	0.012	0.015	0.012
40	0.012	0.011	0.012	0.013	0.013	0.015	0.015
45	0.013	0.012	0.013	0.014	0.016	0.016	0.022
50	0.014	0.013	0.016	0.023	0.040	0.046	0.060

 | Show Table

DownLoad: CSV

3. CONCLUSIONS

The Tikhonov regularization method in dealing with the ill-conditioned problems for the regional gravity filed modeling is investigated. The regional gravity field is modeled by combining the terrestrial, shipboard and airborne gravity data using Poisson wavelets, and the results show that the choice of regularization matrices affects the quality of the solutions. Compared to the solution with zero-order regularization, the application of the first-order regularization derives the better results, where the accuracy of the gravimetric geoid increases by 0.2 cm. In addition, the optimal regularization parameters derived from three approaches, i.e., L-curve, VCE and MSTD, are quite consistent with each other, which validate the reliability of using these three approaches for estimating the regularization parameters. However, the VCE method could estimate the weights for various observation groups as well as regularization parameters in the same time, which is more efficient for regional gravity field modeling based on huge data. Moreover, we design the optimal network of Poisson wavelets based on the first-order regularization and VCE method, the accuracy of the corresponding gravimetric geoid reaches 1.1 cm in the Netherlands. This result also validates the reliability of applying Tikhonov regularization method in tackling the ill- conditioned problems in the regional gravity field modeling. Although the L-curve and VCE approaches work fine in this case study, it is also noticeable that the potential pitfalls exist in these two methods (Xu et al., 2006; Xu, 1998). The future work may involve implementing the bias-corrected VCE approach (Xu et al., 2006) in a numerical efficient way, which could be used in regional gravity field modeling based on millions of heterogeneous observations. It is also worthy to be mentioned that even we derive relative high-quality solutions over the Netherlands, the effects on the solutions caused by various error sources in the data need to be investigated in future, which may contribute to further improve the solutions.

ACKNOWLEDGMENTS

Thanks go to the three anonymous reviewers, who gave constructive comments and beneficial suggestions which help us to improve this manuscript. Thanks also go to Prof. Roland Klees from Delft University of Technology for kindly providing the gravity data. This research was mainly supported by the National Natural Science Foundation of China (Nos. 41374023, 41131067, 41474019), the National 973 Project of China (No. 2013CB733302), the China Postdoctoral Science Foundation (No. 2016M602301), the Key Laboratory of Geospace Environment and Geodesy, Ministry of Education, Wuhan University (No. 15-02-08), and the State Scholarship Fund from Chinese Scholarship Council (No. 201306270014). The final publication is available at Springer via https://doi.org/10.1007/s12583-017-0771-3.