Three-Dimensional Inversion of Multi-Component Semi-Airborne Electromagnetic Data in an Undulating Terrain for Mineral Exploration

Abstract

For the mineral exploration in complex terrain areas, the semi-airborne transient electromagnetic (SATEM) technology is one of the most powerful methods due to its high efficiency and low cost. However, since the mainstream SATEM systems only observe the component dBz/dt and the data are usually processed by simple interpretation or one-dimensional (1D) inversion, their resolutions are too low to accurately decipher the fine underground structures. To overcome these problems, we proposed a novel 3D forward and inversion method for the multi-component SATEM system. We applied unstructured tetrahedron grids to finely discretize the model with complex terrain, subsequently we used the vector finite element method to calculate the SATEM responses and sensitivity information, and finally we used the quasi-Newton method to achieve high-resolution underground structures. Numerical experiments showed that the 3D inversion could accurately recover the location and resistivities of the underground anomalous bodies under the complex terrain. Compared to a single component data, the inversion of the multi-component data was more accurate in describing the vertical boundary of the electrical structures, and preferable for high-resolution imaging of underground minerals.

1. Introduction

Geophysical methods are important tools for mineral resources exploration [1]. The conventional ground geophysical methods require humans to access and operate equipment, which is quite difficult and almost impossible for the work in dangerous mountain areas. In contrast, the airborne electromagnetic (AEM) technology can overcome the inaccessibility of terrain and achieve efficient and high-precision data to detect the underground structures [2]. However, the AEM technology generally requires a manned helicopter as a carrier, and the cost is too high. The semi-airborne transient electromagnetic (SATEM) method, also known as the ground-air transient EM method, uses the grounded electrical source or the loop source as a transmitter, and it uses the unmanned aerial vehicle (UAV) to carry the EM sensors as the receivers. This method not only has the advantages of large exploration depth of ground transient EM systems, but it can also achieve rapid exploration at a lower cost [3,4,5]. The SATEM method is one of the preferred tools for mineral explorations in mountain areas and other complex environments [6,7].

The SATEM method was first proposed by Nabighian [8] in 1988 as a geophysical exploration method. In the 1990s, Elliot developed the FLAIRTEM (Fixed Loop Airborne Transient Electromagnetic) system [9,10]. Elliot used a large loop source as a transmitter and successfully detected a sulfide ore deposit in Papua New Guinea [11]. Smith et al. [12] and Yang and Oldenburg [13] conducted experiments on the ground, including the use of SATEM and AEM systems in a sulfide mining area. Mogi et al. [14] developed the first SATEM system called GREATEM using light helicopters for data acquisition, and they used this system to investigate the geothermal field in the Beppu area of southwestern Japan. Subsequently, Ito et al. [15,16] used this system to conduct deep exploration in the Aso volcanic area of Japan. Ji et al. [17] carried out significant work on the theoretical research, the system development, the data processing, and the interpretations of the SATEM system. Recently, a German consortium, comprising the Federal Institute for Geosciences and Natural Resources (BGR), three universities, two research institutes and two companies, established the Deep Electromagnetic Sounding for Mineral Exploration (DESMEX) project [5,6,7] to develop a SATEM system. Kotowski et al. [3] evaluated the results of the system developed in the DESMEX project.

Compared to the rapid development of the SATEM equipment, the data processing and interpretation are still underdeveloped and mostly use conductivity-depth imaging [18] and 1D inversion [19]. Although these imaging and 1D inversions can deliver a rough image to the underground structures, they cannot as yet achieve accurate results for the SATEM data acquired in areas with complex terrain, where the mineral exploration is often desirable; thus there is a need to develop a more accurate and practical data interpretation technology [5,6,7]. Furthermore, although the SATEM method provides a good solution to the rapid acquisition of data to detect the underground electrical structures in complex terrain areas, it usually only measures the dBz/dt component, which is useful for detecting the horizontal conductive anomalies, but has low resolution for dipping and vertical anomalies, representing typical ore bodies and fracture zones, etc. Taking these problems into account, we present our study about the 3D forward and the inversion method under complex terrain conditions for the multi-component SATEM in this research.

