Reduction of Electromagnetic Reflections in 3D Airborne Transient Electromagnetic Modeling: Application of the CFS-PML in Source-Free Media

Abstract

To solve the problem of electromagnetic reflections caused by the termination of finite-difference time-domain (FDTD) grids, we apply the complex frequency-shifted perfectly matched layer (CFS-PML) to airborne transient electromagnetic (ATEM) modeling in a source-free medium. To implement the CFS-PML, two important aspects are improved. First, our method adopts the source-free Maxwell’s equations as the governing equations and introduces the divergence condition, consequently, the discrete form of Maxwell’s third equation is derived with regard to the CFS-PML form. Second, because our method adopts an inhomogeneous time-step, a recursive formula composed of convolution items based on a nonuniform time-step is proposed. The proposed approach is verified via a calculation of the electromagnetic response using homogeneous half-space models with different conductivities. The results show that the CFS-PML can reduce a 60 dB relative errors in late times. Moreover, this approach is also applied to 3D anomalous models; the results indicate that the proposed method can reduce reflections and substantially improve the identification of anomalous bodies. Consequently, the CFS-PML has good implications for ATEM modeling in a source-free medium.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

The airborne transient electromagnetic (ATEM) system, which is an economic alternative for acquiring electromagnetic data with highly efficient detection capabilities and a substantial depth of investigation, has been widely applied to problems associated with hydrogeological surveying, mineral exploration, and environmental monitoring [1–3]. Furthermore, high-accuracy ATEM modeling can provide a theoretical basis for subsequent data inversion and prospective instrument design. Numerous studies have been conducted using the three-dimensional (3D) TEM modeling approach developed by Oristaglio and Hohmann based on finite-difference time-domain (FDTD) grids [4]. More recently, Commer and Newman improved the accelerated simulation scheme for 3D TEM modeling using geometric multigrid concepts [5–7]. Subsequently, Guan calculated the 3D ATEM response based on the graphics processing unit (GPU) [8], and Sun et al. improved upon the 3D FDTD modeling of TEM in consideration of the ramp time [9]. However, boundary reflections still constitute one of the most important challenges for the accuracy of TEM modeling. The most widely used Dirichlet boundary condition (DBC) exhibits better effects only at earlier times before the diffusion field has arrived at the boundary [10]. With an increase in the computing time, the modeling accuracy is progressively affected by electromagnetic reflections at the boundaries. Therefore, it is necessary to develop a more effective boundary condition. Berenger proposed a perfectly matched layer (PML), which can absorb electromagnetic waves with any incident angle and frequency [11]. Subsequently, several approaches, such as the uniaxial perfectly matched layer (U-PML), multiaxial perfectly matched layer (M-PML), and complex frequency-shifted perfectly matched layer (CFS-PML), have been developed with many absorbing boundary conditions based on PML theory [12–14]. In these improvements, the CFS-PML absorbing boundary condition has shown the best performance with regard to the absorption of low-frequency induction fields and late reflections [15, 16]. Roden and Gedney more efficiently implemented the CFS-PML condition by utilizing it with the FDTD method based on recursive convolution; this condition, known as the convolutional perfectly matched layer (C-PML), is highly absorptive of evanescent modes, and it is independent of the host medium [17]. Drossaert and Giannopoulos effectively reduced the computation time and memory by using the complex frequency-shifted stretching function while applying the CFS-PML with the velocity-stress wave equations. Based on the FDTD technique [18], Giannopoulos proposed an electromagnetic modeling method using the CFS-PML with the complex frequency-shifted stretching function and analyzed the effects of the absorption coefficient in different media [19]. Gedney and Zhao applied the CFS-PML boundary condition to the simulation of two-dimensional (2D) electromagnetic waves based on an auxiliary differential equation (ADE), thereby improving the absorbing properties of the CFS-PML in their research [20]. Moreover, Li and Huang used the CFS-PML to calculate the TEM response in a source medium while discretizing the air [21]. Feng et al. reduced the memory requirements in an implementation of the CFS-PML based on the memory-minimized method (Tri-M) and adopted the discrete Zernike transform (DZT) to guarantee the optimal accuracy [22]. Yu et al. improved the C-PML parameters to absorb low-frequency electromagnetic waves in both the ground and the air and observed good absorption [23]. Furthermore, Hu et al. applied the CFS-PML within a fictitious wave domain and discussed the selection of various CFS parameters [24]. Zhao et al. improved the discretization of Maxwell’s divergence equation in the CFS-PML boundary based on a uniform time-step [25]. Feng et al. proposed a Crank–Nicolson cycle-sweep-uniform FDTD method based on CFS-PML and applied it to 3D low-frequency subsurface electromagnetic sensing problems [26].

In this paper, we focus on the problem of electromagnetic reflections at boundaries in source-free media and apply the CFS-PML boundary condition to ATEM modeling based on a nonuniform time-step. Maxwell’s curl equation with a source was used as the governing equation in previous studies about CFS-PML, as the source functions are difficult to approximate and the analytical solution in the source medium is unavailable. In contrast to previous investigations, our research adopts the source-free Maxwell’s equations as the governing equations and the analytical solutions as the initial conditions to perform the iterative calculations; therefore, we can obtain modeling results with a high accuracy. Because Maxwell’s divergence equation must be incorporated into the governing equations to ensure the uniqueness of the solution [10], one of the fundamental aspects of our approach is the discretization of Maxwell’s divergence equation in the stretched coordinate space. Another crucial component of the proposed method is the derivation of a recursive convolution formula based on nonuniform time-steps in the discretized stretched coordinate formulation. Finally, the proposed method is verified via a calculation of the electromagnetic response using homogeneous half-space models and anomalous models with 3D bodies, thereby effectively demonstrating the performance of the proposed method with regard to the absorption of electromagnetic reflections during the modeling of ATEM.

