**Three-Dimensional Regularized Focusing Migration: A Case Study from the Yucheng Mining Area, Shandong, China**

#### **Yidan Ding 1,\* , Guoqing Ma 1,\*, Shengqing Xiong <sup>2</sup> and Haoran Wang <sup>1</sup>**


Received: 23 March 2020; Accepted: 21 May 2020; Published: 22 May 2020

**Abstract:** Gravity migration is a fast imaging technique based on the migration concept to obtain subsurface density distribution. For higher resolution of migration imaging results, we propose a 3D regularized focusing migration method that implements migration imaging of an entire gravity survey with a focusing stabilizer based on regularization theory. When determining the model parameters, the iterative direction is chosen as the conjugate migration direction, and the step size is selected on the basis of the Wolfe–Powell conditions. The model tests demonstrate that the proposed method can improve the resolution and precision of imaging results, especially for blocky structures. At the same time, the method has high computational efficiency, which allows rapid imaging for large-scale gravity data. It also has high stability in noisy conditions. The developed novel method is applied to interpret gravity data collected from the skarn-type iron deposits in Yucheng, Shandong province. Migration results show that the depth of the buried iron ore in this area is 750–1500 m, which is consistent with the drilling data. We also provide recommendations for further mineral exploration in the survey area. This method can be used to complete rapid global imaging of large mining areas and it provides important technical support for exploration of deep, concealed deposits.

**Keywords:** 3D focusing migration; regularization; conjugate migration direction; fast imaging; skarn-type iron deposits

#### **1. Introduction**

Gravity exploration is an effective method for obtaining information regarding subsurface density distribution; hence, it is an important tool in geophysical exploration and it is widely used to estimate target mineralized zones. Subsurface density distribution can be obtained via inversion of gravity data. However, from the perspective of classical theory, all geophysical inversions, including gravity inversion are ill-posed because their solutions are either non-unique or unstable. Nevertheless, geophysicists can solve this problem to some extent by adopting different types of regularization algorithms. Tikhonov [1] established regularization theory (also see [2–4]). The regularized inversion problem can be solved by finding a solution to the optimization problem that minimizes the objective functional, which can be expressed as [5–7]

$$P(\mathfrak{m}) = P\_d(\mathfrak{m}) + \lambda P\_m(\mathfrak{m}) \tag{1}$$

where *P<sup>d</sup>* (*m*) is the misfit functional, *Pm*(*m*) is the stabilizing functional that ensures the inversion result contains a priori information of the model parameters, and λ is the regularization parameter that balances the effects of *P<sup>d</sup>* (*m*) and *Pm*(*m*).

On the basis of this equation, a variety of 3D regularized inversion methods for gravity data have been developed with different regularization terms [6,8], such as regularized smooth inversion based on the L2 norm (e.g., [9,10]), sparse inversion based on the L1 norm (e.g., [11]), total variation inversion (e.g., [12,13]) and focusing inversion (e.g., [14,15]). The focusing inversion method obtains relatively convergent results with sharp boundaries, which makes it is very suitable for mineral interpretation. Several classical techniques have also been developed to automatically estimate an appropriate regularization parameter [16], such as the generalized cross-validation criteria [17], L-curve criteria [18] and adaptive regularization method [6]. For the adaptive regularization method, the regularization parameter λ successively decreases during the process of the iteration, thus accelerating convergence and improving robustness.

Inversion methods can be inefficient in actual data processing due to the amount and time of calculation they require, and they are not able to complete the overall 3D density imaging of large regions [19]. These inefficiencies have restricted the study of regional structures and deep geological problems. Therefore, there is an urgent need to improve computing efficiency. Zhdanov introduced the concept of potential field migration and demonstrated its application for rapid 2D and 3D imaging of gravity and gradiometry surveys [6,20–22]. For gravity and its gradient fields, migration imaging is based on the direct integral transformation of gravity and its gradient data (which can also be mathematically described as the action of the adjoint operator on the observed data). The potential field migration method is fast and stable, and can be used for real-time imaging. However, the imaging results appear divergent [23]. This method was later extended to iterative migration [24], in which case adjoint operators are used iteratively. The forward model was introduced to estimate the accuracy of the migration results by computing the residual of the observed data and calculated data during each iteration. When the residual is less than a given value, this migration imaging model is adopted as the final result; thus, the imaging accuracy is improved. Multi-component joint iterative migration was developed subsequently and further improved the vertical resolution [25].

