Compact Matrix-Exponential-Based FDTD with Second-Order PML and Direct Z-Transform for Modeling Complex Subsurface Sensing and Imaging Problems

To simulate complex subsurface sensing and imaging problems with both propagating and evanescent waves by the finite-difference time-domain (FDTD) method, the highly-accurate second-order perfectly matched layer (SO-PML) formulations based on the direct Z-transform (DZT) and the matrix exponential (ME) techniques are compactly and efficiently proposed for modeling open-domain problems. During mathematical deductions, several manipulations, for example, convolution computations, formulation reorganizations, or variable substitutions, can be circumvented due to the fact that the ME-based method shows a compact first-order differential matrix form. Besides, any material attributes can be completely circumvented because of using electric and magnetic flux densities, consequently, the proposed DZT-SO-PML could be applied without needing any alteration. Moreover, the DZT-SO-PML method can not only preserve better absorption accuracies, but also attain palpable improvements in computational efficiencies, even if the distance between the DSP-SO-PML truncation and the target becomes closer for modeling 3D open-domain subsurface sensing and imaging problems. Finally, numerical examples have been carried out to illustrate and validate these proposed formulations.


Introduction
In modern communications, the radio wave in the air can transmit information for radio, radar, navigation systems and other applications [1][2][3][4]. However, limitations [5] have been found in shorter-wavelength radio waves whose signals become weaker and weaker when propagating farther and farther, even never transmit through water, and easily are prevented by the rock layer. For example, the global position system (GPS) [5] signals cannot penetrate into water, soil or building walls, or even can never propagate into these objects, So, the submarine or underground activities, such as mineral surveys, cannot adopt GPS signals. Besides, GPS and other radio signals work quite poorly in indoor rooms, or outdoor skyscrapers in between.
To overcome problems depicted above, nonetheless, the very low frequency (VLF) radio with the longer wavelength can be considered as an impactful approach [6]. As is well known to all, VLF waves have advantages of long-distance propagation and strong-power penetration [7], which can propagate and span over hundreds of feet in water and earth, and thousands of miles in air. Nowadays, the VLF signals have been widely applied in the military and civil fields, such as long-distance communication and navigation, wireless heart rate detectors, disaster relief, and submarine communication and geophysics research.
In many applications, the VLF geophysics exploration research is so few mentioned in the open-domain finite-difference time-domain (FDTD) numerical simulation, as investigations on the absorption performance of the perfectly matched layer (PML) for three-dimensional (3D) VLF subsurface sensing and imaging problems are quite rare in the literatures.
As far as we know, the traditional explicit finite-difference time-domain method has hitherto attracted great popularity as a quite efficacious tool for solving the Maxwell's equations, which are the effective and pinpoint solution of electromagnetic wave interaction problems [8]. With the rapid development of FDTD, it can be used for dealing with numerous types of applications, like antennas, waveguides, propagation, scattering, microwave circuits, non-linear and other special materials, and many other applications [9]. Nevertheless, to model open-domain problems, the split-field perfectly matched layer scheme [10], as one of the most popular and accurate absorbing boundary conditions (ABCs), is developed to truncate the physical domain due to the fact that computers can never store an unlimited amount of field data.
Afterward, several unsplit-field-based PML methods have been gradually adopted to replace previous PMLs in that much memory and computational time can be saved [11][12][13][14][15]. As with split-field PMLs above, nevertheless, strong evanescent waves can be ineffectively absorbed by unsplit-field ones. To avoid this case, therefore, the complex-frequency-shifted PML (CFS-PML) [16] is firstly proposed to attenuate low-frequency evanescent waves and reduce late-time reflections, leading to attracting substantial attention so much so that a large number of the literatures appears to apply this kind of technique [17][18][19][20][21].
To further conquer limitations of the above PMLs with/without CFS scheme, several higher-order PMLs [22][23][24][25][26][27] are presented, which can simultaneously absorb both strong propagating and evanescent waves. Moreover, among these higher-order PMLs, the secondorder PML (SO-PML) [28] has been validated and proved to be a better selection, which not only requires less CPU time and memory, but also maintains almost the same absorption accuracies as compared with higher-order PMLs (n > 2).
For engineering applications in electromagnetics, as far as we know, the maximum relative reflection error with below −40 dB can reach engineering requirements. Consequently, to reduce unnecessary consumption in memory and CPU time, the unsplit-field CFS-PML based on the direct Z-transform (DZT) is proposed to achieve this goal [29]. Very recently, a modified CFS-PML implementation [30] using the matrix exponential (ME) technique is presented, which has a higher accuracy than presented CFS-PMLs and maintains fairly good efficiency. For complicated subsurface sensing and imaging cases with smaller FDTD physical domains, nevertheless, one-pole PMLs may not achieve good absorption accuracies so that larger physical domains are required, leading to occupying much more memory and CPU time.
In this paper, to achieve satisfactory absorption accuracies for subsurface sensing and imaging problems and decrease unnecessary FDTD physical regions, an efficient and compact second-order PML with DZT and ME techniques is developed. The proposal is referred to as the DZT-ME-SO-PML. This paper is properly organized as follows: In Section 2, the proposed DZT-ME-SO-PML scheme is developed and introduced. 3D numerical examples, such as the airborne transient electromagnetic (ATEM), are shown in Section 3. In Section 4, for subsurface imaging problems based on measured data detected from the Chang'E-5 Lunar-Exploration Lander, the proposed formulations are used to implement the reverse-time migration (RTM) technique according to the Chang'E-5 Data and help the Chang'E-5 avoid hard rocks when exploring the lunar soil. The conclusions are drawn in Section 4.

