Nonlinear Inversion Using Very Fast Simulated Annealing for Horizontal Electric Dipole Time-Domain Electromagnetic Data

Article information

J Electromagn Eng Sci. 2021;21(5):379-390
Publication date (electronic) : 2021 October 7
doi : https://doi.org/10.26866/jees.2021.5.r.46
1Physics of Earth and Complex System, Physics Department, Faculty of Mathematics and Natural Sciences, Institut Teknologi Bandung, Bandung, Indonesia
2Applied Geophysics and Exploration, Faculty of Mining and Petroleum Engineering, Institut Teknologi Bandung, Bandung, Indonesia
3Institute for Geothermal Sciences, Graduate School of Science, Kyoto University, Kyoto, Japan
*Corresponding Author: Wahyu Srigutomo (e-mail: wahyu@fi.itb.ac.id)
Received 2020 September 22; Revised 2021 January 28; Accepted 2021 April 29.

Abstract

A nonlinear stochastic inversion scheme, called very fast simulated annealing (VFSA), was applied to the time-domain electromagnetic data generated from a horizontal electric dipole. The forward formulation of the vertical magnetic field was expressed in the Laplace domain by applying the Hankel integral transform. Time-domain transformation was performed by applying the inverse Laplace transform using the Gaver–Stehfest algorithm. In this study, for noise-free synthetic data, the VFSA scheme yielded the smallest misfit and an inverted resistivity model that resembled the test model. The addition of 5% random noise to the synthetic data produced the same level of misfit and a model that still mimicked the test model. However, the addition of 10% noise to the synthetic data resulted in a misfit value that was three times that of the first two values and a resistivity model with a large discrepancy with the test model, particularly at large depths. These results indicate the efficacy of the VFSA inversion scheme for inferring the subsurface resistivity structure from the electromagnetic data. This inversion scheme was applied to field data measured in a volcanic environment. The general pattern of the resistivity structure inferred by the VFSA inversion is consistent with the structure obtained previously by using a deterministic inversion scheme.

I. Introduction

The time-domain electromagnetic (TDEM) method, often called the transient electromagnetic (TEM) method, has been widely applied in geophysical explorations. In this method, current is fed into a transmitter for a specific TDEM configuration, such as a transmitter loop or vertical magnetic dipole (VMD), horizontal magnetic dipole (HMD), horizontal electric dipole (HED), and long grounded wire [13]. Groundwater exploration was previously studied using VMD configuration [4]. Meanwhile, a rectangular-loop TDEM was used to characterize the subsurface weathered layer [5] and analyze and interpret the anomalous conductivity and magnetic permeability effects simultaneously [6]. The TDEM method was applied to determine the aquifer potential [7], study subsurface structures in volcanic areas [8], and study saltwater intrusion [9]. It was also used to characterize leachate contamination in landfills [10].

The TDEM forward modeling problem is expressed in the frequency, time, or Laplace domain [11]. Time-domain transformation for solving this problem has been conducted using Fourier or Laplace transform. Fourier transform has often been applied to solve TDEM responses generated by VMD sources [1216]. The TDEM data generated by an HED source were solved using the inverse Laplace transform in the forward modeling part [8]. The inverse Laplace transform was manifested using the Gaver–Stehfest algorithm [17, 18]. In this study, the TDEM forward formulation is expressed in the Laplace domain.

One intriguing challenge in the TDEM method is solving the inversion problem in which the TDEM data are numerically processed to obtain information on the resistivity structure of the subsurface. The inversion methods are grouped into deterministic and non-deterministic methods. The examples of the first group are least-squares inversion [19] and Occam’s inversion [8], which are used to determine the minimum misfit value between the observed and calculated data in order to obtain a reliable subsurface model. The non-deterministic methods, which belong to the global optimization approach, have recently been applied frequently to minimize misfit and avoid solutions trapped in the local minima [20]. The first application of the global optimization approach to TDEM data effectively revealed a subsurface resistivity profile [21]. An efficacy of global optimization methods in solving inversion problems is their procedure that mimics natural phenomena such as cooling temperature based on slow cooling to achieve the ground state, condition with the minimum possible energy, and a physical phenomenon called simulated annealing (SA) [22]. In this study, we elaborate a nonlinear inversion of TDEM data using a global approach in the form of very fast simulated annealing (VFSA), a modified version of SA, for overdetermined problems that contain fewer model parameters than the observed data.

VFSA inversion has been successfully used to infer the subsurface structure underneath Tibet by using the seismic method [23] and to interpret gravity anomaly caused by simple-shaped buried bodies [24]. The advantage of VFSA over other methods is that it can prevent the local minimum from being reached. VFSA inversion ensures the solution’s stability and can be used to make the noise data robust [25]. It was also applied for the interpretation of stacked seismic data [26] and used to obtain a static shift correction of magnetotelluric (MT) data in one-dimensional (1D) MT data interpretation [27], 1D DC resistivity data interpretation [28], analysis of self-potential due to a 2D inclined structure [29], analysis of residual gravity anomaly to define the gross crustal density distribution of a tectonic region [30], and magnetic anomaly inversion caused by a rod-shaped buried structure [31].

