Next Article in Journal
Improving the Spatial Prediction of Soil Organic Carbon Content Using Phenological Factors: A Case Study in the Middle and Upper Reaches of Heihe River Basin, China
Next Article in Special Issue
3D Airborne EM Forward Modeling Based on Finite-Element Method with Goal-Oriented Adaptive Octree Mesh
Previous Article in Journal
Validation of Swarm Langmuir Probes by Incoherent Scatter Radars at High Latitudes
Previous Article in Special Issue
Metallogenic Prediction of Magnetite in the Pandian Area at the Northwest Margin of Luxi Uplift, China: Constraints of Wide-Field Electromagnetic Data
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Three-Dimensional Dual-Mesh Inversions for Sparse Surface-to-Borehole TEM Data

College of Geo-Exploration Sciences and Technology, Jilin University, Changchun 130026, China
*
Author to whom correspondence should be addressed.
Remote Sens. 2023, 15(7), 1845; https://doi.org/10.3390/rs15071845
Submission received: 15 January 2023 / Revised: 18 March 2023 / Accepted: 27 March 2023 / Published: 30 March 2023
(This article belongs to the Special Issue Multi-Scale Remote Sensed Imagery for Mineral Exploration)

Abstract

:
The surface-to-borehole transient electromagnetic (SBTEM) method can provide images at higher resolution for deep earth because its receivers are close to targets. However, as usually the boreholes distribute sparsely, the limited EM data can result in an “equivalent trap” in SBTEM inversions, i.e., the data are well-fitted, but the model is not properly recovered. To overcome this non-unique problem, we propose a dual-mesh three-dimensional (3D) SBTEM inversion scheme. We first use a coarse mesh to obtain a rough resistivity distribution near the borehole, and then we map the coarse mesh attribute into a fine one and capture details from the fine mesh inversion. We test our method on both synthetic data and survey data acquired in Daye, Hubei Province, China. Numerical experiments show that our dual-mesh inversion strategy can better recover the location and resistivity of targets.

Graphical Abstract

1. Introduction

The surface-to-borehole transient electromagnetic (SBTEM) method uses a grounded wire as transmitter, while the receivers are laid in the borehole to measure EM fields [1]. Since the borehole receivers are close to the targets, the anomaly signal is much stronger than that in the conventional ground EM, so the deep targets concealed around the borehole or at the shaft bottom can be detected [2]. In past decades, the SBTEM method has been widely used in mineral explorations. It has played an important role in the detection of massive sulfide, pyrite-pyrrhotite ore, and gold deposits and in the exploration of deep targets like volcanic veins [2,3,4,5,6,7].
Although the SBTEM equipment has now been well-developed, the interpretation of SBTEM data is still underdeveloped, usually assuming a simple half-space or layered-earth model. The mainstream imaging methods simply transform TEM data into some intermediate parameters to predict the rough structures in the underground [8], which cannot guarantee that the recovered model can fit the observed data. Zhang and Xiao [9] proposed a method that runs joint inversions of surface and borehole data to one-dimensional (1D) models and found that they can well-detect both the shallow and deep structures. Chen et al. [10,11] successfully conducted the joint inversions and 1D Occam’s inversions of SBTEM data. However, all these imaging and 1D inversions have difficulty in recovering the complex underground structures, so the development of three-dimensional (3D) inversions for SBTEM data is necessary.
Since, until now, little has been reported on 3D SBTEM inversions, we review the mainstream 3D forward modeling and inversions of conventional TEM methods. At present, the mainstream numerical simulations for TEM methods include the integral equation (IE) [12,13] the finite-difference (FD) [14], and the finite-element (FE) methods [15,16]. As for the TEM inversions, most conventional 3D algorithms are based on structured grids [17,18]. For accurate modeling and inversions, one needs very small grids to fit the undulating topography or complex underground structures. This can result in heavy computational cost. To solve this problem, we take in this study the locally refined tetrahedral grids to discretize the computational domain and use the vector FE method to calculate EM responses and sensitivity information. For the time discretization, we choose the unconditionally stable backward Euler method [16,19,20,21]. Furthermore, we apply the L-BFGS (limited-memory BFGS) method for model updates at two inversion stages. In the SBTEM method, we usually have only sparsely distributed borehole data. To obtain the fine structures close to the borehole, however, the grids around the borehole receivers are generally refined; thus, we actually deal with an inverse problem with sparse data for the fine model. In this case, the borehole data can be easily fitted by a combination of resistivities of fine grids close to the receivers as they usually have large sensitivities. This can bring serious non-uniqueness to 3D EM inversions. To deal with this problem, we propose an inversion strategy for SBTEM that adds the surface EM data to the inversion and applies a dual-mesh scheme. Different from other multi-mesh inversion methods [22,23], we use a coarse mesh to do the first round of inversion to invert the distribution of the background structures, and then, we keep the model and switch to a fine mesh to recover model details and to fit the data.
In the sequence, we first introduce the SBTEM inversion algorithm and then present in detail the dual-mesh inversion strategy. We test our 3D SBTEM inversion scheme via both synthetic and field survey data.

2. Methods

As shown in Figure 1, in the SBTEM method, the transmitter placed at the earth surface transmits the EM waves into the earth and creates the eddy currents in the underground. The eddy currents will be distorted when they travel to the targets with anomalous conductivities, so that anomalies will be observed by the sensors placed at the borehole. From these anomalies, we can infer the anomalous bodies. Since the receivers are generally laid close to the anomalous bodies, one can observe large anomalous responses and, thus, can easily detect the targets.

