Application of Laplace Domain Waveform Inversion to Cross-Hole Radar Data

Abstract

Full waveform inversion (FWI) can yield high resolution images and has been applied in Ground Penetrating Radar (GPR) for around 20 years. However, appropriate selection of the initial models is important in FWI because such an inversion is highly nonlinear. The conventional way to obtain the initial models for GPR FWI is ray-based tomogram inversion which suffers from several inherent shortcomings. In this paper, we develop a Laplace domain waveform inversion to obtain initial models for the time domain FWI. The gradient expression of the Laplace domain waveform inversion is deduced via the derivation of a logarithmic object function. Permittivity and conductivity are updated by using the conjugate gradient method. Using synthetic examples, we found that the value of the damping constant in the inversion cannot be too large or too small compared to the dominant frequency of the radar data. The synthetic examples demonstrate that the Laplace domain waveform inversion provide slightly better initial models for the time domain FWI than the ray-based inversion. Finally, we successfully applied the algorithm to one field data set, and the inverted results of the Laplace-based FWI show more details than that of the ray-based FWI.

1. Introduction

Ground Penetrating Radar (GPR) is extensively applied in archaeological, environmental, hydrological and engineering investigations to detect subsurface structures [1]. Cross-hole radar (normally 20–250 MHz), with transmitting antennas in a borehole and receiving antennas in a second, neighboring borehole, can provide high-resolution tomographic images of permittivity ($\epsilon$) and conductivity ($\sigma$) between the two boreholes.

Traditional radar tomographic inversions, such as travel time inversion and attenuation inversion, are based on geometrical ray theory [2,3,4,5,6,7]. The resolution of ray theory tomography is around the width of the first Fresnel-zone [8,9]. Full-waveform inversion (FWI) was first applied in seismic exploration (e.g., [10,11,12]) and was subsequently introduced to GPR (e.g., [13,14,15,16,17,18]). Under the most ideal conditions, the resolution of FWI can reach up to the scale of half the propagated wavelength [19].

Although FWI has been developed for over a few decades, there are still several challenges to confront to make it more popular [19]. One challenge is building accurate initial models. Conventional initial models in the cross-hole FWI are from ray-based inversions [2,3,4,5,6,7]. However, these ray-based inversions suffer from several shortcomings in inversion processing, including the high-frequency limitation of the ray approximation, the limited angular coverage of the target and the small amount of signal attributes [16].

In exploration seismic, Shin and Cha [20] suggest using a Laplace domain waveform inversion to build an initial velocity model for FWI. By back-propagating the long-wavelength residuals in the Laplace domain, the results of the inversion can provide a smooth reconstruction of the velocity model as an initial model for the subsequent time or frequency domain FWI [21]. Subsequently, Shin and Ha [22] compare the behavior of objective functions in the frequency and Laplace domains and conclude that the objective functions in the Laplace domain are smoother than that in the frequency domain. Waveform inversion in the Laplace domain is insensitive to the frequency content and yields consistent results when lacking low frequency data [23,24]. Park et al. [25] suggest using a power objective function in the Laplace domain waveform inversion for data with a small amplitude difference. In addition, 3D inversions in the Laplace domain have been applied to synthetic and field data to obtain 3D velocity models [26,27,28,29,30].

In this study, we propose an adaptive Laplace domain waveform inversion of cross-hole radar data, which provides an alternative approach for obtaining initial models for the time domain FWI. The gradients are deduced based on the derivation of the objective function with respect to permittivity and conductivity. Permittivity and conductivity distributions are updated by adopting the conjugate gradient method and a line search approach. For the forward modeling in the inversion processing, we solve Maxwell’s equation for the transverse electric (TE) mode by using the 2D finite-difference time-domain method (FDTD) with the Convolutional Perfected Marched Layer (CPML) absorbing boundary condition [31].

The proposed inversion is tested by two synthetic examples. First, Laplace domain waveform inversion at a set of damping constants is performed to choose the optimal value of the damping constant. Second, we compare the results of the time domain FWI that use the Laplace domain waveform and ray-based results as the initial models. Finally, we apply the Laplace domain waveform inversion to one field cross-hole radar data. All images of the Laplace domain waveform inversion are compared with the respective ray-based results.

