** 3 Regional gravity field modeling with the point mass method**

**3.4 Point mass method with fixed positions**

converges rapidly within the smallest number of iterations, and then the NLCG method follows. The main reason can be attributed to the large descents of the objective function at the first several iterations for the LM method, in particular the first iteration in which the Hessian is computed based on Eq. (3.66b) with the given Jacobian matrix, while the Hessians for NLCG, L-BFGS, and L-BFGS-B are the identity matrices. The different behavior of the first iteration between NLCG and L-BFGS/L-BFGS-B is due to the differently implemented line search algorithms in the chosen software. For the case of using the initial point D, the required number of iterations is larger for the LM method than the other three methods. This may be due to the singularity of the Hessian computed by Eq. (3.66b). Pujol (2007) showed that the Gauss-Newton method fails to find the target point in this case. Instead, the trust region based LM method solves this problem by employing a regularization parameter and a regularization matrix, see Eqs (3.71) and (3.72), resulting in smaller search directions in the iterations, especially in the first iteration (see Figs 3.5 and 3.6).

The above discussion is only based on the simple example. At least, we can say that the performance of the four algorithms is quite similar when solving an unconstrained nonlinear problem. Often, the LM method can provide a faster convergence than the other three methods (see the cases of the model points A, B, and C); however, a slower convergence can be observed in the case of the model point D, indicating the higher dependence of the LM method on the chosen initial model point.

If the nonlinear problem is bound constrained, of course, the dedicated L-BFGS-B method should outperform the other three methods. A related numerical comparison between the four algorithms for regional gravity field modeling will be carried out in the next chapter.

**A**=

*D*_{l}*B*^{PM}(r_{1}*,***r**^{0}_{1}) *D*_{l}*B*^{PM}(r_{1}*,***r**^{0}_{2}) · · · *D*_{l}*B*^{PM}(r_{1}*,***r**^{0}* _{K}*)

*D*

_{l}*B*

^{PM}(r

_{2}

*,*

**r**

^{0}

_{1})

*D*

_{l}*B*

^{PM}(r

_{2}

*,*

**r**

^{0}

_{2}) · · ·

*D*

_{l}*B*

^{PM}(r

_{2}

*,*

**r**

^{0}

*)*

_{K}... ... . .. ...

*D*_{l}*B*^{PM}(r_{I}*,***r**^{0}_{1}) *D*_{l}*B*^{PM}(r_{I}*,***r**^{0}_{2}) · · · *D*_{l}*B*^{PM}(r_{I}*,***r**^{0}* _{K}*)

*.* (3.85b)

In the framework of the Gauss-Markov model, the unknown magnitude vector**x**can be estimated by
the linear least-squares procedures described in Section 3.3.3.

In addition to the observation equation system given by Eqs (3.84) and (3.85), (n^{0}+ 1)^{2} additional
equations, which play the role of constraints (see also Eq. (3.11)), are sometimes taken into account
for estimating the magnitudes, especially in the case when the full point mass RBFs are used, but
the input data are bandlimited. In this case, a combined observation equation system is obtained
according to

**l**
0

−

**e**
**e**^{0}

=

**A**
**C**

**x,** (3.86)

where**C**is the (n^{0}+ 1)^{2}×*K* design matrix and **Cx**= 0 is the matrix-vector product of Eq. (3.11).

Similar to the process of Tikhonov regularization described in Section 3.3.4, the least-squares solution of Eq. (3.86) can be obtained by solving the following normal equation

**A**^{T}**P**_{e}**A**

| {z }

**N**1

+λ**C**^{T}**P**_{e}^{0}**C**

| {z }

**N**2

**x**=**A**^{T}**P**_{e}**l**

| {z }

**y**

*,* (3.87)

where*λ*is a scaling factor, which is determined empirically in this thesis: let*m*1 and*m*2 be the mean
value of the diagonal elements of **N**_{1} and **N**_{2}, respectively. The scaling factor is then defined as the
ratio between*m*1 and*m*2, yielding*λ*=*m*1*/m*2. Defining ¯**N**=**N**_{1}+*λN*2, Eq. (3.87) can be rewritten
as

**Nx**¯ =**y.** (3.88)

