3D Numerical Modeling of Induced-Polarization Grounded Electrical-Source Airborne Transient Electromagnetic Response Based on the Fictitious Wave Field Methods

The grounded electrical-source airborne transient electromagnetic (GREATEM) system is widely used in mineral exploration. Meanwhile, the induced polarization (IP) effect, which indicates the polarizability of the earth, is often found. In this paper, the Maxwell equations in the frequency domain are transformed into fictitious wave domain, where Maxwell equations are solved by the time domain finite difference method. Then, an integral transformation method is used to convert the calculation results back to the time domain. A three-dimensional (3D) numerical simulation in a polarizable medium is presented. The accuracy of this method is proven by comparing it with the analytical solution and the existing method, and the calculation efficiency is increased five-fold. The simulation results show that the GREATEM system has a higher response amplitude in the conductive region, while IP effects cannot be identified in the conductive area. The GREATEM system has a higher response amplitude in the low-resistance region, but IP effects cannot be identified in the low-resistance area, and the detection of IP effects is more suitable for the high-resistance area. Therefore, it is necessary to improve the detection ability of the GREATEM system in the low-resistance area.


Introduction
Electromagnetic technology has a wide range of applications in various fields. Nasim (2019) proposed a metamaterial platform for solving integral equations using a monochromatic electromagnetic field [1]. The application of electromagnetic technology in medical diagnosis equipment realizes the measurement of blood glucose and the detection of cancer stage [2,3]. In addition, it is necessary to consider electromagnetic fields in plasma electronics [4,5].
Here, we apply electromagnetic technology to geophysics. Grounded electrical-source airborne transient electromagnetic (GREATEM) is a useful detection geophysical method where the transmitter is set on the ground and the receiver is installed on the aircraft [6]. Combining the ground and airborne electromagnetic systems, GREATEM has the advantages of large detection depth, high resolution, wide range, and fast speed. In particular, it has unique superiority in mountains, forest-covered areas, swamps, and other special areas [7,8]. Ji (2013) used an airship as the carrier and the two-dimensional finite difference method to calculate the time-domain GREATEM response with a long wire source [9]. Ren (2017) investigated a three-dimensional time-domain airborne electromagnetic model based on the finite volume method, and then used the strategy of separating the primary field from the secondary field, which significantly saved on calculation time [10]. For airborne transient electromagnetic measurement, a three-dimensional forward simulation method which considered attitude change was proposed, and the electromagnetic response of a shallow surface was analyzed in Reference [11]. Later, the spectral element method (SE) was used to carry out three-dimensional modeling of GREATEM, effectively improving the modeling accuracy [12].
To interpret GREATEM, researchers often make the assumption that the earth is just conductive or resistive. However, there is an induced polarization (IP)effect in the earth, which is usually reflected as a reverse signal. The existence of the IP effect has a great impact on the electromagnetic response [13,14]. The IP effect is related to conductivity and frequency, and its dispersion characteristics can be expressed by the Debye model and the Cole-Cole model [15,16]. At present, numerical simulations of the polarization effect are mainly conducted in the frequency domain. Most of the early studies simulated the IP effect by converting the electromagnetic field from the frequency domain to the time domain [17,18]. Unfortunately, when broadband frequency content is considered, the methods are not so efficient [19]. In order to simulate the IP effect in the time domain, Kang (2015) proposed a method to generate a three-dimensional pseudo charge distribution using airborne time-domain electromagnetic data and effectively extracted the polarization signal [20]. Wu (2017) used the FDTD (finite difference time-domain) method to simulate the IP effect with a small loop source [21]. Commer (2017) expanded the finite difference time-domain scheme, modeled and analyzed the IP effect, and developed a simple calculation of the Cole-Cole model [22].
However, three-dimensional (3D) numerical simulations of IP effects based on FDTD require a high computation load, when including the air layer. In addition, it is difficult to calculate the source in the FDTD method; thus, the upward continuation method is often used to avoid directly loading the emission source. Mittet (2010) put forward the theory of wave field transformation and successfully applied it to marine CSEM (controllable source electromagnetic method), increasing the calculation efficiency by nearly 10-fold [23]. Mittet (2018) used the least squares fitting method to calculate the IP response of sea water in the fictitious wave field [24].  improved the emission source, applied the wave field transformation method to the ground transient electromagnetic simulation, and realized a three-dimensional transient electromagnetic numerical simulation including the air layer [25]. However, it did not consider the IP effect of the earth. Therefore, here, we perform the wave field transformation in GREATEM, and we realize a three-dimensional numerical simulation in polarizable media using the Cole-Cole model. In addition, the computing efficiency is also improved.
The wave field transformation method is proposed for the grounded electrical-source airborne transient electromagnetic system in this paper. The Maxwell equations in the real diffusion field are transformed into the fictitious wave field by using the transformation relationship between the diffusion field and fictitious wave field. The fractional order Cole-Cole model is introduced into the fictitious wave field, and the Maxwell equations with the Cole-Cole model are solved using the finite difference time-domain method. Finally, an integral transformation method is used to transform the calculation results back to the real diffusion field. In the case of including polarized media, the two-field numerical simulation of the induction field and polarization field is realized. In addition, the computing efficiency is also improved.