2. Methods

Maxwell’s divergence equation is employed as the governing equation in FDTD-based ATEM modeling in a source-free medium to advance the magnetic field of H_z, while the other fields are advanced using Maxwell’s curl equations. Following the method proposed by Wang and Hohmann [10], the initial conditions are calculated based on a homogeneous half-space. The staggered Yee grid is used to discretize the earth model inhomogeneously. In addition, the time-steps are advanced using a modified Du Fort-Frankel method [4, 27]. The ATEM model contains a ground-air boundary and a computational domain boundary; therefore, we adopt an upward continuation at the ground-air boundary and the CFS-PML as the computational domain boundary condition. However, two problems are encountered during the application of the CFS-PML: we must derive the discrete form of Maxwell’s divergence equation in the stretched coordinate, and the recursive convolution formulation based on a nonuniform time-step must be derived. Consequently, these key issues in the application of the CFS-PML boundary condition are discussed in this paper.

2.1. Governing Equations in the Stretched Coordinate Space

In the modeling of ATEM in a source-free medium, Maxwell’s divergence equation must be included within the governing equations during the modeling of ATEM in a source-free medium to ensure the uniqueness and stability of the solution. Therefore, Maxwell’s divergence equation (equation (3)) is adopted to advance the magnetic field $B_z$ in the stretched coordinate space [28–32], while the other electromagnetic fields are advanced using Maxwell’s curl equations ((1) and (2)). $$\begin{array}{}\text{(1)}& {\nabla }_s\times \mathit{H}=\epsilon j\omega \mathit{E}+\sigma \mathit{E},\text{(2)}& {\nabla }_s\times \mathit{E}=-\mu j\omega \mathit{H},\text{(3)}& {\nabla }_s·\mathit{H}=0,\end{array}$$where $E$ is the electric field, $H$ is the magnetic field, $j\omega$ is the complex frequency variable, $\mu$ is the permeability, $\epsilon$ is the permittivity, and $\sigma$ is the electrical conductivity.

In the stretched coordinate space, the Hamilton operator (${\nabla }_s$) can be expressed as follows [33]: $$\begin{array}{}\text{(4)}& {\nabla }_s=\overrightarrow{x}\frac{1}{s_x}\frac{\partial }{\partial x}+\overrightarrow{y}\frac{1}{s_y}\frac{\partial }{\partial y}+\overrightarrow{z}\frac{1}{s_z}\frac{\partial }{\partial z},\end{array}$$where $s_i=k_i+\left({\sigma }_{pi}/\left({\alpha }_i+j\omega \mu \right)\right)$, $k_i$ and ${\alpha }_i$ are positive and less than 1 and ${\sigma }_{pi}$ is the conductivity of the CFS-PML layers ($i=x,y,z$).

The difference scheme of Maxwell’s curl equations in the stretched coordinate space was developed by Roden and Gedney similar to the application of the CFS-PML [17]. Thus, in this paper, we mainly deduce the discrete form of Maxwell’s divergence equation (equation (3)).

Inserting (4) into (3) leads to $$\begin{array}{}\text{(5)}& \frac{1}{s_x}\frac{\partial {\mathit{H}}_x}{\partial x}+\frac{1}{s_y}\frac{\partial {\mathit{H}}_{\mathit{y}}}{\partial y}+\frac{1}{s_z}\frac{\partial {\mathit{H}}_{\mathit{z}}}{\partial z}=0.\end{array}$$

Equation (5) is next transformed into the time domain to obtain the renewal equation of H_z, after which the expression of S_i (i = x, y, z) is inserted into (5), which can consequently be expressed as $$\begin{array}{}\text{(6)}& \frac{1}{k_x}\frac{\partial {\mathit{H}}_x}{\partial x}+\frac{1}{k_y}\frac{\partial {\mathit{H}}_y}{\partial y}+\frac{1}{k_z}\frac{\partial {\mathit{H}}_z}{\partial z}+{\psi }_{hzx}+{\psi }_{hzy}+{\psi }_{hzz}=0,\end{array}$$where ${\psi }_{hzi}$ is the auxiliary expression of a convolution item and is implemented as $$\begin{array}{}\text{(7)}& {\psi }_{hzi}=\left(-\frac{{\sigma }_{pi}}{\epsilon k_i^2}\exp \left(-\left(\frac{{\sigma }_{pi}}{\epsilon k_i}+\frac{{\alpha }_i}{\epsilon }\right)t\right)\right)\ast \frac{\partial {\mathit{H}}_i}{\partial i}.\end{array}$$

The governing equations of other fields in the stretched coordinate can be similarly derived.

2.2. The Discrete Form of Maxwell’s Divergence Equation in the Stretched Coordinate Space

During ATEM modeling, the time-step of an FDTD method usually increases gradually according to the attenuation characteristics of TEM responses (i.e., the TEM field exhibits relatively sharp variations at earlier times and gradually becomes smooth thereafter). Here, the time-step is advanced using a modified Du Fort-Frankel method. The traditional recursive convolution formula is deduced based on a uniform time-step; therefore, we need to modify the derivation procedure based on inhomogeneous time-steps which is shown in the appendix. In the stretched coordinate space, ${\psi }_{hzx}$ at $t_{n+1/2}$ can be expressed as $$\begin{array}{}\text{(8)}& {\psi }_{hzx}\left(n+\frac{1}{2}\right)=\left\{-\frac{{\sigma }_{px}}{\epsilon k_x^2}\exp \left[-\left(\frac{{\sigma }_{px}}{\epsilon k_x}+\frac{{\alpha }_x}{\epsilon }\right)t_{n+1/2}\right]\right\}\ast \frac{\partial {\mathit{H}}_x^{n+1/2}}{\partial x}={\xi }_{hx}\left(t_{n+1/2}\right)\ast \frac{\partial {\mathit{H}}_x^{n+1/2}}{\partial x}.\end{array}$$