ME-Based DZT-SO-PML Formulations
For the frequency domain, the modified Maxwell's curl equations can be developed for the three-dimensional second-order PML truncation, shown as where By applying constitutive relations above, we can obtain and update electric and magnetic fields using D (electric flux density) and B (magnetic flux density) relationship directly due to its full independence of material properties.
As far as we know, the operator ∇ s executed in Equations (1) and (2) can be defined as where Sη (η = x, y, z) is the n-pole complex stretched coordinate metric, shown below It should be noted that we here adopt the SC-PML with the two-pole CFS scheme, and hence the Sη (η = x, y, z) is complex stretched coordinate metrics, set as below where κi,η must be real and ≥1, and σi,η and αi,η are assumed to be positive real. Now, we re-range and tidy up Equation (7), and have where Then we have where ξ −1 where , υ p,η = − exp(υ p,η ∆t), all of them are coefficients.
So far, we have finished converting jω→s→Z. Afterward, we are ready to expand Equations (1) and (2) in the Cartesian coordinates, take y-component for a case Incorporate S η (Z) into Equation (10) and introduce auxiliary variable ψ and δ, we have Tide up, and we have Now, we further tidy up Equation (17), at first, we add ψ Dyz -ψ Dyx and δ Dyz -δ Dyx in the left-hand side of Equations (18)- (21), and then multiply 1/∆t for both side, we have where Now, all of fields' format can be shown as where the symbols ψ and δ are the auxiliary variables, and the symbols Φ and Θ are the substituted electromagnetic fields (E or H). Next, we can get all scalar equations in the following where If the initialize condition t 0 is given, we can obtain below By replacing t 0 with the old timestep n∆t, the above formula in the new timestep (n + 1)∆t can be expressed as (35) Simplify the Equation (25), we can obtain where In FDTD computation, the PML's pseudo code in nth timestep can be expressed in y-direction as follow: Remote Sens. 2021, 13, 94 6 of 12 Step 1: The prestore field on the boundary Step 2: Field iteration for D y D y Step 3: Update the auxiliary quantity We have finished mathematical derivations on the DZT-ME-SO-PML, so far.