At present, there are many numerical simulation methods that can be applied to the SATEM forward modeling, such as the integral equation (IE) method [20], the finite-difference (FD) method [21], the finite-volume (FV) method [22], and the finite-element (FE) method [23], etc. The IE, FD, and FV methods generally use a structured grid to discretize the model, as they have poor adaptability to complex models. In contrast, the unstructured tetrahedral grids are suitable for modeling complex terrain [24,25,26]. In this paper, we use the FE method for the simulation of the SATEM field under complex terrain. The commonly used 3D inversion methods are nonlinear conjugate gradient (NLCG) method [27]. This method avoids the calculation of Hessian matrix but requires multiple linear searches in each iteration to find the optimal step for model updates. The Gauss-Newton (GN) method uses both the information of the first and second derivatives of the objective function, having approximate quadratic convergence and number of required model updates is often less than NLCG [28,29] but this method requires too much time to compute the Hessian matrix. The Limited Memory quasi-Newton method (L-BFGS) [30] uses the gradient and correction of limited iterations to simplify the approximation to the Hessian matrix, and thus requires less memory and converges quickly. Therefore, in this research, we use the L-BFGS method to minimize the objective function for our inversion.

This paper focuses on the detection capability analysis of the multi-component SATEM using 3D forward modeling and inversions by considering complex terrain. We present our work in the following order: firstly, we introduce the basic theory of forward modelling and the inversion algorithm; subsequently, we show the accuracy and the response analysis of the simple plate model; finally, we test on an inclined ore body model with complex terrain and the Ovoid massive sulfide deposit to show the advantages of our method.

2. Method

2.1. Governing Equations

The Maxwell’s equations in time-domain can be written as [31]

$$\nabla \times \mathbf{e}\left(\mathbf{r},t\right)=-\mu \frac{\partial \mathbf{h}\left(\mathbf{r},t\right)}{\partial t},$$

(1)

$$\nabla \times \mathbf{h}\left(\mathbf{r},t\right)=\sigma \mathbf{e}\left(\mathbf{r},t\right)+\epsilon \frac{\partial \mathbf{e}\left(\mathbf{r},t\right)}{\partial t}+{\mathbf{j}}_s(\mathbf{r},t).$$

(2)

Neglecting the displacement current in Equation (2), we can derive the diffusion equation for the electric field as [28]

$$\nabla \times [\frac{1}{\mu }\nabla \times \mathbf{e}\left(\mathbf{r},t\right)]+\sigma \frac{\partial \mathbf{e}(\mathbf{r},t)}{\partial t}+\frac{\partial {\mathbf{j}}_s(\mathbf{r},t)}{\partial t}=0,$$

(3)

where e(r,t) and h(r,t) are the electric and magnetic field at position r and time t, j_s is the imposed source current, while μ, ε, and σ are respectively the magnetic permeability, dielectric permittivity and conductivity.

2.2. Spatial and Time Discretization

We adapt the unstructured finite-element method [32] for our modeling, and discretize the computational domain into tetrahedral elements. In each element, the electric field can be expressed as,

$$\mathbf{e}(\mathbf{r},t)={\displaystyle \sum _{i=1}^6{e}_i^k(t)}{\mathbf{N}}_i^k(\mathbf{r}),$$

(4)

where, N_i^k(r) and e_i^k(t) are, respectively, the vector basis functions [33] and the electric field along the ith edge of the kth element.

From Equation (3), we can define a residual vector, i.e.,

$$R(\mathbf{r},t)=\nabla \times [\frac{1}{\mu }\nabla \times \mathbf{e}(r,t\left)\right]+\sigma \frac{\partial \mathbf{e}(\mathbf{r},t)}{\partial t}+\frac{\partial {\mathbf{j}}_s(\mathbf{r},t)}{\partial t},$$