To obtain the discrete form of the convolution formulation, we assume that $\Delta t_n=t_{n+1/2}-t_{n-1/2}$ ($n=0,1,2,\dots ,\mathrm{n}$), and thus, equation (8) can be rewritten as $$\begin{array}{}\text{(9)}& {\psi }_{hzx}\left(n+\frac{1}{2}\right)={\xi }_{hx}\left(t_{n+1/2}\right)\ast \frac{\partial {\mathit{H}}_x^{n+1/2}}{\partial x}={\int }_0^{t_{n+1/2}}{\xi }_{hx}\left(\tau \right)\frac{\partial {\mathit{H}}_x\left(t_{n+1/2}-\tau \right)}{\partial x}d\tau ,\end{array}$$where ${\xi }_{hx}\left(t\right)=-A_x{e}^{-B_{x}t}$, $A_i=\left({\sigma }_{pi}/\epsilon k_i^2\right)$, $B_i=\left({\sigma }_{pi}/\epsilon k_i\right)+\left({\alpha }_i/\epsilon \right)$ ($i=x,y,z$).

Then ${\psi }_{hzx}\left(n+\left(1/2\right)\right)$ can be performed recursively: $$\begin{array}{}\text{(10)}& {\psi }_{hzx}\left(n+\frac{1}{2}\right)=\frac{1}{2}\left(\frac{A_x}{B_x}{e}^{-B\mathrm{\Delta }{t}_n}-\frac{A_x}{B_x}\right)\left(\frac{\partial H_x^{n+1/2}}{\partial x}+\frac{\partial H_x^{n-1/2}}{\partial x}\right)+e^{-B_x\mathrm{\Delta }{t}_n}{\psi }_{hzx}\left(n-\frac{1}{2}\right).\end{array}$$

We can obtain ${\psi }_{hzz}\left(n+\left(1/2\right)\right)$ and ${\psi }_{hz\mathrm{x}}\left(n+\left(1/2\right)\right)$ in the same manner: $$\begin{array}{}\text{(11)}& {\psi }_{hzz}\left(n+\frac{1}{2}\right)=\frac{1}{2}\left(\frac{A_z}{B_z}{e}^{-B_z\mathrm{\Delta }{t}_n}-\frac{A_z}{B_z}\right)\left(\frac{\partial H_z^{n+1/2}}{\partial z}+\frac{\partial H_z^{n-1/2}}{\partial z}\right)+e^{-B_z\mathrm{\Delta }{t}_n}{\psi }_{hzz}\left(n-\frac{1}{2}\right),\text{(12)}& {\psi }_{hzx}\left(n+\frac{1}{2}\right)=\frac{1}{2}\left(\frac{A_x}{B_x}{e}^{-B_x\mathrm{\Delta }{t}_n}-\frac{A_x}{B_x}\right)\left(\frac{\partial H_x^{n+1/2}}{\partial x}+\frac{\partial H_x^{n-1/2}}{\partial x}\right)+e^{-B_x\mathrm{\Delta }{t}_n}{\psi }_{hzx}\left(n-\frac{1}{2}\right).\end{array}$$

As shown in (6) and (11), ${\psi }_{hzz}$ contains H_z at the current time in the discretization of Maxwell’s divergence equation to obtain the iterative formula of ${\mathit{H}}_{\mathrm{z}}$ at the time $t_{n+1/2}$. Therefore, the iterative formula of ${\mathit{H}}_{\mathrm{z}}$ requires further derivation. Equation (6) can be discretized in both space and time according to the staggered Yee scheme as $$\begin{array}{}\text{(13)}& \frac{1}{k_x}\left(\frac{{\mathit{H}}_x^{n+1/2}\left(i+1,j+1/2,k+1/2\right)-{\mathit{H}}_x^{n+1/2}\left(i,j+1/2,k+1/2\right)}{\Delta x}\right)+\frac{1}{k_y}\left(\frac{{\mathit{H}}_y^{n+1/2}\left(i+1/2,j+1,k+1/2\right)-{\mathit{H}}_y^{n+1/2}\left(i+1/2,j,k+1/2\right)}{\Delta y}\right)+\frac{1}{k_z}\left(\frac{{\mathit{H}}_z^{n-1+1/2}\left(i+1/2,j+1/2,k+1\right)-{\mathit{H}}_z^{n-1+1/2}\left(i+1/2,j+1/2,k\right)}{\Delta z}\right)+{\psi }_{hzx}^{n+1/2}+{\psi }_{hzy}^{n+1/2}+{\psi }_{hzz}^{n+1/2}=0.\end{array}$$

Inserting (10)–(12) into (13), we obtain $$\begin{array}{}\text{(14)}& {\mathit{H}}_z^{n+1/2}\left(i+\frac{1}{2},j+\frac{1}{2},k\right)={\mathit{H}}_z^{n+1/2}\left(i+\frac{1}{2},j+\frac{1}{2},k+1\right)+\frac{C_z{D}_z}{\Delta z\left(i+1/2,j+1/2,k\right)}\left({\mathit{H}}_z^{n-1+1/2}\left(i+\frac{1}{2},j+\frac{1}{2},k+1\right)-{\mathit{H}}_z^{n-1+1/2}\left(i+\frac{1}{2},j+\frac{1}{2},k\right)\right)+\frac{D_z}{k_y\Delta y\left(i+1/2,j,k+1/2\right)}\left({\mathit{H}}_y^{n+1/2}\left(i+\frac{1}{2},j+1,k+\frac{1}{2}\right)-{\mathit{H}}_y^{n+1/2}\left(i+\frac{1}{2},j,k+\frac{1}{2}\right)\right)+\frac{D_z}{k_x\Delta x\left(i,j+1/2,k+1/2\right)}\left({\mathit{H}}_x^{n+1/2}\left(i+1,j+\frac{1}{2},k+\frac{1}{2}\right)-{\mathit{H}}_x^{n+1/2}\left(i,j+\frac{1}{2},k+\frac{1}{2}\right)\right)+D_z{\psi }_{hzx}^{n+1/2}+D_z{\psi }_{hzy}^{n+1/2}+e^{-B_z\Delta t_n}{D}_z{\psi }_{hzz}^{n-1/2},\end{array}$$where $C_z=1/2\left(\left(A_z/B_z\right)e^{-B_z\mathrm{\Delta }{t}_n}-\left(A_z/B_z\right)\right)$, $D_z=\left(k_z/\left(1+k_z{C}_z\right)\right)\Delta z\left(i+\left(1/2\right),j+\left(1/2\right),k\right)$.