Sharma and Biswas [32] successfully conducted TDEM data inversion using VFSA for a coincident loop configuration to detect anomalous conductive bodies in subsurfaces at depths down to 250 m. They used the fast Fourier transform to transform the frequency-domain expression of the forward formulation into the time-domain expression. Li et al. [33] employed VFSA inversion for a large-loop TDEM configuration yielding interpretation depths of up to approximately 1,500 m. They used Laplace transform in the form of the Gaver–Stehfest method in their forward formulation to transform the frequency domain into time domain. Yogi and Widodo [34] applied a hybrid inversion between the conjugate gradient method and the VFSA method for a rectangular transmitting loop and successfully revealed a resistivity profile of down to 1,000 m. Meanwhile, Prabawa and Warsa [35] conducted VFSA inversion for a VMD configuration to infer a carbonate reservoir at depths of 750–1,200 m. They also applied Fourier transform to calculate the electromagnetic responses in the time domain.

This study applies the VFSA method to the data obtained from an HED-TDEM configuration, which can reveal a deeper subsurface resistivity profile of approximately 8–10 km depth [8]. The forward formulation is expressed in the Laplace domain by applying the Hankel transform. Meanwhile, the time-domain transformation is performed by applying the inverse Laplace transform using the Gaver–Stehfest algorithm. The efficacy of this VFSA scheme is tested for a set of synthetic data, both noise-free and noisy. The scheme is also applied to real data obtained from the field to test its applicability in real situations.

II. Forward Modeling

1. TDEM for HED Configuration

Electromagnetic phenomena can fundamentally be represented by relationships among five physical field vectors, namely e (electric field intensity), b (magnetic induction), d (dielectric displacement), h (magnetic field intensity), and j (electric current density). The relations among these vectors are governed by a set of Maxwell’s equations, which are expressed in the form of four first-order differential equations.

(1) ×E+bt=0
(2) ×h+dt=j
(3) ·b=0
(4) ·d=ρ

where ρ is the electric charge density (C/m3).

These equations can be decomposed into two equations by employing the constitutive relations in earth materials, given by

(5) D=ɛE
(6) J=σE
(7) B=μ0H

Here, ɛ is the dielectric permittivity of the material, σ is the electric conductivity, and ω is the angular frequency. The magnetic permeability μ0 is assumed to have the same value as the free-space value.

Applying Fourier transformation to Eqs. (1) and (2) and substituting in the above constitutive relations, a frequency-domain version of Maxwell’s equations is obtained:

(8) ×E+iωμ0H=0
(9) ×H-(σ+iɛω)E=0.

Eqs. (8) and (9) are homogeneous equations, which are only valid in source-free regions. These equations must be substituted by inhomogeneous equations in regions having artificial sources, such as magnetic source current ( Jms) and electric source current ( Jes):

(10) ×E+iωμ0H=-Jms
(11) ×H+(σ+iɛω)E=-Jes.

Using the expressions of potential Schelkunoff A (magnetic term) and F (electric term) for E and H eases the steps required to solve the inhomogeneous differential equation, as the potentials are always parallel to the source fields [1]. The inhomogeneous equations can be solved for infinite homogeneous regions as long as Jes and Jms can be expressed explicitly; hence, the electric and magnetic fields can be regarded as a superposition of electric and magnetic sources

(12a) E=Em+Ee
(12b) H=Hm+He.

Problems that deal with only Em and Hm are associated with cases in which Jes=0, whereas those that deal with only Ee and He are associated with cases in which Jms=0.

In this study, an HED transmitter is used; hence, only the electric source ( Jes) appears in the formulation, simplifying the calculation of Schelkunoff potential to only include A, which by definition is related to He:

(13) He×A.

Substituting Eq. (13) into Eq. (11) yields the inhomogeneous Helmholtz equation:

(14) 2A+k2A=-Jes

where k is a wave number, expressed by

(15) k2=μ0ɛω2-iμ0σω.

Applications of electromagnetic methods in geophysics, such as TDEM, usually consider electromagnetic field propagation in earth materials where the electric conductivity values are within a limited range (100 to 10−4 S/m) and the use of low frequencies (<105 Hz). In return, the conduction current is much greater than the displacement current or μ0σωμ0ɛω2, yielding k2 ≈ −0σω. This condition is often called the quasi-static approximation, where the electromagnetic field propagation in earth materials is predominantly a diffusion phenomenon. In the Laplace domain, k2 = −0σω is expressed as k2 = −0σ, where s = .