2.1. Regularized Inversion

Since 3D SBTEM inversion is nonlinear and underdetermined, we use regularization to solve the problem and construct the objective function based on an L2 norm measure, i.e.,:
φ ( m , d pre ) = W d ( d obs d pre ) 2 2 + λ W m ( m m 0 ) 2 2 ,
where dobs denotes the observed data of the Nd × nchannel dimension, while dpre denotes the predicted data from the inversion model. Nd and nchannel denote the number of receivers and time channels. Here, we use the forward modeling code from Wang et al. [24] to calculate EM responses for the vector dpre. To better introduce the inversion scheme, in Appendix A, we give a simple introduction to the forward modeling algorithm based on the FE method using tetrahedral grids. λ is a regularization factor, Wd is the data variance matrix with its elements being the reciprocal of the standard noise deviation, Wm is the model roughness operator, while m0 and m are the M-dimensional vectors, respectively, for the reference model and the solution.
To obtain the solution mn of nth, we differentiate both sides of Equation (1) with respect to the model parameters mn−1 and obtain the gradient of the objective function, i.e.,:
g n = 2 J T W d T W d r n 1 + 2 λ W m T W m ( m n 1 m 0 ) ,
where J denotes the sensitivity matrix, and r n 1 = d obs d pre n 1 .
The sensitivity matrix in above equation is calculated by the adjoint forward modeling for fast calculation and reduction in memory consumption [25,26]. We assume that the predicted data satisfy d = Le, where L denotes a spatial interpolation operator. For the forward modeling in time-domain, L can be written as L = LtLd, where Ld is the interpolation operator that interpolates the solution of the electric field to the receiver locations, while Lt is the interpolation operator that interpolates the solution of the calculated time channels to the observation ones. From the forward modeling scheme and d = Le, the sensitivity of EM responses with respect to the conductivity of the kth tetrahedral element, namely the element of the kth column in the overall sensitivity matrix, can be written as:
j k = d m k = L e m k = L K 1 ( s TEM m k K m k e ) ,   k = 1 , 2 , , M .
Define a matrix G of the dimension N × M (N = Nedge × nchannel) as:
G = { s TEM m 1 K m 1 e ,   s TEM m 2 K m 2 e ,   , s TEM m M K m M e   } ,
and substitute Equation (4) into (3); we obtain:
J T = G T ( K 1 ) T L T .
Equation (2) shows that to obtain the gradient of the objective function, we need to calculate the product of the transpose of the sensitivity matrix and a vector. Assuming that the multiplication vector is v, then multiplying both sides of Equation (5) with v and defining u = ( K 1 ) T L T v , we obtain the following equations for the adjoint forward modeling, i.e.,:
K T u = L T v ,
where the vector v is the term after the sensitivity matrix J in Equation (2). After solving the adjoint forward equation, we obtain the vector u. Then, the product of the transpose of the sensitivity matrix J and the vector v can be calculated by:
J T v = G T u .
When calculating the matrix G in Equation (4), the term of K m k e is analytical. As for the term of s TEM m k , the derivative of the source s i m k is zero; the rest term can be written as:
A e 0 m k = A e 0 m k + A m k e 0 ,
where A m k e 0 can also be analytically calculated, while the term A e 0 m k can be calculated by adjoint method when executing the forward modeling for the DC field [25].
Since we use the unstructured tetrahedral mesh to discretize the calculation domain, the conventional model roughness operators for structured girds are not suitable. Thus, we modify the model roughness operator proposed by Key [27] for 2D problems to our 3D case [28] and obtain:
W m ( m m 0 ) 2 = i = 1 M V i [ j = 1 N ( i ) ω j ( Δ m i j Δ r i j ) 2 ] ,
where
Δ m i j = m i m j ,
Δ r i j = ( x i x j ) 2 + ( y i y j ) 2 + ( z i z j ) 2 ,
ω j = V j k = 1 N ( i ) V k
In the above equations, Vi denotes the volume of the ith tetrahedron, while N(i) denotes the number of the tetrahedrons that share vertices with the ith tetrahedron. The distance Δrij between tetrahedrons is calculated from the centroid of each element.

2.2. Dual-Mesh Inversion Strategy

In the practical SBTEM survey, the boreholes are sparsely distributed in the working area. In most cases, only one borehole can be used. On the other hand, in 3D EM inversions, the areas close to the borehole receivers are generally refined to achieve accurate results. The sparse data combined with fine meshes can easily result in underdetermined 3D inversions. In this situation, an “equivalent trap” can frequently occur, namely the borehole data are well-fitted, but the model is not properly recovered. Under the condition of limited data, we can alternatively reduce the number of unknowns in the inversion by optimizing the mesh. In this paper, we propose a dual-mesh inversion strategy. First, we use a coarse mesh in the inversion to obtain rough structures in the underground. After that, we map the attributes in the coarse mesh to a fine one and continue the inversion. In this way, the inversion results from the coarse mesh can be taken as a good initial model, so that we can finally recover the underground structures in more detail. This can also help improve the data fitting in the inversions.
Referring to Figure 2, to do the mapping between the coarse and fine meshes, we take four basic situations into account: (1) the centroid of a fine mesh is inside a coarse one; (2) the centroid of a fine mesh is on the surface shared by two coarse meshes; (3) the centroid of a fine mesh is on the edge shared by several coarse meshes; and (4) the centroid of a fine mesh is at the vertex shared by several coarse meshes. For the first case, the attribute of the fine mesh is set to be equal to that of the coarse mesh that contains its centroid, while for all other cases we choose the coarse meshes related to the centroid of the fine mesh for calculating the attribute via volume-weighted average. Then, the attribute of the fine grid is calculated by:
( σ i e ) fine = j = 1 m ( σ i j e ) coarse j = 1 m ( σ i j e V i j e ) coarse ,   i = 1 , 2 , , n ,
where σ i j e denotes the attribute of ith mesh. m is the number of coarse meshes related to the centroid of the fine mesh, n is the total number of fine meshes, and σ i j e and V i j e are the conductivity and volume, respectively, of the jth coarse mesh related to the centroid of ith fine mesh.