2. Methods

2.1. Laplace-Transformed Wavefield

The time domain wavefield is transformed into the Laplace domain by:

$$\tilde{E}\left(s\right)={\int }_0^{\infty }E\left(t\right)e^{-st}dt$$

(1)

where s is a real Laplace damping constant, t is time, $E\left(t\right)$ is the time domain electric fields and $\tilde{E}\left(s\right)$ is the Laplace domain electric fields.

In Equation (1), ${\int }_0^{\infty }E\left(t\right)e^{-st}$ can be viewed as a damped electric field at a given s. The electric field in the Laplace domain is the zero frequency component of the Fourier transform of the damped wavefield at a complex-valued frequency, in which the real part is zero and the imaginary part is related to the damping constant s [20].

The characteristics of the Laplace domain wavefield of a radargram at a damping constant are similar to that of a seismogram [20]. At a low damping constant, the damped wavefield contains most parts of the raw wavefield following the first arrival event. When the damping constant moves to a high value, the proportion of the early part of the wavefield gradually ascends. The amplitude values of the damped wavefields decrease rapidly as the damping constant increases. The higher the damping constant, the lower the value of the Laplace domain wavefield.

The wavefields from cross-hole radargrams are dominant by the first direct wave, which is different from the wavefields in the reflection seismograms [20]. Laplace domain waveform inversion of the cross-hole radar data also provides long-wavelength results because of the smooth features of the virtual source in the processing of gradients [21]. In the next section, we apply the Laplace domain waveform inversion to cross-hole radar data and test its performance.

2.2. Forward Problem

In the time domain, Maxwell’s equations are written as follows:

$$\left[\begin{array}{cc}-\epsilon \left(\mathrm{x}\right){\partial }_t-\sigma \left(\mathrm{x}\right)& \nabla \times \\ \nabla \times & {\mu }_0{\partial }_t\end{array}\right]\left[\begin{array}{c}E\left(\mathrm{x},t\right)\\ H\left(\mathrm{x},t\right)\end{array}\right]=\left[\begin{array}{c}J\left(\mathrm{x},t\right)\\ 0\end{array}\right]$$

(2)

where E and H are the electric and magnetic fields generated by the source J. $\epsilon \left(\mathrm{x}\right)$ and $\sigma \left(\mathrm{x}\right)$ are the permittivity and conductivity at space x. ${\mu }_0$ is the permeability at the free space. According to Equation (1), Equation (2) in the Laplace domain is

$$\left[\begin{array}{cc}-s\cdot \epsilon \left(\mathrm{x}\right)-\sigma \left(\mathrm{x}\right)& \nabla \times \\ \nabla \times & s\cdot {\mu }_0\end{array}\right]\left[\begin{array}{c}\tilde{E}\left(\mathrm{x},s\right)\\ \tilde{H}\left(\mathrm{x},s\right)\end{array}\right]=\left[\begin{array}{c}\tilde{J}\left(\mathrm{x},s\right)\\ 0\end{array}\right]$$

(3)

where $\tilde{E}$, $\tilde{H}$ and $\tilde{J}$ are the electric fields, magnetic fields and source in the Laplace domain, respectively. If eliminating $\tilde{H}$ from Equation (3), we can reformulate Equation (3) by a simple representation:

$$\tilde{\mathrm{M}}\tilde{E}=\tilde{J}$$

(4)

where $\tilde{\mathrm{M}}$ is the linear operator of Maxwell’s equations in the Laplace domain.

$$\tilde{\mathrm{M}}=\left[\begin{array}{cc}-s\cdot \epsilon \left(\mathrm{x}\right)-\sigma \left(\mathrm{x}\right)& \nabla \times \\ \nabla \times & s\cdot {\mu }_0\end{array}\right]$$

(5)

It is important to note that the elimination of $\tilde{H}$ from Equation (4) is only for simplifying the formulation. In the Laplace domain, the electric fields at an arbitrary position can be considered as the multiplication between the source $\tilde{J}$ and $\tilde{\mathrm{G}}$ as in the frequency domain:

$$\tilde{E}=\tilde{\mathrm{G}}\tilde{J},\mathrm{therefore}\tilde{\mathrm{G}}={\tilde{\mathrm{M}}}^{-1}$$