In this study, conductivity is assumed to vary only vertically in the z-direction. The earth is composed of N layers, each of which has a constant conductivity. The HED source with a current moment Idx is located at z = −h, where dx denotes the length of the dipole (Fig. 1). σj and lj are the thickness and conductivity of the i-th layer, respectively, where zi−1 < z < zi, and zj is the depth from the surface to the bottom of the j-th layer. Eq. (13) can now be expressed as

Fig. 1

TDEM setup for 1D case where the earth’s conductivity varies only in the z-direction.

(16) Hi=×Ai.

By applying a proper gauge condition, we obtain

(17) (2+ki2)Ai=0,

where ki2=-μ0sσi. In this problem, we use two axes of symmetry, one in the direction of the dipole’s axis (x-axis) and the other in the direction of conductivity changing (z-axis). Accordingly, we can define that the y-component of the vector potential equals zero:

(18) Ai=(Axi,0,Azi).

The vector potential of the electric dipole, which is in the x-axis direction, can be expressed as

(19) A(r)=Ids4πre-ikrux.

In a layer containing the source, the particular (primary) solution of the inhomogeneous differential equation associated with the presence of the source at z = −h is added to the complementary (secondary) solution. Decomposing the source into transverse magnetic (TM) and transverse electric (TE) modes, the particular solutions above and below the surface can be expressed as

(20) Ap(kx,ky)e-u0|z+h|

as the TM mode and

(21) Fp(kx,ky)e-u0|z+h|

as the TE mode. Both potentials Ap and Fp depend on the particular source.

Merging the primary and secondary solutions above the earth, the potentials as functions of x and y are obtained by applying the inverse Fourier transforms

(22) A=14π2--Ape-u0h(e-u0h+rTMeu0z)ei(kxx+kyy)dkxdky

and

(23) F=14π2--Fpe-u0h(e-u0z+rTEeu0z)ei(kxx+kyy)dkxdky

where rTE and rTM are reflection coefficients for the TE and TM modes, respectively, which take into account the admittance and impedance of the free space and the surface. The admittance and impedance of the surface are calculated recursively, starting from the deepest layer and repeated upward until the surface.

The above solutions can now be seen as a superposition of plane waves having a very wide range of the incidence angle. The TE mode has only a vertical magnetic field, and the coefficient of the TE mode, Fp, is obtained by equating the derivative of Eq. (19) with respect to y with Hz, which yields

(24) Hz=0         (TM mode)

and

(25) Hz=1iωμ(2z2+k2)Fz         (TM mode)

to obtain

(26) Fp=-z^0Ids2u0ikykx2+ky2.

Substituting Eq. (26) into Eqs. (20) and (21), we obtain the expression for TE potentials between the source and the earth:

(27) F(x,y,z)=-z^0Ids8π2--[e-u0(z+h)+rTEeu0(z-h)]ikyu0(kx2+ky2)ei(kxx+kyy)dkxdky

The components of the magnetic field are obtained by substituting Eq. (24) into Eq. (25). The vertical magnetic field that depends only on the potential F of the TE mode at the surface (z = 0) is given by

(28) Hz=-Ids8π2--(1+rTE)eu0zikyu0ei(kxx+kyy)dkxdky,

which after transformation to the Hankel transform becomes

(29) Hz=Ids4πyρ0(1+rTE)eu0zλ2u0J1(λρ)dλ,

where ρ=x2+y2 is the horizontal distance from the center of HED to the receiver and J1 is the Bessel function of order one of the first kind. The variable rTE can be calculated as follows

(30) rTE=Y0-Y^1Y0+Y^,

where

  • Y0=u0iωμ0;

  • Y^1=Y1Y^2+Y1tanh(u1h1)Y1+Y^2tanh(u1h1);

  • Y^n=YnY^n+1+Yntanh(unhn)Yn+Y^n+1tanh(unhn);

  • ŶN = YN;

  • Yn=uniμnω;

  • un=(λ2-kn2)1/2;

  • kn2=ω2μnɛn-iωμnɛn;

  • u0 = (λ2ω2μ0ɛ0);

and ω is the angular frequency.

2. Time-Domain Transformation

The final step of the forward modeling part is to transform the magnetic field previously expressed in the Laplace domain into the time domain:

(31) h(t)=12πic-c+Hz(s)s-stds.

The transformation is manifested by performing numerical integration using the Gaver–Stehfest algorithm [17, 18]:

(32) h(t)=ln2tl=2LKlHz(sl)sl=lln2t,

where

(33) Kl=(-1)l+L2k=l+12min(l,L2)kL2(2k)!(L2-k)!k!(l-1)!(1-k)(2k-l),