Transformation from Real Diffusion Field to Fictitious Wave Field
The Maxwell equation in the frequency domain of the real diffusion field is written as follows: where E and H are the electric and magnetic fields in the real diffusion field, respectively, J denotes the current density, and K denotes the magnetic current density. ω is the angular frequency, µ is permeability tensor, and σ is conductivity, while x is the direction. The IP effect can be described by the Cole-Cole model [16].
where σ ∞ is the conductivity corresponding to infinite frequency, τ is the characteristic time constant, and c and η are the frequency dependence and chargeability, respectively. Applying Equation (3) to the frequency-domain Maxwell equation (Equation (1)), we can get A fictitious dielectric constant is defined by the conductivity tensor (Mittet, 2010).
where ω 0 = 2π × f 0 is the scaling parameter, f 0 = 1Hz, and ε ∞ is the fictitious dielectric permittivity tensor. Substituting Equation (5) into Equation (4), we can get Multiplying both sides of Equation (6) Equation (2) can be rewritten as The relationship between the real diffusion field and fictitious wave field was given by Reference [23].
where ω 0 = 2π is the scaling parameter, ω is the angular frequency in the fictitious wave field, J and K are the current and magnetic current densities in the fictitious wave field, and E and H are the electric and magnetic fields in the fictitious wave domain.
Because the appropriate kernel used to decompose the spectra is closer to a Cole-Cole function with an exponent c of 0.5 [26], then Assuming (16) and (17)) can be obtained through inverse Fourier transform.
By combining Equations (18) and (19), the electric and magnetic field x-components can be found as where Equation (22) can be solved by iteration, as shown in Equation (23).
where t n = n∆t, and the forward staggered time step is t N = (n + 1 2 )∆t.

Electrical Source Loading in Fictitious Wave Field
As shown in Figure 1, Tx is the long wire source and Rx is the receiver. The fictitious emission source is only related to x-and y-components of the electric field, and the z-component is zero.
Here, the first-order Gaussian pulse is selected as the fictitious emitter, as shown in Equation (24).
where β = π f 2 max , f max is the maximum frequency in the fictitious wave field. Here, the first-order Gaussian pulse is selected as the fictitious emitter, as shown in Equation (24).
where 2 max f π β = , max f is the maximum frequency in the fictitious wave field.

The Transformation from Fictitious Wave Field to Real Diffusion Field
It should be noted that the fictitious wave field is only used for the convenience of calculation, and it does not exist in reality. Therefore, we need to transform it back to the real diffusion field. The transformations are as follows [23]: It should be noted that i is the imaginary unit in Maxwell equations, and i, j, and k are the three direction vectors of the coordinate system.