3D Airborne Transient Electromagnetic Problems
To our knowledge, signals in shorter-wavelength radio waves will be weakened to be neglected when propagating farther and farther, and easily are prohibited by the rock layer. To overcome this problem, the extreme low-frequency electromagnetic wave is considered due to its advantages of long-distance propagation and strong penetration. To validate absorption accuracy of the DZT-ME-SO-PML, 3D numerical simulations on complicated subsurface sensing are carried out.
We consider a 3D dispersive numerical simulation for a half-space problem with three ores distributing at different locations, as shown in Figure 1. The whole computational domain contains air, rock (conductivity σ rock = 0.005 S/m), and mines (conductivity σ ore = 4 S/m), respectively.
Remote Sens. 2021, 13, x FOR PEER REVIEW 6 of 13 Step 1: The prestore field on the boundary x Step 2: Field iteration for Dy Step 3: Update the auxiliary quantity We have finished mathematical derivations on the DZT-ME-SO-PML, so far.

D Airborne Transient Electromagnetic Problems
To our knowledge, signals in shorter-wavelength radio waves will be weakened to be neglected when propagating farther and farther, and easily are prohibited by the rock layer. To overcome this problem, the extreme low-frequency electromagnetic wave is considered due to its advantages of long-distance propagation and strong penetration. To validate absorption accuracy of the DZT-ME-SO-PML, 3D numerical simulations on complicated subsurface sensing are carried out.
A bipolar square wave is recognized as the pulse excitation, and then is used with the repetition frequency of 25 Hz (namely, with duration of 0.04 s, rise time 0.0025 s, and fall time 0.0025 s).
It can be observed in Figure 2  interface of the PML and the FDTD computational domains and 1 at the end.
A bipolar square wave is recognized as the pulse excitation, and then is used with the repetition frequency of 25 Hz (namely, with duration of 0.04 s, rise time 0.0025 s, and fall time 0.0025 s).
It can be observed in Figure 2    It is concluded in Figure 2 that these PML methods can achieve good effects for an enough FDTD domain. Furthermore, it is clearly seen that the ME-based second-order PMLs can obtain obvious improvements of the MRREs as compared with others. Especially for the DZT-ME-SO-PML method, better absorption performance can be attained due to combining the DZT method than BZT-and MZT-ME-SO-PMLs.
To further indicate the advantage of DZT-ME-SO-PML for absorbing performances, we bring PML boundaries of x and y directions to very close to targets (gap distances between both are changed to 3 cells). As shown in Figure 3 Figure 3, the relative errors of all PMLs become larger apparently, when physical domains get smaller.
It can be concluded from the numerical result in Figure 3 that the DZT-PML cannot meet the engineering requirement (below −40 dB) due to using smaller domain; the ME-PML becomes much weaker than SO-PMLs based on Z-transform methods, in aspect of absorbing performance, due to existing strong propagating and evanescent waves; the absorption accuracies of Z-transform based SO-PMLs are improved at the price of more complex iteration forms and lower computational efficiency, in spite of obvious reductions for both early-and late-time reflections; To improve this problem, the ME technique is introduced into the BZT-and MZT-SO-PMLs so that higher efficiency is obtained; However, to our knowledge, the DZT approach can achieve the most accuracy as compared with the BZT and MZT, therefore, to sum up, the proposed ME-based DZT-SO-PML is developed to overcome these problems, resulting in achieving satisfactory accuracies and efficiencies for any general problems. It can be concluded from the numerical result in Figure 3 that the DZT-PML cannot meet the engineering requirement (below −40 dB) due to using smaller domain; the ME-PML becomes much weaker than SO-PMLs based on Z-transform methods, in aspect of absorbing performance, due to existing strong propagating and evanescent waves; the absorption accuracies of Z-transform based SO-PMLs are improved at the price of more complex iteration forms and lower computational efficiency, in spite of obvious reductions for both early-and late-time reflections; To improve this problem, the ME technique is introduced into the BZT-and MZT-SO-PMLs so that higher efficiency is obtained; Relative reflection e Next, to illustrate robustness and validity of the proposed method, the 3D subsurface imaging problem for measured data from Chang'E-5 lunar-exploration Lander is considered. To the best of our knowledge, the Chang'E-5 Project are going to be launched to explore a lunar subsurface region and take lunar geological samples back to the earth in the near future. Consequently, it is necessary to implement the high-resolution subsurface imaging firstly in the earth laboratory.