(5)

and by assuming that the weighted residual integral of each element is zero, we have

$${\displaystyle \underset{\mathrm{V}}{\iiint }{\mathbf{N}}_{}^e\cdot {\mathbf{R}}^e(\mathbf{r},t)\mathbf{d}\mathrm{v}=0},$$

(6)

where, R^e is the residual vector of element e, v is the volume of element e. Taking Equation (5) into (6), we obtain

$${\mathbf{A}}^e\frac{d{\mathbf{e}}^e\left(t\right)}{dt}+{\mathbf{B}}^e{\mathbf{e}}^e\left(t\right)+{\mathbf{S}}^e=0.$$

(7)

The (i,j) element of matrix A,B and the ith element of matrix S can be written as

$${\mathbf{A}}_{i,j}^e={\displaystyle \underset{\mathrm{V}}{\iiint }{\sigma }^e{\mathbf{N}}_i^e(\mathbf{r})\cdot {\mathbf{N}}_j^e(\mathbf{r})}\mathbf{d}\mathrm{v},$$

(8)

$${\mathbf{B}}_{i,j}^e=\frac{1}{{\mu }^e}{\displaystyle \underset{\mathrm{V}}{\iiint }\nabla \times {\mathbf{N}}_i^e(\mathbf{r})\cdot \nabla \times {\mathbf{N}}_j^e(\mathbf{r})}\mathbf{d}\mathrm{v},$$

(9)

$${\mathbf{S}}_i^e={\displaystyle \underset{\mathrm{V}}{\iiint }{\mathbf{N}}_i^e(\mathbf{r})\cdot \frac{\partial {\mathbf{j}}_s(\mathbf{r},t)}{\partial t}}\mathbf{d}\mathrm{v}.$$

(10)

For time discretization, we adopt unconditionally the stable backward Euler method [29,34,35]. The second order form can be obtained approximately by a Taylor expansion, i.e.,

$$\frac{d{\mathbf{e}}^{i+2}}{dt}=\frac{1}{2\Delta t}\left(3{\mathbf{e}}^{i+2}-4{\mathbf{e}}^{i+1}+{\mathbf{e}}^i\right).$$

(11)

Taking Equation (11) into (7), we obtain

$$(3\mathbf{A}+2\Delta t\mathbf{B})e^{i+2}=\mathbf{A}(4{\mathbf{e}}^{i+1}-{\mathbf{e}}^i)-2\Delta t{\mathbf{S}}^{i+2}$$

(12)

2.3. Boundary Conditions and Source Term

To ensure that there is a unique solution to the EM field, we assume that the outer boundary of the simulation area is far enough away from the transmitting source, so the electric field satisfies the Dirichlet boundary condition [35], i.e.,

$${(\mathbf{n}\times \mathbf{e})}_{\mathsf{\Gamma }}=\mathbf{0},$$

(13)

where, n is the normal vector of the outside boundary Γ of the calculation domain.

We consider that the SATEM method is generally based on electrical transmitting sources. Therefore, when the transmitting waveform is a step-up or pulse wave, the initial electric field is zero and when the transmitting waveform is a step-off wave, the initial electric field is the solution of a DC problem. The potential generated by the point source satisfies Poisson’s equation. The potential at the boundary meets the Dirichlet boundary condition, so the boundary value problem can be expressed as

$$\nabla \cdot (\sigma \nabla \phi (\mathbf{r}\left)\right)=-\nabla \cdot j_s\left(\mathbf{r}\right),$$

(14)

$$\phi {|}_{\mathsf{\Gamma }}=\mathbf{0},$$

(15)

where φ(r) is the potential at position r. We adopt the unstructured node-based FE method to solve the problem [34].

2.4. Forward Modeling Equation

Assembling Equation (12) for all time channels [28], we obtain the following equation:

$$\left(\begin{array}{cccccc}{\mathbf{C}}_1& & & & & \\ -4\mathbf{A}& {\mathbf{C}}_2& & & & \\ \hspace{0.17em}\mathbf{A}& -4\mathbf{A}& {\mathbf{C}}_3& & & \\ & \hspace{0.17em}\mathbf{A}& -4\mathbf{A}& {\mathbf{C}}_4& & \\ & & \ddots & \ddots & \ddots & \\ & & & \mathbf{A}& -4\mathbf{A}& {\mathbf{C}}_{\mathbf{n}}\end{array}\right)\left(\begin{array}{c}{\mathbf{e}}_1\\ {\mathbf{e}}_2\\ {\mathbf{e}}_3\\ {\mathbf{e}}_4\\ ⋮\\ {\mathbf{e}}_{\mathbf{n}}\end{array}\right)=\left(\begin{array}{c}{\mathbf{S}}_1+3\mathbf{A}{\mathbf{e}}_0\\ {\mathbf{S}}_2-\mathbf{A}{\mathbf{e}}_0\\ {\mathbf{S}}_3\\ {\mathbf{S}}_4\\ ⋮\\ {\mathbf{S}}_{\mathbf{n}}\end{array}\right),$$

(16)

where C_n = 3A+2Δt_nB, s_n = 2Δt_nS_n. Equation (16) can be rewritten as

$$\mathbf{K}\mathbf{e}={\mathbf{S}}_{\mathsf{s}},$$

(17)

where K is the coefficient matrix, e is the electric field, and s_s is the source term. We adopt the direct solver MUMPS [36] to factorize Equation (12) to obtain the electric field at the edges of each tetrahedron element and obtain dBx,y,z/dt at the receiver location via Faraday’s law.

2.5. Regularized Inversion

For the geophysical 3D inversion, the number of model parameters is usually much larger than that of survey data. Mathematically, such problems are underdetermined problems. To overcome the non-unique solutions brought by the underdetermined problems, Tikhonov [37] proposed a regularization method, which reduces the multiple solutions by reducing the size of the solution space. Based on this, we propose the following objective function for 3D time-domain EM inversions, i.e.,

$$\phi (\mathbf{m},\mathbf{d})={\Vert {\mathbf{C}}_{\mathbf{d}}^{}(\mathbf{d}-\mathbf{f}(\mathbf{m}\left)\right)\Vert }_2^2+\lambda {\Vert {\mathbf{C}}_m^{}(m-m_0)\Vert }_2^2,$$

(18)

where, m denotes the parameter vector of the conductivity model, d is the N_d dimensional data vector, C_d is a diagonal matrix whose elements are composed of the covariance of the data error, while f(m) is the forward response calculated by the vector FE method. m₀ is the initial model, λ is trade-off parameter, and C_m is the model roughness operator. By differentiating both sides of Equation (18) with respect to the model parameter m, we obtain the gradient of the objective function [28] as

$$\mathbf{g}=-2{\mathbf{J}}^{\mathrm{T}}{\mathbf{C}}_d^{\mathrm{T}}{\mathbf{C}}_d^{}\mathbf{r}+2\lambda {\mathbf{C}}_m^{\mathrm{T}}{\mathbf{C}}_m^{}(m-m_0),$$

(19)

where r = d $-$ f(m_n), J represents the sensitivity matrix. In this paper, the BFGS method is used to search the step length to obtain the model update Δm, subsequently the new model is obtained by m_n+1 = m_n + Δm.

2.6. Sensitivity Matrix

To obtain the gradient of the objective function, we need to calculate the sensitivity matrix first. We use the adjoint forward method [28] to compute the sensitivity matrix implicitly. This method avoids directly storing the sensitivity matrix and thus can largely reduce the memory consumption. The sensitivity matrix J can be expressed as

$$\mathbf{J}=\frac{\partial d}{\partial m}=L\frac{\partial e}{\partial m},$$

(20)

