Quantitative multi-layer electromagnetic induction inversion and full-waveform inversion of crosshole ground penetrating radar data

0. INTRODUCTION

The electromagnetic tools, EMI and GPR, can be used for a wide range of applications to non-invasively image the subsurface. Due to the fast data acquisition, where the measured data can be directly linked with high resolution GPS systems, it is becoming more and more feasible to map GPR and EMI over large areas. The use of multi-channel EMI and GPR makes it possible to acquire simultaneously EMI and GPR data for different source-receiver offsets that enable an improved subsurface characterization. Ray-based or approximate techniques are often used where only part of the data is exploited. Improved subsurface characterization can be obtained by including advanced modeling tools that are able to calculate the electromagnetic wave propagation with high accuracy using Maxwell's equations. In the following, we will describe recent advancements in the high-resolution imaging of EMI and GPR data.

1. QUANTITATIVE MULTI-CONFIGURATION ULTI-LAYER EMI INVERSION

EMI measurements are often performed with single source-receiver systems to qualitatively map the spatial apparent electrical conductivity (ECa) patterns and to relate them with the soil properties (e.g., Abdu et al., 2008). Recently, new generation multi-configuration EMI systems have become commercially available which consist of multiple coils with different separations and different orientations (Rudolph et al., 2015; Stadler et al., 2015; Saey et al., 2012; Monteiro Santos et al., 2010). Since the measured ECa is a weighted average value over the coil configuration specific penetration depth, this enables a simultaneous sensing over different depth ranges. In principle, these multi-configuration data can be used to obtain electrical conductivity changes with depth. However, for a reliable inversion, the measured ECa values needs to be quantitative, which is often not the case because EMI data can significantly be affected by external influences that alter EMI measurements, such as the presence of the operator or cables near the coils, and the field set-up used to move the EMI system (Nüsch et al., 2010; Robinson et al., 2004). One option to obtain reliable ECa measurements with EMI is to calibrate the data using inverted ERT data measured together with EMI along the same profile. Using this data as input in a Maxwell-based full-solution forward model, the predicted apparent electrical conductivities are used to calibrate the measured data using linear regression parameters (Lavoué et al., 2010).

Instead of applying an inversion algorithm that is based on the low induction number approximation, we invert large-scale multi-configuration EMI data using a Maxwell-based full-solution EM forward model to resolve the 3D multi-layer electrical conductivity structures. Based on the two-layer inversion as discussed by Mester et al. (2011), an improved three-layer inversion scheme was developed (von Hebel et al., 2014) that explores the information content of M=6 different EMI data configurations and minimizes the normalized misfit between the measured magnetic field H_m^mea, and the modeled magnetic field H_m^mod(p)

Here, the modeled magnetic field is calculated using a full-solution horizontally layered EM (electromagnetic) forward model (von Hebel et al., 2014), and p contains the three-layer parameters σ₁, σ₂, σ₃, h₁, and h₂ that are being optimized using an efficient parallelized global minimization algorithm, the shuffled complex evolution (SCE) algorithm (see von Hebel et al. (2014) for more details about the used inversion algorithm).

At the Scheyern test site (Germany), ECa measurements were performed using the CMD-Mini Explorer (GF instruments, Brno, Czech Republic). The sensor consists of three receiver coils separated by 0.32, 0.71, and 1.18 m from the transmitter coil resulting in a theoretical depth of exploration (DOE) of 0.25, 0.5, and 0.9 m in the vertical coplanar (VCP) and 0.5, 1.1, and 1.8 m for the horizontal coplanar (HCP) mode, respectively.

Figure 1 shows interpolated large-scale ECa values measured for a coil distance of 1.18 m and a VCP configuration which clearly indicates lateral changes in apparent electrical conductivity. Due to the influence of the operator and cables close to the EMI sensor, the measured apparent conductivity values can only be qualitatively interpreted. Inverting the multi-configuration EMI data results in unreliable inversion results (not shown). To calibrate the EMI data, reference line P03 was selected that contains lateral changes in conductivity as shown in Fig. 1. Here, ERT measurements were performed over 30 m using electrode distances of 0.25 m. Figure 2 clearly shows the lateral and vertical changes in conductivity obtained by inverting the ERT data. The ERT data were used to calibrate the measured data (Lavoué et al., 2010) and to obtain quantitative EMI data. The calibrated data were inverted using a novel three-layer inversion that uses the shuffled complex evolution (SCE) optimization and a Maxwell-based electromagnetic forward model (von Hebel et al., 2014). The three-layer inversion results that are shown in Fig. 3 approximate well the quantitative values for the different layers of the ERT structures. The dipping 2D structure between 20 and 25 m is less well reconstructed due to the horizontal layer assumption. Note that the EMI measurements were carried out in only a fraction of the time needed to perform the ERT measurements. Recently comma it was shown that the obtained calibration parameters can be applied to large scale data, where for each grid location the experimental EMI data are independently inverted and stitched together to form a quasi-3D image (von Hebel et al., 2014).