where L is the coefficient value that depends on the computer used for calculating the inverse. In this calculation, the value of L is 8 [8]. In the real measurement, we are only interested in the transient response after switching off the current at t = 0, and the magnetic field measured as a function of time can be expressed as

(34) h-(t)=hDC-h(t);         t0,

where hDC is the induced magnetic field at t → ∞ or ω → 0 due to the HED source, calculated by a similar step as h(t) for layered earth, except that the value of the Laplace variable s is set to zero.

III. Global Optimization using VFSA Method

Crystal formation in a material caused by a thermodynamic process is adopted as a fundamental concept in SA inversion. In annealing, a solid substance is heated until it melts into liquid. The temperature then decreases slowly, enabling the cooled liquid to grow large single crystals, a condition that is associated with the global minimum energy state. Analogous to annealing, in the geophysical inversion problem, the temperature cooling process resembles an iteration process of searching a solution. The liquid substance represents a model, and the energy of the system is analogous to an objective function. In the inversion, the process of finding a model that converges toward a solution whose objective function is the minimum is repeated at each iteration.

The SA method incorporates the Boltzmann probability distribution function, which describes the relation between the probability of a model m at temperature T, whose energy is E:

(35) P(m)exp (-E(m)kT)

where k is Boltzmann’s constant. The system configuration is expressed in M number of model parameters (m1, m2, …, mM). In the inversion scheme, E is the objective function to be minimized and T is a fixed control factor, commonly referred to as the iterative process controller [25]. As k is a constant, we can set its value to 1 without loss of generality. The model space is sought after by obtruding perturbations on the estimates of the model parameters being optimized. The perturbations depend on the chosen scheme of temperature decay for each iteration.

VFSA was first introduced by Ingber [31, 36] as a modification of the SA method for two main reasons [37,38]. First, in an NM-dimensional model space, each model parameter has a different range of variations and influences the misfit function differently; hence, each model parameter has a different level of perturbation from their current value. Second, the algorithm of the SA method does not perform sufficiently precise and fast calculation if the Cauchy random numbers are the same as the number of model parameters. The attempt to establish an NM-dimensional Cauchy distribution can be avoided by using the NM-product of 1D Cauchy distributions. This modification permits larger sampling of the model space at high temperature and narrower sampling at low temperature, and each model parameter is allowed to have its own cooling schedule and sampling scheme in the model space.

The implementation of the algorithm is started by assuming a model parameter mi at k-iteration ( mik), which can be expressed as [39]:

(36) miminmikmimax,

where mimin is the minimum and mimax is the maximum value of the parameter mi. The value of mik is perturbed at the (k+1)-th iteration using the following rule:

(37) mik+1=mik+yi(mimax-mimin),

where y ∈ [−1,1] and miminmik+1mimax. The value of yi is produced by the following distribution:

(38) gT(y)=i=1NM12(|yi|+T)ln (1+1Ti)=i=1NMgTi(yi),

whose cumulative probability is presented by

(39) GTi(yi)=12+sgn(yi)2ln (1+|yi|Ti)ln (1+1Ti).

It is equivalent to having a random number ui within the uniform distribution of U[0,1], which can be mapped into the previous distribution expressed as

(40) yi=sgn(ui-12)Ti[(1+1Ti)|2ui-1|-1].

Li et al. [35] showed that for such distribution, the global minima can be randomly determined by applying a cooling schedule, which can be expressed as

(41) Ti(k)=T0iexp(-cik1/NM),

where T0i is the initial temperature of the i-th parameter, and ci is a control parameter for the temperature scheduling. In this paper, we set T0 = 5, ci = 1, and k = 1. Each particle is set to move NV times at each iteration. These movements aim to give chance to each particle to achieve the best misfit randomly as the temperature cools down. The VFSA algorithm is systematically explained by the pseudocode shown in Fig. 2.

Fig. 2

A pseudocode for a VFSA algorithm. Modified from Ingber [31] and Sen and Stoffa [20].

IV. Inversion Result and Discussion

1. Inversion on Synthetic Data

The forward modeling calculation is performed to generate a transient curve of vertical magnetic field Bz (tesla, T). In a real situation, the measurement of Bz instead of its time derivative (∂Bz/∂t) is made possible by using a fluxgate sensor as the receiver [8]. The HED transmitter lies on the surface and has a 1,425-m long cable, both ends of which are grounded, injecting an 8 A current. The transmitter is in parallel with the x-axis and the receiver is located at x = 3 km and y = 4 km from the transmitter’s center; thus, the transmitter–receiver distance is 5 km. The earth is modeled by a three-layer subsurface structure that has a moderately resistive first layer (25 Ωm), followed by a conductive layer (1 Ωm) and a resistive layer (100 Ωm). The magnetic field Bz is calculated at 212 measurement times, starting at 1.58 ms and ending at 4.28 seconds. The resistivity model (test model) and the generated vertical magnetic field Bz (synthetic data) are shown in Fig. 3(a) and 3(b), respectively.