Note that $H_{\mathrm{z}}$ at $t_{n+1/2}$ is related to the value of ${\psi }_{hzz}$ of the previous time-step, and thus, we must refresh ${\psi }_{hzz}$ at the end of the calculation for the next iteration.

3. Numerical Tests

To investigate the feasibility of the proposed method, we first compare the solutions of our method with the numerical solutions calculated by the integral method in homogeneous half-space models [8, 34, 35]. Then, the proposed method is applied to several models with 3D conductors. With the exception of the slanted model, every model has 101 × 101 × 50 grids [36]. The grids are nonuniform with a minimum spacing of 10 m and a maximum spacing of 120 m. The spacing of the CFS-PML boundary is 120 m, and the boundary is not contained within any of the abovementioned grids. A transmitting coil is located at the center of the earth model with a height of 120 m. The magnetic moment is 4π × 10⁻⁷A·m², and the transmitter current is 0.7 × 10⁷A. The receiving coil is situated at a height of 60 m, and it is 130 m away from the transmitting coil in the x-direction.

3.1. Validation with Homogeneous Half-Space Models

The electromagnetic responses of the homogeneous half-space models are calculated, and the FDTD solutions are compared with the numerical solutions to verify the accuracy of our method. To demonstrate the absorbing effect of the proposed method in different host media, the conductivities are set as 0.1 S/m, 0.01 S/m, and 0.005 S/m; the diagrammatic drawing of the homogeneous half-space model is shown in Figure 1.

[figure omitted; refer to PDF]

The FDTD-based solutions (calculated using the DBC and the CFS-PML with 8 layers) are compared with the numerical solutions, as shown in Figure 2. In the models with different conductivities, we can discern that these solutions are in good agreement at earlier times because the diffusion fields have not reached the boundaries yet. With an increase in the time; however, the effectiveness of the CFS-PML gradually becomes clearer, especially in the models with a smaller conductivity (as presented in Figure 2 b~c). The average relative errors of the FDTD solutions using the CFS-PML and DBC are shown in Figure 2(d). In the model with a large conductivity (0.1 S/m) and a slower propagation speed, these boundary conditions both demonstrate a better absorption of electromagnetic waves with average relative errors of 1.61% and 0.79% (corresponding to the DBC and CFS-PML, respectively). When the model conductivity is 0.01 S/m, the average relative errors of the DBC and CFS-PML are 21.47% and 2.63%, respectively, effectively revealing the better ability of the CFS-PML to absorb electromagnetic waves. In addition, when the ground conductivity is 0.005 S/m, the strong electromagnetic wave reflected from the DBC leads to an enormous deviation in the numerical solution, whereas the CFS-PML still exhibits good absorption; in this case, the average relative errors of the DBC and CFS-PML are 42.9% and 4.17%, respectively.

[figures omitted; refer to PDF]

Except for the decay curves at the center of the receiving coil, we also investigate the electromagnetic reflections at the boundaries. Therefore, time slices of the electromagnetic responses from a homogeneous half-space model with a conductivity of 0.01 S/m are presented in Figure 3. As illustrated, the electromagnetic responses calculated by the DBC and CFS-PML boundaries are compared at three separate moments. At earlier times, the electromagnetic responses of the DBC and CFS-PML boundaries display the same variations; meanwhile, at later times, serious distortions are observed at both 4.5 ms and 9 ms when the DBC is used. In contrast, no distortions are evident when the CFS-PML is adopted regardless of the point in time, indicating that the CFS-PML boundary condition demonstrates a better ability to absorb electromagnetic waves.

[figures omitted; refer to PDF]

To study the reflection errors, the time-dependent field at the center of the receiving coil is recorded when the conductivity of the host medium is 0.01 S/m. The error relative to the numerical solution (in dB) is computed using: $$\begin{array}{}\text{(15)}& {\mathrm{E}\mathrm{r}\mathrm{r}\mathrm{o}\mathrm{r}}_{\mathrm{d}\mathrm{B}}=20{\log }_{10}\left|\frac{V_z\left(t\right)-V_{z\mathrm{r}\mathrm{e}\mathrm{f}}\left(t\right)}{V_{z\mathrm{r}\mathrm{e}\mathrm{f}\max }}\right|,\end{array}$$where $V_z\left(t\right)$ represents the voltage received by the receiving coil, V_zref(t) represents the numerical solutions, and V_zrefmax represents the maximum of V_zref(t). The relative errors (in dB) computed via (13) are recorded in Figure 4, from which we can see that the relative errors are similar at earlier times. Moreover, compared with the DBC, the CFS-PML can reduce the relative error by nearly 60 dB at later times.

[figure omitted; refer to PDF]

3.2. Validation with Anomalous Models