The potential field migration method instantaneously provides an estimation of the density distribution in a single step, while iterative migration improves the spatial resolution. However, the introduction of an iterative process reduces the computational efficiency. In order to maintain the high efficiency of migration and improve the resolution and stability, we propose a new method called 3D regularized focusing migration. We introduce an objective functional similar to that of regularized inversion and use the focusing stabilizer, which helps to generate a sharp and focused image of the model. When solving the optimization problem of minimizing the objective functional, the search direction of our algorithm is the conjugate migration direction instead of the conjugate gradient direction that is used for conventional regularized inversion. Good results can be obtained at the first migration and using the conjugate migration direction results in faster convergence than using the conjugate gradient direction in inversion. Thus, we can obtain relatively convergent results with fewer iterations. Unfortunately, the line searches are still performed multiple times while minimizing the objective functional. In order to improve the efficiency of a line search under the premise of ensuring algorithm convergence, the iterative step size is selected based on the inexact line search step size criterion proposed by Wolfe and Powell [26,27]. Adopting the Wolfe–Powell conditions allows a misfit to take less time to reach the given threshold when the model scale is relatively large.

The proposed method is efficient and stable, and can provide high-resolution subsurface density distribution imaging in real-time. It is well suited for evaluating mining targets because our efficient algorithm allows for fast imaging without being constrained by millions of cells from large mining areas. Our regularization parameter selection criterion is applicable to the case where the observed data are contaminated with Gaussian noise of uniform but unknown standard deviation, just like the

real data. The focusing stabilizer can distinguish the boundary of each geological body, especially dense ore bodies.

In the following sections, we briefly introduce the principle of the regularized focusing migration method, and apply it to several synthetic models to verify the rationale and performance of our method. Then it is applied to the interpretation of gravity data collected in Yucheng, Shandong province. This province has multiple potential mining targets for skarn-type iron ore deposits that have no surface expression. We discuss the geological background of this area, interpret recently acquired gravity data, and provide a road map for further mineral exploration.

#### **2. Methodology**

#### *2.1. Forward Modeling of Gravity Anomalies*

The purpose of gravity forward modeling is to calculate the effect of the subsurface density distribution on the observation surface. The modeling domain is discretized into rectilinear cuboids of constant density (Figure 1). Discrete forward modeling for gravity anomalies produced by arbitrary density distribution can be expressed in general matrix notation as:

$$d = A m\_{\prime} \tag{2}$$

where *d* is a vector of the observed data of order *n*; *m* is a vector of the model parameters of order *m*; *A*, which is called the sensitivity matrix, is a matrix of size *n* by *m* formed by the corresponding gravity field kernels. If the subsurface modeling region is divided into *n<sup>x</sup>* × *n<sup>y</sup>* × *n<sup>z</sup>* geometric cuboids, then *n* = *n<sup>x</sup>* × *ny*, *m* = *n<sup>x</sup>* × *n<sup>y</sup>* × *nz*. ݉ ݊ ݉ ݊<sup>௫</sup> × ݊<sup>௬</sup> × ݊<sup>௭</sup> ݊=݊<sup>௫</sup> × ݊௬ ݉ = ݊<sup>௫</sup> × ݊<sup>௬</sup> × ݊<sup>௭</sup>

(ݖ ,ݕ ,ݔ) <sup>ଵ</sup>ߞ ,ଵߟ ,ଵߝ) <sup>ଶ</sup>ߞ ,ଶߟ ,ଶߝ) ( ) **Figure 1.** Representation of the subsurface density distribution. (**a**) A layer of rectilinear cuboids; (**b**) a single cuboid. The observation position is described by coordinate (*x*, *y*, *z*). The coordinates of cuboid corners can be completely described by two vertices (ε1, η1, ζ<sup>1</sup> ) and (ε2, η2, ζ2).

The gravity anomalies generated by any geometric cuboid M at observation point D can be expressed as [28]

$$\mathbf{d} = \gamma \rho\_M \left\{ (\varepsilon - \mathbf{x}) \ln[\mathbf{R} - (\eta - y)] + (\eta - y) \ln[\mathbf{R} - (\varepsilon - \mathbf{x})] - (\zeta - \mathbf{z}) \tan^{-1} \frac{(\varepsilon - \mathbf{x})(\eta - y)}{(\zeta - \mathbf{z})\mathbf{R}} \right\} \begin{vmatrix} \varepsilon\_2 & \eta\_2 & \zeta\_2 \\ \varepsilon\_1 & \eta\_1 & \zeta\_1 \end{vmatrix} \zeta\_2 \tag{3}$$

ଵߝ ଵߟ ଵߞ ெߩ ߛ <sup>ଶ</sup>ߞ <sup>ଵ</sup>ߞ <sup>ଶ</sup>ߟ <sup>ଵ</sup>ߟ <sup>ଶ</sup>ߝ <sup>ଵ</sup>ߝ (ݔ−ߝ)ඥ = ܴ (ݕ−ߟ) + <sup>ଶ</sup> <sup>ଶ</sup>)ݖ − ߞ) + <sup>ଶ</sup> where γ is the gravitational constant; ρ*<sup>M</sup>* is the density of the geometric cuboid; *(x, y, z)* are the coordinates of surface observation point D; and ε1, ε2, η1, η2, ζ1, ζ<sup>2</sup> are coordinates to determine the subdivision cuboid space, *R* = q (ε − *x*) <sup>2</sup> + (η − *y*) <sup>2</sup> + (ζ − *z*) 2 .