Fig. 3

(a) Subsurface resistivity model used for generating the theoretical or synthetic magnetic field Bz. (b) The transient magnetic field Bz, which is generated by performing forward modeling using the subsurface resistivity model described in (a).

To test the pertinency of the proposed VFSA inversion scheme, the synthesis data are inverted to reveal a resistivity model that is sought to be as similar as possible to the test model. The inversion is conducted for several sets of VFSA parameter values to obtain the possible smallest misfit, as listed in Table 1. Within the scope of this study, large values of NV correspond to fairly small misfit values, and relatively smaller NV values produce larger misfit values. In this inversion scheme, the values of the constants C, k, and NM are set to 1, whereas MN is set to 9, which are 5 resistivity values (in Ωm) and 4 thicknesses (in m). The T0 value is varied to observe the effect of the initial temperature in achieving a thermal equilibrium state. The final state of an iteration, which is indicated by its misfit value, becomes an initial condition for the next iteration.

Combination of parameter values in the VFSA global optimization

Table 1 shows that the smallest root mean square (RMS) misfit value of 0.1011% is obtained by a combination of the following parameters: NV = 50, number of iterations = 250, and initial temperature T0 = 5. The inversion results obtained from this combination are shown in Fig. 4. The fitting between the synthetic data and the calculated data resulting from the inversion scheme is remarkably good (Fig. 4(a)). Furthermore, the inverted model of resistivity bears a resemblance to the test or true model, which generates the synthetic data (Fig. 4(b)), with a slight but noticeable discrepancy at lower depths. However, within the scope of our variation of parameter combinations, the relative error between the test and inverted models for this particular combination is the smallest (1.06%), as shown in Table 1.

Fig. 4

(a) Comparison between the noise-free observed Bz data (in this case, the synthetic data) and the calculated data obtained from the inversion scheme. (b) Comparison between the subsurface resistivity test model or the true model and the inverted model for the noise-free data in (a). (c) Curve of the RMS error convergence toward minimum values vs. the number of iterations.

The real measured geophysical data always have a portion of noise in them. In this study, 5% and 10% random noises are added to the synthetic data. As before, these data are inverted using the same combination of parameters. The results obtained for the inversion of the noisy data are shown in Fig. 5. The comparisons between the noisy synthetic data and calculated data for 5% noise and 10% noise are shown in Fig. 5(a) and 5(c), respectively. The RMS misfit of the former is 1.009%, whereas that of the latter is 2.8965%. This shows that until a certain level of noise, the RMS misfit will not be affected much, and the inverted model will somewhat recover the true model optimally. This deduction is supported by the comparison between the test and inverted models for the 5 noisy data, as shown by Fig. 5(b). The inverted model mimics the general feature of the test model optimally, with only slight discrepancies at deeper depths. However, the addition of 10% random noise to the data exhibits an inverted model that resembles the true model at shallow depths only and results in large discrepancies at large depths, as shown in Fig. 5(d). This indicates that a higher level of noise will significantly affect the efficacy of this VFSA inversion scheme to reveal the true model, especially at larger depths, which is associated with the data curve at later times.

Fig. 5

(a) Comparison between the 5% noisy synthetic (observed) data and the calculated data; (b) comparison between the test and inverted models for data in (a). (c) and (d) are the same as (a) and (b), respectively, except that the random noise in the synthetic data is 10%. The parameters in both cases are NV = 50, number of iterations = 250, and initial temperature T0 = 5.

2. Application to the Field Data

The VFSA inversion scheme is also applied to the real data obtained from a TDEM survey conducted at a station called UZ01, which is located around Unzen volcano, Kyushu island, Japan. The profile of the transient data is depicted in Fig. 6. The maximum magnitude of the vertical magnetic field Bz is approximately 105 pT at the early times, and its minimum value is approximately 1 pT. The conditions of NV = 50 and T0 = 5 are set for this inversion, while the maximum number of iterations is set to 400 to give more time for the inversion to reach the smallest misfit possible. An acceptable conformity between the observed and calculated data is achieved at the maximum number of iterations.

Fig. 6

Observed field data and calculated data at Station UZ01 (in pT).

The VFSA inversion successfully reveals the subsurface resistivity model beneath Station UZ01, as shown in Fig. 7(a). The resistivity values start from a moderately resistive layer of approximately 80 Ωm and decrease to a conductive layer of 2 Ωm, and then increase to a highly resistive layer of approximately 250 Ωm. The subsurface resistivity profile exhibits the same general pattern as the profile that resulted from the smoothness-constrained Occam’s inversion (Fig. 7(b)). The latter profile comprises 40 fixed and gradually increased thicknesses from the surface down to the lowermost infinite half-space in contrast to 5 varying thicknesses for VFSA inversion. Both inversion schemes result in the same general resistivity pattern.

