3D DC Resistivity Inversion with Topography Based on Regularized Conjugate Gradient Method

Abstract

During the past decades, we observed a strong interest in 3D DC resistivity inversion and imaging with complex topography. In this paper, we implemented 3D DC resistivity inversion based on regularized conjugate gradient method with FEM. The Fréchet derivative is assembled with the electric potential in order to speed up the inversion process based on the reciprocity theorem. In this study, we also analyzed the sensitivity of the electric potential on the earth’s surface to the conductivity in each cell underground and introduced an optimized weighting function to produce new sensitivity matrix. The synthetic model study shows that this optimized weighting function is helpful to improve the resolution of deep anomaly. By incorporating topography into inversion, the artificial anomaly which is actually caused by topography can be eliminated. As a result, this algorithm potentially can be applied to process the DC resistivity data collected in mountain area. Our synthetic model study also shows that the convergence and computation speed are very stable and fast.

1. Introduction

In DC resistivity exploration method, complex topography can generate artificial anomalies which will cause difficulty for the data interpretation. Based on Qiang and Luo’s work [1] on 3D DC finite element resistivity modeling with complex topography, we conducted 3D regularized inversion and imaging for this complex model.

The efficiency for 3D inversion problem depends primarily on 3 factors: efficient inversion algorithm, method for computing sensitivity matrix and the solver for a large liner system. Tripp et al. [2] introduced a method for calculating the sensitivity matrix, based on the relationship between electric potential and model parameter. This developed method for computing sensitivity matrix was successfully applied to a 2D DC resistivity inversion problem.

Park and Van [3] introduced 3D inversion based on finite difference method. Sasaki [4] also described similar 3D inversion algorithm but based on finite element method. These introduced 3D inversion methods work well to recover shallow resistivity anomalies but fail to produce high-resolution image for the deep anomalous bodies. The synthetic model studies show that the recovered resistivity imaging is quite different from the true model if the anomaly is located deep underground. Zhang et al. [5] conducted the research on 3D DC resistivity inversion using conjugate gradient method. Loke and Barker [6] introduced the E-SCAN (pole-pole array) 3D inversion technique where the Fréchet derivative matrix can be obtained by the halfspace analytical solution. The computation cost can be reduced significantly and it makes the 3D inversion become feasible. Wu and Xu [7–9] and Liu et al. [10] tested the application of E-SCAN method of measurement to 3D inversion problem with complex topography and synthetic study shows that this method works well. However, it is difficult to obtain the raw 3D data in real exploration since this method of measurement is time and money consuming. Moreover, E-SCAN method is based on AM (two electrodes) array and the resolution is lower than gradient array. Papadopoulos et al. [11] improved the efficiency of computing Jacobi matrix which can eliminate some unused parameters for inversion; save memory space, and speed up inversion without losing the accuracy. Tsourlos and Ogilvy [12] studied the 3D inversion of borehole-to-surface resistivity and IP data. Based on the smoothness constrained FEM algorithm, they get reasonable result. Huang et al. [13] increased the accuracy of inversion by introducing the volume factor and switching from global inversion to local inversion. However, this modified inversion process can only be applied to flat surface. Günther et al. [14] achieved 3D resistivity inversion for arbitrary topography based on finite element with unstructured mesh. Oldenborger and Routh [15] also introduced 3D resistivity inversion with the assistance of point spread function.

In this paper, we will focus on the study of fast 3D resistivity inversion problem for complex topography based on weighted regularized conjugate gradient method. The modified sensitivity matrix can be used to recover relatively deep anomalous body.

2. 3D DC Resistivity Forward Modeling and Inversion Theory

2.1. 3𝐷 Resistivity Forward Modeling with Complex Topography