#### *2.2. Principles of of Gravity Field Migration*

 Given that gravity field dataset *d obs* is observed on the observation surface S, the objective is to determine the subsurface density distribution, *m*. Zhdanov proposed a fast migration imaging method on the basis of the Newton method. Its imaging principle is described briefly below [6].

ଵି

−= ଵ ఒ

Within the framework of the Newton method, a solution to the problem of parameter functional minimization can be found with a single step. When both the initial model and the prior model of anomalous density distribution are equal to zero, the solution can be expressed as:

$$\mathcal{m} = -\frac{1}{\lambda} \left( \mathbf{W}\_m^\* \mathbf{W}\_m \right)^{-1} l\_{\mathbf{0}} \tag{4}$$

where λ is the regularization parameter; *W<sup>m</sup>* is a model weighting matrix; *W*<sup>∗</sup> *<sup>m</sup>* is the adjoint matrix of *Wm*; and *l<sup>0</sup>* is the direction of steepest descent that is introduced as a result of applying adjoint gravity operator, *A \** , to the observed gravity data, *d obs*:

*l<sup>0</sup>* = −*A \*d obs* , (5)

where adjoint gravity operator, *A \** , is specified as

$$A^\* = \gamma \iint \frac{d}{\left| \mathbf{r}' - \mathbf{r} \right|^3} (\mathbf{z}' - \mathbf{z}) d\mathbf{s}.\tag{6}$$

Let <sup>1</sup> <sup>λ</sup> = *k*, and substitute Equation (5) into Equation (4), then, migration density can be rewritten as

$$m = k \left(\mathcal{W}\_m^\* \mathcal{W}\_m\right)^{-1} A^\* d^{obs} = k \omega\_z^2 A^\* d^{obs}.\tag{7}$$

Weighting operator, ω*z*, is constructed based on the square root of the integrated sensitivity of the data to the density, *Sz*:

$$
\omega\_{\mathbf{z}} = \sqrt{\mathbf{S}\_{\mathbf{z}}}.\tag{8}
$$

Unknown coefficient, *k*, can be determined by linear line search:

*k* = k *A* <sup>ω</sup>*\*d obs* k 2 k *A* ω*A* <sup>ω</sup>*\*d obs* k 2 , (9)

where *A* <sup>ω</sup> is the weighted sensitivity matrix and *A* ω*\** is the weighted adjoint gravity operator.

Migration gives us a method to analytically determine the model parameters, so rather than running an iterative algorithm, we find the optimal value for *m* in a single step. It is obvious that migration cannot obtain very accurate imaging results, although it has many advantages in terms of efficiency. Therefore, the regularized focusing migration method is introduced in this paper.

#### *2.3. Regularized Focusing Migration*

We can consider *k* in Equation (7) as a step size and ω<sup>2</sup> *zA \** as the migration direction. In this case, the migration is similar to the first iteration of a general inversion problem.

In regularization theory, the solution of linear inverse Equation (2) is substituted by a minimization of Tikhonov objective functional [6,7]:

$$P(m) = \left\| \left| \left. \mathbf{W\_{d}}(A\mathbf{m} - d) \right\| \right|^{2} + \lambda \left\| \left. \mathbf{W\_{m}}(m - m\_{apr}) \right\| \right\|^{2} \to \min,\tag{10}$$

where *mapr* is a prior model and here, *mapr* is equal to zero; *W<sup>d</sup>* is the data weighting matrix; and *W<sup>m</sup>* is the model weighting matrix.

We introduce a focusing weighting matrix *We*, and the reweighted objective functional can be rewritten as:

$$P(m^{\omega}) = \|A^{\omega}m^{\omega} - d^{\omega}\|^2 + \lambda \parallel m^{\omega} \parallel^2,\tag{11}$$

where

$$A^{\omega} = \mathbf{W}\_d \mathbf{A} \mathbf{W}\_m^{-1} \mathbf{W}\_e^{-1} \tag{12}$$

*Minerals* **2020**, *10*, 471

$$d^w = \mathbf{W}\_d d\_\prime \tag{13}$$

$$
\hbar \mathbf{m}^{\omega} = \mathbf{W}\_{\varepsilon} \mathbf{W}\_{m} \mathbf{m}.\tag{14}
$$

There are different options for weighting matrices *We*, *Wd*, *Wm*. Here, we chose the following forms [6]:

$$\mathcal{W}\_{\mathbf{e}} = \operatorname{diag} \frac{1}{\sqrt{(\mathfrak{m})^2 + \mathfrak{e}^2}} \,\,\,\,\tag{15}$$

$$\mathbf{W}\_d = \operatorname{diag} \left( \mathbf{A} \mathbf{A}^T \right)^{\frac{1}{2}},\tag{16}$$