If ¯**N** is singular, the standard Tikhonov regularization will be applied to solve the above normal
equation, resulting in

ˆ

**x*** _{α}* =

^{}

**N**¯ +

*αI*

^{}

^{−1}

**y.**(3.89)

In the following, we call the solution considering the constraints the constrained solution, and the one obtained without using the constraints the unconstrained solution. In addition, the point mass method with fixed positions is abbreviated as PM-FIX in the remaining of this thesis.

**3.4.2 Grid setup**

Besides the estimation step of this method, another vital step is to design the grids by choosing proper grid factors (see Fig. 3.3 or Section 2.7.4 for a more detailed discussion). In this context, only the geographic grid is used, and the grid extent, grid spacing, as well as grid depth will be determined sequentially. The total number of grid points (i.e., total number of point mass RBFs) must not exceed the number of observations. It is dependent on the selected grid extent and grid spacing. Klees et al.

(2008) and Wittwer (2009) found empirically that choosing a grid with an amount of RBFs equals to 1/4−1/3 of the number of observations is advisable. However, in the case of oversampling of the observations, such a grid may be too dense. Often, to avoid possible over-parameterization, the grid spacing is preferred to be chosen in such a way that all signals in the areas with smooth features

can be modeled. Then, further steps can be carried out for representing rough features, such as the local refinement used in Klees and Wittwer (2007), Klees et al. (2008), and Wittwer (2009) or adding more grids. No explicit strategy to define the optimal grid spacing can be given as it is affected by the data distribution, signal variation, as well as spatial bandwidths of the RBFs. In practice, a proper grid spacing is usually chosen by experimenting with various candidates and comparing the obtained solutions. After defining the grid factors associated with the horizontal positions of the point mass RBFs, the grid depth is to be fixed. Because of the high correlation between the grid depth and the grid spacing, some attempts have been made to give the relationship between them.

For example, Needham (1970) suggested that the depth-spacing ratio should be at least 0.8; however, Vermeer (1995) proposed a ratio of 2, and furthermore, the ratios of about 1 were used for the geoid computations in Ihde et al. (1998) and Tenzer et al. (2009).

In this thesis, we use an empirical and a heuristic approach for the grid depth determination.

The principle of this empirical approach is the same as the one for finding the proper regularization
parameter as described in Section 3.3.5. The only change is that the variable of the function becomes
the grid depth *d* below the reference sphere. Given *n* depth candidates, our goal is to find the one
with which the value of the function*RM S*(d) in the form

*RM S*(d) =
r1

*N* kBˆ**x*** _{d}*−

**qk**

^{2}(3.90)

is minimal. For the notations of the other variables in the above equation, see Eq. (3.46). It should be
noted that the linear problem for estimating ˆ**x*** _{d}*may be ill-posed. In this case, Tikhonov regularization
should be applied, even though this will make the computation procedure more complicated. The
successful application of this approach for finding a proper grid depth in local gravity field modeling
can be found in Tenzer and Klees (2008).

The second approach is based on the GCV technique, which is frequently used for finding the proper regularization parameter (see also Section 3.3.5). Considering Eq. (3.50), the GCV function to be minimized in this case is given by

*GCV* (d) = *I*kA_{d}**x**ˆ* _{d}*−

**lk**

^{2}

_{P}*e*

[I−*tr*(Q* ^{d}*)]

^{2}(3.91a)

with

**Q*** ^{d}*=

**A**

_{d}^{}

**A**

^{T}

_{d}**P**

_{e}**A**

_{d}^{}

^{−1}

**A**

^{T}

_{d}**P**

*(3.91b) for the case that the linear problem for estimating ˆ*

_{e}**x**

*is well-posed, or*

_{d}**Q*** ^{d}*=

**A**

_{d}^{}

**A**

^{T}

_{d}**P**

_{e}**A**

*+*

_{d}*α*

_{d}**R**

_{d}^{}

^{−1}

**A**

^{T}

_{d}**P**

*(3.91c) for the case of ill-posedness.*

_{e}**A**

*,*

_{d}*α*

*, and*

_{d}**R**

*in Eq. (3.91) are the design matrix, the regularization parameter, and the regularization matrix associated with the grid depth*

_{d}*d, respectively.*