The Transformation from Fictitious Wave Field to Real Diffusion Field
It should be noted that the fictitious wave field is only used for the convenience of calculation, and it does not exist in reality. Therefore, we need to transform it back to the real diffusion field. The transformations are as follows [23]: where T is the total calculation time in the fictitious wave field, and t is the sample time.

Loading Complex Frequency Shifted Perfectly Matched Layer (CFS-PML) in Fictitious Wave Field
Roden (2000) proposed CFS-PML (complex frequency shifted perfectly matched layer) boundary conditions, which significantly saved on memory [27]. The CFS-PML boundary conditions are introduced into the fictitious wave field to improve the absorption effect of the boundary on the low-frequency wave [28].
In x-components, Maxwell's curl equations are as follows: It should be noted that i is the imaginary unit in Maxwell equations, and i, j, and k are the three direction vectors of the coordinate system. s x , s y , and s z are the coordinate stretching factors, α (x,y,z) +iωε α (x,y,z) > 0, σ (x,y,z) > 0, and k (x,y,z) ≥ 1. By setting α (x,y,z) = σ (x,y,z) = 0, k (x,y,z) = 1, we can get the original stretched coordinate metrics, which is the PML boundary condition.
Taking E x as an example, and using Equations (29) and (30), we can get Appl. Sci. 2020, 10, 1027 6 of 15 Applying Laplace transform to (31) to the time domain, the iterative solution can be obtained.
On the boundary parameters, σ i and k i , which depend on the relative positions of nodes and target areas, are not fixed, whereas Ψ n+1/2 e xy i+1/2, j,k and Ψ n+1/2 e xz i+1/2,j,k are two discrete variables.
where k 0 indicates the intersection of the boundary and the target area boundary, d is the boundary thickness, and m is the polynomial parameter.

Half Space Model with IP Effect
In the field measurement with GREATEM, only the change of the underground secondary field is considered. In the process of numerical simulation, people usually divide the real space into two parts (air and ground), without considering the interface between the ground and air. Therefore, in the research process of this paper, the ground and the interface are not considered. In order to verify the effectiveness of this method, we design a half-space model incorporating the IP effect, as shown in Figure 2. Let the conductivity of the air layer σ air = 10 −6 S/m, conductivity at the infinite frequency of the half-space model σ∞ = 0.1 S/m, chargeability η= 0.6, characteristic time constant τ= 1 s, and frequency dependence c = 0.5 . Tx represents the transmitter, Rx is the receiver, the horizontal distance between Tx and Rx (offset distance) is 300 m, and the flight altitude is h = 0 m. In Figure 3, the calculation results are compared with Reference [29]. We can see that the solution of the fictitious domain finite difference method (FW-FDTD) is in good agreement with that of the transformed frequency domain method. The average relative error of both methods is less than 10%.          From Figures 4-6, it can be seen that characteristic time constant, chargeability, and conductivity at infinite frequency affect the occurrence time of negative response, and all responses are monotonic. Compared with the conductivity at infinite frequency, the effect of characteristic time constant and chargeability on the IP response is lower, mainly in the late stage, and the change trend of the early induction field is basically the same. With the decrease in characteristic time constant, the time of minus response is gradually advanced. With the decrease in chargeability, the time of minus response is delayed. The conductivity at infinite frequency has a great influence on the IP effect. In addition to the minus response at the later stage, it also has a great influence on the early induction field. In the Appl. Sci. 2020, 10, 1027 8 of 15 early stage, the response amplitude of the induction field increases gradually with the increase in conductivity at infinite frequency. In the later period, a smaller conductivity at infinite frequency results in the minus response appearing earlier.   From Figure 4 to Figure 6, it can be seen that characteristic time constant, chargeability, and conductivity at infinite frequency affect the occurrence time of negative response, and all responses are monotonic. Compared with the conductivity at infinite frequency, the effect of characteristic time constant and chargeability on the IP response is lower, mainly in the late stage, and the change trend of the early induction field is basically the same. With the decrease in characteristic time constant, the time of minus response is gradually advanced. With the decrease in chargeability, the time of minus response is delayed. The conductivity at infinite frequency has a great influence on the IP effect. In addition to the minus response at the later stage, it also has a great influence on the early induction field. In the early stage, the response amplitude of the induction field increases gradually    From Figure 4 to Figure 6, it can be seen that characteristic time constant, chargeability, and conductivity at infinite frequency affect the occurrence time of negative response, and all responses are monotonic. Compared with the conductivity at infinite frequency, the effect of characteristic time constant and chargeability on the IP response is lower, mainly in the late stage, and the change trend of the early induction field is basically the same. With the decrease in characteristic time constant, the time of minus response is gradually advanced. With the decrease in chargeability, the time of minus response is delayed. The conductivity at infinite frequency has a great influence on the IP effect. In addition to the minus response at the later stage, it also has a great influence on the early induction field. In the early stage, the response amplitude of the induction field increases gradually