3. Numerical Experiments

In this section, we will test our 3D SBTEM inversion scheme via both synthetic and field survey data. To remain consistent with the survey data, for all models in this section we adopt long wires as transmitting sources and assume the range of the calculation to be between 2.5 × 10−4 s and 1 × 10−2 s.

3.1. Flat Model Inversion

To validate our inversion algorithm, we first design a simple model with a block body buried in a uniform half-space (Figure 3). The two transmitting sources with the same length of 700 m are laid on the surface. Tx 1 is the transmitter corresponding to the ground receivers, while Tx 2 is the transmitter for the receivers in the borehole. The borehole receivers are located from 0 m to 600 m with an interval of 10 m. The borehole in the model is an open-hole well, which is hollow and full of air. The Ez data in the borehole are measured by letting the electrodes touch the rock formation. To overcome the influence of limited data on our 3D inversion, we also put 31 × 31 = 961 survey stations at the earth surface with an interval of 20 m. The abnormal body’s resistivity is 1 Ω·m, while the surrounding rock and the air is 100 Ω·m and 1.0 × 108 Ω·m, respectively. The forward model is discretized into 663,443 grids.
We perform 3D forward modeling to calculate Ex on the surface and Ez in the borehole and then added 3% Gaussian random noise to simulate the real situation. In the inversions, the model is first discretized into 562,888 grids as the coarse mesh, without refining elements close to the borehole and receivers. We use this mesh to do the inversion to the surface data, the borehole data, and the joint surface and borehole data, respectively. Figure 4 shows that the RMS for the surface data inversion is reduced to 0.998, but the other two inversions have difficulty achieving good data fittings. Next, we refine the elements close to the receivers and the borehole and discretize the model into 659,955 grids, and then we map the preliminary inversion results of the borehole data from the coarse mesh to the fine one. This operation will cause an RMS rising, but it does not obviously affect the final convergence. The final RMSs for the inversions of both the borehole data and the joint surface and borehole data reduce to the target level of 1.
Figure 5 shows the inversion results using different strategies for different data. Figure 5a shows that the inversion only with the surface data roughly recover the model, but it has difficulty in recovering the boundaries of abnormal body, the recovered resistivity of the abnormal body is much larger than the true value. When we only invert the borehole data with fine mesh, since the sizes of the meshes near the receivers are very small, it is difficult to observe small meshes with drastic changes in attributes in plots of relatively large inversion areas. Figure 5b shows that the anomalous body is hardly to be observed as the inversion is too ill-posed. In this case, the borehole data can be easily fitted by a combination of resistivities of fine grids close to the receivers as they usually have large sensitivities. In contrast, when we use dual mesh for borehole data inversion, the anomalous body close to the borehole is well-inverted, especially its resistivity is well-recovered (Figure 5c). This succeeds because the coarse grids near the borehole cannot fit the data well, and thus, the inversion tries to find a model to fit the data by seeking the solution in a larger space of fine mesh. Figure 5d shows the joint inversion of surface and borehole data using fine mesh. The abnormal body is roughly recovered; however, the inverted body has a long tail in all sections. Finally, when we use the dual-mesh strategy to invert the joint surface and borehole data, from Figure 5e we can see that the best results are obtained; the boundaries and the resistivity of the abnormal body are well-recovered.
Summarizing the above discussion, we conclude that, by adding the surface data to SBTEM inversions, we can improve the resolution to the target around the borehole. Using a dual-mesh inversion strategy, we can improve the resolvability not only to the location but also to the resistivity of the target body.

3.2. Topographic Model Inversion