**3.4.3 Grid formation**

In this section, the issue of grid formation for PM-FIX will be discussed. It can be subdivided into two groups: (1) a single-grid formation and (2) a multi-grid formation. Fig. 3.7 shows two examples of the point mass RBF locations for the cases of a single-grid formation and a two-grid formation.

In the first case, only one grid at the chosen depth is used for the modeling (see Fig. 3.7a). Such a formation is simple and quite suitable for regularly distributed observations. One main drawback of this formation is that it is not data-adaptive due to the same spatial bandwidths for all point mass

i Topography

Grid

k Sphere

d

**a)** i

Topography

Grid 2

Grid 1

k Sphere

d_{2}

d_{1}

**b)**

**Figure 3.7:**Illustration of the point mass RBF locations in the case of**a)**a single-grid formation;**b)**a two-grid
formation.

RBFs. This makes it less suitable for the case of irregularly distributed data as RBFs with the same density are in the areas with dense and sparse data. A possible way to circumvent the drawback of the single-grid formation is to set up two or more grids at different depths. In this case, the point mass RBFs at different depths seems to be suitable for approximating different signal frequencies.

This can also be interpreted from a geophysical point of view. The anomalous features with small scales but high frequencies are usually caused by anomalous masses near the Earth’s surface, such as coal and metal deposits or oil and gas reservoirs, whereas large-scale and low-frequency features are mainly caused by anomalous masses in the deep Earth, such as the Moho undulations and mantle convection. Often, the deep grid should have a larger grid extent and grid spacing, and the shallower one has a smaller extent and spacing. Regarding the issue of numerical instabilities, it is possible to choose the extents of the grids at different depths as being similar to the data area in practical applications. If the result obtained by using only one grid is sufficiently satisfactory, no more grids are needed. Otherwise, two schemes can be utilized to set up the grid formation. In the first scheme, all grids with different grid factors are defined at one time. In the second scheme, the grids are arranged sequentially from the deepest to the shallowest. In this thesis, we employ the second scheme for the numerical tests in Chapter 4, as the approaches described above for finding a proper grid depth can be easily implemented in this case.

In the light of the above discussion, the general procedures of PM-FIX with the use of a single-grid and a multi-grid formation are described as follows:

*Single-grid formation*

1 Define the grid extent according to the extents of data area and model area. If the data area is larger than the model area, set the grid extent to be similar to the data area. If the model area is the same as the data area, the grid extent is usually preferred to be larger than the data area; however, the larger the grid extent is, the more severe the numerical instabilities become.

2 Give a set of grid spacings as candidates. For each candidate, solve a linear equation system, which is constructed by the point mass RBFs on the grid at a proper depth derived from the empirical or the GCV approach.

3 Assess the quality of the solutions by comparing the predicted values to the observed values on a set of control points. The solution with the smallest RMS (or STD) error is regarded as the final solution.

To reduce the computation complexity, only one grid spacing is usually used in 2 .

*Multi-grid formation*

1 For Grid 1. Define Grid 1 with a given grid extent and grid spacing. Determine the grid depth by the empirical or the GCV approach, and then solve the linear equation system constructed by the point mass RBFs on the nodes of Grid 1 and the initial input data. Subtract the contributions of the estimated point mass RBFs from the initial input data, resulting in the residuals that will be used as input for Grid 2.

2 For Grid *i*(i≥2). If the given input residuals are sufficiently small by satisfying a predefined
criterion, go to 4 . Otherwise, define Grid*i*with a smaller grid spacing than Grid *i*−1. Note
that the extent of Grid *i* can either be the same or smaller than Grid *i*−1. Determine the
grid depth by the empirical or the GCV approach, but with the constraint that Grid *i* must
be shallower than Grid*i*−1. Solve the linear equation system constructed by the point mass
RBFs on the nodes of Grid *i* and the input residuals, and then subtract the contributions of
the estimated point mass RBFs from the input, resulting in new residuals that will be used as
input for the next grid.

3 Let *i*=*i*+ 1 and go to 2 .

4 Solve the new linear equation system constructed by all point mass RBFs located on the chosen grids and the initial input data, yielding the final solution.

Each new grid is added manually. Often, two or three grids are sufficient for practical applications.