The forward modeling method is based on finite method with triangular prism discretization [1] which can simulate variable topography. We approximated the unknown electric potential within each element by dual-linear function and mixed boundary condition [16] is implemented. By assembling each element equation, we can formulate a global linear system of equations as follows 𝐾_{(𝑥,𝑦,𝑧,𝜎)}𝑈_{(𝑥,𝑦,𝑧,𝜎)}=𝑃_{(𝑥𝑠,𝑦𝑠,𝑧𝑠,𝜏,Ω)}, (1) where 𝑥, 𝑦, and 𝑧 are the coordinate for the node; 𝜎 is the conductivity for the element; Ω and 𝜏 indicate the boundary at infinity and the earth’s surface; 𝐾 is the element equation matrix which is related to the location of nodes and element conductivity; 𝑈 is the electric potential, which is a function of node location and conductivity, in the node; 𝑃 contains the information of source location and boundary conditions. Equation (1) can be solved using incomplete Cholesey conjugate gradient method (ICCG). By decomposing the matrix 𝐾 and substituting back to (1), we can compute the electric potential 𝑈 (these values need to be stored to formulate the Jacobi matrix in inversion) at every node corresponding to all point electric sources. Apparent resistivity can be computed for different source-receiver array configuration.

The algorithm is implemented by Fortran language on PC machine. We adopted triangular prism discretization of the subsurface in order to simulate the topography and complex anomalous bodies underground.

2.2. Adaptive Weighed Regularization Inversion Algorithm

In a compact form, the forward modeling process described previously can be written as follows: 𝑑_𝑖=𝐴𝑚_𝑗,𝑖=1,2,…,𝑚,𝑗=1,2,…,𝑛, (2) where 𝑑_𝑖 is the observed data, 𝑚_𝑗 is the model parameter (conductivity in the element), and 𝐴 is a nonlinear operator which is implemented using finite element method. The inverse problem in (2) to recover model parameter from observed data is an ill-posed problem with nonunique solution.

To obtain stable solution for this ill-posed problem, one needs to consider regularization theory. The regularized theory has been applied to 3D DC resistivity inversion for a long time [11, 17–22]. In our study, we consider a minimization of Tikhonov parametric functional as follows 𝑃^𝛼(𝑑,𝑚)=‖𝐴(𝑚)−𝑑‖²+𝛼𝑠𝑚,𝑚_apr=min, (3) where 𝛼 is regularization parameter (Lagrange parameter), 𝑠(𝑚) is a nonnegative stabilizing functional of model parameter with the property of monotonic decreasing, and 𝑚_apr is some a priori information of model parameters.

It was shown by Zhdanov [23] that the regularization parameter can be selected using an adaptive scheme as follows 𝛼_𝑘=𝛼₀𝑞^𝑘,𝑘=1,2,3,…,𝑛,0<𝑞<1, (4) where 𝑞 is a real number smaller than 1 and 𝛼₀ is the initial value for regularization parameter which is selected in such a way to balance the misfit functional and the stabilizer as follows 𝛼₀=‖𝐴(𝑚)−𝑑‖²𝑚−𝑚_apr². (5)

The quality and amplitude of different data vary in geophysical application. Data weighting matrix is preferred in the inversion process in order to increase the contribution of data with good quality without suppressing geophysical information from data with low amplitude. After taking into account of data weighting, the Tikhonov parametric functional in (3) can be modified as following: 𝑃^𝛼(𝑑,𝑚)=𝑊_𝑑𝐴(𝑚)−𝑊_𝑑𝑑^𝑇𝑊_𝑑𝐴(𝑚)−𝑊_𝑑𝑑+𝛼𝑊_𝑚𝑚−𝑊_𝑚𝑚_apr^𝑇𝑊_𝑚𝑚−𝑊_𝑚𝑚_apr=min, (6) where, 𝑊_𝑑 is the data weighting matrix which can be selected according to the amplitude and quality of different data, the superscript 𝑇 indicates the complex transpose, 𝑊_𝑚 is the model weighting matrix which results in practically equal resolution of the inversion with respect to different parameters of the model. It was shown that the model weighting matrix can be computed as follows in such equal resolution criteria: 𝑊_𝑚=diag𝐹^∗𝐹^1/8, (7) where 𝐹 is the Fréchet derivative matrix corresponding to the nonlinear forward modeling operator 𝐴 in (2). 𝐹^∗ is the transpose of 𝐹.