In this section, we further designed a synthetic model with real topography from a field dataset that will be illustrated in the next section. The locations of two grounded wire transmitters and data collection points at the surface are both taken from the field dataset. Referring to Figure 6, transmitting source 1 has a length of 2964 m, while transmitting source 2 has a length of 2791 m. The borehole is located at (−600 m, 0 m) with a depth of 600 m. The borehole receivers are laid at depths from 20 m to 600 m at an interval of 10 m. To remedy the limited data amount, we also take data acquired at the earth surface for our 3D inversion.
Referring to Figure 7, three conductive abnormal bodies with the resistivity of 1 Ω·m are embedded in the underground. The surrounding rock and the air have resistivities of 100 Ω·m and 1.0 × 108 Ω·m, respectively. The model is first discretized into 472,151 coarse grids and, then, into 664,173 fine grids. We perform 3D forward modeling on this synthetic model and compare the responses with those for a uniform half-space. From Figure 8, one sees that the conductive bodies cause obvious anomalies. Considering that the electric fields in the x- and y-directions are not easy to measure in the borehole, we use only the Ez component as the data and add 10% Gaussian random noise to simulate the real situation. In our inversion process, we use a cooling schedule to determine regularization factor λ, namely when the data fitting decreases very slowly, the regularization factor becomes 0.1 times the previous value and remains unchanged after reaching a threshold.
To undertake our dual-mesh inversions, we first use the coarse mesh in Figure 7a to recover the preliminary model. After 25 iterations, we interpolate the inversion results to the fine grids in Figure 7b using the mapping scheme in Figure 2 and continue the inversion with the fine grids until the stop condition is satisfied (i.e., RMS ≤ 0.98). The inversion runs 6 h and 25 min on the coarse mesh with the RMS reduced to 5.05. After mapping to the fine mesh, the inversion runs again 16 h and 28 min on the fine mesh, the RMS is finally reduced to 0.98. In addition, we also run 3D inversions using the conventional single-mesh scheme with the fine mesh. The inversion takes 26 h and 37 min, with the RMS finally reduced to 0.97. Figure 9 compares the inversion parameters between our dual-mesh scheme and the conventional one. For both methods, the RMS and the objective function decrease rapidly at the early time and then get stabilized. Although the parameters for our dual-mesh inversions jump when mapping between two meshes, the overall RMS and objective function attenuate, and the inversion converges.
Figure 10 shows the comparison of inversion results from the two methods. The conventional inversion has difficulty in recovering the conductive anomalous bodies, while our dual-mesh inversion can well-recover the target bodies. The reason for this difference is that when we do the conventional inversion, we invert all underground structures simultaneously. For very limited data, the inversion is heavily dependent on the initial model and can easily be trapped in local minima. At this time, although the data are well-fitted, the recovered model may be far away from the true one (the so-called “equivalent trap”). In contrast, using our dual-mesh scheme, we invert the data first for rough but main underground structures, and then, we map them to the fine mesh as the starting model. Considering that our inversion on the fine mesh keeps the main underground structures unchanged while modifying the local ones, our method cannot be easily trapped in local minima. This inversion scheme follows the basic concept of multiscale inversions. From the data fitting for the two inversions shown in Figure 11, one sees that for the dual-mesh inversion, the data fitting of the inversion with coarse mesh is not good. After the inversion with fine mesh, however, the data fitting becomes as good as the single-mesh inversion. From the data fitting at the earth surface shown in Figure 12, we can also see that the data at the ground surface are well-fitted in both inversions.

3.3. Field Data Inversion

In this section, we take a dataset acquired from the Western Tonglu Mountain copper mine in Daye, China to test the effectiveness of our dual-mesh method. The lithology in the survey area are as follows. From 0 to 727.16 m, it is mainly quartz monzodiorite porphyrite; from 727.16 to 1031.80 m, it is marble with an embedded layer of copper and an iron ore body at the depth of 763.16–778.86 m. This ore body is controlled by the upper contact zone of low resistivity. From 1031.8 to 1036.3 m, it is skarn, while from 1036.3 to 1107.2 m, it is quartz monzodiorite. The copper-bearing hematite and magnetite located at the depth of 1040.4–1046.95 m are controlled by the lower contact zone. Compared to the upper contact zone, the contents of chalcopyrite and magnetite in the lower contact zone are significantly reduced, so that we can see the low-resistivity characteristics in this area.
Referring to Figure 6, the borehole ZK409 was drilled for SBTEM exploration. The depth of ZK409 is 1107.2 m. The positions of transmitting source and survey points are also shown in Figure 6. The geodetic coordinates of the starting and ending points for Source 1 are (3,332,207.8, 385,87,811.2) and (3,329,244.3, 38,587,804.4), while for Source 2 they are (3,332,207.8, 38,587,811.2) and (3,329,845.9, 38,589,298.8). The borehole ZK409 has coordinates of (3,329,616.93, 38,589,691.10).
Figure 13 shows the dual meshes used in our inversions. As we only observe the borehole data within the depth range of 600–900 m, we discretize this section into fine grids. We take the ground Ex and the borehole Ez as the inversion data and assume 5% of amplitudes as the noise floor. The air resistivity is set to 1 × 108 Ω·m, while the background resistivity is assumed to be 500 Ω·m. The regularization factor is selected following the same strategy as in the synthetic examples. The inversion for the coarse mesh runs 12 h and 51 min, and the RMS decreases to 11.8. After mapping to the fine mesh, the inversion runs again 52 h and 28 min, and the RMS decreases to 4.43. Figure 14 and Figure 15 show the data fitting at the earth surface and in the borehole. The surface data are well-fitted, while for the borehole data, the electric field Ez is overall well-fitted, except for the depth range of 620–720 m at early time channels. The good data fitting implies that the inversion using our dual-mesh scheme converges.
Figure 16 shows the inversion results of the survey data from the Western Tonglu Mountain copper mine. The low-resistivity anomalous body at 650–770 m is clearly revealed with a minimum resistivity of about 5 Ω·m that is consistent with the resistivity of the rock sample in Table 1 (mainly the copper, iron ore body). The recovered resistivity of the surrounding rock is from 300 to 1000 Ω·m, which is consistent with the high-resistivity quartz monzodiorite porphyrite with copper mineralization or middle-resistivity dolomitic marble and marble with copper mineralization. The logging data show that there exist two low-resistivity zones, respectively located at the depth ranges 70–340 m and 660–800 m, corresponding well to the locations of the conductive anomalies revealed in Figure 16. It should be pointed out that since, from the borehole survey, we only obtain the data at the depth range of 600–900 m, we are not able to clearly characterize the low-resistivity anomalies outside this range.
In Figure 16f–i, we plot four slices corresponding to the location of rock sampling. The rock sampling and the inversion results show relatively high-resistivity characteristics at z = 596.74 m and z = 802.43 m, while the rock sampling and the inversion results at z = 764.84 m and z = 879.74 m show relatively low-resistivity characteristics. The low-resistivity anomaly at z = 764.84 m is our target.