DZT-ME-SO-PML-Based FDTD Applied to the Reverse-Time Migration Method
The RTM method is the most reliable exploration methods in geological structure prediction. In 2020, the Chang'E-5 is planned for the first time to drill lunar soil in the moon, and then return it to earth, which is the most critical exploration mission in China's lunar exploration project. Due to the inevitable collision of those hard rocks during lunar-soil drilling, the purpose of the RTM method is to predict the distribution information of hard rocks in the soil in advance, so as to prevent the drill bit damage from affecting the soildrilling task of the Chang'E-5. The RTM method must be based on multiple calculations with FDTD as the underlying forward code, so we apply the DZT-ME-SO-PML scheme as its absorbing boundary condition to obtain the subsurface imaging.
In Chang'E-5, to transmit the broadband pulses and to obtain subsurface information, we employed 12 Vavaldi antennas as the basic element on its ground penetrating radar. As shown in Figure 4, the transmitted signals of 12 antennas from Table 1     In the past study [31][32][33][34][35], the property of the lunar soil was initially defined as the permittivity ε r = 3. To model the lunar environment, we adopt the ε r = 3 volcanic ash instead of the lunar soil as the background media. As shown in Figure 5a,b, two single-station laboratory data from an experimental lunar exploration system Chang'E-5 are extracted, and the subsurface laboratory model is filled with the volcanic ash in 2.5 m × 2.0 m × 1.0 m experiment region. In the process of measurement in the earth laboratory, a PEC barrier is placed at 2.5 m to block the unknown electromagnetic echoes that may be transmitted upward. According to the signal transmitted by one antenna and received by 11 antennas at the same time, the received data of 12 × 11 = 132 channels were effectively obtained by transmitting in turn. Those measurement results are shown in Figure 6. In the past study [31][32][33][34][35], the property of the lunar soil was initially defined as the permittivity εr = 3. To model the lunar environment, we adopt the εr = 3 volcanic ash instead of the lunar soil as the background media. As shown in Figure 5a,b, two singlestation laboratory data from an experimental lunar exploration system Chang'E-5 are extracted, and the subsurface laboratory model is filled with the volcanic ash in 2.5 m × 2.0 m × 1.0 m experiment region. In the process of measurement in the earth laboratory, a PEC barrier is placed at 2.5 m to block the unknown electromagnetic echoes that may be transmitted upward. According to the signal transmitted by one antenna and received by 11 antennas at the same time, the received data of 12 × 11 = 132 channels were effectively obtained by transmitting in turn. Those measurement results are shown in Figure 6.
(a) (b) Figure 5. Two single-station laboratory data acquisition systems from an experimental lunar exploration system Chang'E-5. The εr = 3 volcanic ash replaces lunar soil as background for buried Figure 5. Two single-station laboratory data acquisition systems from an experimental lunar exploration system Chang'E-5. The ε r = 3 volcanic ash replaces lunar soil as background for buried objects. In (a), the buried objects are respectively 1 Basalt, 2 Teflon, 3 5 Granite and 4 Metal. In (b), the buried objects are respectively 1 Basalt, 2 Teflon, and 3 4 Granite.
Remote Sens. 2021, 13, x FOR PEER REVIEW 10 of 13 objects. In (a), the buried objects are respectively ① Basalt, ② Teflon, ③⑤ Granite and ④ Metal. In (b), the buried objects are respectively ① Basalt, ② Teflon, and ③④ Granite. It can be seen from Figure 6 that the amplitude of the underground electromagnetic echo decreases with the increase of the receiving time of the antenna echo signal, and the underground electromagnetic signal cannot be understood even after the data is amplified. Hence the transmitting antenna can be considered as an electric dipole in a long dis- It can be seen from Figure 6 that the amplitude of the underground electromagnetic echo decreases with the increase of the receiving time of the antenna echo signal, and the underground electromagnetic signal cannot be understood even after the data is amplified. Hence the transmitting antenna can be considered as an electric dipole in a long distance area, so that the received data can be post processed. For a more intuitive understanding of subsurface information, the relation between echo dataset and propagation time is adjusted below where D meas and D opti are respectively the initial measured data and the optimal measured data. We can obtain the distribution of received data as shown in Figure 7. Obviously, the PEC barrier area buried at 2.5 m underground produces strong electromagnetic reflection and is received by antenna equipment, but the distribution of buried targets is still uncertain. It can be seen from Figure 6 that the amplitude of the underground electromagnetic echo decreases with the increase of the receiving time of the antenna echo signal, and the underground electromagnetic signal cannot be understood even after the data is amplified. Hence the transmitting antenna can be considered as an electric dipole in a long distance area, so that the received data can be post processed. For a more intuitive understanding of subsurface information, the relation between echo dataset and propagation time is adjusted below where Dmeas and Dopti are respectively the initial measured data and the optimal measured data. We can obtain the distribution of received data as shown in Figure 7. Obviously, the PEC barrier area buried at 2.5 m underground produces strong electromagnetic reflection and is received by antenna equipment, but the distribution of buried targets is still uncertain. In the real laboratory measurement, excitation signals can be transmitted from the Chang'E-5's UWBR system and recorded from the channel antenna. After the Chang'E-5 measurement is finished, the DZT-ME-SO-PML method can be applied to implement the subsurface imaging, shown in Figure 8a,b. The aim of reverse-time migration is to find the In the real laboratory measurement, excitation signals can be transmitted from the Chang'E-5's UWBR system and recorded from the channel antenna. After the Chang'E-5 measurement is finished, the DZT-ME-SO-PML method can be applied to implement the subsurface imaging, shown in Figure 8a,b. The aim of reverse-time migration is to find the buried objection in approximate location and help the Chang'E-5 to avoid the hard rocks. The specific shape cannot be restored unless the signal is without any noises. The two different responses for the PEC are occurred because the different electromagnetic echoes come from different paths through the different buried objects and different antenna locations. Consequently, it further corroborates the validity and accuracy of the proposal. buried objection in approximate location and help the Chang'E-5 to avoid the hard rocks. The specific shape cannot be restored unless the signal is without any noises. The two different responses for the PEC are occurred because the different electromagnetic echoes come from different paths through the different buried objects and different antenna locations. Consequently, it further corroborates the validity and accuracy of the proposal.

Conclusions
This study reveals a reliable fact that we can make full use of the DZT-ME-SO-PMLbased FDTD method for implementing the subsurface sensing and imaging problems. We introduce the ME technique for developing a compact first-order differential matrix form,

Conclusions
This study reveals a reliable fact that we can make full use of the DZT-ME-SO-PMLbased FDTD method for implementing the subsurface sensing and imaging problems. We introduce the ME technique for developing a compact first-order differential matrix form, incorporate the DZT method into the PML-based FDTD for achieving higher accuracy, and consider the second-order PML ABC for obtaining better absorption performance, for solving the 3D subsurface sensing and imaging problems.
The findings of this study can be summarized as follows: 1.
The proposed DZT-ME-SO-PML scheme could not only attenuate the strong propagating wave, but also absorb the evanescent wave and reduce late-time reflection, for 3D subsurface sensing.

2.
The DZT-ME-SO-PML-based FDTD can be applied to the RTM method so that the distribution information of hard rocks in the soil can be predicted in advance, in order to prevent the drill bit damage from affecting the soil-drilling task of the Chang'E-5.

3.
For the ATEM problems, we could achieve the secondary-field data from the receiver, analyze the distribution of scattering field, and finally predict the location of objects.