Figure  1.  Interpolated large-scale ECa measurements using the VCP configuration and a source-receiver distance of 1.18 m. Along the P03 transect, where large ranges of ECa values were presented, an ERT profile of 30 m was measured to enable the calibration of the EMI data. The ERT inversion results are shown in Fig. 2.

DownLoad: Full-Size Img PowerPoint

Figure  2.  ERT-inversion results along the line P03 in Fig. 1, where lateral and vertical changes in electrical conductivity are clearly visible.

DownLoad: Full-Size Img PowerPoint

Figure  3.  Inversion results of the calibrated multi-configuration EMI data that show similar quantitative conductivity values for the different layers as indicated by the ERT structures in Fig. 2.

DownLoad: Full-Size Img PowerPoint

2. FULL-WAVEFORM INVERSION OF CROSSHOLE GPR DATA

Within the seismic community the full-waveform inversion (FWI) has found widespread applications (Virieux and Operto, 2009). In contrast, FWI methods for GPR are still fairly new. One of the first implementations of crosshole GPR full-waveform inversion (Ernst et al., 2007) is a conjugate gradient optimization algorithm in time domain. A vectorial approach was implemented that honors the vectorial electromagnetic wave propagation and simultaneously updates the permittivity and conductivity (Meles et al., 2011, 2010). Recently, a modified misfit and corresponding gradient equations have been introduced (Yang et al., 2013a) that do not depend anymore on the number of sources and receivers used. The modified misfit and adapted by Yang et al. (2013a) is defined as

where E_syn^s(ε,σ) and E_obs^s are the synthetic and observed data, respectively, that contain the data for all sources, receivers and observation times according to the Meles et al. (2010) formalism. Here, T denotes the transpose conjugate and the synthetic data are calculated using a 2D finite-difference time-domain (FDTD) solution of Maxwell's equations. The gradients of the misfit function with respect to the permittivity and conductivi-ty $\nabla \bar{S}_{\varepsilon}$, and $\nabla \bar{S}_{\sigma}$ are obtained by cross-correlating the incident wave field emitted from the source with the residual wave fields ΔE^S that are back-propagating using the back-propagation operator $\hat G^T$ from the receiver at all medium locations in the time domain for all sources and receiver combinations as follows (see also Yang et al., 2013b; Meles et al., 2010)

Using these gradient directions, the update directions d_ε and d_σ are computed using the conjugate gradient method whereas the step lengths ζ_ε and ζ_σ are computed using a linear step length calculation and carefully chosen perturbation factors (see also Yang et al., 2013a; Meles et al., 2010). Finally, the model parameters are updated using

where k is the iteration number.

To apply the full-waveform method to experimental data and to invert the data using 2D forward model, a 3D to 2D conversion is required and an effective source wavelet must be estimated (Klotzsche et al., 2010). The starting relative permittivity and conductivity tomograms are usually obtained using a standard ray-based inversion that is employing the first-arrival times and first-cycle amplitudes. Applying the full-waveform inversion on several crosshole GPR data sets acquired in gravel aquifers in Switzerland and the U.S.A. shows that sub-wavelength thickness low-velocity wave guiding layers can be correctly imaged (Klotzsche et al., 2014, 2013, 2012), whereas ray-based inversion techniques are not able to image these thin waveguide layers because they only exploit the first-arrival times and first-cycle amplitudes. Converting the permittivity results into porosity and comparing them with Neutron-Neutron logging data showed a good correspondence (Klotzsche et al., 2014, 2013).

Different versions of the full-waveform inversion were compared and applied to an experimental data set acquired at the Boise Hydrogeophysics Research Site (Yang et al., 2013a). The obtained results show that most improvements result in a reduced final misfit between the measured and synthetic data and a reduced remaining gradient at the final iteration. Regions with relatively high remaining gradient amplitudes indicate less reliable inversion results. Frequency-domain full-waveform approaches have been recently implemented by Yang et al. (2013b) and Lavoué et al. (2014). Although here, selected frequencies can be used for the inversion to reduce the computational demands, selecting the appropriate frequencies is not straightforward.

Recently, many crosshole measurements were carried out at the Krauthausen test site in Germany. For a transect covering 25 m length and 10 m depth consisting of 5 crosshole planes, GPR data were inverted using the ray-based and full-waveform inversion as described by Klotzsche et al.(2013, 2010) and Oberröhrmann et al. (2013). Densely spaced cone penetration tests (CPT), located close to the GPR transect, were used to validate and interpret the obtained images obtained from GPR. A strong correlation was observed between CPT porosity logs and porosity estimates derived from GPR, whereas a less pronounced correlation was observed between electrical conductivity data derived from GPR and CPT (Gueting et al., 2015).