To test the effectiveness of the proposed method, we design three anomalous models. In each model, the ground conductivity is 0.01 S/m, and the thickness of the CFS-PML is 10 cells. The first model is shown in Figure 5; two anomalous bodies with different sizes, depths, and conductivities are selected. An analytical solution exists only in the homogeneous half-space. Therefore, to study the reflection errors due to the boundaries, the mesh is extended by an additional 40 cells in each direction to obtain the reference solution, and the CFS-PML is used to minimize the reflections. To ensure that the electromagnetic waves cannot reach the boundaries during the observation time (10 ms), the diffusion depth is set as $$\begin{array}{}\text{(16)}& d=\sqrt{2\frac{t}{{\sigma \mu }_0}},\end{array}$$where $\sigma$ is the conductivity, $t$ is the observation time, and $d$ is the diffusion depth [37].

[figures omitted; refer to PDF]

Time slices of the anomaly responses are shown in Figure 6. Both at earlier times (3 ms) and at later times (10 ms); reflected electromagnetic waves exist in the model when the DBC is adopted. Moreover, especially at 10 ms, the locations of anomalous bodies cannot be identified accurately. In contrast, the reflections are absorbed well in the model with the CFS-PML, and the results are in agreement with the reference solution. A comparison of $\partial \mathrm{B}/\partial \mathrm{t}$ at the center of the lattice along the x-direction between 3 ms and 10 ms is shown in Figure 7. In earlier times, the DBC solution is in good agreement with the reference solution at the center of the model; however, such good agreement is not observed near the boundary. Furthermore, at 10 ms, the DBC solution severely deviates from the reference solution along the whole boundary. Meanwhile, the responses in the model using the CFS-PML are consistent with the reference solution regardless of the time. The relative errors are shown in Figures 7(c) and 6(d). The maximum relative error is less than 3% at earlier times and no more than 4.6% at later times when the CFS-PML is adopted. However, the results obtained with the DBC are consistent with the reference solution only at the center of the model area at earlier times.

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

For further verification, an anomalous model is designed and presented in Figure 8. In this model, the anomalous body adjoins the boundary. The dimensions of the anomalous body are 480 m × 120 m × 180 m, and the depth of the body is 190 m. Time slices of the anomaly responses are shown in Figure 9. Intense reflections of electromagnetic waves appear when the DBC is utilized, and the information of the anomalous body cannot be retrieved accurately. In contrast, through the implementation of the CFS-PML boundary condition, the reflections are effectively suppressed, and the results are approximately equivalent to the reference solution. Consequently, the information of the anomalous body is presented clearly. The responses at the boundary at the center of the model along the x-direction are shown in Figure 10. The DBC solution evidently coincides with the reference solution at the center of the boundary, although the DBC solution is different near the boundary at earlier times and performs progressively worse with increasing time. Meanwhile, the responses obtained using the CFS-PML are still in good agreement. The maximum relative error of the results acquired using the CFS-PML is less than 5% at later times. Furthermore, the DBC solutions exhibit serious reflections, especially near the anomalous body. From Figures 8 and 9, we can conclude that the CFS-PML is very stable and can be established very close to the anomalous body.

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

Finally, we design an oblique model. The anomalous body is oblique with an angle of 45°, and it is established at the center of the x-y plane, as shown in Figure 11. The dimensions of the anomalous body are 400 m × 610 m × 90 m, and the depth is 350 m. Time slices of the anomaly responses are shown in Figure 12. Reflections of electromagnetic waves evidently occur along the boundary at earlier times when the DBC is adopted, while the responses at the center are basically equivalent to those in the reference solution. At 10 ms, the reflections at the DBC boundaries are so intense that the diffusion is completely distorted. In contrast, the responses with the CFS-PML boundary condition are highly consistent with the reference solution, even at later times. The received ∂B/∂t at the boundary across the center of the receiving coil in the x-direction is shown in Figure 13. It is clear that the DBC solution coincides with the reference solution at the center of the model area at earlier times and seriously deviates from the reference solution at later times. However, the CFS-PML can effectively absorb the reflections from the boundary, even at later times, and the maximum relative error is less than 5%.

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

4. Conclusions and Perspectives

We applied the CFS-PML boundary condition to the modeling of ATEM data in a source-free medium and introduced the source-free Maxwell’s equations as the governing equations. Furthermore, we provided discrete forms of the governing equations in a CFS-PML formulation and deduced the recursive convolution formula with a nonuniform time-step in the stretched coordinate space. The validated results of a homogeneous half-space model demonstrate that the proposed method can absorb incident wave better, and it can reduce the reflections at later times. The CFS-PML is highly absorptive and can considerably reduce the relative error with a 60 dB improvement compared with the traditional boundary condition (i.e., the DBC) at later times.

Finally, we employed the CFS-PML in airborne anomalous models and compared the electromagnetic responses with the reference solution. The results show that the CFS-PML demonstrates better absorptive capabilities and can accurately retrieve the information of the anomalous body regardless of whether the anomalous body is set obliquely or near the boundary of the model area. Moreover, the maximum relative error of the improved method is less than 5% at later times (10 ms). The proposed method also displays a better recognition precision in the presence of multiple or complex 3D conductive bodies, thereby verifying the effectiveness of this method for ATEM modeling. Consequently, the proposed method can provide higher-precision simulation results for deep geological exploration.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

Project supported by the National Natural Science Foundation of China (Grant nos. 41674109, 41604084), the Instrument Developing Project of the Chinese Academy of Sciences (Grant no. ZDYZ2012-1-03), and the Key Laboratory of Electromagnetic Radiation and Detection Technology of the Chinese Academy of Sciences.

Appendix

Appendix