Fig. 7

(a) Inverted subsurface resistivity model from the VFSA inversion scheme for the data in Fig. 6. (b) The resistivity model resulting from the smoothness-constrained Occam’s inversion [8].

For both inversion schemes, the differences at the detail level are attributed to those in the number of layers, the constraints used, and the difference in treatment for measurement error. The resistivity structure comprises a moderately resistive surface layer, a conductive second layer several thousand meters deep, and a highly resistive third layer at 5 km depth. Srigutomo et al. [8] interpreted the conductive layer as a combination of a water-saturated layer and hydrothermally altered rocks.

These general features of resistivity obtained by TDEM inversion are also consistent with the features of the resistivity structure previously interpreted from MT surveys [39].

V. Conclusion

In this study, a global optimization approach for the inversion of TDEM-HED magnetic field data to reveal the subsurface resistivity profile was conducted using the VFSA method. We tested the VFSA inversion scheme with noise-free synthetic data generated by a test subsurface model. Within the scope of this study, the combination of VFSA constants that yielded the optimal results was NV = 50, number of iterations = 250, initial temperature T0 = 5, and other constants C, k, and NM = 1. The calculated data fit the synthetic data with only a small RMS error, and the inverted model could mimic the test model. The inversion scheme was then applied to the data with additional 5% and 10% random noises. The inverted resistivity model successfully recovered the test model, with only slight discrepancies for the first two cases. For the case of data with 10% noise, the discrepancies became significantly larger at greater depths, indicating that a higher level of noise can considerably affect the data at later times. To further evaluate the applicability of the VFSA inversion, a set of vertical magnetic field data measured from a station located in a volcano-geothermal area was inverted to reveal the subsurface resistivity model beneath the station. The model resulting from this TDEM data inversion exhibited the same features as those of the model resulting from a smoothness-constrained deterministic Occam’s inversion. These results indicate the effectiveness and usefulness of the VFSA inversion scheme to be used more broadly for TDEM-HED data.

Acknowledgments

The authors wish to thank the Modeling and Inversion Laboratory, Physics of Earth and Complex System, Faculty of Mathematics and Natural Sciences, Institut Teknologi Bandung, Indonesia, for permission to use the computational facility.

References