Three-Dimensional Model
In addition, we calculated a three-dimensional model with a chargeable anomaly and compared the result with Reference [30]. The dimensions of the cell mesh were 70 × 70 × 70 , and the minimum grid size was 25 m. The modeling volume was 3500 × 3500 × 3500 m 3 . The size of the chargeable anomaly was 100 × 100 × 80 m 3 , while the distance from the top to the bottom of the chargeable anomaly was 40 m. The parameters of the Cole-Cole model were σ ∞ = 0.1 S/m, τ= 0.1 s, η= 0.3, and c = 0.5. The background conductivity was σ ∞ = 0.001 S/m. A vertical point dipole emitter, which was located 30 m above the surface, was used, and the receiver and transmitter were at the same location. The calculation results are shown in Figure 7, from which it can been seen that the two IP effect curves were consistent. The relative error was less than 10%, which further verifies the correctness of the method.
In terms of computation load, the iterations in the fictitious wave field are as follows: where N t represents the number of iteration steps, R max is the maximum transceiver distance, σ max and σ min are the maximum and minimum conductivity, respectively, and ∆x is the minimum grid size. Using Equation (37), the number of iteration steps was 76,681. It should be noted that only one transmitter rather than an array was selected in this work. Compared with Reference [30], the calculation efficiency was improved by about five-fold.
minimum grid size was 25 m . The modeling volume was . A vertical point dipole emitter, which was located 30 m above the surface, was used, and the receiver and transmitter were at the same location. The calculation results are shown in Figure 7, from which it can been seen that the two IP effect curves were consistent. The relative error was less than 10%, which further verifies the correctness of the method. In terms of computation load, the iterations in the fictitious wave field are as follows: where t N represents the number of iteration steps, max R is the maximum transceiver distance, max σ and min σ are the maximum and minimum conductivity, respectively, and is the minimum grid size. Using Equation (37), the number of iteration steps was 76,681. It should be noted that only one transmitter rather than an array was selected in this work. Compared with Reference [30], the calculation efficiency was improved by about five-fold.

Calculation Results of the Three-Dimensional Chargeable Body Model
The phenomenon of the inverse sign (IP effect) often appears in GREATEM. Most researchers focused on the influence of the IP parameters. However, in addition to the IP parameters, the system parameters of GREATEM also have an impact on the IP effect. In this section, we mainly analyze the effect of the system parameters, and we provide a theoretical basis for the field measurements.

Calculation Results of the Three-Dimensional Chargeable Body Model
The phenomenon of the inverse sign (IP effect) often appears in GREATEM. Most researchers focused on the influence of the IP parameters. However, in addition to the IP parameters, the system parameters of GREATEM also have an impact on the IP effect. In this section, we mainly analyze the effect of the system parameters, and we provide a theoretical basis for the field measurements.