$$\mathbf{W\_m} = \operatorname{diag} \left( \mathbf{A^T A} \right)^{\frac{1}{2}},\tag{17}$$

Here e is a focusing parameter and it is selected as 0.08 by trial and error [29]. For the *n* th iteration (*n* = 1, 2, 3...), the objective functional can be written as:

$$P(m\_n^{\omega\_n}) = \left| \|A^{\omega\_n} m\_n^{\omega\_n} - d^{\omega\_n} \|\|^2 + \lambda\_n \|\| m\_n^{\omega\_n} \|\|^2 \right. \tag{18}$$

When determining the model parameters, unlike the conjugate gradient method used in the general inversion problem, we use the conjugate migration direction method, which implies that the search direction of each iteration is along the conjugate migration direction. Then for regularized migration, according to Equation (5), we define the migration direction as

$$I\_{\mathfrak{m}} = I(\mathfrak{m}) = A^\*(A\mathfrak{m} - d) + \lambda \mathfrak{m},\tag{19}$$

and the conjugate migration direction is defined as the linear combination of the migration direction for this step and the conjugate migration direction for the previous step. Hence the conjugate migration direction for the current model parameter *m* ω*n n* is:

$$
\overline{\mathbf{I}\_n^{\omega\_n}} = -\mathbf{I}\_n^{\omega\_n} + \beta\_n \overline{\mathbf{I}\_{n-1}^{\omega\_{n-1}}},\tag{20}
$$

where

$$\mathbf{I}\_{\mathbf{k}}^{\omega\_{\rm n}} = \mathbf{I}^{\omega\_{\rm n}}(m\_{\rm n}^{\omega\_{\rm n}}) = \mathbf{A}^{\omega\_{\rm n}}(\mathbf{A}^{\omega\_{\rm n}}m\_{\rm n}^{\omega\_{\rm n}} - \mathbf{d}^{\omega}) + \lambda\_{\rm n}(m\_{\rm n}^{\omega\_{\rm n}}),\tag{21}$$

$$\beta\_n = \frac{||\mathbf{I}\_n^{\omega\_n}||^2}{||\mathbf{I}\_{n-1}^{\omega\_{n-1}}||^2}. \tag{22}$$

Specifically, the conjugate migration direction is defined as the migration direction when *n* = 0, namelye*I* ω0 <sup>0</sup> = −*I* ω0 0 .

The successive line search in the conjugate migration direction can be written as:

$$
\mathfrak{m}\_{n+1}^{\omega\_n} = \mathfrak{m}\_n^{\omega\_n} + k\_n \overline{I}\_{\alpha\_n}^{\omega\_n}.\tag{23}
$$

Here, *kn* is the iteration step. We take the Wolfe–Powell inexact line search conditions instead of the exact line search criterion shown in Equation (9) of the conventional iterative migration to obtain the step size [30]:

$$P\left(\mathfrak{m}\_{n}^{\omega\_{n}} + k\_{l}\widetilde{I}\_{n}^{\omega\_{n}}\right) \le P\left(\mathfrak{m}\_{n}^{\omega\_{n}}\right) + \delta k\_{l}\mathbf{I}\_{n}^{\omega\_{n}T}\widetilde{I}\_{n}^{\omega\_{n}}.\tag{24}$$

$$\left| I^{\omega\_n} (m\_n^{\omega\_n} + k\_n \overline{I}\_n^{\omega\_n})^T \overline{I}\_n^{\omega\_n} \right| \le \sigma \left| I\_n^{\omega\_n T} \overline{I}\_n^{\omega\_n} \right| \tag{25}$$

where 0 < δ < σ < 1. Furthermore, the sufficient descent property, namely,

*I* ω*n n T*e*I* ω*n <sup>n</sup>* ≤ −*c* k *I* ω*n n* k 2 , (26)

where *c* is a positive constant, is crucial to ensure the global convergence of the optimization problem with the inexact line search techniques.

The actual model parameters can be converted as