We used regularized conjugated gradient method [23] for the minimization of (6): 𝐫_𝑛=𝐴𝑚_𝑛𝐥−𝐝,_𝛼𝑛=𝐹^∗𝑊²_𝑑𝐫_𝑛+𝛼𝑊²_𝑚𝑚_𝑛−𝑚_apr,𝛽^𝛼_𝑛=𝐥^𝛼_𝑛²𝐥^𝛼_𝑛−1²,𝐥^𝛼_𝑛̃=𝐥^𝛼_𝑛+𝛽^𝛼_𝑛𝐥̃^𝛼_𝑛−1,𝑘^𝛼_𝑛=𝐥^𝛼_𝑛̃,𝐥^𝛼_𝑛𝑊_𝑑𝐹_𝑚𝑛𝐥̃_𝑛²+𝛼𝑊_𝑚𝐥̃_𝑛²,𝑚_𝑛+1=𝑚_𝑛−𝑘^𝛼_𝑛𝐥^𝛼_𝑛̃, (8)

where 𝐫_𝑛 is the residual vector between the observed data and predicted data.

Based on the algorithm in (8), we implemented 3D regularized inversion for DC resistivity with complex topography.

2.3. The Formulation of Fréchet Derivative Matrix

We usually fix the size of cells for model parameter in the DC resistivity inversion and assume that the apparent resistivity on the earth’s surface is only function of the conductivity distribution in the cells. As such, the Fréchet derivative matrix is computed as follows: 𝐹=𝜕𝜌_si𝜕𝜎_𝑗, (9) where 𝜌_si is the measured apparent resistivity data and 𝜎_𝑗 is the conductivity value for the 𝑗th cell.

As we know, the measured apparent resistivity on the earth’s surface is related to the electric potential in the observation position. Take the dipole-dipole array as an example: 𝜕𝜌_𝑆^{𝐴𝐵𝑀𝑁}=𝜕𝜎−𝐺^{𝐴𝐵𝑀𝑁}𝐼_𝐴𝐵𝜕𝑈^𝐴_𝑀−𝜕𝜎𝜕𝑈^𝐴_𝑁−𝜕𝜎𝜕𝑈^𝐵_𝑀+𝜕𝜎𝜕𝑈^𝐵_𝑁,𝜕𝜎 (10) where 𝐺^{𝐴𝐵𝑀𝑁} is array coefficient which is a function of electrode spacing, 𝐼_𝐴𝐵 is the injected current, and 𝑈^𝐴_𝑀 is the electric potential at 𝑀 when current is injected from electrode 𝐴. Equation (10) indicates that the derivative of apparent resistivity to the model conductivity can be transformed to the derivative of electric potential to the conductivity in each cell. The electric potential on the earth’s surface can be computed from (1). By considering the first order derivative of both sides of (1) to the model parameter, we can obtain the following formula: 𝐊𝜕𝐔+𝜕𝜎𝜕𝐊𝜕𝜎𝐔=0, (11) which can be written as follows 𝐊𝜕𝐔𝜕𝜎=−𝜕𝐊𝜕𝜎𝐔. (12)

Equation (12) is in a similar to that of as (1). The right side of (12) can be treated as an auxiliary current. We can solve the linear system of equations in (1) to obtain the derivative of the electric potential on the earth’s surface to the conductivity in each cell. This derivative can also be computed by using the following equations: 𝜕𝑈^𝐴_𝑀𝜕𝜎_𝑒=−𝑈^𝐴₁𝑈^𝐴₂𝑈^𝐴₃𝑈^𝐴₄𝑈^𝐴₅𝑈^𝐴₆⋅𝜕𝑘₁₁𝜕𝜎𝜕𝑘₁₂𝜕𝜎𝜕𝑘₁₃𝜕𝜎𝜕𝑘₁₄𝜕𝜎𝜕𝑘₁₅𝜕𝜎𝜕𝑘₁₆𝜕𝜎𝜕𝑘₂₁𝜕𝜎𝜕𝑘₂₂𝜕𝜎𝜕𝑘₂₃𝜕𝜎𝜕𝑘₂₄𝜕𝜎𝜕𝑘₂₅𝜕𝜎𝜕𝑘₂₆𝜕𝜎𝜕𝑘₃₁𝜕𝜎𝜕𝑘₃₂𝜕𝜎𝜕𝑘₃₃𝜕𝜎𝜕𝑘₃₄𝜕𝜎𝜕𝑘₃₅𝜕𝜎𝜕𝑘₃₆𝜕𝜎𝜕𝑘₄₁𝜕𝜎𝜕𝑘₄₂𝜕𝜎𝜕𝑘₄₃𝜕𝜎𝜕𝑘₄₄𝜕𝜎𝜕𝑘₄₅𝜕𝜎𝜕𝑘₄₆𝜕𝜎𝜕𝑘₅₁𝜕𝜎𝜕𝑘₅₂𝜕𝜎𝜕𝑘₅₃𝜕𝜎𝜕𝑘₅₄𝜕𝜎𝜕𝑘₅₅𝜕𝜎𝜕𝑘₅₆𝜕𝜎𝜕𝑘₆₁𝜕𝜎𝜕𝑘₆₂𝜕𝜎𝜕𝑘₆₃𝜕𝜎𝜕𝑘₆₄𝜕𝜎𝜕𝑘₆₅𝜕𝜎𝜕𝑘₆₆⋅𝑈𝜕𝜎_𝐴1𝑈^𝐴₂𝑈^𝐴₃𝑈^𝐴₄𝑈^𝐴₅𝑈^𝐴₆. (13)

As we can see the computation cost in formulating the Jacobi matrix is trivial since the electric potential in all the nodes corresponding to different source location already precomputed at the stage of forward modeling and the stiffness matrix is also computed in forward modeling process. The computation of Jacobi matrix is only a transformation which utilizes the forward modeling result. In real application, we only need the forward modeling result on the nodes which locate on the earth’s surface instead of the electric potential on all nodes. Also, one needs to notice that only one forward modeling is required to obtain Jacobi matrix. As such, the computation cost for inversion is reduced dramatically and it makes the inversion be implemented on a PC machine.

3. Discussions for Inversion

3.1. Stability Issue for Inversion

In 3D resistivity inversion, the misfit for inversion will stay in a relative large level if the initial model is not properly chosen which will cause the background conductivity to be either too small or too large. In this study, we implemented some computation to calculate the maximum and minimum value of apparent resistivity and select reasonable background conductivity from this information to be our initial model for inversion.

Another main issue is that the inversion will be unstable if the variation of inversion parameter is either too large or too small. In order to solve this problem, we will take the logarithm of the model parameter and do the inversion in this new space. As a result, the variation of model parameter is in control. Moreover, the upper and lower limits for model parameter can be set in the inversion by applying this method.

3.2. Speed for Inversion

The modern PC machine is already very powerful now with the development of CPU and increase of memory. However, it can still take hours and even days to run an inversion if the number of nodes increases dramatically especially in 3D problem. As a result, it is very important to choose proper method of data storage and fast solver for linear system of equations. Incomplete Cholesey conjugate gradient method (ICCG) works well for DC resistivity method with many moving electric source. By using this method, we only need to decompose the stiffness matrix once which takes approximately several minutes and the substitution is also very fast. The computation of Jacobi matrix is a time-consuming process. In this paper, we used the method of electric potential combination to form the Jacobi matrix which can reduce the computation time dramatically. We run the inversion in a PC machine and the model discretization is 44 × 22 × 17 × 2 which generates 32912 cells. The number of observation is 1505 with 5 profiles and 155 points to inject currents. The forward modeling takes around 5 minutes and each iteration of inversion takes approximately 6 minutes. The formulation of Jacobi matrix takes around 30 seconds and it takes other 30 seconds to update the model. The total time for inversion is around 1 hour if the number of iteration is set to be 10.

3.3. Sensitivity Analysis

The sensitivity we mentioned here describes how sensitive the measured apparent resistivity to the change of conductivity in each element is. As we know, the electric potential decreases as the distance between the observation point and the source increases. The apparent resistivity also decreases as the observation point moves away from the source since the apparent resistivity is proportional to the electric potential. Figure 1 shows a sensitivity distribution in 𝑋𝑍 vertical section. As we can see from Figure 1, the data is more sensitive to shallow structure than deep structure. Figure 2 shows the sensitivity changes as a function of depth at 𝑋=0. The curve in Figure 2 concaves upward and the sensitivity decreases dramatically as depth increases. As a result, the inverted anomaly will move upward if we use the sensitivity shown in Figures 1 and 2. We can modify the sensitivity distribution to make it decay slower than original sensitivity in order to increase the resolution of deep structure. Figures 3 and 4 show the modified sensitivity distribution. The synthetic model shows that this modified sensitive distribution works well to recover the true model (Figure 5).

Figure 1 

Original sensitivity distribution in 𝑋𝑍 section.

Figure 2 

The original sensitivity changes as a function of 𝑍 at 𝑋=0.

Figure 3 

Modified sensitivity distribution in 𝑋𝑍 section.

Figure 4 

The modified sensitivity changes as a function of 𝑍 at 𝑋=0.

(a)

(b)

Figure 5 

A comparison of inversion between original (a) and modified (b) sensitivity at section 𝑌=0.

4. Numerical Examples

The following synthetic model studies are based on the dipole-dipole array.

4.1. Simple Model

For the simple case, we have two different models. Model 1 (Figure 6(a)) is a 4 m × 4 m × 4 m conductive body buried 3 meters underground with the resistivity of 10 Ωm. The resistivity of wall rock is 100 Ωm. Both the receiver spacing and line spacing are 1 meter. There are 25 points for injecting currents in each profile. The background resistivity is approximated from forward modeling and it is used for the initial model in the inversion.

(a) 3D low resistivity body without topography

(b) 3D low resistivity body model with topography

(c) Apparent resistivity at 𝑌=0 for model 1

(d) Apparent resistivity at 𝑌=0 for model 2

(e) Vertical section of inversion result at 𝑌=0 for model 1

(f) Vertical section of inversion result at 𝑌=0 for model 2

Figure 6 

A comparison of inversion result for the model without and with topography.

Model 2 (Figure 6(b)) is almost the same as the previous one. Based on the previous model, we added some topography to this model. The highest elevation of the topography is 4 meters. The slope in Figure 6(b) is 30 degrees which leads the length of slope to be 8 meters.

The inversion results belong the profile at 𝑌=0 are shown here. Figures 6(c) and 6(d) are the contourline of apparent resistivity at 𝑌=0 for models 1 and 2. Figure 6(c) shows an apparent resistivity distribution for a standard dipole-dipole anomaly. The resistivity of wall rock below and above the anomaly is well recovered and the value is close to 100 Ωm. For Figure 6(d), the apparent resistivity is distorted by the artificial anomaly caused by topography. We can see that the resistivity of wall rock is almost 3 times as the true value and the anomaly caused by the 3D low resistivity structure is unclear.

Figures 6(e) and 6(f) show one of vertical sections of inversion results at 𝑌=0. We can see that the 3D conductive anomaly is well recovered either for the flat surface model or the model with topography without the redundant structures.

4.2. Complex Model

The illustration of this model is shown in Figure 7(a) with two anomalous bodies below some valley. The first anomalous body is a highly resistive body (500 Ωm) buried 2 meters underground with the size of 4 m × 4 m × 4 m. The second anomaly is a conductive body (10 Ωm) burried 4 meters below the surface with the size of 4 m × 4 m × 6 m. The horizontal distance between these two anomalous bodies is 4 meters. The resistivity of the wall rock (background resistivity) is 100 Ωm. The observation profile is oriented to 𝑥 direction with observation spacing of 1 meter. There are 31 points for injecting current along each profile. The survey is conducted along five profiles with the spacing of 1 meter. The initial model is obtained in a similar way to the one we described in Section 4.1.

Figure 7 

The inversion result for a complex model. (a) A model with two bodies under a valley. (b) Slices view of apparent resistivity. (c) Slices view of 3D inversion result.

Figure 7(b) shows some slices of the apparent resistivity distribution for a standard dipole-dipole anomaly. The apparent resistivity is distorted by the artificial anomaly caused by topography.

After 11 iterations which take 4115 seconds, the normalized misfit decreased from 12.63% to 4.97% which is around the noise level of the data. The inversion result is shown in Figure 7(c). From the figure, we can see that the two anomalous bodies are well recovered and the artificial anomaly caused by the topography is eliminated. Figure 8 shows that the predicted data is very similar to the observed data on section 𝑌=0.

(a)

(b)

Figure 8 

A comparison between observed apparent resistivity (a) and predicted apparent resistivity (b) for profile 𝑌=0.

5. Case Study

The primary object of the survey is to estimate the basement depth which can be used as some information for building house. The survey site is on the river bed where the topography is quite flat. We employed the multielectrode resistivity method with dipole-dipole array configuration. Three survey lines are conducted with line distance of 30 m and electrode distance of 1.5 m; totally, we have 60 electrodes. The survey is conducted in the winter season and the frozen soil with thickness of 1 m–1.5 m caused the high resistivity layer near the surface. According to geological information, the surface layer is quaternary alluvium which mostly consisted of loess, silt, and gravel. The lower layer is Jurassic red sandstone and as such it is a relative stable layer.

Based on the observed data for our survey configuration, we discretized our subsurface into 80 × 60 × 40 cells. We choose 100 Ωm as the resistivity of the wall rock. After 10 iterations, the normalized misfit reaches 12.5%. The inversion result is shown in Figure 9. In 55 m along profile 135, drilling shows that the first 5 m consisted of frozen soil, sand, and gravel; from 5–9.8 m, we observed water-bearing ooze; below 9.8 m is red sandstone. Our inversion result shows very good fitting with the information that we get from drilling.

(a) Survey configuration

(b) Inversion result

Figure 9 

The inversion of real data.

6. Conclusion

In this paper, we implemented the 3D DC resistivity inversion based on regularized conjugate gradient method. We used finite element method with triangular prism discretization for forward modeling of apparent resistivity. Synthetic model studies show that the inversion algorithm works well to recover the true model with complex topography. The inversion process is very stable and computation speed is fast.

The electric potential at each node caused by different source is precomputed in the forward modeling and saved for the inversion stage to compute Jacobi matrix. By doing this, the computation cost can be reduced dramatically in the inversion process. Only one forward modeling is required in order to compute the sensitivity matrix.

We conducted a sensitivity analysis and introduced a weighting function to optimize sensitivity matrix. By applying this proper model weighting, the vertical location of the anomalous body can be well recovered.

Acknowledgments

The authors acknowledge National Natural Science Foundation of China, Grants nos. 40974076, 41174104, Chinese Academy Geological Sciences (Sinoprobe-03-02-04), and the Key Laboratory of Metallogenic Prediction of Nonferrous Metals of Ministry of Education for the support of this project.