Here, we compare the obtained ray-based relative permittivity and conductivity inversion results shown in Fig. 4 with the corresponding full-waveform inversion results shown in Fig. 5. The full-waveform inversion results clearly show improved resolution images compared to the ray-based images. Since the gradients have a high sensitivity close to the sources, artefacts can arise close to the boreholes resulting in unreliable inversion results. Therefore, the gradients of 7 cells (7×9 cm= 63 cm) around each borehole are muted during the inversion (see Fig. 5). However, relatively large amplitude conductivity changes are still present close to the boreholes.

Figure  4.  (a) Ray-based relative permittivity and (b) conductivity inversion results of five crosshole planes covering a total range of 25 m length and 10 m depth.

DownLoad: Full-Size Img PowerPoint

Figure  5.  (a) Relative permittivity and (b) conductivity results obtained by full-waveform inversion. The white areas between the planes indicate the region where the gradients are muted.

DownLoad: Full-Size Img PowerPoint

To reduce these artefacts and prevent a relatively large tapering, the gradients $\nabla \bar{S}_{\varepsilon}$ and $\nabla \bar{S}_{\sigma}$ are preconditioned with the maximum values of the forward propagated field and the back-propagated residual field. The preconditioner operator P^k for the updating domain x is defined by (see also Kurzmann et al., 2013).

where

and

Here, a denotes the spatial average of a(x), and C_stab=10 stabilizes the computation of P^k and was manually estimated (see also Kurzmann et al., 2013).

Applying this gradient enables a reduction of the muted domain to 3 cells (3×9 cm= 27 cm), which corresponds approximately to the drilling diameter of about 30 cm. Note that the soil material within this drilling diameter is strongly disturbed due to the drilling, casing, water-filled plastic tubes, etc.. The obtained full-waveform inversion results (see Fig. 6) contain smaller blank areas around the boreholes and the artefacts close to the boreholes are strongly reduced. More research is needed to optimize the parameterization of this preconditioning and to compare the obtained results with ground truth borehole Neutron-Neutron and electrical conductivity logs.

Figure  6.  (a) Relative permittivity and (b) conductivity results obtained by full-waveform inversion using a preconditioned gradient. The white areas between the planes are strongly reduced and fewer artefacts are present close to the boreholes.

DownLoad: Full-Size Img PowerPoint

Figure 7a shows the measured data for receivers present in borehole 30 and sources present in borehole 62 and allow a comparison with the modeled data obtained using the (b) ray-based inversion, (c) full-waveform inversion without gradient preconditioning, and (d) the full-waveform inversion with gradient preconditioning. The modeled data using the ray-based inversion results show similar first arrival times, however, the small-scale amplitude changes present in the measured data are not reconstructed. This is due to the use of only the first arrival time and first cycle amplitudes to obtain the inversion results shown in Fig. 4. In contrast, the modeled data using the full-waveform inversion without and with the preconditioning does explain the measured data well and show similar travel-time and amplitude changes. Note that all the data shown in Fig. 7a is used to obtain the full-waveform inversion results shown in Figs. 5 and 6.

Figure  7.  (a) Measured data for receivers present in borehole 30 and sources present in borehole 62, modeled data using the (b) ray-based inversion (c) full-waveform inversion without gradient preconditioning, and (d) full-waveform inversion with gradient preconditioning.

DownLoad: Full-Size Img PowerPoint

3. CONCLUSION AND OUTLOOK

High-resolution, sub-wavelength information can be obtained from EMI and GPR data by using inversion tools that are based on exact vectorial electromagnetic wave approaches. A novel three-layer EMI inversion algorithm has been used that employs the different sensing depths of multi-configuration EMI data. After proper calibration, the three-layer medium properties are reconstructed using an efficient parallelized global minimization algorithm and a full-solution horizontally layered EM forward model. Inversion of these calibrated multi-configuration EMI data provides lateral and vertical conductivity variations very similar to those observed in an elaborate ERT experiment. Since for every grid location the data are independently inverted and stitched together to form a quasi-3D image, this inversion could be parallelized to invert large scale EMI data. Full-waveform inversion of crosshole GPR data returns significantly higher resolution images than conventional ray-based images and better depicts the medium properties within the subsurface. The use of a preconditioning of the gradient reduces artefacts close to the source locations and improves the obtained results. A good amplitude and phase fit between the measured and modeled data using the full-waveform inversion results was obtained, indicating the reliability of the full-waveform inversion results. Currently the full-waveform inversion is being extended for on-ground GPR data which can then be combined with the multi-configuration EMI inversion and will extend the possible range of applications significantly. The obtained results show novel inversion algorithms that are based on a high resolution modeling program can exploit the full information content present in the data which enables the unprecedented high resolution characterization of the subsurface.