$$\mathbf{m}\_{n+1} = \mathbf{W}\_m^{-1} \mathbf{W}\_{e\_n}^{-1} \mathbf{m}\_{n+1'}^{\omega\_n} \tag{27}$$

where *We<sup>n</sup>* is the focusing weighting matrix updated according to the *m<sup>n</sup>* obtained on the (*n* − 1) st iteration.

Reweight *mn*+*<sup>1</sup>* to get *m* ω*n*+*<sup>1</sup> n*+*1* :

$$
\mathcal{W}\_{n+1}^{\omega\_{n+1}} = \mathcal{W}\_m \mathcal{W}\_{\mathfrak{e}\_{n+1}} \mathcal{m}\_{n+1}.\tag{28}
$$

To ensure that the misfit functional converges to the global minimum, the following procedures are adopted to reduce the regularization factor [31]. First, let

$$p\_n = \frac{||\,m\_{n+1}^{\omega\_{n+1}}||^2}{||\,m\_n^{\omega\_n}||^2},\tag{29}$$

and then update the regularization factor as follows:

$$
\lambda\_{n+1} = \begin{cases}
\lambda\_n & p\_n < 1 \\
\lambda\_n / p\_n & p\_n > 1
\end{cases}
\quad . \tag{30}
$$

If the misfit functional is not shrinking fast enough, namely k *r* ω*n n* k <sup>2</sup> − k *r* ω*n*+*<sup>1</sup> n*+*1* k <sup>2</sup>< 0.01 k *r* ω*n n* k 2 , the value of λ*n*+*<sup>1</sup>* calculated by Equation (30) should be reduced to q times of the original value. We selected *q* = 0.6 by trial and error. It is clear that for the solution obtained in each iteration, if the value of regularization parameter is continuously reduced, then the value of the corresponding objective functional also continues to decrease. This guarantees that the objective functional will converge faster on the premise of obtaining a stable solution.

When the misfit reaches a given range:

$$\left| \text{Misfit} = \right| \left| \frac{\mathbf{A} \mathbf{m}\_{\text{N}} - \mathbf{d}}{\mathbf{d}} \right| \left| \cdot 100\% \le \varepsilon\_{0\prime} \tag{31}$$

the iteration process terminates and outputs *mN*, otherwise it proceeds to the next iteration.

#### **3. Synthetic Examples**

We examined the feasibility and effectiveness of 3D regularized focusing migration with several typical synthetic models. First, we designed a single cube with a density contrast of 1 g/cm<sup>3</sup> , with the center located at a depth of 150 m. The three-dimensional view of the subsurface model and its surface gravity anomalies are shown in Figure 2a,b. The subsurface inversion domain was divided into 50 × 50 × 25 regular cells.

Figure 2c,d show the cross-sections through the density volume resulting from the regularized focusing migration at *y* = 250 m and *z* = 150 m, respectively. The density extremum of the imaging results corresponded well with the center of the anomalous body, and the high-density distribution precisely outlined the anomalous body, which demonstrated the feasibility of the method.

Our second synthetic example includes three anomalous bodies representing a dike complex and a prism (Figure 3a). Figure 3b shows the gravity anomaly computed by that model.

ାଵߣ

=

శ∥ ∥శ࣓ మ

∥ <sup>∥</sup>࣓ మ

1 < ⁄ߣ

ߝ ≥ 100% .

 ∥ ࣓ ∥

<sup>ଶ</sup>−∥ ା ∥ శ࣓

∥ 0.01 ><sup>ଶ</sup>

࣓ ∥ ଶ

1 > ߣ ൜ = ାଵ<sup>ߣ</sup>

‖ି‖

‖‖

= ݐ݂݅ݏ݅ܯ

**Figure 2.** (**a**) Three-dimensional view of a single cube model; (**b**) gravity anomalies at *z* = 0 m generated from the model; (**c**) vertical cross-section of *y* = 250 m obtained from 3D regularized focusing migration; (**d**) horizontal cross-section of *z* = 150 m obtained from 3D regularized focusing migration.

**Figure 3.** (**a**) Three-dimensional view of the subsurface model; (**b**) gravity anomalies at *z* = 0 m **Figure 3.** (**a**) Three-dimensional view of the subsurface model; (**b**) gravity anomalies at *z* = 0 m generated from the model; (**c**) vertical cross-section of *y* = 250 m obtained from conventional iterative migration; (**d**) horizontal cross-section of *z* = 150 m obtained from conventional iterative migration; (**e**) vertical cross-section of *y* = 250 m obtained from 3D regularized focusing migration; and (**f**) horizontal cross-section of *z* = 150 m obtained from 3D regularized focusing migration.

Figure 3c,f show the cross-sections of density distribution obtained by conventional iterative migration and regularized focusing migration. Regularized focusing migration method was able to distinguish the trends in the individual dikes, resulting in better convergence and a more realistic output, while better imaging the deep part of the model. However, the conventional iterative migration method barely distinguished three anomalous bodies and underestimated the dips in the modeled dikes, showing divergent density distribution that was far beyond the actual extent. In contrast, 3D regularized focusing resulted in a better spatial resolution and preserved sharper and finer boundaries.

Next, we investigated how regularized focusing migration performs in noise conditions. Figure 4a shows the three-dimensional diagram for two isolated cubes with a density contrast of 1 g/cm<sup>3</sup> , located at depths of 70 m and 100 m. Figure 4b shows the surface gravity anomalies for this model. We added 5% random Gaussian noise to the observed gravity data, as shown in Figure 4c to simulate the data collected in actual surveys. The imaging results for noiseless and noise-corrupted data are shown in Figure 4d,g. The imaging results for both noiseless and noise-corrupted data were in good agreement with the true model. Although the two bodies were close in distance, which was equal to their side length, our method still distinguished the discrete bodies. The extremum of the derived density distribution corresponded well to the center of each geological body, and the horizontal positions and buried depths could be clearly delineated. However, the recovered density values were somewhat different from the actual conditions, especially for the deeper bodies. This is because the gravity field declines rapidly with depth, as shown in Equation (3). It is clear that the direct imaging of noise-containing data was still able to produce stable and smooth results. The volume images with a density cut-off greater than 0.3 g/cm<sup>3</sup> of noiseless and noise-containing data, are also shown in Figure 4h,i, respectively. The result of the direct migration of noise-containing data is very similar to that of noiseless data. Based on the results of this test, we concluded that the regularized focusing migration method has good noise resistance since the introduction of a regularization term prevents local overfitting. This method shows good performance and is able to obtain the global optimal results.

We also analyzed the efficiency of the proposed method. All tests were run on an Intel Xeon E5-2650 (Intel, Santa Clara, CA, USA) with 2.2 Ghz, 64GB of main memory. The model shown in Figure 2a was adopted for this purpose. Figure 5 shows the misfit curves as a function of the number of iterations for the regularized inversion [32], the conventional iterative migration method and our proposed method. The migration method obtains a result with a lower misfit in the first iteration, which gives rise to convergence and less calculation time. Numerical results show that the migration method can decrease the iterative number by 35% compared with the inversion method and greatly reduce the computation time.

Figure 6 shows the computation time of the regularized focusing migration method as a function of the model scale. The iteration process stopped when the misfit calculated by Equation (31) reached 5%. For models less than 6 × 6 × 6 cells, the computation time for both line search criterion is constant because the model scale is relatively small and most of the time is spent on the initialization of the work. For models larger than 6 × 6 × 6 cells, computation time is linearly dependent on the scale of the model. For computations based on the Wolfe–Powell conditions, the time that is consumed is higher than that of the exact line search with models under 16 × 16 × 16 cells, but this is the opposite for the models exceeding 16 × 16 × 16 cells. This is because the exact line search for each step involves the calculation of *A* ω*A* ω*\** , which means that the time complexity grows roughly as *nx*.*ny* 2 *nx*.*ny*.*nz* 2 , while for the Wolfe–Powell conditions, which calculate *A* <sup>ω</sup>*m* multiple times, each time the step size is determined, the time complexity is approximately *nx* <sup>2</sup>*n<sup>y</sup>* <sup>2</sup>*nz*. The solution to the underdetermined system of equations does not depend on the exact search process. When the model scale is too large, the cost of using the exact line search is relatively high. According to the Wolfe–Powell conditions, the global convergence is expected to be fast. Although the number of steps increases, the convergence speed is accelerated.

**Figure 4.** (**a**) Three-dimensional view of the subsurface model; (**b**) gravity anomalies at *z* = 0 m generated from the model; (**c**) gravity anomalies at *z* = 0 m generated from the model with 5% random Gaussian noise added; (**d**) vertical cross-section of *y* = 250 m for noiseless gravity data; (**e**) horizontal cross-section of *z* = 150 m for noiseless gravity data; (**f**) vertical cross-section of *y* = 250 m for noise-corrupted data; (**g**) vertical cross-section of *z* = 150 m for noise-corrupted data; (**h**) volume image with a density cut-off greater than 0.3 g/cm<sup>3</sup> for noiseless data; (**i**) volume image with a density cut-off greater than 0.3 g/cm<sup>3</sup> for noise-corrupted data.

To process the model with 50 × 50 × 25 cells, the runtime for our method was only 224 s compared to 2278 s for the conventional migration method and 1.8 h for the regularized inversion method. This efficiency is due to the adoption of the conjugate migration direction and the step size is based on the Wolfe–Powell conditions when the model parameters are updated during each iteration. Besides, the continuous decrease in the regularization parameter and the introduced focusing weight can also accelerate the convergence and obtain a solution that minimizes the objective functional more quickly. Therefore, when studying the regional structures and the deep geological problems, regularized focusing migration imaging can be used for global analysis without the need for data thinning. This is especially applicable to the investigation of large-scale mining areas, as in the following case study of the Yucheng mining area, Dezhou, Shandong, China.

**Figure 5.** Misfit of three different methods: regularized focusing migration, conventional iterative migration and regularized inversion.

**Figure 6.** Computation time of regularized focusing migration method for models of different sizes. The blue line represents use of the Wolfe–Powell line search conditions as shown in Equations (24)–(26), and the red line represents use of the exact line search criterion as shown in Equation (9).

#### **4. Field Example**

In 2011, the Shandong Province Mineral Resources Potential Evaluation Project revealed that the Yucheng area had metallogenic conditions that are indicative of skarn-type ore deposits. In the following years, a variety of geological surveys showed that the area had good prospecting potential, and the iron ore bodies were confirmed during drilling. The traditional interpretation of the geophysical surveys in this area was based on the analysis of gravity or magnetic anomalies in contour maps [33–35], which provided the major geophysical characteristics of the area and was used to determine the approximate extent of the deposits. However, they alone could not determine the exact locations and depths of the mineral deposits.

In the following sections, we first summarize the metallogenic geological conditions of this area, then use our high-resolution method to invert the results of the entire airborne gravity data into 3D density models, and finally, we analyze the density imaging results to identify potential mining targets. Our work provides a reference for follow-up studies and further exploration of skarn-type iron ore under a thick overburden.

#### *4.1. Geological Background of the Survey Area*

The survey area is located near the Qihe county and Yucheng county borders to the south of Dezhou city in Shandong province, and belongs to the northwestern Shandong plain. It is located in the Luxi uplift area of the North China tectonic plate (Figure 7) [36]. The area is covered by strata of Cambrian, Ordovician, Carboniferous, Permian, Neogene and Quaternary from the base to the top [33]. The surface is completely covered by thick Quaternary and Neogene strata without any bedrock outcrops [37]. All kinds of deep geological information can be inferred from geophysical data. Due to regional tectonic movements, fault structures are widespread, including NE-striking faults, NNW-striking faults and SN-striking faults, and local developed fold structures (Figure 8).

**Figure 7.** Tectonic background map of Yucheng area. The survey area is detailed in the black rectangle. 1 = the secondary structural unit; 2 = the tertiary structural unit; 3 = the quaternary structural unit; 4 = five-level structural unit; 5 = (latent) sag zone; 6 = (latent) high zone; 7 = toponym; 8 = survey area.

The carbonate rocks of the Ordovician Majiagou group are widely distributed and serve as host rocks for skarn-type iron ores. In addition, some iron-ore bodies also occur in the clastic rocks of the Carboniferous-Permian group. According to available borehole data, the depth of the Ordovician strata in the northwest reaches 3440 m or more, whereas it varies greatly in the southeastern area, ranging from 1165 m and 610–650 m in the east of Litun and in the northwest of Dazhang, respectively. There is no borehole data for the northwest of Pandian, and the depth of the Ordovician rocks is unknown [36]. In the survey area, magmatic activity was intense during the late Yanshanian period, and a significant volume of igneous rock mass was added at that time. The investigation of iron ores indicates that the ore bodies occur near the contact zone between the body mass of the late Yanshanian intermediate-basic intrusive rocks that are dominated by diorite and gabbro and the state of Ordovician limestone. In addition, a small number of ore bodies occur near the intersection of fault structures and in the xenoliths near the rock mass [35].

**Figure 8.** Geological map of the survey area. The figure is from [36]. 1 = Inferred range of concealed diorite rocks; 2 = inferred concealed fault; 3 = Yucheng-Jiyang latent sag; 4 = Lepingpu latent sag; 5 = Qihe latent high; and 6 = toponym.

#### *4.2. Data Acquisition*

The airborne gravity survey that collected gravity data in the Yucheng area, Shandong was performed by the China Aero-Geophysical Survey and Remote Sensing Center for Land and Resources (AGRS, Beijing, China). During the survey, a gravimeter was installed in the cabin of a CESSNA208 aircraft. The flight altitude was fixed at 450 m above sea level. The 45 × 35 km<sup>2</sup> survey area was flown at line spacings of 1500 m. The flight survey lines were designed to be perpendicular to the regional geological strike to ensure effective sampling. Figure 9 shows the gravity survey lines in the Yucheng area and the black rectangle delineates the region where we perform regularized focusing migration.

#### *4.3. Migration Models of Gravity Data*

We performed a terrain correction, latitude correction, and Bouguer correction on the collected gravity data to obtain Bouguer gravity anomalies (Figure 10a). The characteristics of regional gravity anomalies generally reflected the distribution of late Yanshanian intrusive rocks and the Ordovician strata uplift in the region. The iron-ore deposits of interest are shallower than 2 km because the iron ores below that have no economic mining value. The regional field below 2 km was separated from the Bouguer gravity anomaly by matched filtering [38], thus weakening the influence of deep structures on the imaging results for shallow density migration. Figure 10b shows the radial logarithmic power spectrum of the Bouguer gravity anomalies. Long-wavelength anomalies were removed to obtain the shallow local gravity anomalies. Figure 10c shows the local gravity anomalies, where three gravity highs can be seen. Four wells were drilled, including three holes in the Litun area and one hole in the

Yucheng area, and the ZK1 well in the Litun area and the ZK003 well in the Dazhang area discovered bodies of rich magnetite (Tables 1 and 2). According to borehole data, the density of magnetite was significantly higher than other rocks.

**Figure 9.** Gravity survey lines in the Yucheng area. The black rectangle shows the region used for regularized focusing migration.

**Figure 10.** (**a**) Regional gravity anomaly; (**b**) radial logarithmic power spectrum; (**c**) local gravity anomaly.




**Table 2.** Litun ZK1 drilling data.

Based on the characteristics of the stratigraphic, magmatic rocks distribution and the drilling data in the survey area, the high value of the local gravity anomalies were mainly caused by the skarn-type magnetite deposits. The positive density contrast of the migration imaging results could be used as a prospecting criterion. We performed regularized focusing migration imaging of local gravity anomalies with an initial model based on drilling data, and then delineated the spatial distribution of the iron ore under the range of positive density contrast of the imaging results. The local gravity data were inverted in a region designed as 21,500 m × 39,600 m in the X and Y directions and 2 km in the Z direction, with 88 × 98 × 50 small cells. Figure 11 shows the convergence curve of the migration. Figure 12 shows a perspective view of the migration density model with abnormal parts with a density greater than 0.7 g/cm<sup>3</sup> , which gives a clear view of the anomalous bodies. For a better representation, several horizontal cross-sections with different depths are shown in Figure 13.

**Figure 11.** Convergence curve of the migration.

#### *4.4. Geological Interpretation of the Mining Area*

The results of migration imaging suggest that there are considerable mineral resources underneath the Litun, Pandian and Dazhang anomaly areas. Ore bodies derived from our analysis are mainly distributed between a depth of 700 m and 1500 m.

Figure 14 shows the ore-body distribution superimposed on a simplified geological map of the survey area. One can see that the target mineralized zones are mainly concentrated in the vicinity of fault structures since faults can provide channels for magmatic upwelling. The iron ore in Litun (I) is buried at a depth of about 1050–1350 m, that in Pandian (II, III) is about 1150–1550 m deep, and that in Dazhang (IV) is about 750–950 m. It should be noted that the Pandian area (II, III) is rich in iron ore resources underground, but these resources are deeper and less extensive, and drilling would be more difficult than in the Litun (I) and Dazhang areas (IV). The iron ores in the Litun and Dacheung areas are characterized by their great thickness, shallow burial and easy mining, and these should be regarded as key prospecting drilling areas. When these targets are exhausted and more geological information is obtained, exploration can proceed outward, including under the Pandian area along with other nearby high-density distributions that have not yet been measured.

**Figure 12.** Perspective view of the migration density model. The red anomalous bodies show a density distribution with density greater than 0.7 g/cm<sup>3</sup> .

**Figure 13.** Horizontal cross-sections at (**a**) *z* = 650 m, (**b**) *z* = 750 m, (**c**) *z* = 850 m, (**d**) *z* = 950 m, (**e**) *z* = 1050 m, (**f**) *z* = 1150 m, (**g**) *z* = 1250 m, (**h**) *z* = 1350 m, (**i**) *z* = 1450 m and (**j**) *z* = 1550 m of imaging results.

**Figure 14.** Map of iron ore range co-rendered with the corresponding geological map. 1 = scope **Figure 14.** Map of iron ore range co-rendered with the corresponding geological map. 1 <sup>=</sup> scope of prospecting target area; 2 = drill hole; 3 = inferred range of concealed diorite rocks; 4 = inferred concealed fault; 5 = Yucheng-Jiyang latent sag; 6 = Lepingpu latent sag; 7 = Qihe latent high; 8 = toponym.

The recent well (ZK01) in this area found a body of ore that is 54.4 m thick at a depth of 1444.4 m, which is consistent with our interpretation. Therefore, the regularized focusing migration method has the ability to delineate the horizontal range of the ore body, predict its mineral depth and summarize prospecting targets.

#### **5. Conclusions**

Here, we developed 3D regularized focusing migration of large-scale gravity data with higher horizontal and vertical resolution than conventional gravity migration imaging. As validated by synthetic models, this method can determine complicated geological structures and it has high stability in noise conditions. Furthermore, the proposed algorithm based on the conjugate migration direction and the Wolfe–Powell conditions inherits the high efficiency of conventional migration and can produce real-time imaging with high accuracy. This provides a method that can delineate deep mineral resources.

We have shown how this method was efficiently used for the interpretation of the airborne gravity data collected over the Yucheng mining area, Dezhou, Shandong, China. High-resolution imaging results were obtained for the high-density body distribution, which helped to identify several potential

mining targets with no surface expression. Currently, one of the predicted targets outlined in our study is undergoing further drilling exploration based on these airborne gravity analyses and well results so far coincide with our predictions. This method offers a basis for further exploration in this area to map the skarn-type iron-ore bodies under a thick overburden.

Gravity and magnetic combination migration will be used further in this area to reduce the non-uniqueness of the results. Joint migration of multiple geophysical data is the focus of our future research.

**Author Contributions:** Conceptualization, Y.D. and G.M.; data curation, Y.D. and G.M.; formal analysis, Y.D.; funding acquisition, G.M.; investigation, Y.D.; methodology, Y.D.; project administration, G.M.; resources, S.X.; software, Y.D. and H.W.; supervision, G.M.; validation, Y.D.; visualization, Y.D.; writing—original draft, Y.D.; writing—review & editing, G.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Key Research and Development Program of China, grant numbers 2017YFC0602203 and 2017YFC0602000.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*
