Analysis of Forward Model, Data Type, and Prior Information in Probabilistic Inversion of Crosshole GPR Data

: The crosshole ground penetrating radar (GPR) is a widely used tool to map subsurface properties, and inversion methods are used to derive electrical parameters from crosshole GPR data. In this paper, a probabilistic inversion algorithm that uses Markov chain Monte Carlo (MCMC) simulations within the Bayesian framework is implemented to infer the posterior distribution of the relative permittivity of the subsurface medium. Close attention is paid to the critical elements of this method, including the forward model, data type and prior information, and their inﬂuence on the inversion results are investigated. First, a uniform prior distribution is used to reﬂect the lack of prior knowledge of model parameters, and inversions are performed using the straight-ray model with ﬁrst-arrival traveltime data, the ﬁnite-difference time-domain (FDTD) model with ﬁrst-arrival traveltime data, and the FDTD model with waveform data, respectively. The cases using ﬁrst-arrival traveltime data require an unreasonable number of model evaluations to converge, yet are not able to recover the real relative permittivity ﬁeld. In contrast, the inversion using the FDTD model with waveform data successfully infers the correct model parameters. Then, the smooth constraint of model parameters is employed as the prior distribution. The inversion results demonstrate that the prior information barely affects the inversion results using the FDTD model with waveform data, but signiﬁcantly improves the inversion results using ﬁrst-arrival traveltime data by decreasing the computing time and reducing uncertainties of the posterior distribution of model parameters.