(6)

where $\hat{\mathrm{G}}$ is the Green’s operator of $\tilde{\mathrm{M}}$.

2.3. Inversion Problem

In our inversion, we used a logarithmic wavefield in the object function because of the small values of the electric fields in the Laplace domain [20,32]:

$$O\left(\epsilon ,\sigma \right)=\frac{1}{2}{\displaystyle \sum _{ns}{\displaystyle \sum _{nr}{\left[\ln \tilde{E}\left(\epsilon ,\sigma \right)-\ln {\tilde{E}}^{obs}\right]}^2}}$$

(7)

where $\tilde{E}\left(\epsilon ,\sigma \right)$ and ${\tilde{E}}^{obs}$ are the predicted and observed data in the Laplace domain; and ns and nr indicate the number of sources and receivers, respectively. Later, we use $\tilde{E}$ as the simplification of $\tilde{E}\left(\epsilon ,\sigma \right)$.

The gradients are deduced by the derivation of Equation (7) with respect to permittivity ($\epsilon$) and conductivity ($\sigma$).

$$\frac{\partial O}{\partial p}={\displaystyle \sum _{ns}{\displaystyle \sum _{nr}\frac{\partial \tilde{E}}{\partial p}\left(\frac{\ln \tilde{E}-\ln {\tilde{E}}^{obs}}{\tilde{E}}\right)}}$$

(8)

where p represents the model parameters $\epsilon$ and $\sigma$, and $\frac{\partial \tilde{E}}{\partial p}$ is obtained by taking the partial derivative of both sides of Equation (4) with respect to the model parameter p:

$$\frac{\partial \tilde{\mathrm{M}}}{\partial p}\tilde{E}+\tilde{\mathrm{M}}\frac{\partial \tilde{E}}{\partial p}=0,\mathrm{and}\frac{\partial \tilde{E}}{\partial p}={\tilde{\mathrm{M}}}^{-1}v=\tilde{\mathrm{G}}v$$

(9)

where v is the virtual source vector [20,32].

$$v=-\frac{\partial \tilde{\mathrm{M}}}{\partial p}\tilde{E}$$

(10)

Within M in Equation (5), only term $-s\cdot \epsilon \left(\mathrm{x}\right)-\sigma \left(\mathrm{x}\right)$ is related with the model parameters $\epsilon$ and $\sigma$. Therefore,