where, L = L_tL_d. L_d is the interpolation operator that maps the electric field to the data locations, whereas L_t is used to interpolate data from the computational time channels to observed ones.

Differentiating both sides of Equation (17) with respect to model parameters, we obtain

$$\frac{\partial e}{\partial m}=K^{-1}\left(\frac{\partial s_s}{\partial m}-\frac{\partial K}{\partial m}e\right),$$

(21)

then J^Tr_d can be expressed as

$$J^{\mathrm{T}}{r}_d={\left(\frac{\partial s_s}{\partial m}-\frac{\partial K}{\partial m}e\right)}^{\mathrm{T}}{K}^{-\mathrm{T}}{L}^{\mathrm{T}}{r}_d.$$

(22)

Assuming that

$$K^{-\mathrm{T}}{L}^{\mathrm{T}}{r}_d=u$$

(23)

we have

$$K^{\mathrm{T}}u=L^{\mathrm{T}}{r}_d,$$

(24)

since A and C in Equation (16) are symmetric, K^T can be written as

$$K^{\mathrm{T}}=\left(\begin{array}{cccccc}{\mathbf{C}}_1& -4\mathbf{A}& \mathbf{A}& & & \\ & {\mathbf{C}}_2& -4\mathbf{A}& \mathbf{A}& & \\ \hspace{0.17em}& & {\mathbf{C}}_3& -4\mathbf{A}& \ddots & \\ & \hspace{0.17em}& & {\mathbf{C}}_4& \ddots & \mathbf{A}\\ & & & & \ddots & -4\mathbf{A}\\ & & & & & {\mathbf{C}}_{\mathbf{n}}\end{array}\right).$$

(25)

From Equation (25), a reverse recursive process is needed to solve the adjoint forward modeling of u. After the adjoint forward modeling, we substitute u into Equation (22) and obtain

$$J^{\mathrm{T}}{r}_d={\left(\frac{\partial s_s}{\partial m}-\frac{\partial K}{\partial m}e\right)}^{\mathrm{T}}u.$$

(26)

To obtain the product of the transposed sensitivity matrix and the vector, we still needed to calculate (∂K/∂m)e and ∂S_S/∂m. By referring to Equation (21), (∂K/∂m)e can be obtained by a simple calculation. For ∂S_S/∂m, since we use the total field method for forward modeling, s_i in the source term is independent of the conductivity, so we have

$$\frac{\partial s_{\mathrm{i}}}{\partial m}=0,\begin{array}{cc}& \mathrm{i}=1,2,3,\cdots ,n\end{array}.$$

(27)

When the direct method is used to solve the forward problem, we can save the decomposition results of matrix K, which can largely reduce the calculations by directly calling in the accompanying forward process. It must be pointed out that when the step-off wave is used for calculation, the initial electric field contains underground conductivity information, so another item ∂Ae₀/∂m in Equations (16) and (17) needs special treatment. For this purpose, we write ∂Ae₀/∂m as

$$\frac{\partial A{e}_0}{\partial m}=\frac{\partial A}{\partial m}{e}_0+A\frac{\partial e_0}{\partial m},$$

(28)

where the term (∂A/∂m)e₀ can be obtained by simple calculation, but A(∂e₀/∂m) needs to be calculated by adjoint forward modeling in DC field simulation [28].

2.7. Model Roughness Constraint

We use the following model roughness for our 3D time-domain inversion [28]

$${\Vert {\mathbf{C}}_m^{}(m^{}-m_0)\Vert }_2^2={\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{M}}{V}_{\mathrm{i}}\left[{{\displaystyle \sum _{\mathrm{j}=1}^{{\mathrm{N}}_{\mathrm{i}}}{\omega }_{\mathrm{j}}\left(\frac{m_{\mathrm{i}}-m_{\mathrm{j}}}{\Delta r_{\mathrm{ij}}}\right)}}^2\right]},$$

(29)

where,