Forward Modeling Results of 3D Chargeable Body
In order to analyze the influence of system parameters on the IP effect, we build a three-dimensional model with the consideration of a chargeable body, as shown in Figure 8. The ground conductivity was σ ∞ = 0.001 S/m, the air layer conductivity was σ air = 10 −6 S/m, Tx was the long wire source, the length was L Tx = 400 m, Rx represents the receiving point, the height from the ground was h = 25 m, and the offset distance was r = 200 m. The dimensions of the chargeable body were 400 × 400 × 300 m 3 , which was 100 m below the surface. The model parameters of the chargeable body were selected as follows: chargeability η= 0.5, frequency dependence c = 0.5, and characteristic time constant τ= 0.01 s. The conductivities at infinite frequency were σ ∞ = 0.01 S/m, σ∞ = 0.5 S/m, and σ∞ = 0.1 S/m. The calculation results are shown in Figure 9. The snapshots of the induced current system in the fictitious wave field are shown in Figure 10a

Forward Modeling Results of 3D Chargeable Body
In order to analyze the influence of system parameters on the IP effect, we build a three-dimensional model with the consideration of a chargeable body, as shown in Figure 8      It can be seen from Figure 9 that, with the increase in conductivity at the infinite frequency, the attenuation of the early induction field gradually decreased, the attenuation of the late IP response gradually quickened, and the occurrence time of the inverse sign gradually delayed. This is because the influence of the induction field of the low-resistance body was greater than that of the high-resistance body. Therefore, the second field decayed slowly, and the negative response appeared later. Thus, the polarization phenomenon was more easily observed on the high-resistance body. The chargeable body position can be clearly observed from Figure 10, and the method of wave field transformation could effectively identify the polarizable body. Since we set the air layer as a finite high-resistance body, the conductivity was not zero, and the propagation of induced current in the air layer was not zero. In the fictitious wave field, due to the influence of the chargeable body, the propagation mode of the electromagnetic field changed, and it no longer propagated in the form of symmetry.

After setting
, as well as the ground conductivity ground =0.001 S/m σ , ground =0.01 S/m σ , and ground =0.1 S/m σ , which are associated with high-, medium-, and low-resistance areas, the IP effects were calculated as shown in Figure 11. It can be seen from Figure 9 that, with the increase in conductivity at the infinite frequency, the attenuation of the early induction field gradually decreased, the attenuation of the late IP response gradually quickened, and the occurrence time of the inverse sign gradually delayed. This is because the influence of the induction field of the low-resistance body was greater than that of the high-resistance body. Therefore, the second field decayed slowly, and the negative response appeared later. Thus, the polarization phenomenon was more easily observed on the high-resistance body. The chargeable body position can be clearly observed from Figure 10, and the method of wave field transformation could effectively identify the polarizable body. Since we set the air layer as a finite high-resistance body, the conductivity was not zero, and the propagation of induced current in the air layer was not zero. In the fictitious wave field, due to the influence of the chargeable body, the propagation mode of the electromagnetic field changed, and it no longer propagated in the form of symmetry.
After setting σ∞ = 0.1 S/m, as well as the ground conductivity σ ground = 0.001 S/m, σ ground = 0.01 S/m, and σ ground = 0.1 S/m, which are associated with high-, medium-, and low-resistance areas, the IP effects were calculated as shown in Figure 11.
With the increase in ground conductivity, the amplitude of the early induction field and the late polarization field increased. In the high-and middle-resistance regions, the negative responses occurred at 1.7 ms and 3.3 ms. On the other hand, in the low-resistance region, the negative response could not be observed. When the negative value reached the minimum, the curve decayed approximately exponentially. The GREATEM had better resolution in the low-resistance region, but it is more suitable for the high-resistance region when measuring the IP effect.
propagation mode of the electromagnetic field changed, and it no longer propagated in the form of symmetry.
After setting , as well as the ground conductivity ground =0.001 S/m σ , ground =0.01 S/m σ , and ground =0.1 S/m σ , which are associated with high-, medium-, and low-resistance areas, the IP effects were calculated as shown in Figure 11. With the increase in ground conductivity, the amplitude of the early induction field and the late polarization field increased. In the high-and middle-resistance regions, the negative responses occurred at 1.7 ms and 3.3 ms. On the other hand, in the low-resistance region, the negative response could not be observed. When the negative value reached the minimum, the curve decayed approximately exponentially. The GREATEM had better resolution in the low-resistance region, but it is more suitable for the high-resistance region when measuring the IP effect.

