One-dimensional constrained inversion study of TEM and application in coal goafs’ detection

Si Yuanlei; Li Maofei; Liu Yaoning; Guo Weihong

doi:10.1515/geo-2020-0148

Article Open Access

One-dimensional constrained inversion study of TEM and application in coal goafs’ detection

Si Yuanlei , Li Maofei , Liu Yaoning and Guo Weihong

Published/Copyright: November 30, 2020

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information

From the journal Open Geosciences Volume 12 Issue 1

Abstract

Transient electromagnetic method (TEM) is often used in urban underground space exploration and field geological resource detection. Inversion is the most important step in data interpretation. Because of the volume effect of the TEM, the inversion results are usually multi-solvable. To reduce the multi-solvability of inversion, the constrained inversion of TEM has been studied using the least squares method. The inversion trials were performed using two three-layer theoretical geological models and one four-layer theoretical geological model. The results show that one-dimensional least squares constrained inversion is faster and more effective than unconstrained inversion. The induced electromotive force attenuation curves of the inversion model indicate that the same attenuation curve may be used for different geological conditions. Therefore, constrained inversion using known geological information can more accurately reflect the underground geological information.

Keywords: transient electromagnetic method; constrained inversion; the least squares method; 1D TEM forward simulation; induced emf

1 Introduction

Geophysics methods have been widely used to detect the mineral resources [1,2,3,4,5], geological structures [6,7], and urban underground space construction [8,9,10,11]. As the safety of coal mines is increasingly valued by people, it is necessary to detect the water-rich areas around the working face before coal mining, such as goafs and water-conducting structures. As one of the geophysical exploration methods, the transient electromagnetic method (TEM), on the strength of its high efficiency, is widely used in advanced detection of tunnels, accurate positioning of underground pipelines, and detection of holes underground.

At present, the commonly used processing methods are calculating apparent resistivity [12,13,14,15], drawing multi-track curves of induced electromotive force (emf) [16] and one-dimensional (1D) inversion [17,18,19]. The volume effect of TEM makes it difficult for the apparent resistivity and multi-track curves to identify the edge and depth of the target. The general inversion methods without geological information are marked by slow speed, low accuracy, and multi-solvability.

Before proceeding with exploration, we are usually equipped with some underground information, such as the electrical characteristics of the pipeline and the depth of the coal seam. The known information can be used in the transient electromagnetic data processing. Therefore, to improve the inversion accuracy, reduce the multi-solvability, and increase the speed, we add known geological information to the initial model for constrained inversion. Compared with other inversion methods, the least squares inversion is much simpler with a high inversion speed, thus adopted for the study of the TEM 1D constraint inversion.

2 Basic principle of TEM and 1D forward simulation

2.1 Basic principle of TEM

TEM involves applying a current to power supply electrodes or transmitting loops, which causes it to form an electromagnetic field in the surrounding space, and then suddenly cutting off the current. The instantaneous disappearance of the electromagnetic field will generate the induced electromagnetic field in the surrounding conductor [20]. We obtain geological information by analyzing the inductive electromagnetic fields.

2.2 1D forward algorithm

The basis of inversion is forward simulation. The ground surface electromagnetic response of the horizontal stratum which is generated by step current of loop source can be obtained from the Maxwell equations. The time domain electromagnetic response of rectangular loop excitation source is shown below [21]:

(1)∂B∂t=R1+R2,

where

R1=−ln2tI4π∑j=1JKjz0¯∫−WW∫0∞λ(1+rTE)×L−yρLJ1(λρL)+L+yρ−LJ1(λρ−L)dλdx′,

R2=−ln2tI4π∑j=1JKjz0¯∫−LL∫0∞λ(1+rTE)×W−xρWJ1(λρW)+W+xρ−WJ1(λρ−W)dλdy′,

where λ is the space wave number; Kj is the G–S algorithm coefficient (can be obtained from equation (2) [22]); W is the half-width and L is the half length of the rectangular loop in meters; J1 is the first-order Bessel function; z0¯=iωμ0 and μ0=4π×10−7 show the vacuum permeability, ω angular frequency; rTE is the reflection coefficient (can be obtained from equation (3)); x′ and y′ are the coordinates of each point on the rectangular emission loop; I is the current; and t is the time in seconds. Solving these expressions require the Hankel transformation [23,24]. Furthermore, ρL, ρ−L, ρW, and ρ−W are shown below:

ρL=[(x′−x)2+(L−y)2]1/2ρ−L=[(x′−x)2+(−L−y)2]1/2ρW=[(y′−y)2+(W−x)2]1/2ρ−W=[(y′−y)2+(−W−x)2]1/2,

(2)Kj=(−1)j+M×∑k=mmin(j,M)kM(2k)!(M−k)!k!(k−1)!(j−k)!(2k−j)!,

where M=J/2; m is the integer of (j+1)/2; and J is the number of G–S coefficient.

(3)rTE=Y0−Y¯1Y0+Y1¯,

where Y0=u0iωμ0. Y1¯ can be obtained from the following recursive formula:

Yn¯=YnYn+1¯+Yntanh(unhn)Yn+Yn+1¯tanh(unhn)⋮Yn¯=YN,

where Yn=unZn¯; Z¯=iωμn (a total of N layers in the model and n represents the nth layer); un=(λ2+kn2)1/2; kn2=ω2μnεn−iωμnσn; μn, εn, and σn are the magnetic permeability, dielectric constant, and electrical conductivity of the nth layer; and hn is the thickness of the nth layer.

To verify the accuracy of the above algorithm, comparison with the analytical solution of uniform half space is now performed. The analytical formula of uniform half space is as follows [25]:

(4)Bz=Iμ02a3πθae−θ2a2+1−32θ2a2erf(θa),

where a is the radius of circular emission loop; θ=μ04tρ; t is the time; ρ=1σ is the resistivity of uniform half space; and erf(θa)=2π∫0θae−t2dt is the error function.

From equation (4), we can obtain induced emf.

(5)∂Bz∂t=Iρa33erf(θa)−2πθa3+2θ2a2e−θ2a2.

The equivalent relationship between the radius a of circular emission loop and the side length 2W and 2L of rectangular emission loop is shown below.

a=4WLπ.

Set up a uniform half-space model for comparison. The conductivity of uniform half space is 0.01 S m⁻¹, the transmitting current is 1 ampere, and the length of transmitting is 100 m. The analytical solution of the induced emf curve is consistent with the 1D forward simulation result (Figure 1).

Figure 1

Comparison of analytical solution and 1D forward simulation.

3 Constrained inversion by the least squares method

3.1 Basic principle of constrained inversion

The inversion process of the least squares method iteratively calculates the modification by fitting the observation data, continuously obtains new model parameters, and finally obtains a model in which fit between the response data and the observation data meets the given accuracy. Then the model that fits the given accuracy is the inversion result. The least squares inversion method is the most commonly used inversion method. Its objective function is the sum of the squares of fitting error of the data and the sum of the squares of the model parameters. The iterative formula for the inversion is obtained by minimizing the objective function [26].

First, build an objective fitting function of 1D TEM inversion as shown below:

(6)φ=∑i=1M|uai−uci⁎|2=||u−u⁎||22=(u−u⁎)T(u−u⁎),

where M is the number of the measured data; uai is the measured induced emf; uci⁎ is the theoretically calculated induced emf; u is the vector of measured data; and u⁎ is the vector that is calculated by the forward model.

The theoretically calculated induced emf is approximated by Taylor’s formula, omitting quadratic and above terms. Then the observation data and model approximate linear relationship are shown below:

(7)uci⁎≈uci+∑j=1N∂uci∂mjΔmj,

where N is the number of model parameters; uci is the calculated emf from initial model; mj is the jth parameter; and Δmj is the jth model modification. Then, write equation (7) as a matrix and substitute it into equation (6):

(8)φ=(u−u0−J⋅Δm)T(u−u0−J⋅Δm),

where J is the Jacobian matrix; Δm is the model modification vector; and u₀ is the emf vector that is calculated from the initial model.

J=∂uc1∂m1⋯∂uc1∂mj⋯∂uc1∂mN⋮⋱⋮⋰⋮∂uci∂m1⋯∂uci∂mj⋯∂uci∂mN⋮⋰⋮⋱⋮∂ucM∂m1⋯∂ucM∂mj⋯∂ucM∂mN,

Δm=Δm1⋮Δmj⋮ΔmNu0=uc1⋮uci⋮ucMu=ua1⋮uai⋮uaM.

If the jth parameter is known, we can add the known information to the Jacobian matrix. The elements of the jth column of the J matrix are all 0. Therefore, the J matrix can be rewritten as below:

J=∂uc1∂m1⋮⋯⋱∂uc1∂mj−10∂uc1∂mj+1⋮⋯⋰∂uc1∂mN⋮∂uci∂m1⋯∂uci∂mj−10∂uci∂mj+1⋯∂uci∂mN⋮⋰⋮⋱⋮∂ucM∂m1⋯∂ucM∂mj−10∂ucM∂mj+1⋯∂ucM∂mN.

Let (u−u0)=u¯, to minimize the objection (equation (8)), then ∂φ∂Δm=0, and

(9)Δm=(JTJ)−1JTu¯.

The model modification Δm will be obtained by solving the equation (9) by singular value decomposition. Then the new model can be derived with Δm.

Use the relative mean square deviation to measure the accuracy of the fit as below:

(10)RMSD=1M∑i=1Muai−uciuai2.

3.2 Inversion and analysis of theoretical data

To verify the feasibility of the TEM least squares constrained inversion, we use two three-layer (H-type and K-type) geoelectric models, a four-layer (HK-type) geoelectric model, and a five-layer model as examples for inversion trial calculation (Figure 2). Parameters of three models are shown in Table 1 (C is the conductivity and T is the thickness). The length and width of rectangular transmitting loop are all 200 m. The transmitting current is 1 ampere.

Figure 2

The 1D models of TEM forward simulation. (a) H or K model; (b) HK model; and (c) HKH model.

Table 1

The parameters of geoelectric models

	H-type		K-type		HK-type		HKH-type
	C (S m⁻¹)	T (m)	C (S m⁻¹)	T (m)	C (S m⁻¹)	T (m)	C (S m⁻¹)	T (m)
First layer	0.01	100	0.1	100	0.01	100	0.01	50
Second layer	0.1	50	0.01	50	0.1	20	0.1	50
Third layer	0.01		0.1		0.005	30	0.01	50
Fourth layer					0.01		0.1	50
Fifth layer							0.01

Figure 3 shows the inversion result models and the theoretical model of H-type. It can be observed from Figure 3(a) that the constrained inversion and unconstrained inversion result models of H-type geoelectric model are consistent with the theoretical model. Errors of the induced emf curves are all within 1%. However, the unconstrained inversion took 2,564 s and 17 iterations, and the constrained inversion took 1,367 s and 13 iterations.

Figure 3

The least squares inversion result models and theoretical models of H-type. (a) Geoelectric model and (b) the induced emf attenuation curves.

Figure 4 is the inversion result model and the theoretical model of K-type. The K-type model inversion result shows that the result model of constrained inversion is closer to the theoretical model than unconstrained inversion result model, and the constrained inversion took only 950 s and nine iterations, whereas unconstrained inversion took 1,370 s and nine iterations. Moreover, the second layer of unconstrained inversion model of K-type has larger conductivity and is thicker than the theoretical model, even though they have the same induced emf attenuation curves.

Figure 4

The least squares inversion result models and theoretical models of K-type. (a) Geoelectric model and (b) the induced emf attenuation curves.

Figure 5 shows the inversion result models and the theoretical model of HK-type. It still shows that constrained inversion is better than unconstrained inversion. Although their attenuation curves are the same, unconstrained inversion cannot invert the third layer, but constrained inversion can. Constrained inversion took 4,819 s and 26 iterations, whereas unconstrained inversion took 6,412 s and 26 iterations.

Figure 5

The least squares inversion result models and theoretical models of HK-type. (a) Geoelectric model and (b) the induced emf attenuation curves.

Figure 6 shows the inversion result models and the theoretical model of HKH-type. It still shows that constrained inversion is better than unconstrained inversion. Although their attenuation curves are the same, unconstrained inversion can only invert the first high-conductivity layer, but not the second high-conductivity layer, but after constraining the first high-conductivity layer, the second high-conductivity layer can be well reversed.

Figure 6

The least squares inversion result models and theoretical models of HKH-type. (a) Geoelectric model and (b) the induced emf attenuation curves.

4 Instance analysis

In mining areas where the coal seam is shallow, the randomness of small coal kilns gives rise to many unknown goafs. The presence of these goafs poses a great threat to the safety of coal mine production and compromises the stability of ground buildings and ground engineering structures. Usually, goafs and coal seam have the same depth and thickness. Therefore, the thickness of both the cover layer and the target area can be restricted during the inversion.