$$\Delta {\mathbf{r}}_{ij}=\sqrt{{\left(x_i-x_j\right)}^2+{\left(y_i-y_j\right)}^2+{\left(z_i-z_j\right)}^2}$$

(30)

$${\omega }_j=\frac{{\mathbf{V}}_j}{{\displaystyle \sum _{k=1}^{N\left(i\right)}{\mathbf{V}}_k}}$$

(31)

while, V_i is the volume of the ith element, N_i is the number of elements that share the nodes with the ith element, Δr_ij is the distance between the centroids of the current element and the jth adjacent element.

2.8. Termination Condition for the Inversion

In the inversion process, the root mean square (RMS) is used to judge whether to terminate an inversion that is defined by

$$\mathrm{RMS}=\sqrt{\frac{{\displaystyle \sum _{i=1}^{N_d}\left[\left(d_{\mathrm{obs}}\left(i\right)-d_{\mathrm{pre}}{\left(i\right)}^2\right)/e_i^2\right]}}{N_d}}$$

(32)

where e_i is the noise standard deviation of ith data. When the RMS drops to 1, we believe that the data has been well fitted and the inversion is terminated. However, when the inversion falls into a local minimum, the RMS cannot drop to 1. At that time, we use the gradient value to judge whether to terminate the inversion. The judgment conditions are defined as

$${\Vert g\Vert }^2<\epsilon$$

(33)

where, ε is a small real number.

3. Numerical Experiments

In this section, we performed numerical experiments of forward modeling and inversions to analyze the ability of our newly proposed method. All the tests were performed on a workstation equipped with an Inter (R) Xeon (R) 3.1GHz CPU and 128 GB of memory.

3.1. Accuracy Verification

To verify the accuracy of our 3D time-domain FE method, we assumed a uniform half-space and a layered model to compare our 3D time-domain results with 1D semi-analytical solutions. We set up a Cartesian coordinate system with the horizontal surface at z = 0 m and z axis positive downwards. The uniform half-space had a resistivity of 100 Ω∙m and the H-type 1D model had three layers with the resistivity of 100 Ω∙m, 10 Ω∙m and 100 Ω∙m, respectively, and with the thickness of the first two layers being 50 m. An electrical transmitting source of 1000 m long was used as a source with the center located at (0 m, −500 m, 0 m). The transmitting current was 1 A and the waveform was a step-off wave. The receiver was in the air at (0 m, 0 m, −30 m). Figure 1 shows the comparison of the results. Since the transmitting source is parallel to the x-axis, the response in the x-direction of the receiving point is zero, consequently we only compared the EM field components in the other two directions. It can be seen that the 3D forward modeling results in this paper are in good agreement with 1D semi-analytical solutions, with the maximum relative error being less than 5%.

3.2. Detectability Analysis of Multi-Component Responses for Plate Models

To analyze the detectability of each EM component for the anomalous body, we designed two models, one with a horizontal plate and the other with a vertical plate with resistivity of 1 Ω∙m, buried in a uniform half-space of 100 Ω∙m (Figure 2). The size of plates is 200 m × 200 m × 50 m and 50 m × 200 m × 200 m in horizontal and vertical plate models, respectively, and the depth of the top of the plate is 50 m in both the models. An electrical source parallel to the x-axis is used as the transmitter, with the center located at (0 m, −500 m, 0 m). The transmitting current is step-off wave of 50 A. The flight altitude is 30 m with offset between the transmitter and the receivers as 500 m. The survey line is 400 m long along the x-axis and the interval between surveying points is 25 m.

Figure 3 shows the multi-component responses for the two models in Figure 2. It can be seen that for the horizontal plate model, the EM responses in the x direction and z direction can both reflect the existence of the anomalous body, but dBx/dt can better reflect the locations of the vertical boundaries. For the vertical plate model, it is not possible to distinguish the anomalous body only by the response in the z direction, while dBx/dt can still clearly reflect its existence. Therefore, in practical applications, it is necessary to measure multi-component EM fields to better depict the anomalous bodies.