To describe H_x accurately, the integral cannot be divided as usual into [t₀, t₁], [t₁, t₂], … , [t_n, t_n+1/2], because the time-steps are inhomogeneous. Here, we propose a new subdivision, as shown in $$\begin{array}{}\text{(17)}& {\psi }_{hzx}\left(n+\frac{1}{2}\right)={\int }_0^{\Delta t_n}{\xi }_{hx}\left(\tau \right)\frac{\partial {\mathit{H}}_x\left(t_{n+1/2}-\tau \right)}{\partial x}d\tau +{\int }_{\Delta t_n}^{\Delta t_n+\Delta t_{n-1}}{\xi }_{hx}\left(\tau \right)\frac{\partial {\mathit{H}}_{\mathit{x}}\left(t_{n+1/2}-\tau \right)}{\partial x}d\tau +{\int }_{\Delta t_n+\Delta t_{n-1}}^{\Delta t_n+\Delta t_{n-1}+\Delta t_{n-2}}{\xi }_{hx}\left(\tau \right)\frac{\partial {\mathit{H}}_x\left(t_{n+1/2}-\tau \right)}{\partial x}d\tau +\cdots +{\int }_{\Delta t_n+\Delta t_{n-1}+\cdots +\Delta t_2}^{\Delta t_n+\Delta t_{n-1}+\cdots +\Delta t_2+\Delta t_1}{\xi }_{hx}\left(\tau \right)\frac{\partial {\mathit{H}}_x\left(t_{n+1/2}-\tau \right)}{\partial x}d\tau +{\int }_{\left(\Delta t_n+\Delta t_{n-1}+\cdots +\Delta t_1\right)}^{t_{n+1/2}}{\xi }_{hx}\left(\tau \right)\frac{\partial {\mathit{H}}_x\left(t_{n+1/2}-\tau \right)}{\partial x}d.\end{array}$$

The further calculation of (17) leads to $$\begin{array}{}\text{(18)}& {\psi }_{hzx}\left(n+\frac{1}{2}\right)=\frac{1}{2}\left(\frac{A_x}{B_x}{e}^{-B\mathrm{\Delta }{t}_n}-\frac{A_x}{B_x}\right)\left(\frac{\partial H_x^{n+1/2}}{\partial x}+\frac{\partial H_x^{n-1/2}}{\partial x}\right)+\frac{1}{2}\left(\frac{A_x}{B_x}{e}^{-B\left(\mathrm{\Delta }{t}_n+\mathrm{\Delta }{t}_{n-1}\right)}-\frac{A_x}{B_x}{e}^{-B\mathrm{\Delta }{t}_n}\right)\left(\frac{\partial H_x^{n-1/2}}{\partial x}+\frac{\partial H_x^{n-1-1/2}}{\partial x}\right)++\frac{1}{2}\left(\frac{A_x}{B_x}{e}^{-B\left(\mathrm{\Delta }{t}_n+\mathrm{\Delta }{t}_{n-1}+\cdots +\mathrm{\Delta }{t}_1\right)}-\frac{A_x}{B_x}{e}^{-B\left(\mathrm{\Delta }{t}_n+\mathrm{\Delta }{t}_{n-1}+\cdots +\mathrm{\Delta }{t}_2\right)}\right)\left(\frac{\partial H_x^{1/2}}{\partial x}+\frac{\partial H_x^{1+1/2}}{\partial x}\right)+\frac{1}{2}\left(\frac{A_x}{B_x}{e}^{-B\left(t_{n+1/2}\right)}-\frac{A_x}{B_x}{e}^{-B\left(\mathrm{\Delta }{t}_n+\mathrm{\Delta }{t}_{n-1}+\cdots +\mathrm{\Delta }{t}_1\right)}\right)\left(\frac{\partial H_x^0}{\partial x}+\frac{\partial H_x^{1/2}}{\partial x}\right),\end{array}$$where ${\xi }_{hx}\left(t\right)=-A_x{e}^{-B_{x}t}$, $A_i=\left({\sigma }_{pi}/\epsilon k_i^2\right)$, $B_i=\left({\sigma }_{pi}/\epsilon k_i\right)+\left({\alpha }_i/\epsilon \right)$ (i = x, y, z).

From (18), we can obtain ${\psi }_{hzx}\left(n-\left(1/2\right)\right)$ as $$\begin{array}{}\text{(19)}& {\psi }_{hzx}\left(n-\frac{1}{2}\right)=\frac{1}{2}\left(\frac{A_x}{B_x}{e}^{-B\mathrm{\Delta }{t}_{n-1}}-\frac{A_x}{B_x}\right)\left(\frac{\partial H_x^{n-1/2}}{\partial x}+\frac{\partial H_x^{n-1-1/2}}{\partial x}\right)+\frac{1}{2}\left(\frac{A_x}{B_x}{e}^{-B\left(\mathrm{\Delta }{t}_{n-1}+\mathrm{\Delta }{t}_{n-2}\right)}-\frac{A_x}{B_x}{e}^{-B\mathrm{\Delta }{t}_{n-1}}\right)\left(\frac{\partial H_x^{n-1-1/2}}{\partial x}+\frac{\partial H_x^{n-2-1/2}}{\partial x}\right)+\cdots +\frac{1}{2}\left(\frac{A_x}{B_x}{e}^{-B\left(\mathrm{\Delta }{t}_{n-1}+\mathrm{\Delta }{t}_{n-2}+\cdots +\mathrm{\Delta }{t}_1\right)}-\frac{A_x}{B_x}{e}^{-B\left(\mathrm{\Delta }{t}_{n-1}+\mathrm{\Delta }{t}_{n-2}+\cdots +\mathrm{\Delta }{t}_2\right)}\right)\left(\frac{\partial H_x^{1/2}}{\partial x}+\frac{\partial H_x^{1+1/2}}{\partial x}\right)+\frac{1}{2}\left(\frac{A_x}{B_x}{e}^{-B\left(t_{n-1/2}\right)}-\frac{A_x}{B_x}{e}^{-B\left(\mathrm{\Delta }{t}_{n-1}+\mathrm{\Delta }{t}_{n-2}+\cdots +\mathrm{\Delta }{t}_1\right)}\right)\left(\frac{\partial H_x^0}{\partial x}+\frac{\partial H_x^{1/2}}{\partial x}\right).\end{array}$$

