Surface MMR enhanced 3D inversion using model reconstruction strategy

Journal of Applied Geophysics 159 (2018) 193–203

Contents lists available at ScienceDirect

Journal of Applied Geophysics journal homepage: www.elsevier.com/locate/jappgeo

Surface MMR enhanced 3D inversion using model reconstruction strategy Yue Yang, Aihua Weng ⁎, Sirui Li, Dajun Li, Shiwen Li, Jianping Li College of Geo-exploration Sciences and Technology, Jilin University, Changchun, Jilin 130021, China

a r t i c l e

i n f o

Article history: Received 1 February 2018 Received in revised form 20 July 2018 Accepted 28 August 2018 Available online 29 August 2018 Keywords: Magnetometric resistivity 3D inversion Non-linear conjugate gradient inversion Model reconstruction

a b s t r a c t Previous inversions of surface magnetometric resistivity (MMR) could obtain only the ratio, not the absolute value of earth electrical conductivity. In this paper, instead of solving Poisson's equation, the magnetic fields generated by low frequency alternating current(AC) line are looked as MMR data and calculated from electric-type Helmholtz's equation in frequency domain, so that the absolute conductivity could be taken into account. Forward modeling scheme is adapted from staggered-grid finite-difference method. Inversion algorithm is non-linear conjugate gradient (NLCG). The magnetic field perpendicular to line source is selected as observed data. To improve depth resolution, we introduce a model reconstruction strategy based on reference model. In the strategy, the upper portion of the resultant model from previous inversion is replaced by a homogeneous layer to construct a new model to trigger next inversion. The thickness of the homogeneous layer is called reconstruction depth, it is equal to the depth where the models difference between last two iterations after the run of the ith inversion is almost none-zero. The reconstruction could be applied until the model difference is almost zero or the difference of model differences between two adjacent inversions is less than a predefined threshold. Inversion results from synthetic and real data have proved that the proposed multi-reconstruction inversion strategy is able to obtain absolute resistivity, as well as enhance the deep resolution of surface MMR data inversion. © 2018 Elsevier B.V. All rights reserved.

1. Introduction Magnetometric resistivity(MMR) method measures the magnetic fields associated with artificially created, non-inductive (DC or pseudo DC) current flow input into the earth through two electrodes (Chen et al., 2002). Unlike the traditional resistivity method measuring electrical potential difference on the earth surface, only magnetic field is measured in MMR method, and the anomalous part of the magnetic field induced by target body is a function of its geometry and relative resistivity (Bouchedda et al., 2015). Thus, by studying the magnetic field spatial variations, resistivity ratio of the earth could be obtained (Deng, 2005). There are several variants of MMR, including surface magnetometric resistivity (MMR) (Edwards, 1974; Edwards & Nabighian, 1991; Chen et al., 2002; Chen & Oldenburg, 2004), magnetic induced polarization (MIP) (Seigel, 1974; Howland-Rose et al., 1980a, 1980b; Chen & Oldenburg, 2003; Li et al., 2010) and sub-audio magnetic(SAM) (Cattach et al., 1993; Fathianpour & Cattach, 1995; Fathianpour et al., 2005a, 2005b). Unlike the latter two which are mainly used in ground exploration, surface MMR has extended its applications from ground-

⁎ Corresponding author. E-mail address: [email protected] (A. Weng).

https://doi.org/10.1016/j.jappgeo.2018.08.015 0926-9851/© 2018 Elsevier B.V. All rights reserved.

based observations to borehole explorations (Acoata & Worthington, 1983; Bishop et al., 1997; Godber and Bishop, 2007) and even marine explorations (Nobes et al., 1986; Yang, 2005; Weng et al., 2009; Seama et al., 2013). And in addition to its usage as traditional applications like ore prospecting and geological mapping (Edwards & Howell, 1976; Meju, 2002), surface MMR has carved out new places such as engineering investigations and environmental monitoring (Zhu & Yang, 2008; Kofoed, et al., 2011; Maxwell et al., 2013). There may be two factors that hinder the wide application of MMR method. First one is the fact that MMR could obtain only the ratio or the relative difference of the resistivity (Bouchedda et al., 2015). The behind reason is that MMR tried to use DC grounded line as source to create magnetic field, which is in fact static field. This static magnetic field at ground surface is insensitive to 1-D resistivity. When in non 1D case, if the conductivity contrast between targets and their surroundings were same, the magnetic fields would be identical (Edwards et al., 1978; Yang & Tseng, 1992; Chen et al., 2002; Chaladgarn & Yooyuanyong, 2013) (shown in Appendix A). Another factor is the relative slow development of MMR inversion. In the 1980s, researchers used typical curve method (Howland-Rose et al., 1980b) or trial-and-error interpretation scheme (Szarka, 1987; Asten, 1988) for 2D models. Until the beginning of the century, Chen et al. (2002) for the first time had developed conjugate gradient least-

194

Y. Yang et al. / Journal of Applied Geophysics 159 (2018) 193–203

squares (CGLS) and adaptive regularization 3D inversion method. Later, Labrecque et al. (2003) applied Occam's method to jointly inverse ERT and MMR data. In surface MMR inversion, a challenge is poor deep inversion resolution (Chen et al., 2002). Weng et al. (2015) also found only Hy inversion could obtain shallow abnormities in 3D controlled source electromagnetic method. Li & Oldenburg (1996) blamed the insufficient depth resolution of magnetic field on the rapid decaying of kernel function with depth, so a weight function was introduced to compensate the geometrical decaying of the kernel function, thus improving depth resolution. This idea later was used in gravimetrical field data inversion (Li & Oldenburg, 1998). Chen et al. (2002) used almost the same depth weighting methods in 3D surface MMR inversion, but how to select weighting function still relies on empirical estimation or numerical simulation, because it is difficult to separate explicitly the decay factor from kernel function. Consequently, different weighting functions will cause different inversion results (Li & Oldenburg, 1996, 1998; Chen et al., 2002). Portniaguine & Zhdanov have used a minimum support stabilizing functional to generate a sharp, focused inverse image which is applied to induced magnetic field (Portniaguine and Zhdanov, 2000). They used the square root of the integrated sensitivity matrix to determine a weighting matrix of the focusing inversion. This method also improved the effect of the magnetic field data inversion. In order to retrieve the absolute resistivity and meanwhile improve depth resolution, the magnetic fields from a grounded line source at a very low frequency are calculated and regarded as that measured in MMR. Therefore, Helmholtz's equation is solved in frequency domain, thus establishing a direct relationship between magnetic field and resistivity, laying a foundation to obtain conductivity rather than relative resistivity. The staggered-grid finite difference method is used for forward. We also use non-linear conjugate gradient (NLCG) method incorporated with model reconstruction technique (MRT) to improve 3D MMR inversion resolution. MRT has been initially proposed by Ye et al. (2013) in 2D magnetotelluric (MT) inversion to effectively suppress inversion non-uniqueness on reference model. In their work, 2D MT inversions were done in sequential mode. And MRT is used to construct a reference model which is in fact a combination of a homogeneous half-space, called auxiliary model, weighted by previous inversion result, called meta model. The weighted model retains the real information of model in shallow part, while discard the deep part by using a homogeneous electric structure, thus the new model could meet the Dirichlet's boundary condition at bottom of model required by MT forward modeling. Recently, a modified MRT has been applied successfully to AMT 2D inversion (Kong et al., 2015). Contrary to the application of MRT in MT inversion, in MMR inversion, the upper portion of reference model is replaced by a homogeneous layer, forcing the depth resolution increase downward. The thickness of the upper portion, also called reconstruction depth, plays a very important role in MRT. It could be determined as the depth deeper from which the model difference between last two iterations after the ith inversion. The part of inversion model above the reconstruction depth will be replaced by a homogeneous layer. The new model is then used as reference model to run next inversion. MRT will be conducted several times in order to make the inversion model more accurate. This procedure is called multi-construction and terminates when model difference vanishes 0 or the difference of model differences between two adjacent inversions is less than a predefined threshold. In the paper, the forward modeling based on AC is outlined at first. Then, in the second section, brief review of NLCG inversion theory is followed by the concept of reconstruction depth and how to estimate the depth. Thereafter, the method is verified by a synthetic case, and followed by an example done in Australia. Conclusion at last summarizes the papers.

2. Forward modeling 2.1. Data type of MMR In surface MMR method, if x-direction is parallel with the source direction, three possible components of magnetic field Bobs could be observed Bxobs ¼ Bxa þ Bxn ;

ð3aÞ

Byobs ¼ Bya þ Byn ;

ð3bÞ

and Bzobs ¼ Bza þ Bw ;

ð3cÞ

here, Ba is the anomalous magnetic field generated by the current density abnormal due to conductivity disturbance. Bn is the normal magnetic field, which is associated with the current flowing in the background earth. Bw is the magnetic field excited by current flowing wire. Chen et al. (2002) suggested that Byobs is the optimal component for MMR inversion. Therefore, we inverse the y-component of magnetic field as well. There is significance behind this choice. Traditionally, in surface MMR measurement, U-shaped wire is used in order to reduce the influence of Bw (Chen et al., 2002; Fathianpour et al., 2005b; Li et al., 2010). From the Eqs. (3a), (3b), (3c), however, it could be found that Bw only contributes to Bzobs. As a result, when selecting By as observing parameter, the data will be irrelevant to the shape of wire. In such case, we could replace U-shaped wire by a grounded straight line source, which not only facilitate the field work, but also is convenient to do 3D forward modeling. 2.2. Governing equation When the source is an AC wire, the induced magnetic field will relate to absolute resistivity. Therefore, in this paper, the magnetic field in frequency domain is simulated from the electrical type Helmholtz equation
 2 2 2 ∇  ∇  Ea −k Ea ¼ k −kn En

ð1Þ

where Ea is the secondary electric field, and En is the primary electric field. k2n is the background media complex wavenumber. The staggered grid finite-difference method is used to generate discrete form of Eq. (1) (Yee, 1966; Alumbaugh et al., 1996). The total electric fields E = En + Ea, from which the magnetic field B at the center of cell faces is calculated by Faraday's law. The observed magnetic field Bobs is interpolated from B at last (Egbert and Kelbert, 2012). 2.3. Accuracy verification In order to verify the accuracy of the above forward algorithm, we compared the response of Chen's model (2002) with our numerical solution. The model is shown in Fig. 1a. The ground line source is 1200 m long in x direction and locates directly above the anomaly. Frequency is 0.3 Hz. The side length of cube is 400 m, and the depth is 80 m. The resistivity of cube is 10 Ω·m, while the background is 100 Ω·m. The survey area is 800 m × 800 m, with point and line space 40 m. The whole model is discretized into a fine grid of 34 × 34 × 30 cells, with the grid size of each cell 40 m × 40 m × 40 m. The horizontal component of anomalous magnetic field over the area bounded by the coordinates −400 to 400 m in both x and y directions, is shown in Fig. 1. Fig. 1b is the numerical solution calculated by the method above, and Fig. 1c is the numerical solution in Chen et al. (2002). It could be seen that the two solutions are almost the same. To illustrate the difference, we plot

Y. Yang et al. / Journal of Applied Geophysics 159 (2018) 193–203

195

Fig. 1. The numerical solutions of two methods and their difference, the unit in pT. (a) schematic of Chen's model (Chen et al., 2002), A and B present the locations of the source electrodes, respectively; (b) model responses for our method; (c) results from Chen in 2002; (d) response difference between (b) and (c).

the difference between them in Fig. 1d. It can be seen that at the middle of the survey area, the differences are comparatively small; near the ends of wire, the accuracy is lost, but the observations here are usually not allowed. 3. Non-linear conjugate gradient inversion Inversion objective function ϕ(m,d), after Egbert and Kelbert (2012), is defined as T 1 ϕðm; dÞ ¼ ðd−f ðmÞÞ C 1 d ðd− f ðmÞÞ þ λðm−mref Þ C m ðm−mref Þ: T

ð4Þ

here, m is an M-dimensional resistivity model parameter vector, which provides an adequate fit to a data vector d of dimensions Nd. mref is reference model parameter vector, which has a great influence on the inverted model m. f(m) defines the forward mapping. Cd and Cm denote data and model covariance matrix, respectively. Cd is the covariance of data errors, which is always taken to be diagonal, so its operation is just a rescaling of the data predicting error. Although the data error affects inversion (Kelbert et al., 2008), numerical tests showed that an appropriate constant is a candidate of data error. In this paper, according to magnetic field strength, we set it to be 10−8. Cm is model covariance matrix, used to weight or smooth model parameters (Egbert & Kelbert, 2012). Here, in this paper, it is used as smoothing operator (Egbert & Kelbert, 2012). In three directions, they are set as 0.3 (Campanyà et al., 2016). While λ is regularization parameter. Because we use non-linear conjugate gradients method to run inversion, the parameter could be reduced internally during inversion

according to the variation of objective function; therefore any rational value is able to start the inversion. In this paper, we accept 0.8 as the initial value. Newman and Alumbaugh (2000) firstly used NLCG approach to minimize Eq. (4). NLCG inversion tries to modify m through conjugate gradients of the objective function in Eq. (4). It has been applied to magnetotelluric method (Hu et al., 2006; Egbert & Kelbert, 2012; Zhang et al., 2013; Dong et al., 2014), grounded controlled source electromagnetic method (Lin et al., 2012; Weng et al., 2012, 2015) and aerial electromagnetic method (Liu and Yin, 2013). The most advantage of NLCG is that J needs not direct to be computed and stored, instead, a adjoint forward is done to compute the multiplication of J and or JT with data or model vector. Therefore, in the context, NLCG approach is accepted to solve the 3D MMR inversion problem. The flowchart of the algorithm is listed as following. 1) set i = 1, mi = mref and compute ri = − ∇φ(mi), mref is the reference model; 2) stop when ‖ri‖ is sufficiently small. Otherwise, set ui = M‐1 i ri, M is a preconditioning operator, and u is the specified conjugate search direction; 3) find αi that minimizes φ(mi + αui). If αi is too small, reduce the regularization factor λ; when λ is below a threshold value, inversion terminates; 4) set mi+1 = mi + αiui and ri+1 = − ∇φ(mi+1); 5) stop when ‖ri+1‖ is sufficiently small; otherwise go to step (6); 6) set βiþ1 ¼

‐1 T riþ1 M ‐1 iþ1 r iþ1 −r iþ1 M i r i

rTi M ‐1 i ri

;

T 7) set ui+1 = M‐1 i+1ri+1 + βi+1ui; 8) set i = i + 1, and return to step (3).

196

Y. Yang et al. / Journal of Applied Geophysics 159 (2018) 193–203

To verify the feasibility of the inversion method in surface MMR, we design a simple model shown in Fig. 2a–b. A 10 Ω·m cube buried in a uniform half-space of resistivity 400 Ω·m. For simplicity, the cube is located right below the origin of coordinate system. The source is 1200 m along the x-axis at 0.3 Hz. Survey area at the ground extends from −400 to 400 m in both x and y directions, with 21 profiles and 21 sites on each profile. The ground is discretized into 40 × 40 × 30 cells. Reference model is a 400 Ω·m homogeneous half-space. Gaussian random noise at 5% signal level is added to the data. The inversed model after 120 iterations is displayed in Fig. 2c. It can be seen that the resistivity of abnormal body has been recovered to a possible extent, but the abnormal body tends to be concentrated near the surface. Although the Helmholtz equation establishes a direct relationship between magnetic field and resistivity, the depth resolution of the abnormal body remains unsolved (Li & Oldenburg, 1996, 2000; Chen et al., 2002; Zhdanov, 2002). Two possible factors contribute to this problem, one is the fast decay of the field in the body, and the other is the not enough observations and inaccuracy of the observation (Li & Oldenburg, 1998; Zhdanov, 2002; Sobouti et al., 2015). 4. Model reconstruction technique 4.1. Basic idea In order to enhance deep resolution, a model reconstruction technique (MRT) is developed to construct the reference model in Eq. (4). Previously, we run inversion only once, with mref being set to be m0, which is called auxiliary model. But in MRT, inversion will be run one after another for several times. In each inversion, the reference model will be reconstructed from previous inversed model, which is called meta model. In our MRT, the part above a certain depth in the meta model is replaced by a certain homogeneous layer from the auxiliary model to reconstruct reference model. This certain depth is called reconstruction depth. Once an inversion is done, another new hybrid model in the same way is reconstructed as reference model to run next inversion. A complete reconstruction procedure will include several times inversions, each with a gradually modified reference model. This reconstruction procedure is illustrated in Fig. 2. In this way, MRT has a capability to retain the deep resistivity structure from previously inversed model, and transfers previous model as ‘prior information’ to next inversion sequentially in the meantime.

4.3. Model multi-reconstruction In the following, we will use the model in Fig. 2 as an example to discuss in detail how to estimate the reconstruction depth by the model difference. After the first inversion run, we have inversed model in Fig. 2c, which is called meta model. And according to Eq. (5), we get model difference Δmi from the last two iterations and shown in Fig. 2d, from which, the reconstruction depth is determined as plotted by straight line deeper from which model difference almost disappears. Then the new reference model could be reconstructed as illustrated in Fig. 2e–f. In Fig. 2e, a homogeneous half-space, which is called auxiliary model, will be cut down from the reconstruction depth, and replace the portion in the meta model above the reconstruction depth, creating a hybrid new model in Fig. 2f, which will be used as reference model for next run of inversion. We call this procedure model reconstruction. The first reference model can be obtained according to geological information. And during later reconstructions, the auxiliary model stays unchanged. It is also deserved to be pointed out that reconstruction depth may be different for each reconstruction. In order to get more accurate inversion results, above reconstruction could be done several times, and we give this procedure a name ‘multireconstruction’, which could be described mathematically as following. Let mi be the ith inversion model, m1aux is the auxiliary model, it should be noted that the auxiliary models are the same one during the whole procedure, and used as reference model at first inversion. Then the reconstructed model, which is also called reference model of the (i + 1)th inversion mi+1 ref , can be expressed as.
 i 1 miþ1 ref ¼ Ψi m ; maux

The key issue of MRT is how to estimate reconstruction depth. Before continuing, we assume the ith inversion run terminates after j time iterations, and denote the inversed model as mij. Theoretically, at the end of a single inversion, model perturbation corresponding to a reference model should be stable which means the model difference between adjacent iterations should disappear. However, after numerical MMR inversion tests, we find that even a large iteration number is set, MMR inversion always terminates in advance at insufficient convergence due to the regularization parameter smaller than predefined threshold. At this point, we could say that the inversion tries hard to match data by reducing the effect of model constraint on inversion with small regularization parameter; and therefore there is a i small difference between mij and mj−1 , especially in the shallow portion of inversed model, and this model difference is denoted as Δmi, which could be written as ð5Þ

From the view point of inversion, if Δmi is not small enough, it means that there is a possibility to modify the shallow part of the model further after the run of ith inversion. Accordingly, in MMR

ð6Þ

where


i

Ψi m

4.2. Reconstruction depth determination

Δmi ¼ mij −mij−1 :

inversion, we define reconstruction depth as the depth deeper than which model difference between the last two iterations is zero. On the contrary, if Δmi closes to 0, then the inversed models of last two iterations are almost the same, thus more inversion will not refine inversed model significantly, and the inversion catches it real minimal, we terminate the whole inversion.

; m1aux



" ¼

m1auxðabhÞ miðabðc−hÞÞ

# ð7Þ c

Eq. (7) is the mapping function of MRT, and a, b and c are the number of grids in x, y and z directions, respectively, h denotes the reconstruction depth. As the reconstruction goes on, the inversion model will be better, and model difference Δmi will approach 0 (Ye et al., 2013). Through multiple reconstructions, the deep information of previous model could be automatically transferred into reference model. Multireconstruction will continue until encounter a stop criterion. A natural one is that when Δmi is almost 0 or alternatively, when the difference between the Δmi from two adjacent inversions is less than a predefined threshold ε, that is when the reconstruction stops. jΔmiþ1 −Δmi jbε

ð8Þ

The procedure is very computational overload in 3D EM multi-frequency and multi-source inversion (Tang et al., 2007). Fortunately, in MMR method, only one frequency and one source is used, therefore time consuming is acceptable, providing an important prerequisite for multi-reconstruction inversion.

Y. Yang et al. / Journal of Applied Geophysics 159 (2018) 193–203

197

Fig. 2. Sketch map of how to reconstruct a new reference model by MRT after first inversion. Plan view (a) and vertical section in xoz plan (b) of 3D anomaly model; (c) section in xoz plan of the first inversion result, used as meta model; (d) model difference Δmi and the reconstruction depth represented by solid line below which Δm = 0; (e) auxiliary model, the reference model of the first inversion run, and often half space. Auxiliary model is same during the whole MRT procedure; (f) reconstructed model, a hybrid model with the portion shallower than reconstruction depth comes from auxiliary model in (e), and the deeper is from meta model in (c).

5. Synthetic examples 5.1. Single anomaly Fig. 3 illustrates the MRT inversion results of a simple model in Fig. 2. In this example, we reconstruct the reference model seven times and inverse eight times. In the Fig. 3, we use the ‘first inversion’ to represent the inversion result of the homogeneous half-space as the reference model, and ‘reconstruction I’ to represent the inversion result of the first reconstructed model as the reference model. A homogeneous half-space of 400 Ω·m is used as the reference model of the first inversion and the auxiliary model for the next seven reconstructions. To

compare with the inversion models after each reconstruction, the same termination criteria of inversion are accepted. Each inversion takes 120 iterations and consumes 120 min. Fig. 3a–h are the yoz-sections of the each inversion results. Fig. 3i is the xoz-section of the last inversion result. It can be seen that after seven times reconstructions, the abnormal resistivity decreased from 150 Ω·m in the first inversion (Fig. 3a) to about 50 Ω·m (Fig. 3h). It also could be seen clearly, the first inversion result, the abnormal body is at the surface. But after seven reconstructions, the abnormal body is almost the same as the truth model. Because that only y component of magnetic field is inversed, the model resolution in the x direction is not as good as that in the y direction (Fig. 3h). But the size of inversed abnormal shrinks into the

198

Y. Yang et al. / Journal of Applied Geophysics 159 (2018) 193–203

Fig. 3. Inversion results of the first inversion and each reconstruction, the black square box is the position of abnormal body. The ‘first inversion’ represents the inversion result as the reference model is the homogeneous half-space, and ‘reconstruction I' represents the inversion result of the first reconstructed model as the reference model. (a) first inversion, yoz section; (b) reconstruction I,yoz section; (c) reconstruction II,yoz section; (d) reconstruction III, yoz section; (e) reconstruction IV, yoz section; (f) reconstruction V, yoz section; (g) reconstruction VI, yoz section; (h) reconstruction VII, yoz section; (i) reconstruction VII, xoz section.

frame of the actual anomaly. Anyway, the resultant model proves that non-linear conjugate gradient (NLCG) inversion method incorporated with MRT is effective. Reconstruction depth is determined by the model difference between the last two iterations of the inversion, thus it will change after each inversion. Fig. 4 shows how to estimate reconstruction depth using model difference. In brown region, the model difference is nonezero. At the first reconstruction, the reconstruction depth is 100 m (Fig. 4a). The Fig. 4b–d, the brown regions extend to about 120 m, therefore, the reconstruction depth in red line is set 120 m which is deeper than the actual depth of the anomaly. As reconstruction inversion continues, the model differences become smaller, so do the brown regions. The region where the model difference is none-zero shrinks to 100 m (Fig. 4e–f) after 4th reconstruction inversions and 80 m after another two reconstructions (Fig. 4g). After 7th reconstruction inversions, the model differences disappear thus the reconstruction terminates. From

the example, we could see that the reconstruction depth may not be a fixed value, and the best depth should be determined at each reconstruction. 5.2. Stacked anomalies Previous example shows that MRT could enhance the depth resolution of MMR method by pressing shallow abnormal to deep ground. A natural question is that will it lose the shallow resolution? In order to figure it out, we design a stacked model to do further test. There are two near-surface conductive bodies as shown in Fig. 5. The observation mode, frequency, current, grid and noise level settings are similar to those in Fig. 2. We reconstruct the reference models seven times and inverse eight times. The reference model is the homogeneous half-space of 400Ω·m. In order to clearly show the inversion results, we demonstrate results in the slices. As the model is asymmetric, we intercept

Y. Yang et al. / Journal of Applied Geophysics 159 (2018) 193–203

199

Fig. 4. Model differences of the first inversion and seven time thereafter reconstructions in yoz sections (x = 0 m). h is the reconstruction depth. The model shallower than reconstruction depth will be replaced by homogeneous layer. At last reconstruction, in (h), the model difference of the inversion disappeared. (a) first inversion; (b) reconstruction I; (c) reconstruction II; (d) reconstruction III; (e) reconstruction IV; (f) reconstruction V; (g) reconstruction VI; (h) reconstruction VII.

Fig. 5. Anomaly model with two shallow conductors over a deep large abnormal body in half space. AA′ and BB′ are the location of vertical section (yoz) in the Fig. 6. (a) plan view of the model; (b) vertical section of model along line source.

200

Y. Yang et al. / Journal of Applied Geophysics 159 (2018) 193–203

Fig. 6. Inversion results of the first inversion and last reconstruction inversion (construction depth h = 120 m). yoz vertical section at AA′ (a), vertical section yoz at BB′ (b) and xoz vertical section at y = 0 m (c) of the inversed model after first inversion. (d), (e) and (f) are the corresponding vertical sections after reconstruction VII.

the vertical section (yoz) through AA′ and BB′ in Fig. 5. The line AA′ and line BB′ are the lines through the center of conductor A and B, respectively. To be compact, we only show the results of the first inversion and the 7th reconstruction inversion in Fig. 6. The model sections of the first inversion in Fig. 6a–c show that NLCG inversion is able to recover only the near-surface bodies, no sign of deep abnormal body could be found. From the Fig. 6c, it can be seen that the conductor A and B near the surface can be distinguished in the inversion result. Because of the superposition of shallow conductor and deep conductor, the conductor B has a lower resistivity and it is deeper than the true

model, its abnormal shape is more obvious than that of the conductor A. After seven times reconstructions, in Fig. 6d–e, we can distinguish both the near-surface abnormal bodies and deep body clearly, while the two near-surface conductors can be distinguished in the vertical section (xoz) (Fig. 6e). And comparing with the actual model, the inversion result reveals the true model, thus proves that MRT can not only enhance deep model resolution, but also retain shallow model information. Because of the length of an article, only the model differences of the first and the 7th reconstruction inversion are shown in Fig. 7 (vertical section (yoz) at BB′). Due to the existence of shallow conductors, the

Fig. 7. The yoz vertical section through BB′ in model differences of the first inversion (a) and reconstruction VII, or the last inversion (b). (c) The difference defined by Eq. (8) between the reconstruction VI and reconstruction VII, which is a termination condition of the reconstruction. Here, in this example, reconstruction depth is fixed at h = 120 m.

Y. Yang et al. / Journal of Applied Geophysics 159 (2018) 193–203

Fig. 8. MMR response (in percent of background field) over the Mons Cupri deposit, from Chen et al. (2002).

model difference changes little in the deep, so the reconstruction depth remains unchanged. The bottom of none-zero region of the model difference is 120 m, so the reconstruction depth is fixed to 120 m in each reconstruction. After 7th reconstruction inversion, the difference of the model difference between the 6th reconstruction inversion and the 7th reconstruction inversion is less than the threshold (Fig. 7c), so the reconstruction stops. 6. Field data example In order to verify the effectiveness of the proposed algorithm, the 3D ground MMR data from the Mons Cupri mine in Australia is inversed. The typical Cu-Zn-Pb deposit extends from 20 m to 200 m (Chen et al., 2002; Huston, 2006). Fig. 8 shows the MMR measured data in the contours. It should be noted that these measured data are the so-called ‘MMR anomaly’, which is the anomaly magnetic field normalized by the normal field at the center of the survey area (Chen et al., 2002; colton, etc., 2010), i.e. MMR anomaly ¼

Ba B0n

 100%

ð9Þ

In the measurement, the magnetic field can be obtained directly. So we can inverse the y-component of magnetic field at 0.3 Hz from grounded current line directly. The model to be inversed consists of

201

22,500 (30 × 30 × 25) grids. The adaptive regularization strategy is used during the inversion (Commer and Newman, 2008; Weng et al., 2012). Fig. 9 depicts the inversion results after the first inversion and that after the last reconstruction. Fig. 9a and b are the inversion results using a 400 Ω·m homogeneous half-space as reference model, which will be used as auxiliary model to reconstruct a new model for the following inversion. “A” in the figure represents the conductive target mineralization zone. Fig. 9a is the vertical section at x = 1150 m through the conductive target. Fig. 9b is the vertical section at y = 950 m through the conductive target. Fig. 9c and d are the result after three reconstructions. The reconstruction depths of three reconstructions are 60 m, 60 m and 30 m, respectively. From these figures it can be seen that there is an obvious difference before and after reconstruction, i.e. the shallow abnormal which float at the top of surface are pressed down to deep section, coincident to the positions of known geological settings (Chen et al., 2002; Huston, 2006). It also can be seen that the resistivity of abnormal bodies are about 40 Ω·m to 90 Ω·m, and the resistivity difference between abnormal bodies and background become clearer. The effect is not too obvious compared with the synthetic data, it is because that the depth of the conductive target is shallow in this example. 7. Discussions and conclusions We have developed a reference model construction method, and proposed a multi-construction strategy to improve inversion resolution. According to the test of the synthetic data and field data, the following points are the most significant: 1) The alternating current (AC) source, instead of direct current (DC) source, can be used as a source in MMR forward modeling. NLCG method is applied to inverse MMR magnetic field data. 2) To improve deep MMR inversion resolution, we developed a reference model multi-reconstruction strategy. Both the synthetic data and field data inversion results confirm the capability of reconstruction strategy. 3) The reconstruction method or its variants in the paper could also be applied to other kind EM data inversions. Although the method in this paper is able to obtain resistivity instead of resistivity ratio, and can as well improve deep resolution, it

Fig. 9. Inversion results of the first inversion and the last reconstruction inversion. (a) vertical section(yoz) at x = 1150 m and (b) the vertical section (xoz) at y = 950 m of the first inversion; (c) and (d) are the corresponding vertical sections of reconstruction III.

202

Y. Yang et al. / Journal of Applied Geophysics 159 (2018) 193–203

also has some disadvantages. One is the estimation of reconstruction depth. Optimal reconstruction depth is suggested to be at the bottom of the region where the model difference is none-zero, but because that the region is irregular; the depth is difficult to be determined sometimes. Thus, experiences are needed to determine the depth, leading to a difficulty to achieve proceduralization. Another problem is the long inversion time. Comparing with MT and CSAMT, a single MMR inversion time is relatively short. But when trying to improve deep resolution, several times of reconstruction inversion should be conducted, thus the computational efficiency is still needed to be enhanced. In addition, this paper accepts homogeneous half-space as auxiliary model (the reference model of the first inversion). When geological information is available, an auxiliary model can also be established from the geological model.

This paper is financially partly supported by the National Key Scientific Instrument and Equipment Development Projects of China (NO. 2011YQ05006010). The authors wish to thank the two anonymous reviewers whose suggestions and comments are very important to improve the paper. Appendix A. Relationship between magnetic field and resistivity The relationship between magnetic field and resistivity are derived from Maxwell equations of static electric field and the frequency domain, respectively. Maxwell equations for static magnetic field from are ∇  E ¼ 0;

ðA1aÞ

∇  H−σE ¼ J se ;

ðA1bÞ

∇∙ðμH Þ ¼ 0

ðA1cÞ

where E represents electric field strength; H denotes magnetic field strength; Jse is the current density, with σ the electrical conductivity and μ the magnetic permeability. From Eq. (A1a), we define a scalar potential ϕ that holds ðA2Þ

Taking divergence to Eq. (A1b), for static electric field, we can obtain Poisson equation expressed as ∇  ðσ∇ϕÞ ¼ ∇  J se

ðA3Þ

Substituting Eq. (A2) further into Eq. (A1b) yields the equation that H satisfied ∇  H ¼ J se −σ ∇ϕ

ðA4Þ

Eq. (A4) shows that if we keep the multiplication of σ∇φ unchanged, then H will not change as well. Assuming ϕ1 is the potential corresponding to σ1, then for another conductivity model, whose σ2(xk, yk, zk) = χkσ1(xk, yk, zk), where χk is a constant, from Eq. (A3), we could easily get ϕ2 = ϕ1/χk, leading to σ1∇ϕ1 = σ2∇ϕ2. Therefore, in static field, even if we change conductivity by a factor, the source term (i.e Jse − σ∇ϕ) at right-hand side in Eq. (A4) does not change, consequently, the magnetic field strength determined by Eq. (A4) does not depend on absolute conductivity. However, in FD case, by introducing impedance ^z and admittance ^, frequency domain Maxwell equations could be written as rate y ∇  E ¼ −^zH−J sm ;

ðA5aÞ

ðA5bÞ

^ ¼ σ þ iεω, ε represents the dielectric constant, while ω here, ^z ¼ iεω,y represents the angular frequency. Jsm represents the magnetic flux density, which was introduced by Nabighian (1987) for using Schelkunoff potential to calculate Maxwell equations in frequency domain. ^ is represented by σ and In order to compare Eq. (A5b) with (A1), y iεω in Eq. (A5b), so the Eq. (A5b) can be written as ∇  H−σE ¼ J se þ ^zE

ðA5cÞ

The relationship between the magnetic field and determinate conductivity can be obtained by solving Eq. (A5a) and (A5c): 1 ^ F þ ∇ð∇∙F Þ þ ∇  A: H ¼ −y ^z

Acknowledgments

E ¼ −∇ϕ

^E ∇  H ¼ J se þ y

ðA6Þ

^ ¼ σ þ iεω, F and A are the vector potential functions of electric here, y and magnetic fields, respectively. References Acosta, J.E., Worthington, M.H., 1983. A borehole magnetometric resistivity experiment. Geophys.Prospect. 31:800–809. https://doi.org/10.1111/j.1365-2478.1983.tb01087.x. Alumbaugh, D.L., Newman, G.A., Prevost, L., et al., 1996. Three-dimensional wide band electromagnetic modeling on massively parallel computers. Radio Sci. 31 (1), 1–23. Asten, M.W., 1988. The downhole magnetometric resistivity(DHMMR) method. Explor. Geophys. 19 (2), 12–16. Bishop, J., Carroll, N., Asten, M., et al., 1997. Finding sphalerite at Broken Hill with drillhole magnetometric resistivity. Explor. Geophys. 28, 6–10. Bouchedda, A., Giroux, B., Glencore, M.A., 2015. Down-hole Magnetometric Resistivity inversion for zinc and lead lenses localization at Tobermalug, County Limerick, Irland. 2015 SEG New Orleans Annual Meeting, 2007–2011 https://doi.org/10.1190/ segam2015-5849363.1. Campanyà, J., Ogaya, X., Jones, A.G., Rath, V., Vozar, J., Meqbel, N., 2016. The advantages of complementing mt profiles in 3-d environments with geomagnetic transfer function and interstation horizontal magnetic transfer function data: results from a synthetic case study. Geophys. J. Int. 207. Cattach, M.K., Stanley, J.M., Lee, S.J., et al., 1993. Sub-Audio Magnetics (SAM)—a high resolution technique for simultaneously mapping electrical and magnetic properties. Explor. Geophys. 24 (4), 387–400. Chaladgarn, T., Yooyuanyong, S., 2013. Mathematical model of magnetometric resistivity sounding for a conductive host with a bulge overburden. Appl. Math. Sci. 7 (7), 335–348. Chen, J., Oldenburg, D.W., 2003. 3-D inversion of magnetic induced polarization data. Aseg Extended Abstr. Adelaide 2003 (1):245–253. https://doi.org/10.1071/EG06245. Chen, J., Oldenburg, D.W., 2004. Magnetic and electric fields of direct currents in a layered earth (short note). Explor. Geophys. 35 (2), 157–163. Chen, J., Haber, E., Oldenburg, D.W., 2002. Three-dimensional numerical modelling and inversion of magnetometric resistivity data. Geophys. J. Int. 149 (3):679–697. https://doi.org/10.1046/j.1365–246X.2002.01688.x. Commer, M., Newman, G.A., 2008. New advances in three-dimensional controlled-source electromagnetic inversion. Geophys. J. Int. 172 (2), 513–535. Deng, J.W., 2005. Forward Theory Study of Magnetometric Resistivity and Induced Polarization Methods. Dong, H., Wei, W.B., Ye, G.F., et al., 2014. Study of Three-dimensional magnetotelluric inversion including surface topography based on Finite-difference method. Chin. J. Geophys. 57 (3):939–952. https://doi.org/10.6038/cjg20140323 (in Chinese). Edwards, R.N., 1974. The Magnetometric Resistivity Method and its Application to the Mapping of a Fault. Can.J.Earth Sci. 11:1136–1156. https://doi.org/10.1139/e74-108. Edwards, R.N., Howell, E.C., 1976. A field test of the Magnetometric Resistivity (MMR) method. Geophysics 41 (6A(6)), 1170–1183. Edwards, R.N., Nabighian, M.N., 1991. The magnetometric resistivity method. Electromagnetic Methods in Applied Geophysics:Investigations in Geophysics, Part A. vol. 2. Society of Exploration Geophysicists., pp. 47–104. Edwards, R.N., Lee, H., Nabighian, M.N., 1978. On the theory of magnetometric resistivity (MMR) methods. Geophysics 43 (6):1176–1203. https://doi.org/10.1190/1.1440887. Egbert, G.D., Kelbert, A., 2012. Computational recipes for electromagnetic inverse problems. Geophys. J. Int. 189 (1):251–267. https://doi.org/10.1111/j.1365246X.2011.05347.x. Fathianpour, N., Cattach, M.K., 1995. Analytical solutions for the Total Field Magnetometric Resistivity (TFMMR) technique. Explor. Geophys. 26 (3):158–166. https://doi.org/ 10.1071/EG995158. Fathianpour, N., Heinson, G., White, A., 2005a. The total field magnetometric resistivity (TFMMR) method: Part I: theory and 2.5D forward modelling. Explor. Geophys. 36 (2), 181–188. Fathianpour, N., Heinson, G., White, A., 2005b. The total field magnetometric resistivity (TFMMR) method: Part II: 2D resistivity inversion of data from the Flying Doctor Deposit, Broken Hill, Australia. Explor. Geophys. 36 (2), 189–197.

Y. Yang et al. / Journal of Applied Geophysics 159 (2018) 193–203 Godber, K.E., Bishop, J.R., 2007. DHMMR:Coming of age. Exploration in the new millennium. Proceeding of the Fifth Decennial International Conference on Mineral Exploration, pp. 1119–1123. Howland-Rose, A.W., Linford, G., Pitcher, D.H., et al., 1980a. Some recent magnetic induced-polarization developments—Part I: Theory. Geophysics 45 (1), 37–54. Howland-Rose, A.W., Linford, G., Pitcher, D.H., et al., 1980b. Some recent magnetic induced-polarization developments- Part II Survey results. Geophysics 45 (1), 55–74. Hu, Z.Z., Hu, X.Y., He, Z.X., 2006. Pseudo-Three dimensional magetotelluric inversion using nonlinear conjugate gradients. Chin. J. Geophys. 49 (4), 1226–1234 (in Chinese). Huston, D.L., 2006. Mineralization and regional alteration at the Mons Cupri stratiform Cu-Zn-Pb deposit, Pilbara Craton, Western Australia. Mineral. Deposita 41 (1): 17–32. https://doi.org/10.1007/s00126-005-0036-4. Kelbert, A., Egbert, G.D., Schultz, A., 2008. Non-Linear Conjugate Gradient Inversion for Global EM Induction: Resolution Studies. Geophys. J. Int. 173 (2), 365–381. Kofoed, V.O., Jessop, M.L., Wallace, M.J., et al., 2011. Unique applications of MMR to track preferential groundwater flow paths in dams, mines, environmental sites, and leach fields. Leading Edge 30 (2), 192–204. Kong Z Z, Shan K S, Wu Y, et al. 2015. The improvement and applications of the impression method in AMT data processing. Geophys. Geochem. Explor.. (in Chinese), 39 (2): 416–420, doi:org/10.11720/wtyht.2015.2.34. Labrecque, D., Sharpe, R., Dan, C., et al., 2003. Combined electrical and magnetic resistivity tomography: synthetic model study and inverse modeling. J. Environ. Eng. Geophys. 8 (4), 251–262. Li, Y., Oldenburg, D.W., 1996. 3-D inversion of magnetic data. Geophysics 61 (2):394–408. https://doi.org/10.1190/1.1443968. Li, Y., Oldenburg, D.W., 1998. 3-D inversion of gravity data. Geophysics 63 (1):109–119. https://doi.org/10.1190/1.1444302. Li, J.H., Lin, P.R., Guo, P., 2010. A trial study of magnetic induced polarization. Seismol. Geol. 32 (3):492–499. https://doi.org/10.3969/j.issn.0253-4967.2010.03.016 (in Chinese). Lin, C.H., Tan, H.D., Shu, Q., et al., 2012. Three-dimensional conjugate gradient inversion of CSAMT data. Chin. J. Geophys. 55 (11):3829–3839. https://doi.org/10.6038/j. issn.0001-5733.2012.11.030 (in Chinese). Liu, Y.H., Yin, C.C., 2013. 3D inversion for frequency-domain HEM data. Chin. J. Geophys. 56 (12):4278–4287. https://doi.org/10.6038/cjg20131230 (in Chinese). Maxwell, M., Eso, R., Oldenburg, D., 2013. Using 2D and 3D electrical resistivity and magnetometric resistivity techniques for investigating dam and dike soil conditions for leak detection — field examples and forward modelling. Symposium on the Application of Geophysics to Engineering and Environmental Problems 2013. Society of Exploration Geophysicists and Environment and Engineering Geophysical Society: p. 711 https://doi.org/10.4133/sageep2013-264.1. Meju, M.A., 2002. Geoelectromagnetic exploration for natural resources: models, case studies and challenges. Surv. Geophys. 23 (2–3), 133–206. Nabighian, M.N., 1987. Electromagnetic Methods in Applied Geophysics. vol. 1. SEG, Tulsa, Oklahoma. Newman and Alumbaugh, 2000. Three-dimensional magnetotelluric inversion using nonlinear conjugate gradients. Geophys. J. Int. 140:410–424. https://doi.org/10.1046/ j.1365-246x.2000.00007.x.

203

Nobes, D.C., Law, L.K., Edwards, R.N., 1986. The determination of resistivity and porosity of the sediment and fractured basalt layers near the Juan de Fuca Ridge. Geophys. J. R. Astron. Soc. 86 (2), 289–317. Portniaguine, O., Zhdanov, M.S., 2000. 3-D magnetic inversion with data compression and image focusing. Geophysics 67 (5), 1532–1541. Seama, N., Tada, N., Goto, T.N., et al., 2013. A continuously towed vertical bipole source for marine magnetometric resistivity surveying. Earth Planets Space 65 (8):883–891. https://doi.org/10.5047/eps.2013.03.007. Seigel, H.O., 1974. The magnetic induced polarization (MIP) method. Geophysics 39 (3), 321–339. Sobouti, A., Motagh, M., Sharifi, M.A., 2015. Inversion of surface gravity data for 3-D density modeling of geologic structures using total variation regularization. Stud. Geophys. Geod. 60 (1):69–90. https://doi.org/10.1007/s11200-014-0671-2. Szarka, L., 1987. Geophysical mapping by stationary electric and magnetic field components: a combination of potential gradient mapping and magnetometric resistivity methods. Geophys. Prosp. 35, 424–444. Tang, J.T., Ren, Z.Y., Hua, X.R., 2007. The forward modeling and inversion in geophysical electromagnetic field. Progr. Geophys. 22 (4), 1181–1194 (in Chinese). Weng, A.H., Liu, Y.H., Gao, L.J., et al., 2009. Magnetoelectric resistivity reponse characteristics for natural gas hydrate. OGP 44 (Supplement1):158–161. https://doi.org/ 10.13810/j.cnki.issn.1000-7210.2009.sl.037 (in Chinese). Weng, A.H., Liu, Y.H., Jia, D.Y., et al., 2012. Three-dimensional controlled source electromagnetic inversion using non-linear conjugate gradients. Chin. J. Geophys. 55 (10): 3506–3515. https://doi.org/10.6038/j.issn.00053-5733.2012.10.034 (in Chinese). Weng, A.H., Li, D.J., Li, Y.B., et al., 2015. Selection of parameter types in Controlled Source Electromagnetic Method. Chin. J. Geophys. 58 (2):697–708. https://doi.org/10.6038/ cig20150230 (in Chinese). Yang, J., 2005. Geo-electrical responses associated with hydrothermal fluid circulation in oceanic crust: feasibility of magnetometric and electrical resistivity methods in mapping off-axis convection cells. Explor. Geophys. 36 (3), 281–286. Yang, C.H., Tseng, H.W., 1992. Topographic responses in magnetometric resistivity modeling. Geophysics 57 (10), 1409–1478. Ye, T., Chen, X.B., Yan, L.J., 2013. Refined techniques for data processing and two-dimensional inversion in magnetotelluric(III):using the Impressing Method to construct starting model of 2D magnetotelluric inversion. Chin. J. Geophys. 56 (10): 3596–3606. https://doi.org/10.6038/cjg20131034 (in Chinese). Yee, K.S., 1966. Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media. IEEE Trans. Ant. Prop. 14 (3), 302–307. Zhang, K., Dong, H., Yan, J.Y., et al., 2013. A NLCG inversion method of magnetotellurics with parallel structure. Chin. J. Geophys. 56 (11), 3922–3931 (in Chinese). (doi: 10.6038.cjg20131134). Zhdanov, M.S., 2002. Geophysical Inverse Theory and Regularization Problems. Elsevier Academic Press, Amsterdam. Zhu, K., Yang, J., 2008. Time-dependent magnetometric resistivity anomalies of groundwater contaminatin:Synthetic results from computational hydro-geophysical modeling. Appl. Geophys. 5 (4):322–330. https://doi.org/10.1007/s11770-008-0041-3.