$$\{\begin{array}{c}{v}_{\epsilon }=-\frac{\partial \mathrm{M}}{\partial \epsilon }\tilde{E}=s\cdot \tilde{E}\\ v_{\sigma }=-\frac{\partial \mathrm{M}}{\partial \sigma }\tilde{E}=\tilde{E}\end{array}$$

(11)

By combining Equations (8), (9) and (11), the gradients of permittivity and conductivity can be written as

$$\left[\begin{array}{c}\nabla O_{\epsilon }\\ \nabla O_{\sigma }\end{array}\right]={\displaystyle \sum _{ns}\begin{array}{c}\left(s\cdot \tilde{E}\right)\tilde{\mathrm{G}}r\\ \left(\tilde{E}\right)\tilde{\mathrm{G}}r\end{array}},\mathrm{and}r={\displaystyle \sum _{nr}\frac{\ln \tilde{E}-\ln {\tilde{E}}^{obs}}{\tilde{E}}}$$

(12)

where $\tilde{\mathrm{G}}r$ is interpreted as the backward propagation of the residual r at the receiver locations.

To minimize the objective function in Equation (7), the model parameters can be updated as follows:

$$\left[\begin{array}{c}\epsilon {\left(\mathrm{x}\right)}_{k+1}\\ \sigma {\left(\mathrm{x}\right)}_{k+1}\end{array}\right]=\left[\begin{array}{c}\epsilon {\left(\mathrm{x}\right)}_k\\ \sigma {\left(\mathrm{x}\right)}_k\end{array}\right]-\left[\begin{array}{c}{\zeta }_{\epsilon ,k}\\ {\zeta }_{\sigma ,k}\end{array}\right]\cdot \left[\begin{array}{c}{C}_{\epsilon }{\left(\mathrm{x}\right)}_k\\ C_{\sigma }{\left(\mathrm{x}\right)}_k\end{array}\right]$$

(13)

where ${\zeta }_{\epsilon ,k}$ and ${\zeta }_{\sigma ,k}$ are the iterative step lengths, and $C_{\epsilon }{\left(\mathrm{x}\right)}_k$ and $C_{\sigma }{\left(\mathrm{x}\right)}_k$ are the conjugate gradients for the kth iteration.

$$C_{\epsilon }{\left(\mathrm{x}\right)}_k=\nabla O_{\epsilon }{\left(\mathrm{x}\right)}_k+\frac{\nabla O_{\epsilon }{\left(\mathrm{x}\right)}_k\left[\nabla O_{\epsilon }{\left(\mathrm{x}\right)}_k-\nabla O_{\epsilon }{\left(\mathrm{x}\right)}_{k-1}\right]}{\nabla O_{\epsilon }{\left(\mathrm{x}\right)}_{k-1}\nabla O_{\epsilon }{\left(\mathrm{x}\right)}_{k-1}}{C}_{\epsilon }{\left(\mathrm{x}\right)}_{k-1}$$

(14)

$$C_{\sigma }{\left(\mathrm{x}\right)}_k=\nabla O_{\sigma }{\left(\mathrm{x}\right)}_k+\frac{\nabla O_{\sigma }{\left(\mathrm{x}\right)}_k\left[\nabla O_{\sigma }{\left(\mathrm{x}\right)}_k-\nabla O_{\sigma }{\left(\mathrm{x}\right)}_{k-1}\right]}{\nabla O_{\sigma }{\left(\mathrm{x}\right)}_{k-1}\nabla O_{\sigma }{\left(\mathrm{x}\right)}_{k-1}}{C}_{\sigma }{\left(\mathrm{x}\right)}_{k-1}$$

(15)

when k = 1, $C_{\epsilon }{\left(\mathrm{x}\right)}_1=\nabla O_{\epsilon }{\left(\mathrm{x}\right)}_1$ and $C_{\sigma }{\left(\mathrm{x}\right)}_1=\nabla O_{\sigma }{\left(\mathrm{x}\right)}_1$.

The updating of the step lengths follows the approach proposed by Pica et al. [33]:

$${\zeta }_{\epsilon ,k}={\kappa }_{\epsilon }\frac{{\displaystyle \sum _{ns}{\displaystyle \sum _{nr}\left[\frac{{\tilde{E}}_{i,j}\left(\epsilon +{\kappa }_{\epsilon }{C}_{\epsilon ,k},\sigma \right)-{\tilde{E}}_{i,j}\left(\epsilon ,\sigma \right)}{{\tilde{E}}_{i,j}\left(\epsilon ,\sigma \right)}\right]\left[\ln {\tilde{E}}_{i,j}\left(\epsilon ,\sigma \right)-\ln {\tilde{E}}_{i,j}^{obs}\right]}}}{{\displaystyle \sum _{ns}{\displaystyle \sum _{nr}\left[\frac{{\tilde{E}}_{i,j}\left(\epsilon +{\kappa }_{\epsilon }{C}_{\epsilon ,k},\sigma \right)-{\tilde{E}}_{i,j}\left(\epsilon ,\sigma \right)}{{\tilde{E}}_{i,j}\left(\epsilon ,\sigma \right)}\right]\left[\frac{{\tilde{E}}_{i,j}\left(\epsilon +{\kappa }_{\epsilon }{C}_{\epsilon ,k},\sigma \right)-{\tilde{E}}_{i,j}\left(\epsilon ,\sigma \right)}{{\tilde{E}}_{i,j}\left(\epsilon ,\sigma \right)}\right]}}}$$

(16)

$${\zeta }_{\sigma ,k}={\kappa }_{\sigma }\frac{{\displaystyle \sum _{ns}{\displaystyle \sum _{nr}\left[\frac{{\tilde{E}}_{i,j}\left(\epsilon ,\sigma +{\kappa }_{\sigma }{C}_{\sigma ,k}\right)-{\tilde{E}}_{i,j}\left(\epsilon ,\sigma \right)}{{\tilde{E}}_{i,j}\left(\epsilon ,\sigma \right)}\right]\left[\ln {\tilde{E}}_{i,j}\left(\epsilon ,\sigma \right)-\ln {\tilde{E}}_{i,j}^{obs}\right]}}}{{\displaystyle \sum _{ns}{\displaystyle \sum _{nr}\left[\frac{{\tilde{E}}_{i,j}\left(\epsilon ,\sigma +{\kappa }_{\sigma }{C}_{\sigma ,k}\right)-{\tilde{E}}_{i,j}\left(\epsilon ,\sigma \right)}{{\tilde{E}}_{i,j}\left(\epsilon ,\sigma \right)}\right]\left[\frac{{\tilde{E}}_{i,j}\left(\epsilon ,\sigma +{\kappa }_{\sigma }{C}_{\sigma ,k}\right)-{\tilde{E}}_{i,j}\left(\epsilon ,\sigma \right)}{{\tilde{E}}_{i,j}\left(\epsilon ,\sigma \right)}\right]}}}$$

(17)

where ${\kappa }_{\epsilon }$ and ${\kappa }_{\sigma }$ are small stability factors. In this paper, we applied an automatic strategy to choose the values of the stability factors, which are equal to 1% of the model parameter (relative permittivity or conductivity) divided by the middle value of its gradient.

In our inversion, the forward modeling is performed in the time domain, whereas the residuals between the observed and predicted data are computed in the Laplace domain. However, in the back-propagation processing, the real-valued Laplace domain residuals cannot be transformed to the time domain because of the ill-posed problem [34]. We solve the problem by applying the method of Ha et al. [29], in which a predefined time domain source wavelet with the amplitude of the residuals is used.

3. Results

We first tested the performance of the Laplace domain waveform inversion at different damping constant (s) values through a simple synthetic example. Then, a complex synthetic example was used to compare the inverted results of the time domain FWI, whose initial models are from the ray-based inversion and Laplace domain waveform inversion, respectively. Finally, the Laplace domain waveform inversion was successfully applied to field data and provides initials models for the time domain FWI.

3.1. Synthetic Data 1: A Simple Model

Figure 1a,c show the permittivity and conductivity models for testing the performance of our Laplace domain waveform inversion. The model size was 6 m in both the horizontal and vertical directions. Two cylindrical anomalies with diameters of 0.5 m were placed in the homogeneous medium of relative permittivity ${\epsilon }_r=5.5$ and conductivity $\sigma =0.005$ S/m. The upper left anomaly was located at the position of (2 m, 2 m), and had a relative permittivity of 7 and conductivity of 0.008 S/m; the relative permittivity and conductivity of the lower right anomaly at the position of (4 m, 4 m) were 4 and 0.003 S/m, respectively. The two 6 m deep boreholes were 6 m apart. Thirteen transmitter positions (circle indicated) were spaced at 0.5 m intervals in the left borehole and thirteen receiver positions (crosses indicated) were spaced at 0.5 m in the right borehole (Figure 1).

The initial models in the ray-based inversion and Laplace domain waveform inversion have the same relative permittivity and conductivity values as the background medium of the true models. To weaken the near-field effect in the forward and backward modeling, a median filter [35] was applied near the source regions [18]. We used a stepped parameter updating approach because of the large gradient differences between the permittivity and conductivity. The permittivity was inverted with the fixed conductivity and then the conductivity was inverted with the fixed permittivity [14]. In the inversion process, the relative permittivity and conductivity are updated in the logarithmic domain to ensure their positive values and improve convergence [36].

Figure 1b,d are the inverted results of the ray-based inversion in which the two cylindrical anomalies cannot be distinguished because the diameter of the anomalies is smaller than the resolution of the ray-based method [6,7]. Figure 2 shows the results of Laplace domain waveform inversion at six different damping constants after 50 iterations. The two anomalies are well reconstructed when the Laplace domain waveform inversion is performed. However, the inverted results are smoother than the true models due to the long-wavelength features in the Laplace domain. The sizes of the reconstructed anomalies in Figure 2 are larger than those of results from the time domain FWI [18].

In Figure 2a,b,g,h, the results of Laplace domain waveform inversion at the low damping constants of s = 50 × 10⁶ and s = 100 × 10⁶ fail to reconstruct the shape of the cylindrical anomalies. This is due to the inaccuracy of the numerical Laplace transform of the time domain wavefield at a low damping constant [20]. When the damping constant is 300 × 10⁶, both the results of relativity permittivity and conductivity reconstruct the two anomalies (Figure 2c,i). In the results of the Laplace domain waveform inversion at s = 500 × 10⁶, 700 × 10⁶ and 900 × 10⁶ (Figure 2d–f,j–l), the two cylindrical anomalies are well reconstructed. The oscillations in Figure 2 are caused by the high wave-number updates in the models, and the high wave-number components dominate the wave-number information as the s values are increasing.

Figure 3 shows cross sections of the relative permittivity and conductivity through the results obtained by the ray-based inversion, Laplace domain waveform inversion, and the true models. In the results of the ray-based inversion and Laplace domain waveform inversion at damping constants of s = 50 × 10⁶ and 100 × 10⁶, the positions and sizes of the anomalies are not well reconstructed. When the damping constants are 300 × 10⁶, 500 × 10⁶, 700 × 10⁶ and 900 × 10⁶, the permittivity tomograms from the Laplace domain waveform inversion successfully reconstruct the positions of anomalies at a smooth scale (Figure 3a). Similar to the relative permittivity, the conductivity tomograms of the Laplace domain waveform inversion at high damping constants (s = 300 × 10⁶, 500 × 10⁶, 700 × 10⁶ and 900 × 10⁶) depict the positions of the anomalies well (Figure 3b). As the damping constant increases, the conductivity values of the anomalies are more accurate. This may be because the Laplace domain wavefield at high damping constants has higher contrast gradients than that at low damping constants [20]. However, the higher the Laplace damping constant, the less the damped wavefield contains the time domain wavefield. In addition, the damped wavefield at high damping constants is very easily disturbed by noise [20]. Therefore, the damping constant in the Laplace waveform inversion of cross-hole radar data should be chosen in a moderate range. We propose that approximately five times the value of the dominant frequency of the radar data is appropriate.

Laplace domain wavefields at different damping constants are also the zero frequency component of the damped wavefield, but at different values. The inverted results of relative permittivity are approximately the same. Since the first direct event always dominates the damped wavefield of the cross-hole radar data, multiple damping constants in the inversion would not help if an optimal damping constant is chosen. In this study, we used a single damping constant in our Laplace domain waveform inversion.

3.2. Synthetic Data 2: A Complex Model

To further check the performance of the Laplace domain waveform inversion, we built a three-layered model containing three cylindrical anomalies (Figure 4a,d). The horizontal and vertical sizes of the model were 6 m. The three anomalies with diameters of 0.5 m were placed at the positions of (1 m, 3 m), (3 m, 3 m) and (5 m, 3 m). The top layer and the bottom layer had the same relative permittivity of 5 and conductivity of 0.001 S/m. The three cylindrical anomalies of ${\epsilon }_r=7$ and $\sigma =0.008$ S/m were buried in the middle layer of ${\epsilon }_r=5.5$ and $\sigma =0.0028$ S/m. The initial models in the Laplace domain waveform inversion were the same as the middle layer. The locations of transmitters and receivers were the same as the simple example. In this example, the damping constant was 500 × 10⁶ for the Laplace domain waveform inversion.

Figure 4 shows the relative permittivity and conductivity tomograms that result from the ray-based inversion and Laplace domain waveform inversion. The reconstructed anomalies in both relative permittivity and conductivity tomograms from the ray-based results (Figure 4b,e) are mixed together, whereas the layered medium is delineated roughly. In the results of the Laplace domain waveform inversion, the layer interfaces are well reconstructed (Figure 4c,f). However, the three anomalies in Figure 4c,f cannot be distinguished in the horizontal direction because of the small interval between the anomalies. To further compare the results of the ray-based inversion and Laplace domain waveform inversion, we apply their results as the initial models in the time domain FWI [18], and then compare their performance.

Figure 5a,c show the results of the time domain FWI based on the ray-based results, and Figure 5b,d are the results of the time domain FWI based on the Laplace domain waveform results. For both relative permittivity and conductivity, the layered medium and three cylindrical anomalies are rebuilt well. The interface between the middle layer and the bottom layer in Figure 5 agrees with that in the true model.

Figure 6a displays the curves of object function values during the 50 iterations of the time domain FWI in Figure 5. We note that the time domain FWI based on the Laplace domain waveform inversion has slightly large errors in the first five iterations, after which the error curve (red line in Figure 6a) decreases quickly and finally reach convergence at a lower value than the error curve (blue line in Figure 6a) of the time domain FWI based on the ray-based inversion. Figure 6b,c show cross-sections of relative permittivity and conductivity through the tomograms in Figure 4 and Figure 5. In the results of the time domain FWI, the positions and values of the cylindrical anomalies are clearly depicted. The relative permittivity values from the Laplace-based FWI are more accurate than those from the ray-based FWI (Figure 6b). In Figure 6c, FWI based on the ray-based inversion provides better values for the left and right anomalies. For the middle anomaly in Figure 6c, we note that the FWI based on the Laplace domain waveform inversion performs better. Similar to the results shown in Figure 3, inverted results of conductivity perform better than those of relative permittivity, as shown in previous studies [15,16,18]. In this example, Laplace domain waveform inversion is slightly better than the ray-based inversion to provide initial models for the time domain FWI.

3.3. Field Data Measured at Guizhou

Field data were collected in the Guizhou Province using radar systems (MALA Geoscience) with 100 MHz borehole antennas. The horizontal distance between the two boreholes was 18 m. Forty source positions with 1 m intervals were placed between depths of 9 m and 48 m. The receiver positions corresponding to each source position were different (

Table A1. Structures of field data.

Source Index	Source Position (m)	Traces Number	First Receiver Position (m)	Final Receiver Position (m)
1	9	43	8	29
2	10	45	7.5	29.5
3	11	46	7	29.5
4	12	45	7	29
5	13	46	6.5	29
6	14	48	6.5	30
7	15	47	6.5	29.5
8	16	46	7	29.5
9	17	47	6.5	29.5
10	18	46	7	29.5
11	19	47	7	30
12	20	47	7	30
13	21	46	7.5	30
14	22	44	8.5	30
15	23	43	8	29
16	24	41	9	29
17	25	39	10	29
18	26	38	10.5	29
19	27	38	12	30.5
20	28	36	13	30.5
21	29	61	14	44
22	30	37	26	44
23	31	41	26.5	46.5
24	32	47	23.5	46.5
25	33	49	23	47
26	34	40	27	46.5
27	35	40	27	46.5
28	36	42	25.5	46
29	37	43	25.5	46.5
30	38	40	27	46.5
31	39	43	25	46
32	40	38	27.5	46
33	41	39	27	46
34	42	36	29	46.5
35	43	37	28.5	46.5
36	44	33	30.5	46.5
37	45	33	30.5	46.5
38	46	33	30.5	46.5
39	47	31	31.5	46.5
40	48	22	36	46.5

). The interval of receiver positions was 0.5 m. The raw data measured at the site is shown in Figure 7.

The estimation of source wavelet takes an important role in the waveform inversion of field data because the source wavelet must be known for forward modeling. However, the source wavelet is commonly unknown. We followed the deconvolution method proposed by Ernst et al. [15] to estimate the source wavelet, which has been presented in previous studies [15,18].

Before starting Laplace domain waveform inversion, the field data was subjected to careful preprocessing to improve the signal-to-noise ratio (SNR), including removal of the direct-current (DC) and band-pass (BP) filtering. In addition, careful muting of the signals before the first arrival event is needed due to the characteristics of the Laplace domain wavefield [20]. The dominant frequency of the field data is approximately 50 MHz due to the water-saturated environment. Based on the conclusion presented in Section 3.1, we chose 200 × 10⁶, which is approximately four times larger than the dominant frequency of the field data, as the damping constant in the Laplace domain waveform inversion.

Figure 8a,c show the inverted results of the ray-based inversion and Figure 8b,d show the results of the Laplace domain waveform inversion after 83 iterations. From Figure 8a,b, we find that the relative permittivity tomograms from the ray-based inversion and Laplace domain waveform inversion are similar. A tunnel or cave exists at a depth of ~25 m and runs through the two boreholes. The upper right and bottom left regions have high permittivity values, whereas the permittivity of the upper left and bottom right regions is low (Figure 8a,b). However, the conductivity from the ray-based inversion (Figure 8c) shows low values (<10⁻³), which are probably unrealistic due to the existence of underground water. The reason is that the ray-based inversion can only provide relative spatial distributions of conductivity [7]. Conductivity results from the Laplace domain waveform inversion are significantly improved compared with that from the ray-based inversion.

Figure 9 shows the results of the time domain FWI by using the initial models from both ray-based inversion and Laplace domain waveform inversion. Figure 10 shows the root mean square (RMS) error curves of the field inversions with the ray-based (blue) and Laplace-based (red) results as the initial models. The convergence curves of the ray-based FWI and Laplace-based FWI stop after 43 and 69 iterations, respectively. The increase at the 12th iteration in the RMS curve of the ray-based FWI is caused by an additional source estimation during the inversion. We find that the Laplace-based FWI converges to a little lower value than the ray-based FWI. Compared to the initial models, we find obvious improvements in the results of the time domain FWI in Figure 9. At the depth of ~20 m, ~36 m and ~45 m, conductivity tomograms show horizontal layers with high values. We note that the tomograms from the Laplace-based FWI reveal more details than that from the ray-based FWI, especially for the conductivity tomograms. Figure 11 shows the predicted receiver gathers based on the inverted results for the 36th source position. Figure 11a,c are the predicted waveforms using the ray-based and Laplace domain waveform results. Though different initial models are used, waveforms of radar traces generated from the two FWI results agree well with the field waveforms (Figure 11b,d). The comparisons of waveforms indicate that inversions using both ray-based and Laplace domain waveform results are successful in the time domain. Initial models from the Laplace domain waveform inversion improve their corresponding results of FWI with more details (Figure 9b,d).

4. Discussion

FWI has been applied to interpret GPR data for approximately 20 years [13,14,15,16,17,18]. However, building accurate initial models is a challenge in the development of FWI. This is because the gradient-based FWI easily falls into local minima due to high nonlinearity. Accurate initial models are needed for FWI to overcome local minima to find the global minima. The traditional method to obtain initial models for GPR FWI is ray-based inversion, such as the first arrival time method [6] and the centroid frequency downshift method [7]. The first arrival time method can provide smooth distributions of relative permittivity, which requires accurate picking of first arrival times. However, the centroid frequency downshift method can only provide relative spatial distributions of conductivity [7]. Waveform inversion in the Laplace domain show a small advantage on permittivity tomograms because the travel time inversion is a well-established technology. More importantly, the Laplace domain waveform inversion can provide an improved distribution of conductivity, which contributes to the improvement in the time domain FWI.

Due to the characteristic of the Laplace domain wavefields, noise is amplified with the damping constants. Fortunately for cross-hole radar data, Laplace domain waveform inversion at a moderate damping constant, which is less sensitive to noise, is enough to provide initial models for the conventional FWI. Waveform inversion in the Laplace domain is robust because only a small portion of information about the data is used in the inversion. In addition, there is no need to extract first arrival times and amplitudes in the Laplace domain waveform inversion, making it less artificial intervention than the ray-based inversion. Therefore, the Laplace domain waveform algorithm has potential in the development of GPR FWI.

5. Conclusions

We successfully applied Laplace domain waveform inversion to synthetic and field cross-hole radar data. For the Laplace domain waveform inversion, the algorithm was derived via the derivation of Maxwell’s equation and a logarithmic object function. In the inversion, permittivity and conductivity can be obtained together by using a stepped updating strategy. The choice of damping constant in the inversion is investigated using a simple synthetic model. It is shown that the optimum value of the damping constant is approximately five times the value of the dominant frequency of radar data. Results of the time domain FWI based initial models from the ray-based inversion and Laplace domain waveform inversion indicate that the Laplace domain waveform inversion performs well. In the end, we applied the algorithm to a cross-hole field data. In the results of the field data, the permittivity and conductivity distributions from the time domain FWI using both ray-based and Laplace domain waveform results are comparable. Waveform inversion in the Laplace domain is an alternative algorithm to provide accurate initial models for GPR FWI.