2. School of Geodesy and Geomatics, Wuhan University, Wuhan 430079, China;
3. Key Laboratory of Geospace Environment and Geodesy, Ministry of Education, Wuhan University, Wuhan 430079, China
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 stateoftheart satellite gravity missions, especially the dedicated GRACE/GOCE ones, the gravity field at longwavelength up to tens or hundreds kilometers has been significantly improved. The incorporation of more local highquality heterogeneous gravity related data also contributes to further improvement of the regional gravity field, especially in the shortwavelength 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 illconditioned, where regularization is mandatory for deriving the reliable results (Wittwer, 2010).
Typically, several approaches could be applied for illposed 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. Lcurve (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 Lcurve and GCV methods for estimating the optimal regularization parameter in the global gravity field modeling based on simulated GOCEderived gradients. Koch and Kusche (2002) proposed the VCE approach in dealing with the illconditioned 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 secondorder 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 Lcurve 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 biascorrected 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 Lcurve and VCE approaches are the two typical methods that are most extensively used in the illconditioned problems in physic geodesy, we mainly investigate the performances of these two approaches in dealing with illconditioned 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 IglewskaNowak, 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 locallydistributed 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 illconditioned problem.
1 METHODOLOGY 1.1 Functional Model and Parameters EstimationBased on the removecomputerestore (RCR) method, the residual disturbing potential T could be expressed as the linear combination of Poisson wavelets (Klees et al., 2008)
$T\left(\mathit{\boldsymbol{z}} \right) = \sum\limits_{i = 1}^K {{\beta _i}W_{{\mathit{\boldsymbol{y}}_i}}^{ext, d}(\mathit{\boldsymbol{z}})} $  (1) 
where z and y are the vectors that include 3D 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 IglewskaNowak, 2007), which is expressed as
$W_\mathit{\boldsymbol{y}}^{ext, d}(\mathit{\boldsymbol{z}}) = \sum\limits_{l = 0}^{d + 1} {(2\alpha _l^{d + 1} + \alpha _l^d)\frac{{\rho l!{{\left \mathit{\boldsymbol{z}} \right}^l}{P_l}\left({\mathit{\boldsymbol{\hat w}} \cdot \mathit{\boldsymbol{\hat y}}} \right)}}{{{{\left {\mathit{\boldsymbol{z}}  \mathit{\boldsymbol{y}}} \right}^{l + 1}}}}} $  (2) 
with w=z–y is the vector difference between z and y,
$\begin{array}{l} \alpha _l^d = \alpha _{l  1}^{d  1} + l\alpha _l^{d  1}, \;\;d \ge 1\\ \alpha _k^0 = {\delta _{k, 0}} \end{array} $  (3) 
where
${\delta _{k, 0}} = \left\{ \begin{array}{l} 1\;\;\;{{\rm{if}}}\;k = 0\\ 0\;\;{{\rm{if}}}\;k \ne 0 \end{array} \right.\;\; $  (4) 
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
$\Delta g\left(\mathit{\boldsymbol{z}} \right) \cong  \frac{2}{{\left \mathit{\boldsymbol{z}} \right}}T\left(\mathit{\boldsymbol{z}} \right)  \frac{{\partial T\left(\mathit{\boldsymbol{z}} \right)}}{{\partial \left \mathit{\boldsymbol{z}} \right}} $  (5) 
$\delta g\left(\mathit{\boldsymbol{z}} \right) \cong  \frac{{\partial T\left(\mathit{\boldsymbol{z}} \right)}}{{\partial \left \mathit{\boldsymbol{z}} \right}} $  (6) 
The geoidal height N can also be related to residual disturbing potential through Bruns formula
$N\left(\mathit{\boldsymbol{z}} \right) = \frac{{T\left(\mathit{\boldsymbol{z}} \right)}}{{\gamma \left(\mathit{\boldsymbol{z}} \right)}} $  (7) 
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
${\mathit{\boldsymbol{L}}_p} + {\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_p} = {\mathit{\boldsymbol{f}}_p}T\left(\mathit{\boldsymbol{z}} \right) = \sum\limits_{i = 1}^K {{\beta _i}{\mathit{\boldsymbol{f}}_p}W_{{\mathit{\boldsymbol{y}}_i}}^{ext, d}(\mathit{\boldsymbol{z}})}, \;p = 1, 2..., J $  (8) 
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 vectormatrix notation
${\mathit{\boldsymbol{l}}_p}  {\mathit{\boldsymbol{e}}_p} = {\mathit{\boldsymbol{A}}_p}\mathit{\boldsymbol{x}} $  (9) 
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
$\mathit{\boldsymbol{l}}  \mathit{\boldsymbol{e}} = \mathit{\boldsymbol{Ax}} $  (10) 
Then
$\mathit{\boldsymbol{l}} = \left(\begin{array}{l} {\mathit{\boldsymbol{l}}_1}\\ \; \vdots \\ {\mathit{\boldsymbol{l}}_p}\\ \; \vdots \\ {\mathit{\boldsymbol{l}}_J} \end{array} \right), \;\mathit{\boldsymbol{e}} = \left(\begin{array}{l} {\mathit{\boldsymbol{e}}_1}\\ \; \vdots \\ {\mathit{\boldsymbol{e}}_p}\\ \; \vdots \\ {\mathit{\boldsymbol{e}}_J} \end{array} \right), \mathit{\boldsymbol{A}} = \left(\begin{array}{l} {\mathit{\boldsymbol{A}}_1}\\ \; \vdots \\ {\mathit{\boldsymbol{A}}_p}\\ \; \vdots \\ {\mathit{\boldsymbol{A}}_J} \end{array} \right)\; $  (11) 
The observations in different groups are assumed to be uncorrelated to each other, then the error variancecovariance matrix of the observations is
$\mathit{\boldsymbol{C}} = \left(\begin{array}{l} {\mathit{\boldsymbol{C}}_1}\;\;0\;\;0\;\; \cdots \;\;0\\ 0\;\;\;{\mathit{\boldsymbol{C}}_2}\;0\;\; \cdots \;\;0\\ \; \vdots \;\;\;\; \vdots \;\;\; \vdots \;\;\; \cdots \;\; \vdots \\ 0\;\;\;0\;\;\;0\;\; \cdots \;{\mathit{\boldsymbol{C}}_J}\; \end{array} \right) $  (12) 
with
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 squaresderived equations are usually illconditioned, 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)
$\Phi \left(\mathit{\boldsymbol{x}} \right) = \left\ \mathit{\boldsymbol{e}} \right\_{{\mathit{\boldsymbol{C}}^{  1}}}^2 + \alpha \Gamma (T) $  (13) 
and
$\left\ \mathit{\boldsymbol{e}} \right\_{{\mathit{\boldsymbol{C}}^{  1}}}^2 = \sum\limits_{p = 1}^J {\left\ {{\mathit{\boldsymbol{e}}_p}} \right\_{\mathit{\boldsymbol{C}}_{^p}^{  1}}^2} $  (14) 
$\Gamma (T) = \left\ \mathit{\boldsymbol{x}} \right\_\mathit{\boldsymbol{R}}^2 $  (15) 
where R is the regularization matrix, then the minimum of the quadratic functional of
$\mathit{\boldsymbol{\hat x}} = {\mathit{\boldsymbol{N}}^{  1}}\mathit{\boldsymbol{W}} $  (16) 
with
$\begin{array}{l} \mathit{\boldsymbol{N}} = \sum\limits_{p = 1}^J {{\mathit{\boldsymbol{N}}_p}} + \alpha \mathit{\boldsymbol{R}}, \;{\mathit{\boldsymbol{N}}_p} = \mathit{\boldsymbol{A}}_p^T\mathit{\boldsymbol{C}}_p^{  1}{{\rm{N}}_p}\\ \mathit{\boldsymbol{W}} = \sum\limits_{p = 1}^J {{\mathit{\boldsymbol{W}}_p}}, \;{\mathit{\boldsymbol{W}}_p} = \mathit{\boldsymbol{A}}_p^T\mathit{\boldsymbol{C}}_p^{  1}{\mathit{\boldsymbol{l}}_p} \end{array} $  (17) 
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
$\hat \sigma _p^2 = \frac{{\mathit{\boldsymbol{v}}_p^T\mathit{\boldsymbol{Q}}_p^{  1}{\mathit{\boldsymbol{v}}_p}}}{{{r_p}}} $  (18) 
and
${r_p} = {m_p}  \frac{1}{{\hat \sigma _p^2}}tr\left({\mathit{\boldsymbol{A}}_p^T\mathit{\boldsymbol{Q}}_p^{  1}{\mathit{\boldsymbol{A}}_p}{\mathit{\boldsymbol{N}}^{  1}}} \right) $  (19) 
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 computerestore 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.
For the selection of the optimal regularization matrix, the performances of zero and firstorder Tikhonov regularization matrices are investigated.
1.2.1 Zeroorder Tikhonov regularizationUsing the addition theorem of spherical harmonics, the Poisson wavelets is rewritten as
$\begin{array}{l} W_\mathit{\boldsymbol{y}}^{ext, d}(\mathit{\boldsymbol{z}}) = \frac{\rho }{{\left \mathit{\boldsymbol{z}} \right}}{\sum\limits_{l = 0}^\infty {{l^d}\left({\frac{{\left \mathit{\boldsymbol{y}} \right}}{{\left \mathit{\boldsymbol{z}} \right}}} \right)} ^l}(2l + 1){P_l}\left({{{\mathit{\boldsymbol{\hat z}}}^T}\mathit{\boldsymbol{\hat y}}} \right)\\ \;\;\;\;\;\;\;\;\;\;\;\;\; = \sum\limits_{l = 0}^\infty {\sum\limits_{m =  l}^l {{\psi _l}{{\bar Y}_{lm}}\left({\mathit{\boldsymbol{\hat z}}} \right)} {{\bar Y}_{lm}}\left({\mathit{\boldsymbol{\hat y}}} \right)}, \;\;{\psi _l} = {l^d}{\left({\frac{{\left \rho \right}}{{\left \mathit{\boldsymbol{z}} \right}}} \right)^{l + 1}} \end{array} $  (20) 
and
We assume the target function as the Poisson wavelets on the sphere, i.e.,
$\mathit{\boldsymbol{R}}_{P, Q}^0 = \left\langle {{s_P}, {s_Q}} \right\rangle = \frac{1}{{4\pi {\rho ^2}}}\int_{{\Omega _\rho }} {W_{{\mathit{\boldsymbol{y}}_P}}^{ext, d}(\mathit{\boldsymbol{z}})} W_{{\mathit{\boldsymbol{y}}_Q}}^{ext, d}(\mathit{\boldsymbol{z}})d{\Omega _\rho } $  (21) 
Based on the research by Holschneider and Iglewska Nowak (2007), Eq. (21) is fast synthesized by
$\left\langle {W_{{\mathit{\boldsymbol{y}}_P}}^{ext, d}, W_{{\mathit{\boldsymbol{y}}_Q}}^{ext, d}} \right\rangle = \frac{\rho }{{\left {{\mathit{\boldsymbol{y}}_P}} \right}}W_{{\mathit{\boldsymbol{y}}_Q}}^{ext, 2d}\left({\mathit{\boldsymbol{y}}_P^*} \right) = \frac{\rho }{{\left {{\mathit{\boldsymbol{y}}_Q}} \right}}W_{{\mathit{\boldsymbol{y}}_P}}^{ext, 2d}\left({\mathit{\boldsymbol{y}}_Q^*} \right) $  (22) 
where
The zeroorder Tikhonov regularization derived above is suitable for functions that restrict to a sphere where σ_{U} belongs to L_{2}(σ_{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 firstorder derivative of Poisson wavelets, then the inner product of different target functions could be used to compute the elements of the firstorder Tikhonov regularization matrix.
Assume the target function is the firstorder radial derivative of Poisson wavelets, then elements of radial constrained regularization matrix
$\mathit{\boldsymbol{R}}_{P, Q}^r = \frac{1}{{4\pi {\rho ^2}}}\int_{{\Omega _\rho }} {\frac{{\partial W_{{y_\mathit{\boldsymbol{P}}}}^{ext, d}(\mathit{\boldsymbol{z}})}}{{\partial \left \mathit{\boldsymbol{z}} \right}}} \frac{{\partial W_{{\mathit{\boldsymbol{y}}_Q}}^{ext, d}(\mathit{\boldsymbol{z}})}}{{\partial \left \mathit{\boldsymbol{z}} \right}}d{\Omega _\rho } $  (23) 
As
$\mathit{\boldsymbol{R}}_{P, Q}^r = {\sum\limits_{l = 0}^\infty {{\psi _l}^2\left({\frac{{l + 1}}{{\left \mathit{\boldsymbol{z}} \right}}} \right)} ^2}(2l + 1){P_l}\left({{{\mathit{\boldsymbol{\hat y}}}_P}^T{{\mathit{\boldsymbol{\hat y}}}_Q}} \right) $  (24) 
According to Holschneider and IglewskaNowak (2007), if the operator
$\Lambda :{\bar Y_{lm}} \to {l^f}{\bar Y_{lm}} $  (25) 
then
$\Lambda W_\mathit{\boldsymbol{y}}^{ext, d}(\mathit{\boldsymbol{z}}) = W_\mathit{\boldsymbol{y}}^{ext, d + f}(\mathit{\boldsymbol{z}}) $  (26) 
Combining Eqs. (22), (24) and (26),
$\begin{array}{l} \mathit{\boldsymbol{R}}_{P, Q}^r = \frac{\rho }{{{{\left \mathit{\boldsymbol{z}} \right}^2}\left {{\mathit{\boldsymbol{y}}_P}} \right}}\\ \quad \quad \;\;\left({W_{{\mathit{\boldsymbol{y}}_Q}}^{ext, 2d + 2}\left({\mathit{\boldsymbol{y}}_P^*} \right) + 2W_{{\mathit{\boldsymbol{y}}_Q}}^{ext, 2d + 1}\left({\mathit{\boldsymbol{y}}_P^*} \right) + W_{{\mathit{\boldsymbol{y}}_Q}}^{ext, 2d}\left({\mathit{\boldsymbol{y}}_\mathit{\boldsymbol{P}}^\mathit{\boldsymbol{*}}} \right)} \right) \end{array} $  (27) 
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 = \sum\limits_{l = 0}^\infty {{\psi _l}^2\frac{{l\left({l + 1} \right)}}{{{{\left x \right}^2}}}} (2l + 1){P_l}\left({{{\mathit{\boldsymbol{\hat y}}}_P}^T{{\mathit{\boldsymbol{\hat y}}}_Q}} \right) $  (28) 
Similarly, Eq. (28) is computed through
$\mathit{\boldsymbol{R}}_{P, Q}^h = \frac{\rho }{{{{\left \mathit{\boldsymbol{z}} \right}^2}\left {{\mathit{\boldsymbol{y}}_P}} \right}}\left({W_{{\mathit{\boldsymbol{y}}_Q}}^{ext, 2d + 2}\left({\mathit{\boldsymbol{y}}_P^*} \right) + W_{{\mathit{\boldsymbol{y}}_Q}}^{ext, 2d + 1}\left({\mathit{\boldsymbol{y}}_P^*} \right)} \right) $  (29) 
Lcurve is the plot of the residual norm
$\left({\mu (\alpha), \lambda (\alpha)} \right) = \left({\log \left({{{\left\ {\mathit{\boldsymbol{A}}{{\hat x}_\alpha }  \mathit{\boldsymbol{l}}} \right\}_P}} \right), \log \left({{{\left\ {{{\mathit{\boldsymbol{\hat x}}}_\alpha }} \right\}_K}} \right)} \right) $  (30) 
The corner point of Lcurve 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 Lcurve is derived as (Hansen et al., 1993)
$k\left(\alpha \right) = \frac{{\mu '\lambda ''  \mu ''\lambda '}}{{{{\left({{{\left({\mu '} \right)}^2} + \left({\lambda '} \right)} \right)}^{3/2}}}} $  (31) 
where μʹ, λʹ, μʺ and λʺ are the first and secondorder derivatives of μ(α) and λ(α), respectively.
It is noticeable that the Lcurve 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 Lcurve 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 estimationVCE approach could also be used to properly determine the regularization parameter (Koch and Kusche, 2002). We rewrite Eq. (13) as (Klees et al., 2008)
$\Phi \left(\mathit{\boldsymbol{x}} \right) = \sum\limits_{p = 1}^J {\left\ {{\mathit{\boldsymbol{e}}_p}} \right\_{\mathit{\boldsymbol{C}}_p^{  1}}^2} + \left\ {{\mathit{\boldsymbol{e}}_{J + 1}}} \right\_{C_{J + 1}^{  1}}^2 $  (32) 
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 variancecovariance matrix of observations for this group are expressed as (Klees et al., 2008)
${\mathit{\boldsymbol{e}}_{J + 1}} = {\mathit{\boldsymbol{I}}_K}\mathit{\boldsymbol{x}}, \;\;D\left({{\mathit{\boldsymbol{e}}_{J + 1}}} \right) = {\mathit{\boldsymbol{C}}_{J + 1}} = \frac{1}{\alpha }{\mathit{\boldsymbol{R}}^{  1}} $  (33) 
where I_{K} is the K × K identity matrix, σ^{2}_{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
$\mathit{\boldsymbol{\hat x}} = {\mathit{\boldsymbol{N}}^{  1}}\mathit{\boldsymbol{W}} $  (34) 
where the normal equation matrix and righthand side vector are
$\begin{array}{l} \mathit{\boldsymbol{N}} = \sum\limits_{p = 1}^{J + 1} {\mathit{\boldsymbol{A}}_p^T\mathit{\boldsymbol{C}}_p^{  1}{\mathit{\boldsymbol{A}}_p} = \sum\limits_{p = 1}^J {\mathit{\boldsymbol{A}}_p^T\mathit{\boldsymbol{C}}_p^{  1}{\mathit{\boldsymbol{A}}_p}} + \alpha \mathit{\boldsymbol{R}}}, \;\;\mathit{\boldsymbol{W}}\\ \quad \, = \sum\limits_{p = 1}^J {\mathit{\boldsymbol{A}}_p^T\mathit{\boldsymbol{C}}_p^{  1}{\mathit{\boldsymbol{l}}_p}} \end{array} $  (35) 
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)
$\mathit{\boldsymbol{\hat \sigma }}_{J + 1}^2 = \frac{{\mathit{\boldsymbol{v}}_{J + 1}^T\mathit{\boldsymbol{Q}}_{J + 1}^{  1}{\mathit{\boldsymbol{v}}_{J + 1}}}}{{{r_{J + 1}}}} $  (36) 
And regularization parameter α is estimated as
$\alpha {\rm{ = 1/ }}\hat \sigma _{J + 1}^2 = 1/\left({\frac{{\mathit{\boldsymbol{v}}_{J + 1}^T\mathit{\boldsymbol{Q}}_{J + 1}^{  1}{\mathit{\boldsymbol{v}}_{J + 1}}}}{{{r_{J + 1}}}}} \right) $  (37) 
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 illconditioned 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 approachLcurve 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 Lcurve 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 socalled 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_{res}^i = N_{GPSL}^i  N_{SRBF}^i $  (38) 
And the standard deviation of these differences is computed by (Wittwer, 2010)
${N_{std}} = \sqrt {\frac{{\sum\limits_{i = 1}^m {{{\left({N_{res}^i  \overline {{N_{res}}} } \right)}^2}} }}{{m  1}}} $  (39) 
and
Thus, the regularization parameter that minimizes Eq. (39) is regarded as the optimal one, which could be used for cross validation when Lcurve 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 leaveoneout 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 groundbased control data (e.g., GPS/leveling data) are used instead of the leaving out observation. In this way, the highquality 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 highquality 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 modelBased on RCR methodology, the long and shortwavelength 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 highresolution 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^{3} and ρ_{sea}=1.64 g/cm^{3}, and the RTM corrections were computed based on tesseroid integrals (Heck and Seitz, 2006).
Download:


Terrestrial, shipboard and airborne gravity data are combined for regional gravity field modeling. Pointwise 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 lowpass filter is used for reducing the effect of highfrequency noise in these two data sets. Moreover, DTU10GRA 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 longwavelength parts from the original data, which are shown in Fig. 2. The corresponding statistics are displayed in Table 1.
Download:


To choose the optimal regularization matrix in dealing with the illconditioned 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 illconditioned. The performance of various Tikhonov regularization matrices is investigated, and the Lcurve is applied to estimate the optimal regularization parameters.
Figure 2 shows the Lcurve plots for various regularization matrices, where the optimal regularization parameter with zeroorder regularization is between 10^{11} and 10^{10}. While, for the firstorder 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 diagonaldominated (the regularization matrices are shown in log scale), there are still differences between the zero and firstorder regularization matrix, which lead to the difference between the estimated regularization parameters. To further demonstrate the corner point of Lcurve, 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 firstorder 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 firstorder 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 firstorder regularization, we find the application of firstorder regularization derives better results, and the accuracy of the corresponding geoid reaches 1.3 cm. For the solutions computed from zeroorder 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 firstorder regularization matrix may be more preferable in the regional gravity field modeling based on Poisson wavelets.
Download:


Download:


The Lcurve, 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 firstorder 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 Lcurve and MSTD approaches, see Fig. 3, Tables 2 and 3. Compared to the Lcurve 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 Lcurve 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.
Download:


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 tradeoff 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 illconditioned problem, firstorder 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 illcondition 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 illconditioned 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 overregularization. Figure 7 shows the evaluation results of the gravimetric geoid computed from the optimal Poisson wavelets' network in the Netherlands.
Download:


Download:


The Tikhonov regularization method in dealing with the illconditioned 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 zeroorder regularization, the application of the firstorder 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., Lcurve, 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 firstorder 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 Lcurve 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 biascorrected 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 highquality 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.
ACKNOWLEDGMENTSThanks 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. 150208), and the State Scholarship Fund from Chinese Scholarship Council (No. 201306270014). The final publication is available at Springer via https://doi.org/10.1007/s1258301707713.
REFERENCES CITED
Albertella, A., Sansò, F., Sneeuw, N., 1999. BandLimited Functions on a Bounded Spherical Domain:The Slepian Problem on the Sphere. Journal of Geodesy, 73(9): 436447. DOI:10.1007/PL00003999 
Andersen, O. B., 2010. The DTU10 Gravity Field and Mean Sea Surface. Second International Symposium of the Gravity Field of the Earth (IGFS2). Fairbanks, Alaska

Chambodut, A., Panet, I., Mandea, M., et al., 2005. Wavelet Frames:An Alternative to Spherical Harmonic Representation of Potential Fields. Geophysical Journal International, 163(3): 875899. DOI:10.1111/gji.2005.163.issue3 
Girard, A., 1989. A Fast nMonteCarlo CrossValidationo Procedure for Large Least Squares Problems with Noisy Data. Numerische Mathematik, 56(1): 123. 
Guo, D. M., Bao, L. F., Xu, H. Z., 2015. Tectonic Characteristics of the Tibetan Plateau Based on EIGEN6C2 Gravity Field Model. Earth ScienceJournal of China University of Geosciences, 40(10): 16431652. DOI:10.3799/dqkx.2015.148 
Hansen, P. C., Jensen, T. K., Rodriguez, G., 2007. An Adaptive Pruning Algorithm for the Discrete LCurve Criterion. Journal of Computational and Applied Mathematics, 198(2): 483492. DOI:10.1016/j.cam.2005.09.026 
Hansen, P. C., O'Leary, D. P., 1993. The Use of the LCurve in the Regularization of Discrete IllPosed Problems. SIAM Journal on Scientific Computing, 14(6): 14871503. DOI:10.1137/0914086 
Hashemi, Farahani H., Ditmar, P., Klees, R., et al., 2013. The Static Gravity Field Model DGM1S from GRACE and GOCE Data:Computation, Validation and an Analysis of GOCE Mission's Added Value. Journal of Geodesy, 87(9): 843867. DOI:10.1007/s0019001306503 
Heck, B., Seitz, K., 2006. A Comparison of the Tesseroid, Prism and PointMass Approaches for Mass Reductions in Gravity Field Modelling. Journal of Geodesy, 81(2): 121136. 
Heiskanen, W. A., Moritz H., 1967. Physical Geodesy. WH Freeman and Co., San Francisco

Hirt, C., 2013. RTM Gravity ForwardModeling Using Topography/Bathymetry Data to Improve HighDegree Global Geopotential Models in the Coastal Zone. Marine Geodesy, 36(2): 183202. 
Hoerl, A., Kennard, R., 1970. Ridge Regression:Biased Estimation for Nonorthogonal Problems. Technometrics, 42(1): 8086. 
Holschneider, M., IglewskaNowak, I., 2007. Poisson Wavelets on the Sphere. Journal of Fourier Analysis and Applications, 13(4): 405419. DOI:10.1007/s0004100669099 
Johnston, P. R., Gulrajani, R. M., 2000. Selecting the Corner in the LCurve Approach to Tikhonov Regularization. IEEE Transactions on Biomedical Engineering, 47(9): 12931296. DOI:10.1109/10.867966 
Klees, R., Tenzer, R., Prutkin, I., et al., 2008. A DataDriven Approach to Local Gravity Field Modelling Using Spherical Radial Basis Functions. Journal of Geodesy, 82(8): 457471. DOI:10.1007/s0019000701963 
Koch, K. R., 1987. Bayesian Inference for Variance Components. Manuscr. Geod., 12: 309313. 
Koch, K. R., Kusche, J., 2002. Regularization of Geopotential Determination from Satellite Data by Variance Components. Journal of Geodesy, 76(5): 259268. DOI:10.1007/s001900020245x 
Kusche, J., Klees, R., 2002. Regularization of Gravity Field Estimation from Satellite Gravity Gradients. Journal of Geodesy, 76(6/7): 359368. 
Luthcke, S. B., Sabaka, T. J., Loomis, B. D., et al., 2013. Antarctica, Greenland and Gulf of Alaska LandIce Evolution from an Iterated GRACE Global Mascon Solution. Journal of Glaciology, 59(216): 613631. DOI:10.3189/2013JoG12J147 
Rummel, R., Schwarz, K. P., Gerstl, M., 1979. Least Squares Collocation and Regularization. Bulletin Géodésique, 53(4): 343361. DOI:10.1007/BF02522276 
Simons, F. J., Dahlen, F. A., 2006. Spherical Slepian Functions and the Polar Gap in Geodesy. Geophysical Journal International, 166(3): 10391061. DOI:10.1111/gji.2006.166.issue3 
Tenzer, R., Klees, R., 2008. The Choice of the Spherical Radial Basis Functions in Local Gravity Field Modeling. Studia Geophysica et Geodaetica, 52(3): 287304. DOI:10.1007/s1120000800222 
Tikhonov, A. N., 1963. Regularization of Incorrectly Posed Problems. Sov. Math. Dokl., 4(1): 16241627. 
Wittwer, T., 2010. Regional Gravity Field Modelling with Radial Basis Functions: [Dissertation]. Delft University of Technology, Delft

Wu, Y. H., Luo, Z. C., 2016. The Approach of Regional Geoid Refinement Based on Combining MultiSatellite Altimetry Observations and Heterogeneous Gravity Data Sets. Chinese J. Geophys., 59(5): 15961607. 
Wu, Y. H., Luo, Z. C., Chen, W., et al., 2017. HighResolution Regional Gravity Field Recovery from Poisson Wavelets Using Heterogeneous Observational Techniques. Earth, Planets and Space, 69(34): 115. 
Wu, Y. H., Luo, Z. C., Zhou, B. Y., 2016. Regional Gravity Modelling Based on Heterogeneous Data Sets by Using Poisson Wavelets Radial Basis Functions. Chinese J. Geophys., 59(3): 852864. 
Xu, P. L., 1992. The Value of Minimum Norm Estimation of Geopotential Fields. Geophysical Journal International, 111(1): 170178. DOI:10.1111/gji.1992.111.issue1 
Xu, P. L., 1998. Truncated SVD Methods for Discrete Linear IllPosed Problems. Geophysical Journal International, 135(2): 505514. DOI:10.1046/j.1365246X.1998.00652.x 
Xu, P. L., 2009. Iterative Generalized CrossValidation for Fusing Heteroscedastic Data of Inverse IllPosed Problems. Geophysical Journal International, 179(1): 182200. DOI:10.1111/gji.2009.179.issue1 
Xu, P. L., Rummel, R., 1994. Generalized Ridge Regression with Applications in Determination of Potential Fields. Manuscr. Geod., 20: 820. 
Xu, P. L., Shen, Y. Z., Fukuda, Y., et al., 2006. Variance Component Estimation in Linear Inverse IllPosed Models. Journal of Geodesy, 80(2): 6981. DOI:10.1007/s0019000600321 
Xu, S. F., Chen, C., Du, J. S., et al., 2015. Characteristics and Tectonic Implications of Lithospheric Density Structures beneath Western Junggar and Its Surroundings. Earth ScienceJournal of China University of Geosciences, 40(9): 15561565. DOI:10.3799/dqkx.2015.140 