The Influence of System Parameters on IP Response
In the field measurements of GREATEM, due to the influence of airflow and environment, the flight height was not a constant. Therefore, it was necessary to analyze the influence of flight altitude on the IP effect. We take the 3D model ( Figure 8) as an example. The air layer conductivity was set to σ air = 10 −6 S/m and the earth conductivity was set to 0.001 S/m, while the size of the chargeable body was 400 × 400 × 300 m. The model parameters of the chargeable body were as follows: chargeability  Figure 14.
Based on Figure 12, flight altitude mainly affected the amplitude of the response. The effect on the time of sign reversal was not obvious, and the sign reversal phenomenon was mainly concentrated in the range 1.5-2 ms. With the increase in flight altitude, the received electromagnetic response decreased gradually. It can be seen from Figure 13 that the length of the emission source mainly affected the amplitude of the response. The length of emission source increased gradually, and the amplitude of response increased gradually, which also shows that a longer emission source resulted in a higher resolution of the polarizable body. In Figure 14, with the decrease in offset distance, the time of the sign reversal appeared earlier, and the polarization response became easier to observe at a short offset distance.
Next, we set the air layer conductivity to σ = 10 −6 S/m and earth conductivity to 0.001 S/m. The size of the chargeable body was 400 × 400 × 100 m 3 . The model parameters of the chargeable body were chosen as η= 0.5, c = 0.5, τ= 0.01 s, σ ∞ = 0.1 S/m, offset distance r = 200 m, emission source length L Tx = 400 m, and flight altitude h = 25 m. The chargeable body was buried 25 m and 50 m below the surface. The response curves of the chargeable body at different depths are shown in Figure 15. Moreover, we changed the size of the chargeable body to 300 × 300 × 200 m 3 and the chargeable body depth to 100 m, while keeping other parameters the same, and setting the distance from the left boundary of the chargeable body to the emission source as 0 m, 50 m, and 100 m. The response curves at different source positions are shown in Figure 16.
GREATEM on the IP response, other parameters remained unchanged, while the flight height was set to h = 25 m，and the emission source length was Tx L =50 m , Tx L =100 m , and Tx L =200 m . The responses of different emission source lengths are shown in Figure 13. Taking Tx L =400 m and the offset distance r=100 m , r=200 m , and r=400 m , the responses of different offsets are shown in Figure 14.   GREATEM on the IP response, other parameters remained unchanged, while the flight height was set to h = 25 m，and the emission source length was Tx L =50 m , Tx L =100 m , and Tx L =200 m . The responses of different emission source lengths are shown in Figure 13. Taking Tx L =400 m and the offset distance r=100 m , r=200 m , and r=400 m , the responses of different offsets are shown in Figure 14.   Based on Figure 12, flight altitude mainly affected the amplitude of the response. The effect on the time of sign reversal was not obvious, and the sign reversal phenomenon was mainly concentrated in the range 1.5-2 ms. With the increase in flight altitude, the received electromagnetic response decreased gradually. It can be seen from Figure 13 that the length of the emission source mainly affected the amplitude of the response. The length of emission source increased gradually, and the amplitude of response increased gradually, which also shows that a longer emission source resulted in a higher resolution of the polarizable body. In Figure 14, with the decrease in offset distance, the time of the sign reversal appeared earlier, and the polarization response became easier to observe at a short offset distance.
Next, we set the air layer conductivity to   From Figures 15 and 16, it can be seen that all the amplitudes decreased, while the depth of chargeable body increased. However, the time of the reversal sign did not change. In the early stage, as the distance between the emitter source and the chargeable body decreased, the amplitude of the induced field increased, the attenuation of the curve was delayed. In the later stage, the time of sign reversal appeared later when the distance between the emitter and the chargeable body decreased. 50 m below the surface. The response curves of the chargeable body at different depths are shown in Figure 15. Moreover, we changed the size of the chargeable body to 3 300 300 200 m × × and the chargeable body depth to 100 m , while keeping other parameters the same, and setting the distance from the left boundary of the chargeable body to the emission source as 0 m , 5 0 m , and 100 m . The response curves at different source positions are shown in Figure 16.  From Figures 15 and 16, it can be seen that all the amplitudes decreased, while the depth of chargeable body increased. However, the time of the reversal sign did not change. In the early stage, as the distance between the emitter source and the chargeable body decreased, the amplitude of the induced field increased, the attenuation of the curve was delayed. In the later stage, the time of sign reversal appeared later when the distance between the emitter and the chargeable body decreased.