To verify the effectiveness of our inversion algorithm, we inverted forward modelled data from the plate models with 3% Gaussian noise to single component and multi-component data, respectively. Nine survey lines parallel to the x-axis were planned between −200 m and 200 m, with 50 m line spacing along the y-axis. There were 41 measurement points on each line, with 30 m flight altitude. The same transmitting source as shown in Figure 2 was used to generate the synthetic data. From Figure 4A,B, we can see that the single and the multi-component inversions can both recover the position and the resistivity of the anomalous body, but the multi-component inversion has more capability to depict the vertical boundaries of the plate (Figure 4B,D).

Figure 4C,D show the inversion results for the vertical plate model. It can be seen that for the vertical plate model, since the response of its vertical component can barely see the existence of the anomalous body (see Figure 3F), the inversion results from the vertical component alone cannot recover the anomalous body. In contrast, the multi-component inversion can better recover the location and resistivity of the anomalous vertical body though it fails to recover the correct depth of the bottom of the vertical body. This confirms that the horizontal components are very important in the detection of inclined structures.

3.3. Inclined Ore Body under Complex Terrain

Considering that the mineralization is often found in mountainous areas with complex terrain, the data interpretation and the inversion must consider the effect of the terrain. In this section, we designed an inclined ore body below an undulating terrain (Figure 5) to verify the practicability of our algorithm. The size of the ore body was 750 m × 300 m × 300 m and the depth of the center point of the ore body was 200 m below the horizontal surface. The resistivity of the ore body and background was 1 and 100 Ω∙m, respectively. For this inversion test, we designed 12 survey lines between 700 m and 1800 m in the x direction with a line spacing of 100 m. There were 61 survey points on each line at an interval of 25 m between 0 and 1500 m in the y direction with a flight altitude of 30 m. The transmitting source was 1 km long aligned in the y-axis direction with a 500 m offset to the nearest measurement point. Figure 5A,B show the locations of the source, survey points and the inclined ore body. Figure 5C,D show the grids used for forward modeling and inversions, respectively. The forward model was discretized into 642338 tetrahedral elements, and the inverse model was discretized into 608834 elements. All the synthetic data was contaminated by 3% Gaussian noise.

To analyze the characteristics of the response of the inclined ore body model, we displayed the response along the central surveying line passing over the ore body in Figure 6. Due to the existence of complex terrain, the response curves had obvious fluctuations in the early time channels. Since the transmitting source in this example was parallel to the y axis (different from that in Section 3.2), the component of dBy/dt was obviously affected by the anomaly body (see Figure 6B), so that it had the potential to improve the inversion result. We note that even under the influence of the terrain, the dBz/dt component was still sensitive to the steep anomaly body (Figure 6C), thus it could recover the inclined ore body to some extent.

Figure 7 shows the inversion results for the inclined ore body model. It can be seen that the inversions of the vertical component alone could not well depict the horizontal location and extent of the anomalous body (Figure 7A), while the multi-component inversions could better recover the location, dipping nature and the resistivity of the target. This example shows the practicability of the multi-component inversions together with complex terrain.

Figure 8 shows the parameters of the SATEM inversion for an inclined ore body through various iterations. It can be seen that the RMS and objective function values for the two inversions decreased rapidly at early stage. When the RMS approached 1, the model variations tended to be stable, which proved that the inversions worked well. Meanwhile, compared with the single component inversion, the multi-component inversions proposed in this paper needed fewer iterations to converge. Figure 9 shows the comparison between the noise-added synthetic data and the predicted data for the multi-component inversion at t = 0.92 ms for an inclined ore body. The distribution of each component of the predicted data was very similar to the observed ones, and the noise in the synthetic data was smoothed out. The RMS for most of the survey points was in the range between 0 and 2 (other than a few points which were greater than 2), which means that most of the data have been fitted well.

3.4. Application to the Ovoid Massive Sulfide Deposit