From (18) and (19) we can get $$\begin{array}{}\text{(20)}& {\psi }_{hzx}\left(n+\frac{1}{2}\right)=\frac{1}{2}\left(\frac{A_x}{B_x}{e}^{-B\mathrm{\Delta }{t}_n}-\frac{A_x}{B_x}\right)\left(\frac{\partial H_x^{n+1/2}}{\partial x}+\frac{\partial H_x^{n-1/2}}{\partial x}\right)+e^{-B_x\mathrm{\Delta }{t}_n}{\psi }_{hzx}\left(n-\frac{1}{2}\right).\end{array}$$

Note that within the integral interval $\tau \in \left(0,\Delta t_n\right)$ in equation (17), ${\mathit{H}}_x\left(t_{n+1/2}-\tau \right)$ has two values at ${\mathit{H}}_x^{n+1/2}$ and ${\mathit{H}}_x^{n-1/2}$, and we approximate $\left(\partial {\mathit{H}}_x\left(t_{n+1/2}-\tau \right)\right)/\partial x$ as $1/2\left(\left(\partial {\mathit{H}}_x^{n+1/2}/\partial x\right)+\left(\partial {\mathit{H}}_x^{n-1/2}/\partial x\right)\right)$.

Data Availability

The (MATLAB.mat) data used to support the findings of this study are available from the corresponding author upon request.

References

[1] Y. Ji, H. Luan, S. Y. Li, L. Wan, Y. Wang, Y. C. Xu, L. Li, J. Lin, "Resolution of full-waveform airborne TEM," Journal of Jilin University (Earth Science Edition), vol. 41 no. 3, pp. 885-891, 2011.

[2] J. M. Legault, "Airborne electromagnetic systems—state of the art and future directions," CSEG Recorder, vol. 40 no. 6, pp. 38-49, 2015.

[3] Y. Qi, L. Huang, X. Wang, G. Fang, G. Yu, "Airborne transient electromagnetic modeling and inversion under full attitude change," IEEE Geoscience and Remote Sensing Letters, vol. 14 no. 9, pp. 1575-1579, DOI: 10.1109/LGRS.2017.2724558, 2017.

[4] M. L. Oristaglio, G. W. Hohmann, "Diffusion of electromagnetic fields into a two-dimensional earth: a finite-difference approach," Geophysics, vol. 49 no. 7, pp. 870-894, DOI: 10.1190/1.1441733, 1984.

[5] M. Commer, G. Newman, "A parallel finite-difference approach for 3D transient electromagnetic modeling with galvanic sources," Geophysics, vol. 69 no. 5, pp. 1192-1202, DOI: 10.1190/1.1801936, 2004.

[6] M. Commer, G. A. Newman, "An accelerated time domain finite difference simulation scheme for three-dimensional transient electromagnetic modeling using geometric multigrid concepts," Radio Science, vol. 41 no. 3,DOI: 10.1029/2005rs003413, 2006.

[7] M. Commer, G. M. Hoversten, E. S. Um, "Transient-electromagnetic finite-difference time-domain earth modeling over steel infrastructure," Geophysics, vol. 80 no. 2, pp. E147-E162, DOI: 10.1190/geo2014-0324.1, 2015.

[8] S. Guan, Parallel calculation of 3D FDTD helicopter transient electromagnetic response based on the GPU, [Doctoral dissertation], 2012.

[9] H. Sun, X. Li, S. Li, Z. Qi, W. Wan, M. SU, Y. Xue, B. Liu, "Three-dimensional FDTD modeling of TEM excited by a loop source considering ramp time," Chinese Journal of Geophysics, vol. 56 no. 3, pp. 1049-1064, 2013.

[10] T. Wang, G. W. Hohmann, "A finite-difference, time-domain solution for three-dimensional electromagnetic modeling," Geophysics, vol. 58 no. 6, pp. 797-809, DOI: 10.1190/1.1443465, 1993.

[11] J. P. Berenger, "A perfectly matched layer for the absorption of electromagnetic waves," Journal of Computational Physics, vol. 114 no. 2, pp. 185-200, DOI: 10.1006/jcph.1994.1159, 1994.

[12] S. D. Gedney, "An anisotropic perfectly matched layer-absorbing medium for the truncation of FDTD lattices," IEEE Transactions on Antennas and Propagation, vol. 44 no. 12, pp. 1630-1639, DOI: 10.1109/8.546249, 1996.

[13] K. C. Meza-Fajardo, A. S. Papageorgiou, "A nonconvolutional split-field perfectly matched layer for wave propagation in isotropic and anisotropic elastic media stability analysis," Bulletin of the Seismological Society of America, vol. 98 no. 4, pp. 1811-1836, DOI: 10.1785/0120070223, 2008.

[14] R. Martin, D. Komatitsch, S. D. Gedney, "A variational formulation of a stabilized unsplit convolutional perfectly matched layer for the isotropic or anisotropic seismic wave equation," Computer Modeling in Engineering and Sciences, vol. 37 no. 3, pp. 274-304, 2008.

[15] J. P. Berenger, "An optimized CFS-PML for wave–structure interaction problems," IEEE Transactions on Electromagnetic Compatibility, vol. 54 no. 2, pp. 351-358, DOI: 10.1109/TEMC.2011.2178852, 2012.

[16] A. S. Moura, R. R. Saldanha, E. J. Silva, M. F. Pantoja, A. C. Lisboa, W. G. Facco, "Discretization of the CFS-PML for computational electromagnetics using discrete differential forms," Microwave and Optical Technology Letters, vol. 55 no. 2, pp. 351-357, DOI: 10.1002/mop.27298, 2013.