1. Nabighian MN. Electromagnetic Methods in Applied Geophysics Tulsa, OK: Society of Exploration Geophysicists; 1988.
2. Nabighian MN, Macnae JC. Time domain electromagnetic prospecting methods. Electromagnetic Methods in Applied Geophysics (Volume 2, Application) Tulsa, OK: Society of Exploration Geophysicists; 1991. p. 427–509.
3. Kaufman AA, Keller GV. Frequency and Transient Soundings Amsterdam, Netherlands: Elsevier; 1983.
4. Fitterman DV, Stewart MT. Transient electromagnetic sounding for groundwater. Geophysics 51(4):995–1005. 1986;
5. Farquharson CG, Oldenburg DW, Li Y. An approximate inversion algorithm for time-domain electromagnetic surveys. Journal of Applied Geophysics 42(2):71–80. 1999;
6. Pavlov DA, Zhdanov MS. Analysis and interpretation of anomalous conductivity and magnetic permeability effects in time domain electromagnetic data. Part I: Numerical modeling. Journal of Applied Geophysics 46(4):217–233. 2001;
7. Krivochieva S, Chouteau M. Integrating TDEM and MT methods for characterization and delineation of the Santa Catarina aquifer (Chalco Sub-Basin, Mexico). Journal of Applied Geophysics 52(1):23–43. 2003;
8. Srigutomo W, Kagiyama T, Kanda W, Munekane H, Hashimoto T, Tanaka Y, Utada H, Utsugi M. Resistivity structure of Unzen Volcano derived from time domain electromagnetic (TDEM) survey. Journal of Volcanology and Geothermal Research 175(1–2):231–240. 2008;
9. El-Kaliouby H, Abdalla O. Application of time-domain electromagnetic method in mapping saltwater intrusion of a coastal alluvial aquifer, North Oman. Journal of Applied Geophysics 115:59–64. 2015;
10. Baawain MS, Al-Futaisi AM, Ebrahimi A, Omidvarborna H. Characterizing leachate contamination in a landfill site using Time Domain Electromagnetic (TDEM) imaging. Journal of Applied Geophysics 151:73–81. 2018;
11. Ward SH, Hohmann GW. Electromagnetic theory for geophysical applications. Electromagnetic Methods in Applied Geophysics (Volume 1, Theory) Tulsa, OK: Society of Exploration Geophysicists; 1988. p. 130–311.
12. Newman GA, Hohmann GW, Anderson WL. Transient electromagnetic response of a three-dimensional body in a layered earth. Geophysics 51(8):1608–1627. 1986;
13. Huang H, Palacky GJ. Damped least-squares inversion of time-domain airborne EM data based on singular value decomposition. Geophysical Prospecting 39(6):827–844. 1991;
14. Farquharson CG, Oldenburg DW. Inversion of time-domain electromagnetic data for a horizontally layered earth. Geophysical Journal International 114(3):433–442. 1993;
15. Zhdanov MS, Traynin PN, Portniaguine O. Resistivity imaging by time domain electromagnetic migration (TDEMM). Exploration Geophysics 26(2–3):186–194. 1995;
16. Auken E. 1D time domain electromagnetic interpretations over 2D/3D structures. In : Proceedings of the Symposium on the Application of Geophysics to Engineering and Environmental Problems. Orlando, FL; 1995; p. 329–337.
17. Knight JH, Raiche AP. Transient electromagnetic calculations using the Gaver-Stehfest inverse Laplace transform method. Geophysics 47(1):47–50. 1982;
18. Villinger H. Solving cylindrical geothermal problems using the Gaver-Stehfest inverse Laplace transform. Geophysics 50(10):1581–1587. 1985;
19. Barnett CT. Simple inversion of time-domain electromagnetic data. Geophysics 49(7):925–933. 1984;
20. Sen MK, Stoffa PL. Global Optimization Methods in Geophysical Inversion 2nd edth ed. Cambridge, UK: Cambridge University Press; 2013.
21. Sharma SP, Kaikkonen P. Global optimisation of time domain electromagnetic data using very fast simulated annealing. Pure and Applied Geophysics 155(1):149–168. 1999;
22. Kirkpatrick S, Gelatt CD, Vecchi MP. Optimization by simulated annealing. Science 220(4598):671–680. 1983;
23. Zhao LS, Sen MK, Stoffa P, Frohlich C. Application of very fast simulated annealing to the determination of the crustal structure beneath Tibet. Geophysical Journal International 125(2):355–370. 1996;
24. Biswas A. Interpretation of residual gravity anomaly caused by simple shaped bodies using very fast simulated annealing global optimization. Geoscience Frontiers 6(6):875–893. 2015;
25. Srigutomo W, Heriyanto M, Aufa MH. Gravity inversion of Talwani model using very fast simulated annealing. Journal of Mathematical and Fundamental Sciences 51(2):177–190. 2019;
26. Srivastava RP, Sen MK. Fractal-based stochastic inversion of poststack seismic data using very fast simulated annealing. Journal of Geophysics and Engineering 6(4):412–425. 2019;
27. Sharma SP. VFSARES: a very fast simulated annealing FORTRAN program for interpretation of 1-D DC resistivity sounding data from various electrode arrays. Computers & Geosciences 42:177–188. 2012;
28. Sharma SP, Biswas A. Interpretation of self-potential anomaly over a 2D inclined structure using very fast simulated-annealing global optimization: an insight about ambiguity. Geophysics 78(3):WB3–WB15. 2013;
29. Netto A, Dunbar J. 3-D constrained inversion of gravimetric data to map the tectonic terranes of south-eastern Laurentia using simulated annealing. Earth and Planetary Science Letters 513:12–19. 2019;
30. Biswas A, Acharya T. A very fast simulated annealing method for inversion of magnetic anomaly over semi-infinite vertical rod-type structure. Modeling Earth Systems and Environment 2(4):1–10. 2016;
31. Ingber L. Very fast simulated re-annealing. Mathematical and Computer Modelling 12(8):967–973. 1989;
32. Sharma SP, Biswas A. Global nonlinear optimization for the estimation of static shift and interpretation of 1-D magnetotelluric sounding data. Annals of Geophysics 54(3):249–264. 2011;
33. Li JH, Zhu ZQ, Feng DS, Xiao JP, Peng LX. Calculation of all-time apparent resistivity of large loop transient electromagnetic method with very fast simulated annealing. Journal of Central South University of Technology 18(4):1235–1239. 2011;
34. Yogi IBS. Implementation of hybrid optimization for time domain electromagnetic 1D inversion. AIP Conference Proceedings 1861article no. 030035. 2017; https://doi.org/10.1063/1.4990922 .
35. Prabawa RS. Comparison of very fast simulated annealing and modified particle swarm optimizaton inversion method for 1-D TDEM data modelling. IOP Conference Series: Earth and Environmental Science 318article no. 012045. 2019; https://doi.org/10.1088/1755-1315/318/1/012045 .
36. Ingber L. Simulated annealing: practice versus theory. Mathematical and Computer Modelling 18(11):29–57. 1993;
37. Stoffa PL, Sen MK. Nonlinear multiparameter optimization using genetic algorithms: inversion of plane-wave seismograms. Geophysics 56(11):1794–1810. 1991;
38. Li X, Koike T, Pathmathevan M. A very fast simulated re-annealing (VFSA) approach for land data assimilation. Computers & Geosciences 30(3):239–248. 2004;
39. Kagiyama T, Utada H, Yamamoto T. Magma ascent beneath Unzen Volcano, SW Japan, deduced from the electrical resistivity structure. Journal of Volcanology and Geothermal Research 89(1–4):35–42. 1999;