We presented a realistic model (Figure 10) as the second example. This model is an Ovoid massive sulfide deposit located at Voisey’s Bay, Labrador, Canada. The ore body is flat-lying and composed of 70 % massive sulfide. It is located under approximately 20 m of overburden [38,39,40,41] and the resistivity of the body and background is 0.01 and 100 Ω∙m, respectively. For this inversion test, we designed 10 survey lines between −750 to 150 m in the y direction with a line spacing of 100 m. There were 41 survey points on each line between −500 m and 500 m in the x direction with the interval of 25 m, and the flight altitude was 30 m. The transmitting source was 1 km long aligned in the x-axis direction, and the offset to the nearest survey point was 500 m. Figure 10A shows the source and the receivers locations with the ore body in the xy plane. Figure 10B,C show the ore body and the model discretization in the two clips from the forward model at x = −200 m and y = 0 m, respectively. The forward model was discretized into 814573 tetrahedral elements, and the inverse model was discretized into 537721 elements. All the synthetic data was contaminated by 3% Gaussian noise.

Figure 11 shows the EM responses of the sulfide deposit model at x = 50 m. Since the model was mainly horizontal, the dBz/dt component had the largest anomaly response in the early time channels (Figure 11C). However, we noticed that dBy/dt also had an obvious anomaly response that could provide information additional to the dBz/dt in the inversion.

Figure 12 shows the results of the single- and multi-component inversion for the sulfide ore body model. We can see that the results of the multi-component inversion method effectively reduced the false anomalies in the inversion, and depicted the anomaly boundary more finely. This test demonstrates again that the multi-component data has obvious advantages on the improvement of resolution.

Figure 13 shows the parameters versus the iterations during the inversion progress. We can draw the same conclusion that for both inversions, the RMS, the objective function and the norm of model variation were all changing smoothly during the iteration and converged at the end of inversion, indicating the stability of the inversion algorithm. Meanwhile, compared with the single component algorithm, the multi-component inversions could converge faster because it contained more useful information from dBx/dt and dBy/dt components. Figure 14 shows the comparison between the noise-added synthetic data and the predicted data for the multi-component inversion at t = 0.92 ms. The distribution of each component of the predicted data agreed well with the observed ones, and the RMS at each survey point was distributed in the range of [0, 5] for this time channel, showing that most of the data have been fitted.

4. Conclusions

In this paper, we proposed a 3D forward and inverse algorithm for a multi-component SATEM and successfully tested it on synthetic mineral exploration models with complex terrain. The forward modeling based on the vector FE method for unstructured grids could accurately simulate arbitrary complex terrain and geological structures. After verifying the accuracy of the forward algorithm, we analyzed the detectability of the multi-component SATEM using plate models. We found that the horizontal components were more sensitive to the vertical boundaries of the anomalous body than the vertical component alone. Since there is no available SATEM field data at present, we carried out synthetic data inversion for several typical mineral models. The synthetic data for inversion was obtained by adding 3 % Gaussian noise to the forward modelled data. For all inversion examples, the initial starting model was a uniform half-space model with a resistivity of 100 Ω∙m which was chosen by trial-and-error for fast convergence. The trade-off parameter and its decay in each iteration needed to be selected appropriately by assessing various values. We observed that the inversion converges well and quickly when appropriate inversion parameters are used. The numerical experiments on two mineral exploration models showed that the inversion of the noise-added synthetic data from significant terrain recovers the anomalous bodies well. Therefore, our inversion algorithm can be used for field data interpretation. Comparing the results of multi-component inversions with those with a single component, we confirmed that the multi-component data inversion has obvious advantages to provide a more accurate and real model. The multi-component method proposed in this paper has the potential for detecting fine structures for mineral explorations in complex terrain areas. The measurement of all the components of the magnetic field decay instead of only the vertical component should be emphasized in the future development of SATEM receivers. The 3D inversion of multi-component SATEM field data will be our next target once such data are available.