[17] J. A. Roden, S. D. Gedney, "Convolution PML (CPML): an efficient FDTD implementation of the CFS-PML for arbitrary media," Microwave and Optical Technology Letters, vol. 27 no. 5, pp. 334-339, DOI: 10.1002/1098-2760(20001205)27:5<334::AID-MOP14>3.0.CO;2-A, 2000.

[18] F. H. Drossaert, A. Giannopoulos, "Complex frequency shifted convolution PML for FDTD modelling of elastic waves," Wave Motion, vol. 44 no. 7-8, pp. 593-604, DOI: 10.1016/j.wavemoti.2007.03.003, 2007.

[19] A. Giannopoulos, "An improved new implementation of complex frequency shifted PML for the FDTD method," IEEE Transactions on Antennas and Propagation, vol. 56 no. 9, pp. 2995-3000, DOI: 10.1109/TAP.2008.928789, 2008.

[20] S. D. Gedney, B. Zhao, "An auxiliary differential equation formulation for the complex-frequency shifted PML," IEEE Transactions on Antennas and Propagation, vol. 58 no. 3, pp. 838-847, DOI: 10.1109/TAP.2009.2037765, 2010.

[21] Z. H. Li, Q. H. Huang, "Application of the complex frequency shifted perfectly matched layer absorbing boundary conditions in transient electromagnetic method modeling," Chinese Journal of Geophysics-Chinese, vol. 57 no. 4, pp. 1292-1299, 2014.

[22] N. Feng, Y. Yue, Q. H. Liu, "Direct Z-transform implementation of the CFS-PML based on memory-minimized method," IEEE Transactions on Microwave Theory and Techniques, vol. 63 no. 3, pp. 877-882, DOI: 10.1109/TMTT.2015.2389218, 2015.

[23] X. Yu, X. Wang, X. Li, X. Lin, F. Yang, M. Tang, "Three-dimensional finite difference forward modeling of the transient electromagnetic method in the time domain," Chinese Journal of Geophysics-Chinese Edition, vol. 60 no. 2, pp. 810-819, 2017.

[24] Y. Hu, G. Egbert, Y. Ji, G. Fang, "A novel CFS-PML boundary condition for transient electromagnetic simulation using a fictitious wave domain method," Radio Science, vol. 52 no. 1, pp. 118-131, DOI: 10.1002/2016RS006160, 2017.

[25] X. Zhao, Y. Ji, S. Qiu, S. Guan, Y. Wu, "The CFS-PML in numerical simulation of ATEM," Proceedings of SPIE 10322 Seventh International Conference on Electronics and Information Engineering,DOI: 10.1117/12.2265965, .

[26] N. Feng, Y. Zhang, Q. Sun, J. Zhu, W. T. Joines, Q. H. Liu, "An accurate 3-D CFS-PML based crank–Nicolson FDTD method and its applications in low-frequency subsurface sensing," IEEE Transactions on Antennas and Propagation, vol. 66 no. 6, pp. 2967-2975, DOI: 10.1109/TAP.2018.2816788, 2018.

[27] J. I. Adhidjaja, G. W. Hohmann, "A finite-difference algorithm for the transient electromagnetic response of a three-dimensional body," Geophysical Journal International, vol. 98 no. 2, pp. 233-242, DOI: 10.1111/j.1365-246X.1989.tb03348.x, 1989.

[28] W. C. Chew, W. H. Weedon, "A 3D perfectly matched medium from modified Maxwell's equations with stretched coordinates," Microwave and Optical Technology Letters, vol. 7 no. 13, pp. 599-604, DOI: 10.1002/mop.4650071304, 1994.

[29] W. C. Chew, W. H. Weedon, A. Sezginer, "A 3-D perfectly matched medium by coordinate stretching and its absorption of static fields," Applied Computational Electromagnetics Symposium Digest, vol. 1, pp. 482-489, 1995.

[30] F. L. Teixeira, W. C. Chew, "A general approach to extend Berenger's absorbing boundary condition to anisotropic and dispersive media," IEEE Transactions on Antennas and Propagation, vol. 46 no. 9, pp. 1386-1387, DOI: 10.1109/8.719984, 1998.

[31] F. L. Teixeira, W. C. Chew, "General closed-form PML constitutive tensors to match arbitrary bianisotropic and dispersive linear media," IEEE Microwave and Guided Wave Letters, vol. 8 no. 6, pp. 223-225, DOI: 10.1109/75.678571, 1998.

[32] F. L. Teixeira, W. C. Chew, "Unified analysis of perfectly matched layers using differential forms," Microwave and Optical Technology Letters, vol. 20 no. 2, pp. 124-126, DOI: 10.1002/(SICI)1098-2760(19990120)20:2<124::AID-MOP12>3.0.CO;2-N, 1999.

[33] M. Kuzuoglu, R. Mittra, "Frequency dependence of the constitutive parameters of causal perfectly matched anisotropic absorbers," IEEE Microwave and Guided Wave Letters, vol. 6 no. 12, pp. 447-449, DOI: 10.1109/75.544545, 1996.

[34] M. N. Nabighian, "Electromagnetic methods in applied geophysics," Society of Exploration Geophysicists, vol. 2, pp. 211-216, 1988.

[35] D. Guptasarma, "Computation of the time-domain response of a polarizable ground," Geophysics, vol. 47 no. 11, pp. 1574-1576, DOI: 10.1190/1.1441307, 1982.

[36] K. Yee, "Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media," IEEE Transactions on Antennas and Propagation, vol. 14 no. 3, pp. 302-307, DOI: 10.1109/tap.1966.1138693, 1966.

[37] Z. L. Niu, Principle of Electromagnetic Method in Time Domain, 2007.