We used TEM to detect shallow goafs in a mining area in Inner Mongolia, China. The side length of the transmitting loop is 20 × 20 m, and the emission current is 18.6 ampere. The instrument is a transient electromagnetic instrument named Tsikl-7. The observation data are normalized by receiving coil area and current. According to known geological data, the average thickness of the coal seam in the locality is 15.6 m, and the average depth is 66.4 m. Therefore, we limited the burial depth and thickness of the target layer during the inversion. However, because it is not known whether water is available, there is no restriction on resistivity.

Figure 7 is the inversion result of one of the observation points. In Figure 7(a), there are two conductive layers in the result, the first being a shallow aquifer at 10–20 m underground, and the second 66.4–80 m underground, which is caused by water accumulation in the goaf. Because the coal mine goaf is rich in sulfides and metal ions, the conductivity of the water-containing goaf is much greater than that of the shallow aquifer. Figure 7(b) represents the comparison between the attenuation curves of the inversion result data and measured data. Because of various noise in actual field work, the inversion curve cannot fit the measured curve as well as the theoretical model.

Figure 7

Inversion results of actual observation data. (a) Model of inversion result and (b) attenuation curves of the inversion model and measured data.

5 Conclusion

From the present research work, we can arrive at the following conclusions:

The 1D least squares constrained inversion result model of TEM is better than the unconstrained inversion on the grounds that it takes less time and iterations, and the inversion results can better reflect the real situation underground.
The induced emf attenuation curves of the constrained inversion result model and unconstrained inversion result model are consistent with the theoretical model’s attenuation curve, and the errors are all within 1%, indicating that different inversion models have the same attenuation curves. Therefore, constrained inversion under known geological information can reduce multi-solvability.
The actual observation data show that when the depth and thickness of a certain layer are constrained, the resistivity and shallow electrical information of the modified position can still be effectively reversed under the influence of noise.

Acknowledgments

The authors thank the reviewers and editors for their valuable advice and help. This research was funded by the Nature Science Foundation of China (Grant No. 41674133) and National Science and Technology Support Program of China (Grant No. 2013BAK06B01).

References

[1] Songda J. Comprehensive application of shallow seismic reflection method and gravity to exploring goaf in gypsum mine. Chinese J Eng Geophys. 2019;16(5):620–5. 10.3969/j.issn.1672-7940.2019.05.011.Search in Google Scholar

[2] Du LZ, Qiu JH, Zhang Q, Zhang XP. Development and application of a high-fidelity and high-resolution telemetry seismic data acquisition system. Chinese J Geophys. 2019;62(10):3964–73. 10.6038/cig2019M0483.Search in Google Scholar

[3] Li M, Liu S, Su B, Ma Y, Sun Q. Study on transient electromagnetic response of high resistivity goafs and its application. Elektron Elektrotechn. 2019;25(1):31–5. 10.5755/j01.eie.25.1.22733.Search in Google Scholar

[4] Jiang Z, Liu L, Liu S, Yue J. Surface-to-underground transient electromagnetic detection of water-bearing goaves. IEEE Trans Geosci Remote Sens. 2019;57(8):5303–18. 10.1109/TGRS.2019.2898904.Search in Google Scholar

[5] Shaohua L. Study on the depth of coal seam floor failure detected by direct current resistivity in borehole. Xuzhou: China University of Mining and Technology; 2019.10.1007/s12665-019-8731-0Search in Google Scholar

[6] Dongning L, Yueqiang Q, Qing H, Qiuliang W, Song L, Xue L. Quaternary activity and its seismo-geological implication of eastern segment of Danjiang fault. J Jilin Univ (Earth Sci Ed). 2019;49(5):1362–275. 10.13278/j.cnki.jjuese.20170318.Search in Google Scholar

[7] Pengjun C, Zhibo L, Xinxue X, Xueming Y, Dong D, Xing L, et al. Shallow seismic exploration and new understanding of active faults in central Tianjin under the background of neotectonics movement. Progr Geophys (in Chinese). 2020;35(3):1128–34. 10.6038/pg2020DD0289.Search in Google Scholar

[8] Loera HL, Leal JAR, Harris PD, Gaytan DET, Ruiz VJM, Gogichaishvili A. Geophysical exploration of fractured-media aquifers at the Mexican Mesa Central: Satellite City, San Luis Potosí, Mexico. Surv Geophys. 2015;36(1):167–84. 10.1007/s10712-014-9302-2.Search in Google Scholar

[9] Dina E, Mohamed A, Abdou A, Farouk S. Monitoring temporal variations in groundwater levels in urban areas using ground penetrating radar. Sci Total Environ. 2019;703:1–10. 10.1016/j.scitotenv.2019.134986.Search in Google Scholar PubMed