Introduction
The crosshole ground penetrating radar (GPR) method has been widely applied for geological, environmental, and engineering investigations to characterize subsurface properties [1][2][3][4][5][6]. This technique uses high-frequency electromagnetic (EM) waves and transmission measurements through two boreholes to retrieve the spatial distribution of electrical properties, including mainly the dielectric permittivity and electrical conductivity, which determine the geometrical and physical features of the subsurface medium [7].
The inversion of crosshole GPR data involves estimating the dielectric parameters from the information provided by crosshole measurements [8]. In order to search for model parameters, deterministic inversion algorithms are extensively used to find an optimal solution that best fits the observed data [9][10][11]. In contrast, probabilistic inversion methods are capable of providing an ensemble of solutions that are statistically acceptable and quantifying parameter uncertainties [12][13][14][15][16][17]. Probabilistic inversion treats parameters and data as random variables, and Bayes theorem is often implemented to derive the posterior distribution of model parameters from the prior distribution and measured data [18]. The prior distribution reflects the knowledge of model parameters before observations, and an informative prior helps to reduce inversion uncertainties and improve the parameter searching efficiency [19,20]. The likelihood function measures how well the model fits data by performing forward simulations and computing the distance between the simulated and measured data. According to the forward assumptions, the ray-based forward kernels use the first-arrival traveltime or first-cycle amplitude data, while the wave-based simulators employ the waveform data. The ray-based forward models are usually computational efficient, yet utilize only a small portion of the measured data that may bias the inversion result [21]. Instead, the wave-based forward models (e.g., the finite-difference time-domain (FDTD) simulator of Maxwell's equations) take advantage of the full-waveform data for a better characterization of model parameters at the price of considerable computational resources [22][23][24].
In our previous work, a Bayesian inversion method was developed to infer the relative permittivity values of underground structures from crosshole GPR data [25]. This method used the FDTD forward simulator with waveform data and achieved better accuracy than the traditional ray tomography. Yet one main concern of this approach is that it requires thousands to millions of forward model evaluations for the parameters to converge, leading to dozens of hours of inversion time, which limits the practicability of this method. In a later study, a straight-ray forward model was adopted and the computational efficiency was significantly improved [26]. In probabilistic inversion, the prior information is also an important aspect which is capable of reducing ambiguity by imposing constraints on the model parameters [27,28], yet improper prior distribution may distort the inversion result.
Based on the previous research, this paper focuses on the key elements of probabilistic inversion of crosshole GPR data, and systematically investigates the impact of the forward model, data type, and prior information on the inversion results. The remainder of this paper is structured as follows. First, the formulation of the probabilistic inversion method is briefly introduced. Next, the inversion results using the straight-ray model with firstarrival traveltime data, FDTD model with first-arrival traveltime data, and FDTD model with waveform data are investigated. Then, the influence of the prior information on the inversion results are analyzed, and followed by a discussion of the strengths and weaknesses of this method. Finally, this paper is concluded with a summary of the main findings.

Formulation of Probabilistic Inversion
The inversion of crosshole GPR data aims at estimating model parameters through measured data. In a Bayesian sense, the inverse problem can be formulated as where p(m) and p(m|d) are the prior and posterior probability density of model parameters m, respectively, L(m|d) ≡ p(d|m) is the likelihood function, and p(d) is the probability density of data d that serves as a normalizing constant. In this study, the model parameters refer to the relative permittivity (ε r ) values of the subsurface structure, while the data are the first-arrival traveltimes or waveforms obtained by crosshole GPR measurements. The Bayesian approach treats all information as random variables and infers the posterior distribution of model parameters from the combination of prior knowledge and measured data.

Prior
The prior distribution summarizes all the model parameter information before collecting crosshole GPR data. In this work, a uniform distribution was used because of the lack of knowledge about the model parameters. In practice, the changes in the spatial distribution of underground media occur gradually. Hence, a smooth constraint of the model parameters with the following form is feasible to be included in the prior information to confine the model structure [27].
where W x and W z are the difference operators in x and z directions with rank R x and R z , respectively, and λ is the standard deviation of the model gradient.

Likelihood
The likelihood function describes the distance between the simulated and measured crosshole GPR data. Assuming that the measurement errors are independent and identically distributed with a normal distribution of zero mean and standard deviation of σ, the likelihood function can be written as where f i (m) and d i are the i-th simulated and measured crosshole GPR data, respectively. n denotes the total number of data. To prevent the value of Equation (3) from approaching zero when using a large number of data, it is feasible the following logarithmic form (4)

Forward Model
Evaluation of Equation (4) requires a forward model f (m) to generate simulated data. The crosshole GPR uses EM waves for subsurface characterization, and the wellknown Maxwell's equations govern this process. In two-dimensional space, the Maxwell's equations can be simplified in transverse electric (TE) mode (5) where E x and E z are the electric fields in the x and z directions, and H y signifies the magnetic field in the y-direction. ε, σ e and µ are the dielectric permittivity, electric conductivity, and magnetic permeability of the subsurface medium, respectively.
As analytical solutions of Equation (5) are always tricky, numerical methods are usually used to resolve the electromagnetic fields. In this work, the FDTD algorithm is implemented for numerical solutions [29]. This approach uses central difference to discretize the Maxwell's equations in both the space and time domains, and is capable of simulating the electromagnetic field at each time step for a given time duration with high precision. Therefore, the first-arrival traveltimes, first-cycle amplitudes, and full-waveforms, which are used for the crosshole GPR inversion, can be obtained with an FDTD solver. However, within the Bayesian inversion framework, thousands to millions of forward model evaluations are needed for the model parameters to converge to their target distribution, and this CPU-intensive forward kernel makes the inversion a time-consuming task.
To improve the computational efficiency of the probabilistic inversion method, a straight-ray based forward approach is also considered [30]. This ray-based approach assumes that the propagation of the EM wave from the source to receiver is along a straight raypath l, and computes the first-arrival traveltime t by where i denotes the i-th measurement, and s signifies the slowness of the EM wave along the raypath. s is the reciprocal of velocity v derived by v = c/ √ ε r , in which c is the EM wave velocity in vacuum. Equation (6) simplifies the EM wave propagation into a straight ray and runs much faster than the FDTD model. However, with the straight-ray model, only the first-arrival traveltime data can be simulated [31].

Posterior
MCMC simulation is adopted to derive the posterior distribution of model parameters. This method uses a Markov chain that generates a trial move from the current state, m (i) to a new position, m * , and decides whether to accept this jump using the Metropolis ratio [32] where p a (·) is the acceptance probability. In this paper, the DREAM (ZS) package [33], which is an adaptive MCMC algorithm designed for high-dimensional problems, was used to explore the posterior distribution of model parameters [34].

Reference Model and Data
In this section, a synthetic model is created for simulating crosshole GPR detection of a defect in an underground concrete structure. As shown in Figure 1a, a 1.0 m × 1.0 m model was built with ε r = 9 to simulate the concrete structure (background), and a 0.2 m × 0.2 m higher relative permittivity (ε r = 12) area was created as the defect (target). Crosshole GPR measurements were simulated using 11 transmitters on the left side and 26 receivers on the right side of the model.
In order to obtain crosshole GPR data, the gprMax software [35,36], which is an FDTD solver of the Maxwell's equations, was used to perform numerical simulations. The model depicted in Figure 1a is discretized, with a grid size of 0.02 m × 0.02 m, leading to a total number of 50 × 50 = 2500 unknown relative permittivity values, subject to estimate. The FDTD simulations use the Ricker wavelet with a center frequency of 500 MHz as the source, and collect a 20 ns GPR waveform for each transmitter-receiver pair. At each transmitter position, GPR waveforms were collected at all locations of the 26 receivers. A total of 11 × 26 = 286 GPR traces were obtained as waveform data, and the first breaks of the GPR traces were extracted as the first-arrival traveltime data. Figure 1b,c display the synthetic traveltime and waveform data, respectively. Then, these data were perturbed with white noise and served as the measurement data. For the first-arrival traveltime data, a zero mean and standard deviation of 3% of the mean value (0.24 ns) Gaussian noise was used as the measurement error. For the waveform data, the signal-tonoise ratio (SNR) was used to control the noise level, and a Gaussian random noise with SNR = 20 was generated to simulate the measurement error, which was equivalent to a standard deviation of 2.06.

Effect of Forward Model and Data Type
In this section, the impact of the forward model and data type on the inversion results is investigated. The straight-ray model computes only first-arrival traveltime data, whereas the FDTD model simulates the EM fields, from which both first-arrival traveltime and waveform data can be obtained. Therefore, three cases are discussed in this section, which are the straight-ray model with first-arrival traveltime data, the FDTD model with first-arrival traveltime data, and the FDTD model with waveform data.
The forward model and discretized grid size uses the same values of the reference model shown in Figure 1a, thus a total number of 2500 relative permittivity values need to be estimated. The discrete cosine transform (DCT) strategy was adopted for dimensionality reduction [37,38], and 256 DCT-coefficients were used to represent the full parameter space. The simulated data were computed using source and receiver positions identical to those used in the synthetic measurements, which produced 11 × 26 = 286 first-arrival traveltimes for the straight-ray model, and 286 GPR traces for the FDTD model. For each GPR trace, the waveform from 10 to 20 ns, which contains most of the information, was used as the simulated waveform data. Each waveform consists of 212 data points, and a total of 286 × 212 = 60, 632 data points constitute the waveform data for the FDTD model.
MCMC simulation with the DREAM (ZS) algorithm was deployed to sample the posterior distribution of model parameters. A uniform distribution of model parameters was used in this section as non-informative prior, and Equation (4) was adopted as the likelihood function. The setup of the probabilistic inversions with different forward models and data types is summarized in Table 1. The inversion results are discussed in the following subsections.

Straight-Ray Model with Traveltime Data
First, the straight-ray forward model with first-arrival traveltime data was used to infer the relative permittivity field. After 1,440,000 model evaluations, all parameters converged to their target distribution, and the last 50% of the parameters after convergence were considered as posterior samples. Figure 2 displays the histograms of the prior (blue bars) and posterior (red bars) distribution of three randomly selected parameters. Before inversion, the parameters were uniformly distributed in the prior range. After inversion, the parameters moved toward the true values marked with black crosses. However, the mean values were estimated with large errors, and the posterior parameters were distributed in a wide range with large variance. Figure 3a-c present three randomly selected posterior realizations. They show significant differences in model structure, and none of them resembles the true relative permittivity model. Figure 3d,e illustrate the maximum a-posteriori (MAP) and posterior mean realizations, respectively, which fail to reconstruct the correct model structure again. Figure 3f shows a large variance of posterior model parameters, which means the inverted parameters have large uncertainties. Therefore, with the straight-ray forward model and first-arrival traveltime data, the model parameters converge to a broad range, and the posterior models do not recover the relative permittivity field correctly.    Figure 4 shows the data fit of the inverted models. It can be seen that the posterior simulated data are in good agreement with the measured data, and the root mean squared error (RMSE) values of the posterior models are distributed around the measurement error (0.24 ns), which means the posterior models fit the data very well. However, the posterior realizations differ significantly from the actual model (see Figure 3), indicating that a large number of different model parameters exists that explain the measured data very well. Therefore, the information provided by the first-arrival traveltime data is insufficient to retrieve the correct relative permittivity field.

FDTD Model with Traveltime Data
In the second experiment, the FDTD forward model with first-arrival traveltime data was evaluated. A total of 1,638,000 iterations were required for all parameters to converge to their posterior distribution. Figure 5 summarizes all posterior samples and illustrates the prior and posterior distribution of three randomly selected parameters. After Bayesian inversion, the model parameters were sampled around their actual values. Compared with the results using the straight-ray forward model, the posterior model parameters derived from the FDTD model converged to a smaller range with reduced variances. The posterior reconstructed relative permittivity models are displayed in Figure 6. Similar to the realizations obtained from the inversion using the straight-ray forward model, the FDTD model with first-arrival traveltime data also fails to resolve the reference model correctly. Although the FDTD forward model reduces parameter uncertainties, it is not sufficient to produce relative permittivity models that resemble the true one.
The simulated first-arrival traveltimes of the posterior parameters using the FDTD forward model are shown in Figure 7a. Unlike the case using the straight-ray forward model, the posterior simulated data using the FDTD model are distributed in a much smaller range close to the measured data, and the RMSE values of the posterior models are slightly higher than the measurement error (see Figure 7b). In this case, the posterior models using the FDTD forward simulations with first-arrival traveltime data explain the observed data very well with smaller uncertainties, yet they still failed to converge to the actual relative permittivity field.

FDTD Model with Waveform Data
Next, the use of the FDTD forward model with waveform data to infer the relative permittivity field was investigated. All the parameters converged to their target distribution after 156,000 model evaluations, much less than the cases using the first-arrival traveltime data. The histograms of the prior and posterior distribution of three randomly selected parameters are presented in Figure 8, in which the posterior parameters are distributed in a very narrow range and focus on the true values. The posterior realizations shown in Figure 9 successfully reconstruct the proper relative permittivity field. The MAP and posterior mean models present similar patterns and are visually identical to the randomly selected posterior realizations. Small variances are observed, indicating that the waveform data provide sufficient information for the inversion algorithm to retrieve the correct model parameters.   Figure 10b displays the distribution of the RMSE values of the posterior solutions, which is close to the real value of the measurement error. Therefore, the FDTD forward model with waveform data is capable of characterizing the relative permittivity field. The inversion results of the three cases using uniform prior distribution are summarized as follows. With first-arrival traveltime data, neither the straight-ray model nor the FDTD model reconstructs the relative permittivity field correctly. When using the waveform data, the model parameters converge to the real values, and the model structure is recovered correctly. This is because there are 2500 unknown model parameters, yet the first-arrival traveltime data contain only 286 data points, far less than the number of unknowns. Without additional knowledge about model parameters, a large number of permittivity models exist that fit the data well. Therefore, the information provided by the first-arrival traveltime data is not enough to lead the model parameters to their correct values. On the contrary, the waveform data, which contain 60,632 data points, are much more informative to infer the model parameters correctly.

Effect of Prior Information
In this section, the informative prior defined in Equation (2) is considered, and the impact of prior information on the inversion results is discussed. Three inversion cases including the straight-ray model with first-arrival traveltime data, the FDTD model with first-arrival traveltime data, and the FDTD model with waveform data are analyzed. Except for the prior distribution, the same inversion setup was used as that in the previous section. Table 2 summarizes the inversion parameters of the three cases.
When the prior information was accounted for, it cost 144,000, 236,800, and 156,000 model evaluations for the three cases to converge to their posterior distribution. For the cases using the first-arrival traveltime data, much fewer iterations were needed compared with their flat-prior counterparts. However, when the waveform data were used, the prior distribution had little influence on the convergence rate. Figure 11 illustrates the distribution of a randomly selected parameter for each case. It is obvious that for all three cases, the parameter converges to a small range and is distributed around the actual value.
The prior information reduces the parameter uncertainty substantially for the cases using the first-arrival traveltime data, yet has little impact on the parameter distribution using waveform data.  Figure 11. Histograms of the prior (blue bars) and posterior (red bars) distribution of three randomly selected parameters derived from Bayesian inversion with smooth model constraint using (a) the straight-ray forward model with first-arrival traveltime data, (b) FDTD model with first-arrival traveltime data, and (c) FDTD model with waveform data. The true values of the parameters are denoted with black crosses. Figure 12 displays the posterior realizations using the straight-ray model with firstarrival traveltime data, the FDTD model with first-arrival traveltime data, and the FDTD model with waveform data, respectively. For the two cases using the first-arrival traveltime data, the use of prior information improves inversion results significantly so that the higher permittivity area can be identified in the MAP and posterior mean realizations clearly. The MAP models exhibit better spatial resolution, yet suffers from more noise. The posterior mean models, on the other hand, present a smoother pattern. The posterior model variance also becomes much smaller compared with the results using the uniform prior distribution. For the case using the FDTD model with waveform data, the result is little affected by the prior information.
In Figure 13, the posterior simulated data of the three cases are in good agreement with the measured data. The use of prior information refines the convergence of parameters to the correct values, but does not affect the fitting capacity. The RMSE values of the three cases are close to the correct values and distributed in a small range compared with the uniform prior cases.
Taken together, Table 3 summarizes the performance of the above numerical experiments involving different forward models, data types, and prior distribution. The number of iterations and corresponding inversion time measure the computational efficiencies, the RMSE values of the posterior mean simulations evaluate how well the simulated models fit the measured data, and the peak signal-to-noise ratios (PSNR) quantify the similarity between the posterior mean models and the real relative permittivity model [39]. With uniform prior, the inversions using first-arrival traveltime data (Cases 1 and 2) require more than one million iterations for the parameters to converge, yet the realizations do not resemble the true relative permittivity field with low PSNR values of 17.22 and 19.30, respectively. In contrast, the inversion using the FDTD model with waveform data cost only 120,000 model evaluations (13 times less than those of Case 2) and achieved a better in-version accuracy, with a PSNR value of 28.95. The results indicate that, without additional knowledge, the first-arrival traveltime data are unable to provide sufficient information to infer the model parameters. Although the simulated models nicely fit the measurement data with RMSE values close to the measurement error (0.24 ns), the inverted model structures differ significantly from the true model. When the smooth constraint was taken into consideration as prior information, the inversion results using first-arrival traveltime data (Cases 4 and 5) improved considerably by reducing the number of iterations to 144,000 and 236,800, and increasing the PSNR values to 26.83 and 26.59, respectively. Yet, the inversion results using the FDTD model with waveform data (Cases 3 and 6) were barely affected by the prior information, indicating that the waveform data contain sufficient information to infer the model parameters. The choice of forward model also influenced the computational efficiency greatly. By comparing Cases 4 and 5, the straight-ray model inversion achieved almost the same inversion accuracy (PSNR = 26.83) as that of the FDTD model (PSNR = 26.59), but runs more than 120 times faster than the FDTD model inversion.

Discussion
The probabilistic inversion with Bayesian inference is able to explicitly treat different sources of errors, providing a set of statistically acceptable solutions [25]. Compared with deterministic inversion methods that offer a single realization fitting the observed data, the probabilistic inversion results are delivered in the form of a posterior distribution of model parameters that quantifies the uncertainties of the estimation. The probabilistic inversion method summarizes the information of measured data and prior knowledge, and explores the posterior distribution of parameters using global optimization approaches (e.g., the DREAM (ZS) algorithm in this work). This process usually requires a large number of iterations (thousands to millions, depending on the number of parameters) to converge, and the results are affected by the prior distribution and data fitting. Therefore, the choice of forward model, data type, and prior distribution should be carefully considered. Although the inversion using the FDTD model with waveform data has the best accuracy, the low computational efficiency caused by the CPU-intensive FDTD forward solver is unacceptable for practical use. The straight-ray forward model, which is much more computationally appealing, reduces the inversion time remarkably, yet it uses only a small portion of the information of GPR data and leads to biased posterior distribution. An informative prior (smooth model constraint) was therefore used to confine the model structure, thus deriving an improved inversion result. The use of the smooth constraint reduces the parameter ambiguity, yet it sacrifices the inversion accuracy by forcing the parameters to vary gradually. It turns out that there is a trade-off between computational efficiency and inversion accuracy. For our cases, the inversion using the straight-ray forward model and first-arrival traveltime data with smooth model constraint achieves very high computational efficiency and acceptable inversion accuracy. Note that when the inversion problem involves larger parameter differences, greater target sizes or multiple objects, the applicability of this approach should be further investigated.

Conclusions
In this paper, the key components of the probabilistic inversion for the relative permittivity parameters with crosshole GPR data were investigated. The impact of the forward model, data type, and prior distribution on the inversion efficiency and accuracy was analyzed. The inversion results on a synthetic example demonstrate that, (1) the waveform data are much more informative than the first-arrival traveltime data, and are able to infer the model parameters accurately without additional information; (2) prior knowledge in the form of model constraint is necessary for the inversions using the first-arrival traveltime data to reconstruct the relative permittivity field correctly; (3) the straight-ray forward model is much more CPU-friendly than the FDTD model (more than 120 times faster), and better facilitates the probabilistic inversion, which requires a large number of model evaluations. Future work should involve analysis of the applicability to more complex model structures and real-world cases.

Data Availability Statement:
The data used to support the findings of this study are included within the article.