4. Discussion

The advantage of the SBTEM method is that the receivers are located in the underground near the abnormal body, so it can provide much higher resolution for a deep target than the conventional method. However, due to limited boreholes, when we use traditional dense local meshes (near the borehole) for inversion, the ambiguity is very strong, as can be seen from Figure 16b,c. In contrast, our dual-mesh inversion has distinctive advantages: (1) we use the coarse mesh to invert the data that can outline the main structures in the underground; (2) then, we use the fine mesh to invert the fine structure and obtain good data fitting. In this way, we can get good recovery to the underground structure and the good data fitting, and the uniqueness of the inversion can be improved.
Although the dual-mesh strategy overcomes ambiguity for the SBTEM data inversion to some extent, the limited resolution of the borehole data still cannot be solved. Only the targets close to the borehole receivers can be recovered, while the others cannot be restored. By jointly inverting the borehole and surface data, as can be seen from Figure 16c,e, the inversion result can be improved for the targets close to the borehole, but it does not work for the targets far away from the borehole.
It needs to be pointed out that the induced polarization (IP) effect is often observed when TEM with grounded lines is applied in the field, especially for massive sulfide or graphite exploration. However, from the rock-physics experiments on the samples collected from the borehole in our survey area, we did not see an obvious large IP effect; thus, we did not consider it in our 3D inversions. Considering that an IP effect can seriously affect the inversion results of TEM data, special attention needs to be paid in areas with an obviously strong IP effect.

5. Conclusions

Based on the unstructured finite-element method and back Euler discretization, we successfully developed a dual-mesh scheme for SBTEM data inversions. Numerical experiments for synthetic models showed that by adding the surface data to SBTEM inversions, we can get better results than when only the borehole data are used in the inversion. By taking the dual-mesh strategy, we can further improve the resolutions of both the target location and the resistivity. More inversions of synthetic data with multiple bodies under topographic earth showed that in contrast to the conventional methods that invert the data once, our dual-mesh method inverts the data in a multiscale way that can suppress the non-uniqueness of the SBTEM inversions to some extent. The inversion of the field dataset from the Western Tonglu Mountain copper mine also demonstrated that our dual-mesh method can deliver good results. The inversions are in good agreement with the resistivity measurements of the borehole rock samples and well logging data.

Author Contributions

Conceptualization, L.W., Y.L. and C.Y.; methodology, L.W., C.Y. and Y.L.; software, L.W. and Y.L.; formal analysis, L.W. and Y.L.; investigation, L.W., C.Y., Y.L., B.Z., X.R. and Y.S.; writing—original draft preparation, L.W.; writing—review and editing, L.W. and C.Y.; funding acquisition, C.Y., Y.L., B.Z. and X.R. All authors have read and agreed to the published version of the manuscript.

Funding

This paper is financially supported by the National Key Research and Development Program of China (2021YFB3202104) and the National Natural Science Foundation of China (42030806, 42074120, 42174167, 42274093).

Data Availability Statement

Data associated with this research are available and can be obtained by contacting the corresponding author.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

In the appendix, we will introduce how to calculate the SBTEM responses using the finite-element method. The governing equation for SBTEM forward modeling can be obtained from the following curl–curl equations, i.e.,:
× [ 1 μ 0 × E ( r , t ) ] + σ E ( r , t ) t + j s ( r , t ) t = 0 ,
where E(r,t) is the electric field, μ0 is the permeability of the vacuum, σ is the conductivity, and js(r,t) is the source current.
We first use the open-source code TetGen [29] to discretize the model domain into unstructured tetrahedral elements. Referring to Jin [30], for our forward modeling based on the finite-element method, we write the electric field at time t and location r in an interpolation format as E e ( r , t ) = j = 1 6 e j e ( t ) n j e ( r ) , where eje(t) denotes the electric field at the jth edge, while nje denotes the vector interpolation basis functions. Inserting this interpolation into Equation (A1) and multiplying both sides by a weighting function (it usually takes the basis function) and integrating over the element, we get:
A e d e e ( t ) d t + B e e e ( t ) + s e = 0 ,
where e e = [ e 1 e , e 2 e , , e n e ] are the electric fields at the element edges. The mass matrix Ae, the stiffness matrix Be, and the source term se can be written as:
A i j e = V e σ e n i e ( r ) n j e ( r ) dV ,
B i j e = V e × n i e ( r ) × n j e ( r ) dV ,
s i e = V e n i e ( r ) j s ( r , t ) t dV .
To solve Equation (A2), we adopt the implicit backward Euler scheme for time discretization, which is unconditionally stable [19,31]. For that purpose, we divide the time channels into several ranges and assume a small step for the early time but a large step for the later time. The initial time step is chosen to be 1/100 of the first time channel, and the rest time step is twice the previous time step.
Assembling the element matrix (A3) and (A4) into a global one and substituting it into Equation (A2), we obtain the following recursive formulation for the electric field, i.e.,:
C e n + 2 = A ( 4 e n + 1 e n ) 2 Δ t s n + 2 ,
where C = 3 A + 2 Δ t B . The final equation for the solution of the electric field can be written as:
[ C 1 4 A C 2 A 4 A C 3 A 4 A C 4 A 4 A C n ] [ e 1 e 2 e 3 e 4 e n ] = [ 2 Δ t 1 s 1 + 3 A e 0 2 Δ t 2 s 2 A e 0 2 Δ t 3 s 3 2 Δ t 4 s 4 2 Δ t n s n ] ,
where e0 represents the initial electric field.
Since we assume a grounded wire as the source to transmit a step-off current into the earth, we need to solve a Poisson’s equation for the initial DC field. We use the nodal finite-element method to calculate potentials at the nodes and then calculate e0 by interpolations and spatial derivatives. As the distance from the source increases, the grid size gets bigger. We need to expand the calculation domain far away from the source, so that we can apply the homogeneous Dirichlet boundary condition by setting the electric field on the outer boundary to zero. Equation (A7) can be simplified to:
K e = s TEM ,
where K is the coefficients matrix, and sTEM is the source term. The direct solver MUMPS [32] is used to solve Equation (A8) for the electric field at each edge. Then, the electric field at the receiver locations can be calculated by interpolations, while the magnetic field can be calculated by Faraday’s law.