Conclusions
Three-dimensional numerical simulations are the basis of data interpretation and imaging. In this paper, a three-dimensional numerical simulation of underground polarized media was realized, which achieved simulation results closer to the actual mineral resources. In the actual measurement process with GREATEM, the geodetic model could be established using this method, and the optimal configuration structure of the equipment was obtained using the numerical simulation method. This not only provides the basis for the exploration of underground mineral resources, but also provides theoretical guidance for the design of detection instruments. Wave field transformation was proposed to calculate the response of the GREATEM system using the Cole-Cole model. The electromagnetic wave was absorbed by the boundary condition of CFS-PML in the fictitious wave field. Compared with the analytical solution and the transformed FD method, the relative error was less than 10%, satisfying the standard calculation requirements. Compared with the existing method, the calculation time was reduced five-fold. According to the simulations, the flight altitude, electrical source length, electrical source position, and the offset distance all affected the electromagnetic response monotonically. In the homogeneous half-space model, when ) were too large, the negative response was likely overwhelmed by the noise; thus, the negative response could not be observed. Therefore, the occurrence time of the negative response could be enhanced by reducing the offset distance. Low-altitude flights can facilitate the recognition ability of the chargeable body, but they are often limited by geological conditions. The response of the GREATEM

Conclusions
Three-dimensional numerical simulations are the basis of data interpretation and imaging. In this paper, a three-dimensional numerical simulation of underground polarized media was realized, which achieved simulation results closer to the actual mineral resources. In the actual measurement process with GREATEM, the geodetic model could be established using this method, and the optimal configuration structure of the equipment was obtained using the numerical simulation method. This not only provides the basis for the exploration of underground mineral resources, but also provides theoretical guidance for the design of detection instruments. Wave field transformation was proposed to calculate the response of the GREATEM system using the Cole-Cole model. The electromagnetic wave was absorbed by the boundary condition of CFS-PML in the fictitious wave field. Compared with the analytical solution and the transformed FD method, the relative error was less than 10%, satisfying the standard calculation requirements. Compared with the existing method, the calculation time was reduced five-fold. According to the simulations, the flight altitude, electrical source length, electrical source position, and the offset distance all affected the electromagnetic response monotonically. In the homogeneous half-space model, when the characteristic time constant (τ > 1 s) and conductivity at the infinite frequency (σ ∞ > 1 S/m) were too large, the negative response was likely overwhelmed by the noise; thus, the negative response could not be observed. Therefore, the occurrence time of the negative response could be enhanced by reducing the offset distance. Low-altitude flights can facilitate the recognition ability of the chargeable body, but they are often limited by geological conditions. The response of the GREATEM system had good recognition ability in the low-resistance area; however, the chargeable body could hardly be detected. Therefore, it is necessary to improve the detection accuracy of the GREATEM in low-resistance areas.
This paper realized a three-dimensional numerical simulation of polarized media based on the Cole-Cole model, and improved the calculation efficiency. However, in the actual measurement process, the underground medium is diverse. It is not comprehensive to only use the Cole-Cole model to describe the characteristics of polarized medium; the existence of the GEMTIP model in particular is a major problem in the field of geophysics. In future work, we need to consider the GEMTIP model in the numerical simulation, so as to improve the universality of this method.