[10] Ding X, Hu Z, Fang R, Ma W. Research on accurate identification of underground cavity detected by ground penetrating radar. Bull Surv Map. 2019;S2:188–92. 10.13474/j.cnki.11-2246.2019.0621.Search in Google Scholar

[11] Hao J, Li M, Liao Y, Fan C, Hu D. Detecting and analyzing leak location of high resistivity pipeline by DC resistivity method. Progr Geophys. 2018;33(2):0803–7. 10.6038/pg2018BB0101.Search in Google Scholar

[12] Zhang W, Zhen Q, Di Q. Study on MTEM frequency-domain ratio method and apparent resistivity calculation. Chinese J Geophys. 2018;61(10):4171–81. 10.6038/cjg2018L0744.Search in Google Scholar

[13] Cui J, Xue G, Chen W, Deng J. Calculation of all-time apparent resistivity for B field due to electrical source TEM. Progr Geophys. 2015;30(6):2690–7. 10.6038/pg20150629.Search in Google Scholar

[14] Zhao Y, Li X, Wang W. Full-domain apparent resistivity definition for large-loop TEM. Progr Geophys. 2015;30(6):1856–63. 10.6038/pg20150446.Search in Google Scholar

[15] Raiche AP, Spies BR. Coincident loop transient electromagnetic master curves for interpretation of two-layer earths. Geophysics. 1981;46(1):53–64. 10.1190/1.1441139.Search in Google Scholar

[16] Yaoning L, Shucai L, Maofei L, Xinming L, Weihong G. A study on transient electromagnetic interpretation method based on the seismic wave impedance inversion model. Open Geosci. 2019;11:572–80. 10.1515/geo-2019-0047.Search in Google Scholar

[17] Li F, Yang H, Deng J, Tang H, Liu S. One-dimensional smoke ring inversion technology of ground transient electromagnetic method. Progr Geophys. 2016;31(2):0688–94. 10.6038/pg20160225.Search in Google Scholar

[18] Zhang Z. Research in the application of 1D transient electromagnetic method. Jilin: Jilin University; 2017.Search in Google Scholar

[19] Li Y, Qiang J, Tang J. A research on 1-D forward and inverse airborne transient electromagnetic method. Chinese J Geophys. 2010;53(3):751–9. 10.3969/j.issn.0001-5733.2010.03.031.Search in Google Scholar

[20] Niu Z. Principle of electromagnetic method in time domain. Changsha: Central South University Press; 2007.Search in Google Scholar

[21] Li J, Zhu Z, Liu S, Liu X, Zeng S. Transient electromagnetic field excited by rectangular transmitting loop based on Gaver–Stehfest algorithm. Petrol Geophys Explor. 2011;46(3):489–92. 10.13810/j.cnki.issn.1000-7210.2011.03.016.Search in Google Scholar

[22] Knight JH, Raiche AP. Transient electromagnetic calculation using the Gaver-Stehfest inverse Laplace transform method. Geophysics. 1982;47(1):47–50. 10.1190/1.1441280.Search in Google Scholar

[23] Anderson WL. A hybrid fast Hankel transform algorithm for electromagnetic modeling. Geophysics. 1989;54(2):263–6. 10.1190/1.1442650.Search in Google Scholar

[24] Guptasarma D, Singh B. New digital linear filters for Hankel J0 and J1 transforms. Geophys Prospect. 1997;45:745–62. 10.1046/j.1365-2478.1997.500292.x.Search in Google Scholar

[25] Denghai B, Meju MA, Jian L, Lifeng W, Zhaohai H. Numerical calculation of all-time apparent resistivity for the central loop transient electromagnetic method. Chinese J Geophys. 2003;46(5):697–704. 10.3321/j.issn:0001-5733.2003.05.018.Search in Google Scholar

[26] Ou D. Geophysical inversion tutorial. Beijing: Geological Publishing House; 2015.Search in Google Scholar

Received: 2020-05-19

Revised: 2020-10-29

Accepted: 2020-11-05

Published Online: 2020-11-30

This work is licensed under the Creative Commons Attribution 4.0 International License.

Articles in the same Issue

https://doi.org/10.1515/geo-2020-0148

Keywords for this article

transient electromagnetic method; constrained inversion; the least squares method; 1D TEM forward simulation; induced emf

Creative Commons

BY 4.0