References

  1. Li, J. Inversion Algorithms for Electromagnetic Problems in Well: Numerical Examples in Reservoir Target Detection. Ph.D. Thesis, China University of Geosciences, Wuhan, China, 2015. [Google Scholar]
  2. Augustin, A.M.; Kennedy, W.D.; Morrison, H.F.; Lee, K.H. A theoretical study of surface-to-borehole electromagnetic logging in cased holes. Geophysics 1989, 54, 90–99. [Google Scholar] [CrossRef]
  3. Malmqvist, L.; Pantze, R. Directional EM-measurements in boreholes. Geoexploration 1981, 19, 149–150. [Google Scholar] [CrossRef]
  4. Lane, R.J.L. The downhole EM response of an intersected massive sulphide deposit, South Australia. Explor. Geophys. 1987, 18, 313–318. [Google Scholar] [CrossRef]
  5. Silic, J.; Eadie, E.T. DHEM: The Que-Hellyer Volcanics experience. Explor. Geophys. 1989, 20, 65–69. [Google Scholar] [CrossRef]
  6. Vella, L. Petrophysical characteristics of BIF-hosted gold deposits and the application of downhole EM to their exploration, with examples from Hill 50 gold mine, Mt Magnet, Western Australia. Explor. Geophys. 1995, 26, 106–115. [Google Scholar] [CrossRef]
  7. Paggi, J.; Macklin, D. Discovery of the Eureka volcanogenic massive sulphide lens using downhole electromagnetics. Explor. Geophys. 2016, 47, 248–257. [Google Scholar] [CrossRef]
  8. Wei, B.; Zhang, G. Imaging method in surface-to-borehole electromagnetic system. Shi You Da Xue Xue Bao 1999, 23, 24–29. [Google Scholar]
  9. Zhang, Z.; Xiao, J. Inversions of surface and borehole data from a large-loop transient electromagnetic system over a 1-D earth. Geophysics 2001, 66, 1090–1096. [Google Scholar] [CrossRef]
  10. Chen, W.; Han, S.; Khan, M.Y.; Chen, W.; He, Y.; Zhang, L.; Hou, D.; Xue, G. A Surface-to-Borehole TEM System Based on Grounded-wire Sources: Synthetic Modeling and Data Inversion. Pure Appl. Geophys. 2020, 177, 4207–4216. [Google Scholar] [CrossRef]
  11. Chen, W.; Han, S.; Xue, G. Analysis on the full-component response and detectability of electric source surface-to-borehole TEM method. Acta Geophys. Sin. 2019, 62, 1969–1980. [Google Scholar] [CrossRef]
  12. Lajoie, J.J.; West, G.F. The electromagnetic response of a conductive inhomogeneity in a layered earth. Geophysics 1976, 41, 1133–1156. [Google Scholar] [CrossRef]
  13. West, R.C.; Ward, S.H. The borehole transient electromagnetic response of a three-dimensional fracture zone in a conductive half-space. Geophysics 1988, 53, 1469–1478. [Google Scholar] [CrossRef]
  14. Eaton, P.A.; Hohmann, G.W. The influence of a conductive host on two-dimensional borehole transient electromagnetic responses. Geophysics 1984, 49, 861–869. [Google Scholar] [CrossRef] [Green Version]
  15. Christoph, S.; Ralph-Uwe, B.; Klaus, S. Three-dimensional adaptive higher order finite element simulation for geo-electromagnetics—A marine CSEM example. Geophys. J. Int. 2011, 187, 63–74. [Google Scholar]
  16. Yin, C.; Qi, Y.; Liu, Y.; Cai, J. 3D time-domain airborne EM forward modeling with topography. J. Appl. Geophys. 2016, 134, 11–22. [Google Scholar] [CrossRef]
  17. Gu, G.; Li, T. Three-dimensional magnetotelluric inversion with surface topography based on the vector finite element method. Acta Geophys. Sin. 2020, 63, 2449–2465. [Google Scholar] [CrossRef]
  18. Kordy, M.; Wannamaker, P.; Maris, V.; Cherkaev, E.; Hill, G. 3-D magnetotelluric inversion including topography using deformed hexahedral edge finite elements and direct solvers parallelized on SMP computers; Part I, Forward problem and parameter Jacobians. Geophys. J. Int. 2016, 204, 74–93. [Google Scholar] [CrossRef]
  19. Um, E.S. Three-Dimensional Finite-Element Time-Domain Modeling of the Marine Controlled-Source Electromagnetic Method. Ph.D. Thesis, Stanford University, Stanford, CA, USA, 2011. [Google Scholar]
  20. Oldenburg, D.W.; Haber, E.; Shekhtman, R. Three dimensional inversion of multisource time domain electromagnetic data. Geophysics 2013, 78, E47–E57. [Google Scholar] [CrossRef] [Green Version]
  21. Li, J.; Lu, X.; Colin, G.F.; Hu, X. A finite-element time-domain forward solver for electromagnetic methods with complex-shaped loop sources. Geophysics 2018, 83, E117–E132. [Google Scholar] [CrossRef]
  22. Yang, D.; Fournier, D.; Kang, S.; Oldenburg, D.W. Deep mineral exploration using multi-scale electromagnetic geophysics; the Lalor massive sulphide deposit case study. Can. J. Earth Sci. 2019, 56, 544–555. [Google Scholar] [CrossRef]
  23. Yang, D.; Oldenburg, D.W. Three-dimensional inversion of airborne time-domain electromagnetic data with applications to a porphyry deposit. Geophysics 2012, 77, B23–B34. [Google Scholar] [CrossRef] [Green Version]
  24. Wang, L.; Yin, C.; Liu, Y.; Su, Y.; Ren, X.; Hui, Z.; Zhang, B.; Xiong, B. Three-dimensional forward modeling for the SBTEM method using an unstructured finite-element method. Appl. Geophys. 2021, 18, 101–116. [Google Scholar] [CrossRef]
  25. Hui, Z.; Yin, C.; Liu, Y.; Zhang, B.; Ren, X.; Wang, C. 3D inversions of time-domain marine EM data based on unstructured finite-element method. Acta Geophys. Sin. 2020, 63, 3167–3179. [Google Scholar] [CrossRef]
  26. Qi, Y.; Zhi, Q.; Li, X.; Jing, X.; Qi, Z.; Sun, N.; Zhou, J.; Liu, W. Three-dimensional ground TEM inversion over a topographic earth considering ramp time. Acta Geophys. Sin. 2021, 64, 2566–2577. [Google Scholar] [CrossRef]
  27. Key, K. MARE2DEM: A 2-D inversion code for controlled-source electromagnetic and magnetotelluric data. Geophys. J. Int. 2016, 207, 571–588. [Google Scholar] [CrossRef] [Green Version]
  28. Liu, Y.; Yin, C.; Qiu, C.; Hui, Z.; Zhang, B.; Ren, X.; Weng, A. 3-D inversion of transient EM data with topography using unstructured tetrahedral grids. Geophys. J. Int. 2019, 217, 301–318. [Google Scholar] [CrossRef]
  29. Si, H. TetGen, a Delaunay-Based Quality Tetrahedral Mesh Generator. ACM Trans. Math. Softw. 2015, 41, 1–36. [Google Scholar] [CrossRef]
  30. Jin, J.-M. The Finite Element Method in Electromagnetics, 2nd ed.; Wiley: New York, NY, USA, 2002. [Google Scholar]
  31. Um, E.S.; Harris, J.M.; Alumbaugh, D.L. 3D time-domain simulation of electromagnetic diffusion phenomena; a finite-element electric-field approach. Geophysics 2010, 75, F115–F126. [Google Scholar] [CrossRef]
  32. Amestoy, P.R.; Duff, I.S.; L’Excellent, J.-Y.; Koster, J. A Fully Asynchronous Multifrontal Solver Using Distributed Dynamic Scheduling. SIAM J. Matrix Anal. Appl. 2001, 23, 15–27. [Google Scholar] [CrossRef] [Green Version]
