- Email: [email protected]
Inversion of arbitrary segmented loop source TEM data over a layered earth
Journal of Applied Geophysics 128 (2016) 87–95
Contents lists available at ScienceDirect
Journal of Applied Geophysics journal homepage: www.elsevier.com/locate/jappgeo
Inversion of arbitrary segmented loop source TEM data over a layered earth Hai Li a,b,⁎, Guo-qiang Xue a, pan Zhao a,b, Nan-nan Zhou a, Hua-sen Zhong a,b a b
Key Laboratory of Mineral Resources, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing, 100029, China University of Chinese Academy of Sciences, China
a r t i c l e
i n f o
Article history: Received 29 October 2015 Received in revised form 17 February 2016 Accepted 23 March 2016 Available online 29 March 2016 Keywords: OCCAM inversion Loop source TEM method Mined out area Segmented loop Loop geometry
a b s t r a c t The loop source TEM method has been widely used in the detection of a mined out area in China. In the cases that the laying of traditional rectangle or square transmitting loop is limited due to the presence of obstacle on the path of the loop, the changing of the shape of the transmitting loop to bypass the obstacle is a labor saving solution. A numeric integration scheme is proposed to calculate the response and Jacobian of the segmented loop source from that of an electric dipole source. The comparison of forward response between the segmented loop and square loop shows the effect of loop geometry on the decay curves. In order to interpret the data from an irregular source loop, this paper presents an inversion scheme that incorporate the effect of loop geometry. The proposed inversion scheme is validated on the synthetic data, and then applied to the field data. The result reveals that the developed inversion scheme is capable of interpreting the segmented loop source TEM field data. © 2016 Elsevier B.V. All rights reserved.
1. Introduction The time domain electromagnetic method (TDEM of TEM) is one of the major exploration tools for electromagnetic (EM) methods. The loop source TEM method which is highly efficient and gives an excellent resolution of conductive layer at depth, has been widely used in hydrogeophysical investigation and mineral exploration (Danielsen et al., 2003; Fitterman and Stewart, 1986; Xue et al., 2012b; Zhou and Xue, 2014). The theory and interpretation techniques of the TEM method have drawn a dramatic surge of interest since the early 1950s. Commonly, TEM data can be interpreted by using fast imaging techniques (Tartaras et al., 2000) or 1-D inversion base on fitting the theoretical decay curves from a 1-D structure to the real data in the least square sense (Farquharson and Oldenburg, 1996). The interpretation skills mentioned above are generally based on the square or rectangle transmitting loop. For the airborne TEM method which covers a large portion of TEM surveys, the source loop is wound around the aircraft for fixed-wing transient EM system (Smith et al., 2003) or mounted on a lightweight wooden lattice frame for the helicopter electromagnetic system (Sørensen and Auken, 2004; Smith et al., 2009; Witherly et al., 2004). The geometry of transmitting loop is almost stationary during the data collection procedure so that there is need for considering the variation of the shape of the source loop in process of data interpretation. However, for the land cases, in which ⁎ Corresponding author at: Institute of Geology and Geophysics, Chinese Academy of Sciences, Beitucheng west road, Chaoyang district, Beijing, China. E-mail address: [email protected] (H. Li).
http://dx.doi.org/10.1016/j.jappgeo.2016.03.017 0926-9851/© 2016 Elsevier B.V. All rights reserved.
the laying of the square loop or rectangle loop source is restricted, such as when there exists a scarp, pool, river or dwellings on the path of the transmitting line, we have to redesign the path of the transmitting loop, which is not easy to implement and labor intensive. Otherwise, it is often the case that we simply change the shape of the transmitting loop to bypass these obstacles, which is labor-saving. The adoption of existing interpretation techniques to directly processing these data is likely to get an unreliable result. Zhou and Xue (2012) compared the response of rectangular loop and the circular loop at different location inside the loop and pointed out that the responses are affected by the loop geometry at almost every position inside the loop. Qi et al. (2014) proved that the neglect of the loop geometry will affect the recovered resistivity model. However, the interpretation techniques considering the geometry of the source loop are hardly seen in the literature. In this paper, we propose a 1D inversion scheme that incorporated the geometry of the transmitting loop. Although several papers have been published on the effect of 3D structure on the loop source TEM method (Goldman et al., 1994; Newman et al., 1987; Rabinovich, 1995), the 1D inversion is also commonly used in literatures to invert field data (Danielsen et al., 2003; Nekut, 1987; Tezkan and Siemon, 2014). Besides, the full 3D subsurface models can still be obtained from datasets which consist of many individual sounding by using techniques, such as laterally constrained 1D-inversion (Auken et al., 2004) or spatially constrained inversion 1D inversion(Viezzoli et al., 2008). Indeed, it is limited to apply the 1D inversion to invert multidimensional structures, and we do not expect the recovery of multidimensional resistivity structure with 1D interpretation techniques. Instead, our goal is to develop a 1D inversion scheme especially for the segmented loop
88
H. Li et al. / Journal of Applied Geophysics 128 (2016) 87–95
source TEM data and appraise its feasibility. We first analyze the effect of loop geometry by considering two designed segmented transmitting loop by comparing their responses with the square loop responses. In a further step, the developed inversion scheme is validated on synthetic data, and then applied to field data. A case study shows that the developed inversion scheme is capable of interpreting the field data excited by segmented loop source. 2. Basic theory This section gives the methodology of forward and inverse modeling used for the inversion of arbitrary segmented loop source TEM data. There is no analytical solution for the responses of a loop source with arbitrary shape. The responses presented in this paper are calculated using numerical integration. The integration method and reciprocity theorem have been used in the literatures to calculate the responses excited by a loop source. Poddar (1982) calculated the responses of a rectangular loop using the concept of reciprocity and the known solution for the electric field responses of a vertical dipole source placed over a two-layered earth half space. Qi et al. (2014) extended this method and obtained the solutions for a horizontal loop of arbitrary shape over layered earth. Here, we propose an alternate method to directly calculate the responses for arbitrary segmented loop and derive the formulas for calculating the sensitivity without using the reciprocity theorem.
of the symmetrical property, the vector potential A has two components which are in the direction of the dipole and vertical direction respectively. To reduce Eq. (3) into an ordinary differential equation in the variable z, we apply the 2D Fourier transformation to the x-variable and yvariable of Eq. (3) and get the Helmholtz equation in the wavenumber domain. ̂ S ∂ A 2 2 2 − 2 þ kx þ ky −k A ¼ Ĵ e ∂z 2 ̂
ð4Þ
where kx and ky are the spatial wavenumber variables of the Fourier transform, and A represents the vector potential in the wave number domain. After Eq. (4) is solved, the solution for the vector potential in the spatial domain can be obtained by making use of the relation between the 2D Fourier transformation and the Hankel transformation. Aðx; y; z; ωÞ ¼
1 2π
Z 0
∞
2 2 Â kx þ ky ; z; λ J 0 ðλrÞdλ:
ð5Þ
Hence, the frequency domain EM response for arbitrary oriented dipole source can be obtained by substituting Eq. (5) into Eqs. (1) and (2). Next, we derive the formulations used to calculate the EM field due to arbitrary segmented loop source by integrating the responses of a horizontal electric dipole around the transmitting loop H ðωÞ ¼ ∮ F ðωÞ dl
ð6Þ
L
2.1. Calculation of forward responses Firstly, we consider a 1D model shown in Fig. 1, where the horizontal dipole source is located on the surface of the ground. The algorithm for the calculation of the responses for the horizontal dipole sources have been well studied in the literature. The one adopted in this paper follows the works of Ward and Hohmann (1988) and Key (2009). By using the Schelkunoff vector potential A, the magnetic field B and the electric field E can be specified as B¼∇A
ð1Þ
1 ∇ð∇ AÞ: E ¼ iωA þ μσ
ð2Þ
By assuming a time dependence of eiωt and neglecting the displacement current, the Helmholtz equation for the Schelkunoff potential in the spatial domain can be derived from the Maxwell equation and the constitutive relations.
where F(ω) and H(ω) represent the responses of an arbitrarily oriented dipole source and the arbitrary segmented loop source respectively. L represents the integration path which is along the transmitting loop. For the segmented loop considered here, the source loop can be discretized as H ðωÞ ¼
N Z X n¼1
2
ð3Þ
H ðωÞ ¼
N Z X
yi
Z F x ðωÞdx þ
yiþ1
yi
N Z X F y ðωÞdy ¼ n¼1
xi
xiþ1
F ðωÞ cosαdx
ð8Þ
F ðωÞ cosβ dy
where α and β donate the angle between the directions of the dipole and the x-axis or y-axis respectively. Hence, the EM field excited by arbitrary segmented loop source can be obtained with the integration equation. In the final step, the time domain responses are evaluated using the digital filtering techniques (Anderson, 1979) hðt Þ ¼ −
Fig. 1. The M-layered 1D model and the coordinate system.
xiþ1
xi yiþ1
Z þ
JSe = Iδ(r − r0)
represent the electric dipole where k = − iωμσ, and source imposed at position r0 with the source dipole moment I. Because
ð7Þ
N is the number of segments. The field excited by arbitrarily oriented dipoles can be computed using the vector superposition of the fields produced by the horizontal x-oriented and y-oriented dipoles (Key, 2009). Thus the integral of the arbitrarily oriented dipole source can be calculated using the integral of the responses of the horizontal xoriented and y-oriented dipoles,
n¼1
∇2 A þ k A ¼ − JSe
F ðωÞdl:
Li
2 π
Z 0
∞
HðωÞ Im cosðωtÞdω: ω
ð9Þ
In order to verify the effectiveness of our algorithm for the layered earth, the responses for a 40 × 40 m rectangle central loop configuration on the three layer model have been validated by comparing the results from the Leroi code which is developed by CSIRO and has been integrated on the commercial Maxwell software. The model used for comparison is the typical H-type model and K-type model. The models consist of 100 Ω·m half space with a 100 m thick layer whose resistivity is 10 Ω·m for the H-type model and 1000 Ω·m for the K-type model,
H. Li et al. / Journal of Applied Geophysics 128 (2016) 87–95
and is centered at a 150 m depth. Fig. 2 shows the comparison results which verify the correctness of our algorithm. 2.2. Calculation of sensitivities To conduct the inversion routine, we need to calculate the sensitivities which quantify the effects of the changes in the model parameters on the responses. The derivation of the sensitivities for the arbitrary segmented loop source is based on the sensitivities for the dipole source and the integration around the transmitting loop in the same way as the calculation of the forward responses shown in Eq. (8). For the layered model considered in this paper, the model vector can be described in terms of a linear combination of conductivity functions,
σ ðzÞ ¼
N X
2.3. The inversion algorithm A popular OCCAM 1D inversion algorithm is incorporate into our segmented loop inversion scheme (Constable et al., 1987). For the OCCAM inversion, the model is generally discretized into several layers with fixed thickness. The model vector can be written as m = (m1, …,mN)T. The forward operator and the Jacobian are computed based on the formula derived above. The basic idea behind OCCAM's inversion is an iteratively applied local linearization. The linearization of forward operator with the Taylor theorem leads to the objective function for the inversion 2 min WF mk þ WJ mk Δm−Wd −χ 20 2 2 2 k þ α L m þ Δm ; 2
σ j δ j ðzÞ
ð10Þ
j¼1
where N specifies the number of layers in the model, the function δj(z) is equal to unity in the layer j and zero in another layer. To calculate the sensitivities, we start by specifying the change in vector potential in the wavenumber domain due to the changes in the model parameter. Then the transformations introduced in Eqs. (5), (1) and (2) were adopted to get the sensitivities. Using the model vector defined in Eq. (10) and differentiating Eq. (4), we obtain the Helmholtz equation in the wavenumber domain 2
0 13 2 N X ∂ ∂A 2 2 4− þ @kx þ ky þ σ j δ j ðzÞA5 ¼ iωμ 0 δ j ðzÞA: ∂σ j ∂z2 j¼1
ð11Þ
Eq. (11) is an inhomogeneous ordinary differential equation for ∂A=∂ σ j with boundary condition ∂A=∂σ j →0, as z→ ±∞. Hence, the solution for the boundary value problem here can be obtained with the adjoint Green function method (Dettman, 2013). 0
89
∂A iωμ ðkx ; kx ; 0; ωÞ ¼ − Ids ∂σ j
Z
j
z¼ j−1
h
i2 A j ðkx ; kx ; 0; ωÞ dz
ð12Þ
After the values of ∂A=∂σ j have been obtained, the transformations which are used previously to get the forward response are adopted here to get the sensitivities for the arbitrary segmented loop source.
Fig. 2. The comparison of the responses computed with the algorithm here and the responses calculated with the Leroi routine.
ð12Þ
where F is the forward operator, J is the Jacobian, d is the data vector, χ0 is the target misfit value, and α is the Lagrange multiplier which is used to balance the misfit and the roughness term. L is the roughness matrix to implement the Tikhonov regularization. W is the data weighting matrix to incorporate the standard deviations
−1 −1 W ¼ diag σ −1 1 ; σ 2 ; ; σ M
ð13Þ
where the σi represents an estimation of the error in data. By minimizing the defined objective function (Eq. (12)), we obtain the smoothest resistivity model for a given target misfit. To minimize Eq. (12) with respect to the model vector mk, we rearrange Eq. (12) and let d mk ¼ d− F mk þ J mk mk :
ð14Þ
Then we use the Gauss–Newton method and find the iterative solution, mkþ1 ¼ mk þ Δm ¼
−1 T T J mk J mk þ α 2 LT L J mk d mk ð15Þ
At each iteration, we pick the largest α2 that keeps the misfit value of the solution from exceeding the boundary value χ0 specified in Eq. (12). At the end of the iteration procedure, we will have a solution whose misfit is equal to χ20. 3. The analysis of forward responses To study the effect of the shape of the transmitting loop on the responses and demonstrate that it is necessary to develop the inversion scheme for arbitrary segmented loop sources, we calculate and compare the responses excited by three different loop sources, as is shown in Fig. 3. The segmented loop source 1 shown in Fig. 3a, the segmented loop source 2 shown in Fig. 3b and the rectangle loop source considered as a reference in Fig. 3a and b. There are three types of loop source configurations which have been commonly used in engineering investigation and near surface geophysical exploration, including central loop configuration, fixed loop configuration and the modified central loop configuration (Xue et al., 2012a). Accordingly, we design four recording stations which are shown in Fig. 3 and calculated the time derivative of the vertical magnetic induction at each station to get our comparison results. The geological model considered here consists of 10 Ω·m half space with a 100 m thick and 100 Ω·m resistive layer centered at a depth of −130 m. Firstly, the responses excited by the square loop source and the segmented loop source 1 for the three layer model are calculated using the algorithm introduced previously. The decay curves at recording stations ①–④ are plotted in Fig. 4a–d respectively. In this case, the responses of the square loop source are larger than those of the segmented loop source 1. The recording stations ① and ③ are designed for the modified
90
H. Li et al. / Journal of Applied Geophysics 128 (2016) 87–95
Fig. 3. The geometry of the loop source used for the forward modeling. ①, ②, ③, and ④ denote the position of the recording station. (a) The segmented loop source 1 (consists of the solid line and the red solid line) and the rectangle loop source (consists of the solid line and the dash line); (b) the segmented loop source 2 (consists of the solid line and the red solid line) and the rectangle loop source (consists of the solid line and the dashed line). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Fig. 4. The comparison of responses between a conventional rectangular loop source and the segmented loop source 1. The decay curves in (a) to (d) correspond to the recording stations ① to ④ shown in Fig. 3.
H. Li et al. / Journal of Applied Geophysics 128 (2016) 87–95
91
Fig. 5. The comparison of responses between a conventional rectangular loop source and the segmented loop source 2. The decay curves in (a) to (d) correspond to the recording stations ① to ④ shown in Fig. 3.
central loop source configuration. The decay curves for the segmented loop source at these stations show a different character. The early time response at station ① is affected by the loop geometry while the response at station ③ appears largely unaffected. The station ② is designed for the central loop source, and the effect of the geometry on the response at this station is basically impregnable. The station ④ is designed for the fixed loop source configuration, as is shown in Fig. 4d, the transition period of the response appears affected while the other time range is slightly affected. Besides, the late time responses are slightly offset from the square loop responses due to the difference of the area of the transmitting loop. Similarly, the responses excited by the square loop source and the segmented loop source 2 are compared, as is shown in Fig. 5. Unlike in the case of the segmented loop source 1, the response of the segmented loop source is larger than that of the square loop. For the recording station located at the center of the loop, where the response is basically impregnable for the segmented loop source 1, the response is affected by the loop geometry. The early time response is deviated from the response of the square loop source. From the comparison of the forward response, we conclude that the responses excited by the loop source are influenced by the geometry of the transmitting loop. These effects may be introduced by the different transmitting moments, or irregular sides of the loop. For the recording station inside the loop, the deviation of the response from the square loop mainly appears at an early time, depending on the location of the
station. In the case of the recording station outside of the loop, the responses at the transition period, where the time derivatives of the magnetic induction change its sign, are most likely to be affected. Hence, the use of inversion scheme for square loop source TEM data to invert the arbitrary segmented loop source data may lead to an inappropriate result, which makes it necessary to develop the inversion scheme for the arbitrary segmented loop source. The differences of the calculated responses will have an effect on the inverted model. Generally, the root mean square (RMS) misfit is considered in the inversion scheme
χ RMS
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi " #ffi u i 2 n u1 X F ðmÞi1 −d2 t ; ¼ n i¼1 σ i1
ð16Þ
where n is the number of data and σi is the standard error for the measured data, m is the model vector, and F is the forward operator. For the modeled responses considered here, we adopt the relative Table 1 The relative error at each recording station between the segmented transmitting loop and square loop. Recording station
1–①
1–②
1–③
1–④
2–①
2–②
2–③
2–④
Relative error
0.99
0.51
0.43
1.55
1.47
1.22
0.65
2.00
92
H. Li et al. / Journal of Applied Geophysics 128 (2016) 87–95
Fig. 6. Inversion results for recording station ① of segmented loop source 1. (a) Synthetic data and their fitting for square loop source, (b) synthetic data and their fitting for segmented loop source and (c) the inverted model.
error to evaluate the difference introduced by the change of the geometry of the transmitting loop
4. 1D inversion of synthetic segmented loop source TEM data
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi " # u i i 2 n u1 X d1 −d2 : E¼t i n i¼1 d1
In this section, the inversion scheme is validated on synthetic segmented loop source TEM data. For comparison, these data are also inverted with the square loop inversion routine. A five layer model (resistivity: 100, 70, 30, 60, 20 Ω·m; thickness: 50, 75, 50, 75 m) is considered in the inversion process for it has a combination of resistivity model such as ρ1 N ρ2 N ρ3 b ρ4 b ρ5, the conventional Q- and K-type models (Bhattacharya, 2012; Tezkan and Siemon, 2014). Besides, the inversion result for the layer at different depth can be compared in this model. The synthetic data excited by the segmented loop source in Fig. 3 in which the recording stations ① and ② are considered. The Gaussian noise generated with the noise model proposed by (Munkholm and Auken, 1996) is added into the dataset. The added Gaussian noise whose standard deviation is 3% is log-gated and gated stacked, and the resulting log-gated noise gets different amplitudes at
ð17Þ
where di1 and di2 represent the modeled response of the segmented loop source and the square loop source respectively. Table 1 shows the results calculated by Eq. (17). The relative error introduced by the segmented loop 2 is larger for than that of the segmented loop 1. For the recording station considered, the responses at station ① and station ④ are more largely affected than that of station ② and station ③. The comparison of the relative error can give us some insight to the effect on the models recovered, which may be introduced by the loop geometry.
Fig. 7. Inversion result for the data collected at recording stations ① and ② of segmented loop source 2. (For interpretation of the reference to color in this figure legend, the reader is referred to the web version of this article.)
H. Li et al. / Journal of Applied Geophysics 128 (2016) 87–95
93
Fig. 8. Geographic location of the selected profile.
each time gate to fit the decaying properties of time-domain EM data. The misfit is calculated with Eq. (16), and the target misfit is 1.00. Firstly, we inverted the data excited by the segmented loop source shown in Fig. 3a, and the recording station is located at position ①. The synthetic data is inverted using both the inversion scheme of square loop source and arbitrary segmented loop source. Fig. 6a and b show synthetic data and their fitting for the segmented loop source inversion scheme and the square loop source inversion scheme respectively. Fig. 6c shows the inverted models of segmented loop source 1. We can see that the model results (square loop inversion and segmented loop
inversion) were not consistent for the shallow layers but consistent in the relatively deeper layer. The resistivity at each layer for this model is not well resolved, but depth of each layer is recovered quite well for the segmented loop inversion. Fig. 7 shows the models obtained by inverting the data excited by segmented loop source 2 at the recording stations ① and ②. The inverted models for the segmented loop inversion at these two stations are approximately the same, and the obtained resistivity values are well resolved. This illustrated that with the developed inversion scheme, we can get a reliable result from the segmented loop source data. As a contrast, when we use the square loop inversion scheme to invert the segmented loop source data, the deviation of results from the true model is non-negligible, especially for the shallow layer.
Fig. 9. The segmented transmitting loop at the beginning of the survey line. There are four recording stations in this loop. The coordinated system was built for the inversion.
Fig. 10. The decay curve collected at the station where profile coordinate is 40 m with the error bar of each data.
94
H. Li et al. / Journal of Applied Geophysics 128 (2016) 87–95
Fig. 11. The core from drill hole showing the peat at a depth of −57 m which is partly weathered.
From the inversion results shown in the examples above, we conclude that the square loop inversion scheme will lead to a result deviated from the true model when the geometry of the source loop is changed. For the considered cases, the deviation appears at a shallow depth. On the other hand, using the developed arbitrary segmented loop inversion scheme can well resolve the five layer model by including the loop geometry into the inversion algorithm. Therefore, it is necessary to consider the change of loop geometry during the processing of the loop source TEM data.
5. 1D inversion of segmented loop source TEM field data After the successful validation of the segmented loop source inversion on the synthetic data, it has been applied to field data. In a survey aimed at investigating the mined out area of an iron ore mine in Shan Xi province, China, we collected the field data excited by a segmented loop source data at the first few stations of the survey line. Fig. 8 shows the geographic location of the selected profile. The survey line is designed based on the survey and design of an inter-city railway. The mined-out area in this iron mine is caused by illegal excavation among folks, and is one of the major threats to the safety of the intercity railway. Due to the existence of a scarp at the beginning of the survey line, the shape of the first transmitting loop had to be changed as shown in Fig. 9. The first segmented loop was designed to test our algorithm, but at the same time, the geo-electrical information underground at these recording station is of interest to the mining engineer who is responsible for the mining process in this iron ore mine. The vertex of the transmitting loop was recorded by the GPS so that we can get the shape of the transmitting loop, as is shown in Fig. 9. We built the coordinate system to calculate and invert the response of the segmented transmitting loop which contains four recording stations. The data at the rest of the 800 m survey line was collected with the modified square loop source TEM configurations (Xue et al., 2012a),
and the transmitting loop size was 100 × 100 m. Both the first segmented loop source TEM data and the square loop source data are inverted using our inversion scheme. Hence, the inverted cross-section of the modified central loop TEM data can also be used to validate our result in return. The terraTEM time domain EM survey system designed by Monash GeoScope was used for the measurement. The system was powered by an external 24 V battery power package, and the source waveform was a bipolar rectangle current with 50% duty circle. The receiver was a vertical component magnetic core probe with 10,000 m2 effective area, which is the same as that introduced by Xue et al. (2012b). Fig. 10 shows the field data collected at the profile coordinated 40 m. The measurements are made from 10−4 s to 10−2 s. The error bar represents the standard deviation of each data point. It can be seen in Fig. 10 that the last few data points are submerged in the noise. It is illustrated previously that the data stand deviations have been incorporated into the solution so that we can still get an appropriate result at the presence of noise, which has been verified by the synthetic data test. The program of the survey and design of an inter-city railway carried out several drillings to get some geological knowledge of this area. The basic layered rock model was constructed. This area is overlain by several meters of Quaternary loose sediment. The underlying rocks generally consist of Tertiary sandstone, Permian sandstones and mudstones embedded with peat. Fig. 11 shows the core from drill hole showing the peat at a depth of −57 m. The inverted cross-section of the survey line is shown in Fig. 12. The inversion result of the segmented loop source TEM data has been incorporated into this section, which is delineated by the red rectangle in Fig. 12. It is inferred from Fig. 12 that the inverted models at profile coordinate 0 to 60 m are consistent with the rest of the survey line. This indicates that our inversion scheme will lead to concurrent results regardless of the shape of the transmitting loop. The conductive overburden in Fig. 12 can be interpreted as Quaternary loose sediments. Typically the resistivity of the water-filled mined out area is lower than 50 Ω·m and the air-filled one is larger than 500 Ω·m. Accordingly, we deduced four mined-out areas at profile positions 205 m, 340 m, 480 m, and 640 m at depths of 80 m to 130 m. The deduced minedout area at profile position 340 m is consistent with the information collected among the local residents. The geo-electric information provided by the inversion result will give reference for the design of the path of the railway, and other geophysical tools such as seismic reflection method will be used to confirm this result in the near future. 6. Conclusion The loop source TEM sounding is one of the most popular geophysical tools for the detection of a near surface conductor, especially for the detection of a mined-out area in China. In the cases of the presence of complex terrain, lake and river, or dwellings, laying the commonly used square transmitting loop or rectangle loop is limited, especially when we adopt a large transmitting loop, and changing the shape of the transmitting loop is one of the labor saving solutions.
Fig. 12. The inverted cross-section of the survey line. The mined-out areas are deduced at profile positions 205 m, 340 m, 480 m, and 640 m at depths of 80 m to 130 m.
H. Li et al. / Journal of Applied Geophysics 128 (2016) 87–95
This paper presents the inversion scheme for the segmented loop source TEM data. When changing the shape of the transmitting loop, we can always find a segmented one. Hence, the development of an inversion routine to interpret these data has its practical meaning. We first give the formulas for the calculation of forward response and the Jacobian. The analysis of the forward responses shows that the effect of loop geometry cannot be neglected. An OCCAM inversion scheme which incorporates the loop geometry is developed, and the synthetic data test of this routing demonstrates that it is necessary to consider the change of loop geometry during the processing of the loop source TEM data. This case study illustrates the capability of the inversion scheme of interpreting the field data. The inversion scheme has been applied to the segmented irregular loop source data and the square loop source data. The inverted cross-section is consistent with the basic geological model in the survey area. One of the deduced mined-out areas is in good agreement with the information collected among the local residents. Acknowledgments This research was supported by the R&D of Key Instruments and Technologies for Deep Resources Prospecting (the National R&D Projects for Key Scientific Instruments), Grant No. ZDYZ2012-1-05-04, the State Major Basic Research Program of the People's Republic of China (2012CB416605) and the Natural Science Foundation of China (41174090 and 41474095). References Anderson, W.L., 1979. Numerical integration of related Hankel transforms of orders 0 and 1 by adaptive digital filtering. Geophysics 44, 1287–1305. Auken, E., Christiansen, A., Jacobsen, L., Sørensen, K., 2004. Laterally constrained 1Dinversion of 3D TEM data. 10th European Meeting of Environmental and Engineering Geophysics. Bhattacharya, P., 2012. Direct Current Geoelectric Sounding: Principles and Interpretation. Elsevier. Constable, S.C., Parker, R.L., Constable, C.G., 1987. Occam's inversion: a practical algorithm for generating smooth models from electromagnetic sounding data. Geophysics 52, 289–300. Danielsen, J.E., Auken, E., Jørgensen, F., Søndergaard, V., Sørensen, K.I., 2003. The application of the transient electromagnetic method in hydrogeophysical surveys. J. Appl. Geophys. 53, 181–198. Dettman, J.W., 2013. Mathematical Methods in Physics and Engineering. Courier Corporation.
95
Farquharson, C., Oldenburg, D., 1996. Approximate sensitivities for the electromagnetic inverse problem. Geophys. J. Int. 126, 235–252. Fitterman, D.V., Stewart, M.T., 1986. Transient electromagnetic sounding for groundwater. Geophysics 51, 995–1005. Goldman, M., Tabarovsky, L., Rabinovich, M., 1994. On the influence of 3-D structures in the interpretation of transient electromagnetic sounding data. Geophysics 59, 889–901. Key, K., 2009. 1D inversion of multicomponent, multifrequency marine CSEM data: methodology and synthetic studies for resolving thin resistive layers. Geophysics 74, F9–F20. Munkholm, M.S., Auken, E., 1996. Electromagnetic noise contamination on transient electromagnetic soundings in culturally disturbed environments. J. Environ. Eng. Geophys. 1, 119–127. Nekut, A., 1987. Direct inversion of time-domain electromagnetic data. Geophysics 52, 1431–1435. Newman, G.A., Anderson, W.L., Hohmann, G.W., 1987. Interpretation of transient electromagnetic soundings over three-dimensional structures for the central-loop configuration. Geophys. J. Int. 89, 889–914. Poddar, M., 1982. A rectangular loop source of current on a two-layered earth. Geophys. Prospect. 30, 101–114. Qi, Y., Huang, L., Wu, X., Fang, G., Yu, G., 2014. Effect of loop geometry on TEM response over layered earth. Pure Appl. Geophys. 1–9. Rabinovich, M., 1995. Errors of 1-D interpretation of 3-D TDEM data in the application of mapping saltwater/freshwater contact. J. Appl. Geophys. 34, 23–34. Smith, R., Fountain, D., Allard, M., 2003. The MEGATEM fixed-wing transient EM system applied to mineral exploration: a discovery case history. First Break 21. Smith, R.S., Hodges, G., Lemieux, J., 2009. Case histories illustrating the characteristics of the HeliGEOTEM system. Explor. Geophys. 40, 246–256. Sørensen, K.I., Auken, E., 2004. SkyTEM—a new high-resolution helicopter transient electromagnetic system. Explor. Geophys. 35, 194–202. Tartaras, E., Zhdanov, M.S., Wada, K., Saito, A., Hara, T., 2000. Fast imaging of TDEM data based on S-inversion. J. Appl. Geophys. 43, 15–32. Tezkan, B., Siemon, B., 2014. Appraisal of a new 1D weighted joint inversion of ground based and helicopter-borne electromagnetic data. Geophys. Prospect. 62, 597–614. Viezzoli, A., Christiansen, A.V., Auken, E., Sørensen, K., 2008. Quasi-3D modeling of airborne TEM data by spatially constrained inversion. Geophysics 73, F105–F113. Ward, S.H., Hohmann, G.W., 1988. Electromagnetic theory for geophysical applications. Electromagn. Methods Appl. Geophys. 1, 131–311. Witherly, K., Irvine, R., Morrison, E., 2004. The Geotech VTEM time domain helicopter EM system. ASEG Extended Abstracts 2004, pp. 1–4. Xue, G.q., Bai, C.y., Yan, S., Greenhalgh, S., Li, M.-f., Zhou, N.-n., 2012a. Deep sounding TEM investigation method based on a modified fixed central-loop system. J. Appl. Geophys. 76, 23–32. Xue, G.Q., Qin, K.Z., Li, X., Li, G.-M., Qi, Z.-P., Zhou, N.-N., 2012b. Discovery of a large-scale porphyry molybdenum deposit in Tibet through a modified TEM exploration method. J. Environ. Eng. Geophys. 17, 19–25. Zhou, N., Xue, G., 2012. Simulation Error Of Rectangular Loop Circular Loop in transient electromagnetic Survey. Geophys. Geochem. Explor. 36 (S1), 57–60. Zhou, N., Xue, G., 2014. The ratio apparent resistivity definition of rectangular-loop TEM. J. Appl. Geophys. 103, 152–160.
Contents lists available at ScienceDirect
Journal of Applied Geophysics journal homepage: www.elsevier.com/locate/jappgeo
Inversion of arbitrary segmented loop source TEM data over a layered earth Hai Li a,b,⁎, Guo-qiang Xue a, pan Zhao a,b, Nan-nan Zhou a, Hua-sen Zhong a,b a b
Key Laboratory of Mineral Resources, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing, 100029, China University of Chinese Academy of Sciences, China
a r t i c l e
i n f o
Article history: Received 29 October 2015 Received in revised form 17 February 2016 Accepted 23 March 2016 Available online 29 March 2016 Keywords: OCCAM inversion Loop source TEM method Mined out area Segmented loop Loop geometry
a b s t r a c t The loop source TEM method has been widely used in the detection of a mined out area in China. In the cases that the laying of traditional rectangle or square transmitting loop is limited due to the presence of obstacle on the path of the loop, the changing of the shape of the transmitting loop to bypass the obstacle is a labor saving solution. A numeric integration scheme is proposed to calculate the response and Jacobian of the segmented loop source from that of an electric dipole source. The comparison of forward response between the segmented loop and square loop shows the effect of loop geometry on the decay curves. In order to interpret the data from an irregular source loop, this paper presents an inversion scheme that incorporate the effect of loop geometry. The proposed inversion scheme is validated on the synthetic data, and then applied to the field data. The result reveals that the developed inversion scheme is capable of interpreting the segmented loop source TEM field data. © 2016 Elsevier B.V. All rights reserved.
1. Introduction The time domain electromagnetic method (TDEM of TEM) is one of the major exploration tools for electromagnetic (EM) methods. The loop source TEM method which is highly efficient and gives an excellent resolution of conductive layer at depth, has been widely used in hydrogeophysical investigation and mineral exploration (Danielsen et al., 2003; Fitterman and Stewart, 1986; Xue et al., 2012b; Zhou and Xue, 2014). The theory and interpretation techniques of the TEM method have drawn a dramatic surge of interest since the early 1950s. Commonly, TEM data can be interpreted by using fast imaging techniques (Tartaras et al., 2000) or 1-D inversion base on fitting the theoretical decay curves from a 1-D structure to the real data in the least square sense (Farquharson and Oldenburg, 1996). The interpretation skills mentioned above are generally based on the square or rectangle transmitting loop. For the airborne TEM method which covers a large portion of TEM surveys, the source loop is wound around the aircraft for fixed-wing transient EM system (Smith et al., 2003) or mounted on a lightweight wooden lattice frame for the helicopter electromagnetic system (Sørensen and Auken, 2004; Smith et al., 2009; Witherly et al., 2004). The geometry of transmitting loop is almost stationary during the data collection procedure so that there is need for considering the variation of the shape of the source loop in process of data interpretation. However, for the land cases, in which ⁎ Corresponding author at: Institute of Geology and Geophysics, Chinese Academy of Sciences, Beitucheng west road, Chaoyang district, Beijing, China. E-mail address: [email protected] (H. Li).
http://dx.doi.org/10.1016/j.jappgeo.2016.03.017 0926-9851/© 2016 Elsevier B.V. All rights reserved.
the laying of the square loop or rectangle loop source is restricted, such as when there exists a scarp, pool, river or dwellings on the path of the transmitting line, we have to redesign the path of the transmitting loop, which is not easy to implement and labor intensive. Otherwise, it is often the case that we simply change the shape of the transmitting loop to bypass these obstacles, which is labor-saving. The adoption of existing interpretation techniques to directly processing these data is likely to get an unreliable result. Zhou and Xue (2012) compared the response of rectangular loop and the circular loop at different location inside the loop and pointed out that the responses are affected by the loop geometry at almost every position inside the loop. Qi et al. (2014) proved that the neglect of the loop geometry will affect the recovered resistivity model. However, the interpretation techniques considering the geometry of the source loop are hardly seen in the literature. In this paper, we propose a 1D inversion scheme that incorporated the geometry of the transmitting loop. Although several papers have been published on the effect of 3D structure on the loop source TEM method (Goldman et al., 1994; Newman et al., 1987; Rabinovich, 1995), the 1D inversion is also commonly used in literatures to invert field data (Danielsen et al., 2003; Nekut, 1987; Tezkan and Siemon, 2014). Besides, the full 3D subsurface models can still be obtained from datasets which consist of many individual sounding by using techniques, such as laterally constrained 1D-inversion (Auken et al., 2004) or spatially constrained inversion 1D inversion(Viezzoli et al., 2008). Indeed, it is limited to apply the 1D inversion to invert multidimensional structures, and we do not expect the recovery of multidimensional resistivity structure with 1D interpretation techniques. Instead, our goal is to develop a 1D inversion scheme especially for the segmented loop
88
H. Li et al. / Journal of Applied Geophysics 128 (2016) 87–95
source TEM data and appraise its feasibility. We first analyze the effect of loop geometry by considering two designed segmented transmitting loop by comparing their responses with the square loop responses. In a further step, the developed inversion scheme is validated on synthetic data, and then applied to field data. A case study shows that the developed inversion scheme is capable of interpreting the field data excited by segmented loop source. 2. Basic theory This section gives the methodology of forward and inverse modeling used for the inversion of arbitrary segmented loop source TEM data. There is no analytical solution for the responses of a loop source with arbitrary shape. The responses presented in this paper are calculated using numerical integration. The integration method and reciprocity theorem have been used in the literatures to calculate the responses excited by a loop source. Poddar (1982) calculated the responses of a rectangular loop using the concept of reciprocity and the known solution for the electric field responses of a vertical dipole source placed over a two-layered earth half space. Qi et al. (2014) extended this method and obtained the solutions for a horizontal loop of arbitrary shape over layered earth. Here, we propose an alternate method to directly calculate the responses for arbitrary segmented loop and derive the formulas for calculating the sensitivity without using the reciprocity theorem.
of the symmetrical property, the vector potential A has two components which are in the direction of the dipole and vertical direction respectively. To reduce Eq. (3) into an ordinary differential equation in the variable z, we apply the 2D Fourier transformation to the x-variable and yvariable of Eq. (3) and get the Helmholtz equation in the wavenumber domain. ̂ S ∂ A 2 2 2 − 2 þ kx þ ky −k A ¼ Ĵ e ∂z 2 ̂
ð4Þ
where kx and ky are the spatial wavenumber variables of the Fourier transform, and A represents the vector potential in the wave number domain. After Eq. (4) is solved, the solution for the vector potential in the spatial domain can be obtained by making use of the relation between the 2D Fourier transformation and the Hankel transformation. Aðx; y; z; ωÞ ¼
1 2π
Z 0
∞
2 2 Â kx þ ky ; z; λ J 0 ðλrÞdλ:
ð5Þ
Hence, the frequency domain EM response for arbitrary oriented dipole source can be obtained by substituting Eq. (5) into Eqs. (1) and (2). Next, we derive the formulations used to calculate the EM field due to arbitrary segmented loop source by integrating the responses of a horizontal electric dipole around the transmitting loop H ðωÞ ¼ ∮ F ðωÞ dl
ð6Þ
L
2.1. Calculation of forward responses Firstly, we consider a 1D model shown in Fig. 1, where the horizontal dipole source is located on the surface of the ground. The algorithm for the calculation of the responses for the horizontal dipole sources have been well studied in the literature. The one adopted in this paper follows the works of Ward and Hohmann (1988) and Key (2009). By using the Schelkunoff vector potential A, the magnetic field B and the electric field E can be specified as B¼∇A
ð1Þ
1 ∇ð∇ AÞ: E ¼ iωA þ μσ
ð2Þ
By assuming a time dependence of eiωt and neglecting the displacement current, the Helmholtz equation for the Schelkunoff potential in the spatial domain can be derived from the Maxwell equation and the constitutive relations.
where F(ω) and H(ω) represent the responses of an arbitrarily oriented dipole source and the arbitrary segmented loop source respectively. L represents the integration path which is along the transmitting loop. For the segmented loop considered here, the source loop can be discretized as H ðωÞ ¼
N Z X n¼1
2
ð3Þ
H ðωÞ ¼
N Z X
yi
Z F x ðωÞdx þ
yiþ1
yi
N Z X F y ðωÞdy ¼ n¼1
xi
xiþ1
F ðωÞ cosαdx
ð8Þ
F ðωÞ cosβ dy
where α and β donate the angle between the directions of the dipole and the x-axis or y-axis respectively. Hence, the EM field excited by arbitrary segmented loop source can be obtained with the integration equation. In the final step, the time domain responses are evaluated using the digital filtering techniques (Anderson, 1979) hðt Þ ¼ −
Fig. 1. The M-layered 1D model and the coordinate system.
xiþ1
xi yiþ1
Z þ
JSe = Iδ(r − r0)
represent the electric dipole where k = − iωμσ, and source imposed at position r0 with the source dipole moment I. Because
ð7Þ
N is the number of segments. The field excited by arbitrarily oriented dipoles can be computed using the vector superposition of the fields produced by the horizontal x-oriented and y-oriented dipoles (Key, 2009). Thus the integral of the arbitrarily oriented dipole source can be calculated using the integral of the responses of the horizontal xoriented and y-oriented dipoles,
n¼1
∇2 A þ k A ¼ − JSe
F ðωÞdl:
Li
2 π
Z 0
∞
HðωÞ Im cosðωtÞdω: ω
ð9Þ
In order to verify the effectiveness of our algorithm for the layered earth, the responses for a 40 × 40 m rectangle central loop configuration on the three layer model have been validated by comparing the results from the Leroi code which is developed by CSIRO and has been integrated on the commercial Maxwell software. The model used for comparison is the typical H-type model and K-type model. The models consist of 100 Ω·m half space with a 100 m thick layer whose resistivity is 10 Ω·m for the H-type model and 1000 Ω·m for the K-type model,
H. Li et al. / Journal of Applied Geophysics 128 (2016) 87–95
and is centered at a 150 m depth. Fig. 2 shows the comparison results which verify the correctness of our algorithm. 2.2. Calculation of sensitivities To conduct the inversion routine, we need to calculate the sensitivities which quantify the effects of the changes in the model parameters on the responses. The derivation of the sensitivities for the arbitrary segmented loop source is based on the sensitivities for the dipole source and the integration around the transmitting loop in the same way as the calculation of the forward responses shown in Eq. (8). For the layered model considered in this paper, the model vector can be described in terms of a linear combination of conductivity functions,
σ ðzÞ ¼
N X
2.3. The inversion algorithm A popular OCCAM 1D inversion algorithm is incorporate into our segmented loop inversion scheme (Constable et al., 1987). For the OCCAM inversion, the model is generally discretized into several layers with fixed thickness. The model vector can be written as m = (m1, …,mN)T. The forward operator and the Jacobian are computed based on the formula derived above. The basic idea behind OCCAM's inversion is an iteratively applied local linearization. The linearization of forward operator with the Taylor theorem leads to the objective function for the inversion 2 min WF mk þ WJ mk Δm−Wd −χ 20 2 2 2 k þ α L m þ Δm ; 2
σ j δ j ðzÞ
ð10Þ
j¼1
where N specifies the number of layers in the model, the function δj(z) is equal to unity in the layer j and zero in another layer. To calculate the sensitivities, we start by specifying the change in vector potential in the wavenumber domain due to the changes in the model parameter. Then the transformations introduced in Eqs. (5), (1) and (2) were adopted to get the sensitivities. Using the model vector defined in Eq. (10) and differentiating Eq. (4), we obtain the Helmholtz equation in the wavenumber domain 2
0 13 2 N X ∂ ∂A 2 2 4− þ @kx þ ky þ σ j δ j ðzÞA5 ¼ iωμ 0 δ j ðzÞA: ∂σ j ∂z2 j¼1
ð11Þ
Eq. (11) is an inhomogeneous ordinary differential equation for ∂A=∂ σ j with boundary condition ∂A=∂σ j →0, as z→ ±∞. Hence, the solution for the boundary value problem here can be obtained with the adjoint Green function method (Dettman, 2013). 0
89
∂A iωμ ðkx ; kx ; 0; ωÞ ¼ − Ids ∂σ j
Z
j
z¼ j−1
h
i2 A j ðkx ; kx ; 0; ωÞ dz
ð12Þ
After the values of ∂A=∂σ j have been obtained, the transformations which are used previously to get the forward response are adopted here to get the sensitivities for the arbitrary segmented loop source.
Fig. 2. The comparison of the responses computed with the algorithm here and the responses calculated with the Leroi routine.
ð12Þ
where F is the forward operator, J is the Jacobian, d is the data vector, χ0 is the target misfit value, and α is the Lagrange multiplier which is used to balance the misfit and the roughness term. L is the roughness matrix to implement the Tikhonov regularization. W is the data weighting matrix to incorporate the standard deviations
−1 −1 W ¼ diag σ −1 1 ; σ 2 ; ; σ M
ð13Þ
where the σi represents an estimation of the error in data. By minimizing the defined objective function (Eq. (12)), we obtain the smoothest resistivity model for a given target misfit. To minimize Eq. (12) with respect to the model vector mk, we rearrange Eq. (12) and let d mk ¼ d− F mk þ J mk mk :
ð14Þ
Then we use the Gauss–Newton method and find the iterative solution, mkþ1 ¼ mk þ Δm ¼
−1 T T J mk J mk þ α 2 LT L J mk d mk ð15Þ
At each iteration, we pick the largest α2 that keeps the misfit value of the solution from exceeding the boundary value χ0 specified in Eq. (12). At the end of the iteration procedure, we will have a solution whose misfit is equal to χ20. 3. The analysis of forward responses To study the effect of the shape of the transmitting loop on the responses and demonstrate that it is necessary to develop the inversion scheme for arbitrary segmented loop sources, we calculate and compare the responses excited by three different loop sources, as is shown in Fig. 3. The segmented loop source 1 shown in Fig. 3a, the segmented loop source 2 shown in Fig. 3b and the rectangle loop source considered as a reference in Fig. 3a and b. There are three types of loop source configurations which have been commonly used in engineering investigation and near surface geophysical exploration, including central loop configuration, fixed loop configuration and the modified central loop configuration (Xue et al., 2012a). Accordingly, we design four recording stations which are shown in Fig. 3 and calculated the time derivative of the vertical magnetic induction at each station to get our comparison results. The geological model considered here consists of 10 Ω·m half space with a 100 m thick and 100 Ω·m resistive layer centered at a depth of −130 m. Firstly, the responses excited by the square loop source and the segmented loop source 1 for the three layer model are calculated using the algorithm introduced previously. The decay curves at recording stations ①–④ are plotted in Fig. 4a–d respectively. In this case, the responses of the square loop source are larger than those of the segmented loop source 1. The recording stations ① and ③ are designed for the modified
90
H. Li et al. / Journal of Applied Geophysics 128 (2016) 87–95
Fig. 3. The geometry of the loop source used for the forward modeling. ①, ②, ③, and ④ denote the position of the recording station. (a) The segmented loop source 1 (consists of the solid line and the red solid line) and the rectangle loop source (consists of the solid line and the dash line); (b) the segmented loop source 2 (consists of the solid line and the red solid line) and the rectangle loop source (consists of the solid line and the dashed line). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Fig. 4. The comparison of responses between a conventional rectangular loop source and the segmented loop source 1. The decay curves in (a) to (d) correspond to the recording stations ① to ④ shown in Fig. 3.
H. Li et al. / Journal of Applied Geophysics 128 (2016) 87–95
91
Fig. 5. The comparison of responses between a conventional rectangular loop source and the segmented loop source 2. The decay curves in (a) to (d) correspond to the recording stations ① to ④ shown in Fig. 3.
central loop source configuration. The decay curves for the segmented loop source at these stations show a different character. The early time response at station ① is affected by the loop geometry while the response at station ③ appears largely unaffected. The station ② is designed for the central loop source, and the effect of the geometry on the response at this station is basically impregnable. The station ④ is designed for the fixed loop source configuration, as is shown in Fig. 4d, the transition period of the response appears affected while the other time range is slightly affected. Besides, the late time responses are slightly offset from the square loop responses due to the difference of the area of the transmitting loop. Similarly, the responses excited by the square loop source and the segmented loop source 2 are compared, as is shown in Fig. 5. Unlike in the case of the segmented loop source 1, the response of the segmented loop source is larger than that of the square loop. For the recording station located at the center of the loop, where the response is basically impregnable for the segmented loop source 1, the response is affected by the loop geometry. The early time response is deviated from the response of the square loop source. From the comparison of the forward response, we conclude that the responses excited by the loop source are influenced by the geometry of the transmitting loop. These effects may be introduced by the different transmitting moments, or irregular sides of the loop. For the recording station inside the loop, the deviation of the response from the square loop mainly appears at an early time, depending on the location of the
station. In the case of the recording station outside of the loop, the responses at the transition period, where the time derivatives of the magnetic induction change its sign, are most likely to be affected. Hence, the use of inversion scheme for square loop source TEM data to invert the arbitrary segmented loop source data may lead to an inappropriate result, which makes it necessary to develop the inversion scheme for the arbitrary segmented loop source. The differences of the calculated responses will have an effect on the inverted model. Generally, the root mean square (RMS) misfit is considered in the inversion scheme
χ RMS
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi " #ffi u i 2 n u1 X F ðmÞi1 −d2 t ; ¼ n i¼1 σ i1
ð16Þ
where n is the number of data and σi is the standard error for the measured data, m is the model vector, and F is the forward operator. For the modeled responses considered here, we adopt the relative Table 1 The relative error at each recording station between the segmented transmitting loop and square loop. Recording station
1–①
1–②
1–③
1–④
2–①
2–②
2–③
2–④
Relative error
0.99
0.51
0.43
1.55
1.47
1.22
0.65
2.00
92
H. Li et al. / Journal of Applied Geophysics 128 (2016) 87–95
Fig. 6. Inversion results for recording station ① of segmented loop source 1. (a) Synthetic data and their fitting for square loop source, (b) synthetic data and their fitting for segmented loop source and (c) the inverted model.
error to evaluate the difference introduced by the change of the geometry of the transmitting loop
4. 1D inversion of synthetic segmented loop source TEM data
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi " # u i i 2 n u1 X d1 −d2 : E¼t i n i¼1 d1
In this section, the inversion scheme is validated on synthetic segmented loop source TEM data. For comparison, these data are also inverted with the square loop inversion routine. A five layer model (resistivity: 100, 70, 30, 60, 20 Ω·m; thickness: 50, 75, 50, 75 m) is considered in the inversion process for it has a combination of resistivity model such as ρ1 N ρ2 N ρ3 b ρ4 b ρ5, the conventional Q- and K-type models (Bhattacharya, 2012; Tezkan and Siemon, 2014). Besides, the inversion result for the layer at different depth can be compared in this model. The synthetic data excited by the segmented loop source in Fig. 3 in which the recording stations ① and ② are considered. The Gaussian noise generated with the noise model proposed by (Munkholm and Auken, 1996) is added into the dataset. The added Gaussian noise whose standard deviation is 3% is log-gated and gated stacked, and the resulting log-gated noise gets different amplitudes at
ð17Þ
where di1 and di2 represent the modeled response of the segmented loop source and the square loop source respectively. Table 1 shows the results calculated by Eq. (17). The relative error introduced by the segmented loop 2 is larger for than that of the segmented loop 1. For the recording station considered, the responses at station ① and station ④ are more largely affected than that of station ② and station ③. The comparison of the relative error can give us some insight to the effect on the models recovered, which may be introduced by the loop geometry.
Fig. 7. Inversion result for the data collected at recording stations ① and ② of segmented loop source 2. (For interpretation of the reference to color in this figure legend, the reader is referred to the web version of this article.)
H. Li et al. / Journal of Applied Geophysics 128 (2016) 87–95
93
Fig. 8. Geographic location of the selected profile.
each time gate to fit the decaying properties of time-domain EM data. The misfit is calculated with Eq. (16), and the target misfit is 1.00. Firstly, we inverted the data excited by the segmented loop source shown in Fig. 3a, and the recording station is located at position ①. The synthetic data is inverted using both the inversion scheme of square loop source and arbitrary segmented loop source. Fig. 6a and b show synthetic data and their fitting for the segmented loop source inversion scheme and the square loop source inversion scheme respectively. Fig. 6c shows the inverted models of segmented loop source 1. We can see that the model results (square loop inversion and segmented loop
inversion) were not consistent for the shallow layers but consistent in the relatively deeper layer. The resistivity at each layer for this model is not well resolved, but depth of each layer is recovered quite well for the segmented loop inversion. Fig. 7 shows the models obtained by inverting the data excited by segmented loop source 2 at the recording stations ① and ②. The inverted models for the segmented loop inversion at these two stations are approximately the same, and the obtained resistivity values are well resolved. This illustrated that with the developed inversion scheme, we can get a reliable result from the segmented loop source data. As a contrast, when we use the square loop inversion scheme to invert the segmented loop source data, the deviation of results from the true model is non-negligible, especially for the shallow layer.
Fig. 9. The segmented transmitting loop at the beginning of the survey line. There are four recording stations in this loop. The coordinated system was built for the inversion.
Fig. 10. The decay curve collected at the station where profile coordinate is 40 m with the error bar of each data.
94
H. Li et al. / Journal of Applied Geophysics 128 (2016) 87–95
Fig. 11. The core from drill hole showing the peat at a depth of −57 m which is partly weathered.
From the inversion results shown in the examples above, we conclude that the square loop inversion scheme will lead to a result deviated from the true model when the geometry of the source loop is changed. For the considered cases, the deviation appears at a shallow depth. On the other hand, using the developed arbitrary segmented loop inversion scheme can well resolve the five layer model by including the loop geometry into the inversion algorithm. Therefore, it is necessary to consider the change of loop geometry during the processing of the loop source TEM data.
5. 1D inversion of segmented loop source TEM field data After the successful validation of the segmented loop source inversion on the synthetic data, it has been applied to field data. In a survey aimed at investigating the mined out area of an iron ore mine in Shan Xi province, China, we collected the field data excited by a segmented loop source data at the first few stations of the survey line. Fig. 8 shows the geographic location of the selected profile. The survey line is designed based on the survey and design of an inter-city railway. The mined-out area in this iron mine is caused by illegal excavation among folks, and is one of the major threats to the safety of the intercity railway. Due to the existence of a scarp at the beginning of the survey line, the shape of the first transmitting loop had to be changed as shown in Fig. 9. The first segmented loop was designed to test our algorithm, but at the same time, the geo-electrical information underground at these recording station is of interest to the mining engineer who is responsible for the mining process in this iron ore mine. The vertex of the transmitting loop was recorded by the GPS so that we can get the shape of the transmitting loop, as is shown in Fig. 9. We built the coordinate system to calculate and invert the response of the segmented transmitting loop which contains four recording stations. The data at the rest of the 800 m survey line was collected with the modified square loop source TEM configurations (Xue et al., 2012a),
and the transmitting loop size was 100 × 100 m. Both the first segmented loop source TEM data and the square loop source data are inverted using our inversion scheme. Hence, the inverted cross-section of the modified central loop TEM data can also be used to validate our result in return. The terraTEM time domain EM survey system designed by Monash GeoScope was used for the measurement. The system was powered by an external 24 V battery power package, and the source waveform was a bipolar rectangle current with 50% duty circle. The receiver was a vertical component magnetic core probe with 10,000 m2 effective area, which is the same as that introduced by Xue et al. (2012b). Fig. 10 shows the field data collected at the profile coordinated 40 m. The measurements are made from 10−4 s to 10−2 s. The error bar represents the standard deviation of each data point. It can be seen in Fig. 10 that the last few data points are submerged in the noise. It is illustrated previously that the data stand deviations have been incorporated into the solution so that we can still get an appropriate result at the presence of noise, which has been verified by the synthetic data test. The program of the survey and design of an inter-city railway carried out several drillings to get some geological knowledge of this area. The basic layered rock model was constructed. This area is overlain by several meters of Quaternary loose sediment. The underlying rocks generally consist of Tertiary sandstone, Permian sandstones and mudstones embedded with peat. Fig. 11 shows the core from drill hole showing the peat at a depth of −57 m. The inverted cross-section of the survey line is shown in Fig. 12. The inversion result of the segmented loop source TEM data has been incorporated into this section, which is delineated by the red rectangle in Fig. 12. It is inferred from Fig. 12 that the inverted models at profile coordinate 0 to 60 m are consistent with the rest of the survey line. This indicates that our inversion scheme will lead to concurrent results regardless of the shape of the transmitting loop. The conductive overburden in Fig. 12 can be interpreted as Quaternary loose sediments. Typically the resistivity of the water-filled mined out area is lower than 50 Ω·m and the air-filled one is larger than 500 Ω·m. Accordingly, we deduced four mined-out areas at profile positions 205 m, 340 m, 480 m, and 640 m at depths of 80 m to 130 m. The deduced minedout area at profile position 340 m is consistent with the information collected among the local residents. The geo-electric information provided by the inversion result will give reference for the design of the path of the railway, and other geophysical tools such as seismic reflection method will be used to confirm this result in the near future. 6. Conclusion The loop source TEM sounding is one of the most popular geophysical tools for the detection of a near surface conductor, especially for the detection of a mined-out area in China. In the cases of the presence of complex terrain, lake and river, or dwellings, laying the commonly used square transmitting loop or rectangle loop is limited, especially when we adopt a large transmitting loop, and changing the shape of the transmitting loop is one of the labor saving solutions.
Fig. 12. The inverted cross-section of the survey line. The mined-out areas are deduced at profile positions 205 m, 340 m, 480 m, and 640 m at depths of 80 m to 130 m.
H. Li et al. / Journal of Applied Geophysics 128 (2016) 87–95
This paper presents the inversion scheme for the segmented loop source TEM data. When changing the shape of the transmitting loop, we can always find a segmented one. Hence, the development of an inversion routine to interpret these data has its practical meaning. We first give the formulas for the calculation of forward response and the Jacobian. The analysis of the forward responses shows that the effect of loop geometry cannot be neglected. An OCCAM inversion scheme which incorporates the loop geometry is developed, and the synthetic data test of this routing demonstrates that it is necessary to consider the change of loop geometry during the processing of the loop source TEM data. This case study illustrates the capability of the inversion scheme of interpreting the field data. The inversion scheme has been applied to the segmented irregular loop source data and the square loop source data. The inverted cross-section is consistent with the basic geological model in the survey area. One of the deduced mined-out areas is in good agreement with the information collected among the local residents. Acknowledgments This research was supported by the R&D of Key Instruments and Technologies for Deep Resources Prospecting (the National R&D Projects for Key Scientific Instruments), Grant No. ZDYZ2012-1-05-04, the State Major Basic Research Program of the People's Republic of China (2012CB416605) and the Natural Science Foundation of China (41174090 and 41474095). References Anderson, W.L., 1979. Numerical integration of related Hankel transforms of orders 0 and 1 by adaptive digital filtering. Geophysics 44, 1287–1305. Auken, E., Christiansen, A., Jacobsen, L., Sørensen, K., 2004. Laterally constrained 1Dinversion of 3D TEM data. 10th European Meeting of Environmental and Engineering Geophysics. Bhattacharya, P., 2012. Direct Current Geoelectric Sounding: Principles and Interpretation. Elsevier. Constable, S.C., Parker, R.L., Constable, C.G., 1987. Occam's inversion: a practical algorithm for generating smooth models from electromagnetic sounding data. Geophysics 52, 289–300. Danielsen, J.E., Auken, E., Jørgensen, F., Søndergaard, V., Sørensen, K.I., 2003. The application of the transient electromagnetic method in hydrogeophysical surveys. J. Appl. Geophys. 53, 181–198. Dettman, J.W., 2013. Mathematical Methods in Physics and Engineering. Courier Corporation.
95
Farquharson, C., Oldenburg, D., 1996. Approximate sensitivities for the electromagnetic inverse problem. Geophys. J. Int. 126, 235–252. Fitterman, D.V., Stewart, M.T., 1986. Transient electromagnetic sounding for groundwater. Geophysics 51, 995–1005. Goldman, M., Tabarovsky, L., Rabinovich, M., 1994. On the influence of 3-D structures in the interpretation of transient electromagnetic sounding data. Geophysics 59, 889–901. Key, K., 2009. 1D inversion of multicomponent, multifrequency marine CSEM data: methodology and synthetic studies for resolving thin resistive layers. Geophysics 74, F9–F20. Munkholm, M.S., Auken, E., 1996. Electromagnetic noise contamination on transient electromagnetic soundings in culturally disturbed environments. J. Environ. Eng. Geophys. 1, 119–127. Nekut, A., 1987. Direct inversion of time-domain electromagnetic data. Geophysics 52, 1431–1435. Newman, G.A., Anderson, W.L., Hohmann, G.W., 1987. Interpretation of transient electromagnetic soundings over three-dimensional structures for the central-loop configuration. Geophys. J. Int. 89, 889–914. Poddar, M., 1982. A rectangular loop source of current on a two-layered earth. Geophys. Prospect. 30, 101–114. Qi, Y., Huang, L., Wu, X., Fang, G., Yu, G., 2014. Effect of loop geometry on TEM response over layered earth. Pure Appl. Geophys. 1–9. Rabinovich, M., 1995. Errors of 1-D interpretation of 3-D TDEM data in the application of mapping saltwater/freshwater contact. J. Appl. Geophys. 34, 23–34. Smith, R., Fountain, D., Allard, M., 2003. The MEGATEM fixed-wing transient EM system applied to mineral exploration: a discovery case history. First Break 21. Smith, R.S., Hodges, G., Lemieux, J., 2009. Case histories illustrating the characteristics of the HeliGEOTEM system. Explor. Geophys. 40, 246–256. Sørensen, K.I., Auken, E., 2004. SkyTEM—a new high-resolution helicopter transient electromagnetic system. Explor. Geophys. 35, 194–202. Tartaras, E., Zhdanov, M.S., Wada, K., Saito, A., Hara, T., 2000. Fast imaging of TDEM data based on S-inversion. J. Appl. Geophys. 43, 15–32. Tezkan, B., Siemon, B., 2014. Appraisal of a new 1D weighted joint inversion of ground based and helicopter-borne electromagnetic data. Geophys. Prospect. 62, 597–614. Viezzoli, A., Christiansen, A.V., Auken, E., Sørensen, K., 2008. Quasi-3D modeling of airborne TEM data by spatially constrained inversion. Geophysics 73, F105–F113. Ward, S.H., Hohmann, G.W., 1988. Electromagnetic theory for geophysical applications. Electromagn. Methods Appl. Geophys. 1, 131–311. Witherly, K., Irvine, R., Morrison, E., 2004. The Geotech VTEM time domain helicopter EM system. ASEG Extended Abstracts 2004, pp. 1–4. Xue, G.q., Bai, C.y., Yan, S., Greenhalgh, S., Li, M.-f., Zhou, N.-n., 2012a. Deep sounding TEM investigation method based on a modified fixed central-loop system. J. Appl. Geophys. 76, 23–32. Xue, G.Q., Qin, K.Z., Li, X., Li, G.-M., Qi, Z.-P., Zhou, N.-N., 2012b. Discovery of a large-scale porphyry molybdenum deposit in Tibet through a modified TEM exploration method. J. Environ. Eng. Geophys. 17, 19–25. Zhou, N., Xue, G., 2012. Simulation Error Of Rectangular Loop Circular Loop in transient electromagnetic Survey. Geophys. Geochem. Explor. 36 (S1), 57–60. Zhou, N., Xue, G., 2014. The ratio apparent resistivity definition of rectangular-loop TEM. J. Appl. Geophys. 103, 152–160.