Biography

Wahyu Srigutomo is a professor at Physics of Earth and Complex Systems, Physics Department, Faculty of Mathematics and Natural Sciences, Institut Teknologi Bandung, Indonesia. He received his Ph.D. in time-domain electromagnetic modeling and inversion from the Department of Earth and Planetary Science, School of Science, The University of Tokyo, in 2002. His research topics are mainly related to electromagnetic modeling and inversion, geophysical exploration methods, and geothermal and volcanology.

Cahyo Aji Hapsoro received his doctoral degree at the Physics of Earth and Complex Systems, Physics Department, Faculty of Mathematics and Natural Sciences, Institut Teknologi Bandung, in 2020. He mainly focuses on numerical modeling in geophysics. His current research is electromagnetic induction, especially time-domain electromagnetic methods.

Acep Purqon is a lecturer and researcher at the Physics Department, Faculty of Mathematics and Natural Sciences, Institut Teknologi Bandung, Indonesia. He received his Ph.D. in nonlinear physics and complex system modeling at Kanazawa University in 2008. He is currently continuing his research related to the modeling of fluids and other nonlinear physics phenomena, including econophysics.

Warsa is a lecturer and researcher in Applied Geophysics and Exploration, Geophysical Engineering, Faculty of Mining and Petroleum Engineering, Institut Teknologi Bandung, Indonesia. He finished his doctoral degree in the surface magnetic nuclear resonance (SNMR) method at the same affiliation in 2014. Currently, his research topics are related to geoelectromagnetic modeling and inversion.

Doddy Sutarno is a professor of Physics of Earth and Complex Systems at the Physics Department, Faculty of Mathematics and Natural Sciences, Institut Teknologi Bandung (ITB), Indonesia. He received his PhD in electromagnetic induction at Macquarie University in Australia in 1992. He has been working at the Institut Teknologi Bandung, Indonesia, since 1979. His main research topics are the robust processing of electromagnetic data as well as the numerical modeling of electromagnetic methods.

Tsuneomi Kagiyama is currently a professor Emeritus at the Institute for Geothermal Sciences, Graduate School of Science, Kyoto University, Japan. Previously, he worked as Associate Professor at the Earthquake Research Institute (ERI), the University of Tokyo. Most of his research topics deal with volcano eruption prediction and mechanism and physical structure of volcano based on electromagnetic data as well as other geophysical data.

Article information Continued

Fig. 1

TDEM setup for 1D case where the earth’s conductivity varies only in the z-direction.

Fig. 2

A pseudocode for a VFSA algorithm. Modified from Ingber [31] and Sen and Stoffa [20].

Fig. 3

(a) Subsurface resistivity model used for generating the theoretical or synthetic magnetic field Bz. (b) The transient magnetic field Bz, which is generated by performing forward modeling using the subsurface resistivity model described in (a).

Fig. 4

(a) Comparison between the noise-free observed Bz data (in this case, the synthetic data) and the calculated data obtained from the inversion scheme. (b) Comparison between the subsurface resistivity test model or the true model and the inverted model for the noise-free data in (a). (c) Curve of the RMS error convergence toward minimum values vs. the number of iterations.

Fig. 5

(a) Comparison between the 5% noisy synthetic (observed) data and the calculated data; (b) comparison between the test and inverted models for data in (a). (c) and (d) are the same as (a) and (b), respectively, except that the random noise in the synthetic data is 10%. The parameters in both cases are NV = 50, number of iterations = 250, and initial temperature T0 = 5.

Fig. 6

Observed field data and calculated data at Station UZ01 (in pT).

Fig. 7

(a) Inverted subsurface resistivity model from the VFSA inversion scheme for the data in Fig. 6. (b) The resistivity model resulting from the smoothness-constrained Occam’s inversion [8].

Table 1

Combination of parameter values in the VFSA global optimization

NV Iteration T0 RMS misfit (%) Model error (%)
20 100 10 0.3048 54.50
100 5 0.3460 4.50
200 1 0.5035 25.90
250 1 0.3541 38.90
400 5 0.3492 24.90
25 125 5 0.4945 48.90
50 50 1 0.4408 13.30
100 1 0.6080 97.60
150 1 0.1996 4.03
200 1 0.1798 3.70
250 1 0.1629 3.50
150 5 0.1257 1.86
200 5 0.1221 3.28
250 5 0.1011 1.06
300 1 0.5530 62.90