Figure 1. SBTEM survey configuration.
Figure 1. SBTEM survey configuration.
Remotesensing 15 01845 g001
Figure 2. Schematic diagram of attribute mapping. The coarse mesh is in black, while the fine mesh is in blue. The centroid of the fine mesh is (a) within a coarse mesh, (b) on the surface shared by two coarse meshes, (c) at the edge shared by several coarse meshes, and (d) at the vertex shared by several coarse meshes.
Figure 2. Schematic diagram of attribute mapping. The coarse mesh is in black, while the fine mesh is in blue. The centroid of the fine mesh is (a) within a coarse mesh, (b) on the surface shared by two coarse meshes, (c) at the edge shared by several coarse meshes, and (d) at the vertex shared by several coarse meshes.
Remotesensing 15 01845 g002
Figure 3. A synthetic model. (a) Plane view at z = 250 m; (b) cross section at x = 0 m. The orange blocks show the anomalous body. The two red lines denote the transmitting sources, while the blue dot at the left figure and the blue line at the right figure denote the borehole. The black dots at the left and right figures denote the receiver locations.
Figure 3. A synthetic model. (a) Plane view at z = 250 m; (b) cross section at x = 0 m. The orange blocks show the anomalous body. The two red lines denote the transmitting sources, while the blue dot at the left figure and the blue line at the right figure denote the borehole. The black dots at the left and right figures denote the receiver locations.
Remotesensing 15 01845 g003
Figure 4. Inversion parameters versus iterations for the model in Figure 3. (a) RMS; (b) Objective functional; (c) Model; (d) Trade-off parameter.
Figure 4. Inversion parameters versus iterations for the model in Figure 3. (a) RMS; (b) Objective functional; (c) Model; (d) Trade-off parameter.
Remotesensing 15 01845 g004
Figure 5. Inversions of the surface data, borehole data, and joint surface and borehole data for the model in Figure 3.
Figure 5. Inversions of the surface data, borehole data, and joint surface and borehole data for the model in Figure 3.
Remotesensing 15 01845 g005
Figure 6. Source locations and survey points.
Figure 6. Source locations and survey points.
Remotesensing 15 01845 g006
Figure 7. Three 3D abnormal bodies embedded in a topographic half-space.
Figure 7. Three 3D abnormal bodies embedded in a topographic half-space.
Remotesensing 15 01845 g007
Figure 8. Responses of receivers in the borehole in different times for the synthetic model in Figure 7 with and without anomalous bodies.
Figure 8. Responses of receivers in the borehole in different times for the synthetic model in Figure 7 with and without anomalous bodies.
Remotesensing 15 01845 g008
Figure 9. Inversion parameters versus iterations for the dual-mesh and conventional methods. (a) RMS; (b) Objective functional; (c) Model; (d) Trade-off parameter.
Figure 9. Inversion parameters versus iterations for the dual-mesh and conventional methods. (a) RMS; (b) Objective functional; (c) Model; (d) Trade-off parameter.
Remotesensing 15 01845 g009
Figure 10. True model and inversion results from the dual-mesh and conventional methods: (a,b) true model; (c,d) conventional inversions; (e,f) dual-mesh inversions.
Figure 10. True model and inversion results from the dual-mesh and conventional methods: (a,b) true model; (c,d) conventional inversions; (e,f) dual-mesh inversions.
Remotesensing 15 01845 g010
Figure 11. Borehole data fitting for dual-mesh and conventional methods.
Figure 11. Borehole data fitting for dual-mesh and conventional methods.
Remotesensing 15 01845 g011
Figure 12. Data fitting at the earth surface for dual-mesh and conventional methods.
Figure 12. Data fitting at the earth surface for dual-mesh and conventional methods.
Remotesensing 15 01845 g012
Figure 13. Meshes for our dual-mesh inversion of survey data from Daye, China.
Figure 13. Meshes for our dual-mesh inversion of survey data from Daye, China.
Remotesensing 15 01845 g013
Figure 14. Surface data fitting for the field data inversion.
Figure 14. Surface data fitting for the field data inversion.
Remotesensing 15 01845 g014
Figure 15. Borehole data fitting for the field data inversion.
Figure 15. Borehole data fitting for the field data inversion.
Remotesensing 15 01845 g015
Figure 16. Dual-mesh inversions of the survey data from the Western Tonglu Mountain copper mine in Daye, China. (a) 3D resistivity distribution; (b) slice at x = 0 m; (c) slice at y = 0 m; (d) ore body in borehole at x = 0 m; (e) ore body in borehole at y = 0 m; (fi) x,y slices at rock sample positions.
Figure 16. Dual-mesh inversions of the survey data from the Western Tonglu Mountain copper mine in Daye, China. (a) 3D resistivity distribution; (b) slice at x = 0 m; (c) slice at y = 0 m; (d) ore body in borehole at x = 0 m; (e) ore body in borehole at y = 0 m; (fi) x,y slices at rock sample positions.
Remotesensing 15 01845 g016
Table 1. Resistivities of rock samples from borehole ZK409.
Table 1. Resistivities of rock samples from borehole ZK409.
No.Borehole Depth (m)Core Propertiesρ0 (Ω·m)
wx-1596.74–596.84porphyritic quartz monzodiorite7178
wx-2764.84–764.94copper-bearing hematite magnetite ore8.7
wx-3802.43–802.54marble rock1125
wx-4879.75–879.85marble rock374
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Wang, L.; Liu, Y.; Yin, C.; Su, Y.; Ren, X.; Zhang, B. Three-Dimensional Dual-Mesh Inversions for Sparse Surface-to-Borehole TEM Data. Remote Sens. 2023, 15, 1845. https://doi.org/10.3390/rs15071845

AMA Style

Wang L, Liu Y, Yin C, Su Y, Ren X, Zhang B. Three-Dimensional Dual-Mesh Inversions for Sparse Surface-to-Borehole TEM Data. Remote Sensing. 2023; 15(7):1845. https://doi.org/10.3390/rs15071845

Chicago/Turabian Style

Wang, Luyuan, Yunhe Liu, Changchun Yin, Yang Su, Xiuyan Ren, and Bo Zhang. 2023. "Three-Dimensional Dual-Mesh Inversions for Sparse Surface-to-Borehole TEM Data" Remote Sensing 15, no. 7: 1845. https://doi.org/10.3390/rs15071845

APA Style

Wang, L., Liu, Y., Yin, C., Su, Y., Ren, X., & Zhang, B. (2023). Three-Dimensional Dual-Mesh Inversions for Sparse Surface-to-Borehole TEM Data. Remote Sensing, 15(7), 1845. https://doi.org/10.3390/rs15071845

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop