Waveform Inversion of Crosshole Georadar Data - IEEE Xplore

4610

IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 50, NO. 11, NOVEMBER 2012

Waveform Inversion of Crosshole Georadar Data: Influence of Source Wavelet Variability and the Suitability of a Single Wavelet Assumption Florian A. Belina, James Irving, Jacques R. Ernst, Associate Member, IEEE, and Klaus Holliger

Abstract—Waveform-based tomographic imaging of crosshole georadar data is a powerful method to investigate the shallow subsurface because of its ability to provide images of electrical properties in near-surface environments with unprecedented spatial resolution. A critical issue with waveform inversion is the a priori unknown source signal. Indeed, the estimation of the source pulse is notoriously difficult but essential for the effective application of this method. Here, we explore the viability and robustness of a recently proposed deconvolution-based procedure to estimate the source pulse during waveform inversion of crosshole georadar data, where changes in wavelet shape with location as a result of varying near-field conditions and differences in antenna coupling may be significant. Specifically, we examine whether a single, average estimated source current function can adequately represent the pulses radiated at all transmitter locations during a crosshole georadar survey, or whether a separate source wavelet estimation should be performed for each transmitter gather. Tests with synthetic and field data indicate that remarkably good tomographic reconstructions can be obtained using a single estimated source pulse when moderate to strong variability exists in the true source signal with antenna location. Only in the case of very strong variability in the true source pulse are tomographic reconstructions clearly improved by estimating a different source wavelet for each transmitter location. Index Terms—Convolution, geophysical inverse problems, geophysical tomography, radar imaging, wavelet.

I. I NTRODUCTION

T

HERE is a strong need for high-resolution georadar-based imaging of the shallow subsurface in many environmental, hydrological, and engineering problems [1]–[5]. To this end, waveform-based tomographic imaging of crosshole georadar data has recently gathered significant interest because of its potential for providing highly resolved images of complex nearsurface structure in terms of pertinent petrophysical parameters [6], [7]. Indeed, this technique may provide subsurface images with a resolution that is comparable to that of borehole logs

Manuscript received March 16, 2011; revised August 18, 2011 and March 10, 2012; accepted March 11, 2012. Date of publication May 15, 2012; date of current version October 24, 2012. This work was supported by a grant from the Swiss National Science Foundation. F. A. Belina, J. Irving, and K. Holliger are with the Institute of Geophysics, University of Lausanne, 1015 Lausanne, Switzerland (e-mail: florian.belina@ gmail.com; [email protected]; [email protected]). J. R. Ernst is with EGL, 8953 Dietikon, Switzerland (e-mail: jacques.ernst@ gmail.com). Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TGRS.2012.2194154

and thus roughly one order-of-magnitude higher than that of conventional ray-based tomographic images [8]. Over the past three decades, significant developments have been made in the field of seismic waveform inversion (e.g., [9]–[21]), whereas comparable efforts towards crosshole georadar waveform inversion have only been recently undertaken [6], [7], [22]–[24]. The latter algorithms are all based on conjugate-gradient optimization schemes in conjunction with finite-difference solutions of Maxwell’s equations in 2-D Cartesian coordinates. This methodology has shown much promise for reconstructing the distribution of both the dielectric permittivity and electrical conductivity under challenging real-world conditions, and has seen the successful application to field data. Whereas Ernst et al. [6], [7] and Kuroda et al. [22], [23] invert only for the vertical component of the electric field, Meles et al. [24] have built on this work to allow use of both the vertical and horizontal components. Moreover, the algorithm of Meles et al. [24] allows for a simultaneous inversion of the dielectric permittivity and the electrical conductivity, whereas the algorithms of Ernst et al. [6], [7] and Kuroda et al. [22], [23] require a cascaded, sequential approach. A crucial prerequisite for exploiting the full potential of waveform tomographic inversion schemes is having an accurate estimation of the source wavelet (e.g., [9]). This information is used to generate predicted data, which are then compared with the observed data to determine how the underlying model parameters should be updated. In general, the source wavelet is unknown a priori and, because of the complex interactions that occur between a source and its immediate surroundings, potentially variable as a function of the source location. As a result, it must be estimated in some manner from the recorded data. From a conceptual point of view, the most elegant way to estimate the source wavelet in waveform inversion is to augment the inverse problem with an additional set of model parameters (e.g., [9], [25]). In practice, however, this approach tends to be highly sensitive to the values of the starting model and hence is strongly susceptible to local minima [26]. Alternatively, the source wavelet can be estimated using a deconvolution-based approach (e.g., [18], [27]–[31]). In this regard, Ernst et al. [7] proposed a relatively simple iterative deconvolution technique for crosshole georadar, whereby the source wavelet is estimated prior to waveform inversion and can be updated during the inversion procedure if convergence proves difficult. In previous related studies, we have explored the sensitivity of this source wavelet estimation approach with regard to

0196-2892/$31.00 © 2012 IEEE

BELINA et al.: INFLUENCE OF SOURCE WAVELET VARIABILITY

realistic heterogeneity in the dielectric permittivity and electrical conductivity as well as with regard to the presence of ambient noise and unaccounted for dispersive behavior [26], [32]. The results of these extensive tests with synthetic and field data indicate that strong heterogeneity in the dielectric permittivity and/or the electrical conductivity, the presence of significant amounts of ambient noise in the data, and errors and uncertainties in the initial estimate of the source wavelet have remarkably little effect on the robustness and adequacy of the estimated source wavelet [26]. In the presence of dispersive behavior, this approach provides an estimate of an “effective” source wavelet, which partially accounts for the distortion of the source signal due to frequency-dependent velocity and attenuation along its propagation path. For weakly to moderately dispersive media, as characterized by frequency-independent quality factors Q of 40 and 20, respectively, such effective source wavelet estimates were found provide adequate waveform tomographic reconstructions which are vastly superior to their ray-based counterparts. Conversely, this approach starts to break down for strongly dispersive media, as characterized by Q < 20 [32]. In this context, it is, however, important to note that, for such media, it is very unlikely anyway that one could record crosshole georadar data with sufficiently high signal-tonoise ratios to be amenable to any kind of tomographic imaging, be it ray or waveform based. It is, however, well known that significant variability in the georadar source pulse as a function of the transmitter location may result from changes in near-field conditions. Indeed, roughness along the borehole wall and changes in the electrical properties of the borehole fluid and surrounding materials can significantly affect the coupling of an antenna to its surroundings (e.g., [33]–[35]). Hence, an important issue, which we explore in this study, is whether a single, average source pulse will allow for adequate waveform tomographic reconstructions when the true radiated source pulses vary with antenna position, or whether the source signal should be estimated for each transmitter location separately. To investigate this problem, we first invert a number of synthetic data sets exhibiting different degrees of source pulse variability with transmitter location using 1) the true source wavelets, 2) a single average estimated source wavelet, and 3) multiple source wavelets estimated for each transmitter position separately. To evaluate the results, the estimated source wavelets and waveform tomograms are compared to the true source wavelets and the corresponding subsurface model, respectively. Next, we consider a field data set collected at the Boise Hydrogeophysical Research Site near Boise, Idaho, USA, where we compare neutron porosity log measurements with porosities inferred along the boreholes from waveform-based tomographic reconstructions obtained using a single source wavelet and multiple estimated source wavelets.

II. M ETHODOLOGICAL BACKGROUND A. Forward Modeling The core of the waveform inversion algorithm considered in this study is a finite-difference solution of Maxwell’s equations

4611

in 2-D Cartesian coordinates [36], [37]. In a typical borehole georadar experiment, emitters and receivers correspond to dipole-type antennas that are aligned with the borehole axis, which in turn often coincides approximately with the z-axis of the local coordinate system (e.g., [38]). Therefore, borehole georadar surveys are concerned primarily with the vertical component of the transmitted electric field that is parallel to the transmitting and receiving antennas, such that in 2-D, the so-called transverse electric set of Maxwell’s equations is most appropriate for the purpose of modeling and inversion (e.g., [6], [7]) 
1 ∂Hy ∂Ex = − σEx (1a) − ∂t ε ∂z 
1 ∂Hy ∂Ez = − σEz (1b) ∂t ε ∂x
 1 ∂Ez ∂Ex ∂Hy = − (1c) ∂t µ ∂x ∂z where ε is the dielectric permittivity, µ is the magnetic permeability, σ is the electric conductivity, and t is time. E and H denote the electric and magnetic field components, respectively, and x, y, and z refer to the three spatial directions of a 3-D Cartesian coordinate system. In most practically relevant cases, the magnetic permeability µ can be assumed to be constant and equal to its free-space value µ0 (e.g., [39]–[43]). Indeed, this assumption is essentially universal in georadar studies and can also be expected to be fully valid for the field example from the Boise Hydrogeophysical Research Site considered in this study [45] (J.H. Bradford, personal communication, 2012, W. Barrash, personal communication, 2012). Equations (1a)–(1c) can be solved through a staggered-grid finite-difference time-domain (FDTD) approach that is secondorder accurate in both space and time [6], [7], [35], [36], [44], [46]. It can be shown that for the FDTD algorithm considered in this study, at least 10 to 20 grid points per minimum wavelength are required to limit the detrimental effects of grid dispersion (e.g., [37]). The transmitting and receiving georadar antennas are approximated as infinitesimal dipoles, which is a reasonable approximation for most modern commercial crosshole georadar systems [6], [7], [35], [46], [47]. Artificial reflections at the edges of the model space are minimized through the use of generalized perfectly matched layer absorbing boundaries [48], [49]. B. Waveform Inversion The crosshole georadar waveform inversion scheme that we consider in this study is based on the algorithm of Ernst et al. [6], [7] and corresponds to an electromagnetic adaptation of the method proposed by [9] for seismic data. The overall goal of the waveform inversion procedure is to find the spatial distribution of ε and σ that minimizes the misfit function S=

1 E syn (xtrn , xrec , t, εsyn (x), σ syn (x)) 2 z   2 −Ezobs xtrn , xrec , t, εtrue (x), σ true (x) 

(2)

4612

IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 50, NO. 11, NOVEMBER 2012

where Ezsyn and Ezobs are the synthetic and observed vertical components of the electric field at the receiver locations, and x denotes the location in the model space. The model permittivity and conductivity distributions to be estimated are given by εsyn and σ syn , whereas εtrue and σ true represent the true subsurface parameters. xtrn and xrec denote the transmitter and receiver locations. Please note that bolded characters refer to vectors and matrices. To ensure positive values for the physical properties and to aid convergence, we use logarithmically scaled versions of the model parameters, as given by
syn 
syn  ε (x) σ (x)  syn (x) = log  and σ εsyn (x) = log ε0 σ0 (3)

residual, between the observed and synthetic wavefields at the receiver locations. Essentially, Ez res back can be regarded as the “missing” part of the wavefield that needs to be generated in the next iteration of the inversion to eliminate the difference between observed and modeled wavefields. The integrands in (6a) and (6b) can be interpreted as zero-lag cross correlations (e.g., [12]). Note that the additional time derivative of the forward modeled electric field Ez syn f wd in (6a) and (6b) serves to transform the hard source used by our forward simulator into the corresponding soft current source required by the algorithm used to evaluate the gradients (e.g., [23], [24], [51]–[53]). The complete crosshole georadar waveform inversion scheme then proceeds as follows [6].

where ε0 is the free-space dielectric permittivity and σ0 is set to 1 S/m. To minimize the misfit function in (2), a conjugate-gradient scheme is employed [50]. The model parameters are updated as follows:

1) Choose an initial electrical property model,  εsyn (x)k and syn  (x)k for k = 1. In general, this is obtained through σ a ray-based inversion of the first-arrival traveltimes and maximum first-cycle amplitudes of the observed georadar data Ezobs . 2) Compute the corresponding synthetic electric field Efsyn wd (x, t)k . 3) Determine the difference between the observed and synthetic vertical electric field at each receiver location, and use these residuals as source signals and the receiver positions as source locations to perform a corresponding FDTD simulation backward in time to obtain the backres (x, t)k . propagated wavefield Eback (x)k and γ (x)k using 4) Compute the gradient directions γ ε σ (6a) and (6b). and ζ using the approach 5) Compute the step lengths ζ εk σk of Pica et al. [15]. 6) Update the model using (4a) and (4b). 7) Repeat steps 2 through 6 until convergence has been achieved. Practically, this is assessed by the rate at which the objective function S, as quantified by (2), diminishes between subsequent iterations. In the applications presented in this study, we assume to reach convergence when S changes by less than 1%.

 εsyn (x)k − ζ cε (x)k , εsyn (x)k+1 =  εk   σ

syn

 (x)k+1 = σ

syn

(x)k − ζ cσ (x)k σk 

(4a) (4b)

(x)k and c (x)k are the conjugate gradient update where c ε σ and ζ are the cordirections for the kth iteration and ζ εk σk and ζ , we follow responding step lengths. To calculate ζ  εk σk closely the approach of Pica et al. [15]. The update directions are obtained using the following formulas:   (x)k −γ  (x)k−1 γ T (x)k γ  ε ε ε c (x)k =γ  (x)k +  c (x)k−1 , ε ε ε γ T (x)k−1 γ  (x)k−1 ε  ε (x)1 = γ  (x)1 , (5a) with c ε ε   γ T (x)k γ  (x)k −γ  (x)k−1 σ σ  σ c (x) =γ (x) + c (x)k−1 , k k T  σ σ σ γ (x)k−1 γ  (x)k−1 σ  σ with c (x) = γ (x)1 (5b) 1  σ σ where the superscript T denotes the transpose operator, and and γ  are the gradient directions, for which Ernst et al. γ ε σ [6] derive the following expressions in terms of the original nonscaled electrical material parameters (x)k = ε γ ε

syn

Tmax   ∂Ez syn f wd (x, t)k (x)k dt ∂t trn 0

∂Ez res back (x, t)k × ∂t γ σˆ (x)k = − σ

syn

(6a)

Tmax   ∂Ez syn f wd (x, t)k (x)k dt ∂t trn 0

×

Ez res back (x, t)k .

(6b)

Here, Tmax denotes the maximum observation/modeling time, res Ez syn f wd is the forward propagated wavefield, and Ez back corresponds to the backpropagation in time of the difference, or

Based on the above description it is clear that this algorithm, and indeed backpropagation-based waveform inversion algorithms in general, are not amenable to any standard regularization approaches, such as, for example, damping and smoothing, which are commonly used for linearized-least-squares-type inversions [54], [55]. The waveform inversion algorithm used in this study requires three FDTD simulations for each iteration as well as the storage of the corresponding vertical electrical field components for every cell of the model space. While no obvious shortcuts are possible with regard to the computational load of the forward and backward FDTD simulations, the memory requirements can be dramatically reduced by averaging several modeling cells into a single inversion cell. In general, we therefore average either four or nine FDTD cells into one inversion cell before calculating the conjugate gradient. After evaluating the conjugate gradient for each cell, these larger inversion cells are again subdivided into the corresponding smaller FDTD cells for the subsequent calculations outlined above. Our experience

BELINA et al.: INFLUENCE OF SOURCE WAVELET VARIABILITY

4613

indicates that this procedure allows for reducing the memory requirements by roughly one order of magnitude [6]–[8]. The reason why no obvious shortcuts are possible with regard to the computational load of the forward and backward FDTD simulations is that the spatial discretization of the model required to control the accumulation of numerical errors caused by grid dispersion is much finer than the spatial resolution of the considered wavefield. C. Source Wavelet Estimation In the absence of pronounced nonlinear phenomena, such as multiple scattering or strong dispersion, a crosshole georadar measurement can be effectively approximated as the convolution of the Earth’s impulse response with the georadar source wavelet. As a result, deconvolving this impulse response from the measured data should yield an estimate of the source wavelet, and performing the deconvolution using a number of different measurements should increase the robustness of that estimate. The true impulse response of the earth is a priori unknown, but Ernst et al. [7] showed that a reasonable approximation can be obtained by forward modeling through an approximate subsurface property model, such as that obtained using standard ray-based tomography. In other words, we can first calculate synthetic georadar data through our best guess of the subsurface model using a known “starting” source wavelet and then deconvolve this known wavelet from the synthetic data in order to get an estimate of the Earth’s impulse response. In the frequency domain, this is expressed as follows for a given transmitterreceiver configuration

−1 ˜ syn (xtrn , xrec , ω) S˜ini (ω) ˜ (xtrn , xrec , ω) = E (7) M ˜ the approximate imwhere ω dentoes angular frequency, M syn ˜ the synthetic georadar data, pulse response of the earth, E S˜ini the known initial wavelet, and xtrn and xrec the transmitter and receiver locations, respectively. In our approach, S˜ini is obtained by averaging the first cycle of all waveforms recorded when the corresponding transmitter and receiver antennas are horizontally aligned [56]. However, any other realistic estimate of the source wavelet with comparable amplitude and phase spectra could be used [7], [26]. The next step is to estimate the ˜ from unknown field source wavelet S˜est by deconvolving M obs ˜ the observed field georadar data E ˜ obs (xtrn , xrec , ω) S˜est (xtrn , xrec , ω) = E

−1 ˜ (xtrn , xrec , ω) × M .

(8)

Please note that, in principle, the Earth’s impulse response ˜ and the corresponding estimation of the source wavelet S˜est M need to be evaluated separately for each transmitter-receiver configuration considered [(7) and (8)]. In practice, however, a single average source wavelet is commonly employed [7]. Given the overdetermined nature of the problem described by (8), as characterized by far more data than unknowns, the globally best-fitting single source wavelet can then be estimated

through a linearized least-squares optimization procedure. Arguably, the key objective and the major potential merit of this approach are to reduce the influence of noise and thus to increase robustness of the estimated average source wavelet [7]. Once the source wavelet has been estimated, the previously described crosshole georadar waveform inversion algorithm can be run until convergence is achieved. If eventual convergence seems unlikely after a certain number of iterations, in our experience approximately ten, then an iterative source wavelet refinement may be required. This involves using the current subsurface model to obtain an improved estimate the Earth’s impulse response, which in turn is used to obtain an enhanced estimate of the source wavelet. Through extensive testing of this algorithm in a variety of realistic scenarios, we have found that one iteration of the source wavelet estimation procedure described by (7) and (8) is in most cases sufficient [26]. Only very complex subsurface structures and/or inadequate starting models tend to require an iterative refinement of the estimated source wavelet. Evidently, the source wavelet estimation approach described above is in principle only valid when the variability in the source pulse between different transmitter locations is benign, which is indeed commonly not the case. In the following, we therefore explore the limits of this classical approach as well as the potential benefits to waveform inversion that can be gained from estimating a separate source wavelet for each transmitter position.

III. R ESULTS A. Synthetic Tests To investigate whether a single estimated source wavelet can adequately reproduce the pulses radiated during a crosshole georadar survey in the case where the true source pulses vary with antenna location, we apply the wavelet estimation and waveform inversion algorithms outlined above to three synthetic data sets. All of these data sets are based on the same underlying model of subsurface electrical properties, but each has a different degree of variability in the true source signal as a function of the transmitter location. That is, the data sets simulate different amounts of variability in the wavelet radiated by the transmitter antenna with regard to its position in the borehole. Such source wavelet variability is indeed likely to occur in practice because of changing near-field conditions [33], [47]. We carry out the waveform inversion of the three data sets using 1) the true source wavelets, 2) a single, or average, estimated source wavelet based on the data from all transmitter gathers together, and 3) multiple estimated source wavelets based on the data from each transmitter gather separately. We then compare the source wavelet estimates and resulting waveform tomograms in order to assess the benefits and drawbacks of using a single versus multiple estimated source wavelets in each case. Please note that, for simplicity, we only consider heterogeneity in the dielectric permittivity in this study. The magnetic permeability and electrical conductivity are assumed to be known and constant, and equal to the permeability of free space µ0 and 2 mS/m, respectively. The assumption of a

4614

IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 50, NO. 11, NOVEMBER 2012

Fig. 1. (a) Dielectric permittivity model used to generate the synthetic data in this study. (b) Corresponding ray-based traveltime inversion tomogram that was used as an initial model for the source wavelet estimation and waveform inversions.

nonmagnetic subsurface is, as outlined before, reasonable and adequate for the overwhelming majority of practically relevant situations [39]–[41], [43]. The effects of heterogeneity in the electrical conductivity were explored in a previous related study [26], which demonstrated that this only affects the amplitudes, but not the phase, of the received signal and hence can be ignored when working with normalized or relative amplitudes as we do here. Fig. 1(a) shows the dielectric permittivity model upon which the three synthetic data sets are based. The model exhibits a realistic, highly heterogeneous layered structure that is characterized by band-limited scale invariance. To generate the model, we used a spectral geostatistical simulation technique [8], [57] based on a von Kármán autocovariance function having horizontal and vertical correlation lengths of 2 m and 0.2 m, respectively, and a v-value of 0.05. These values allow us to emulate the seemingly universal and ubiquitous “flicker noise” character of petrophysical data [58]–[60]. The permittivity fluctuations have a standard deviation of ∼15% of the mean value of 10ε0 , which corresponds to a standard deviation of ∼7.5% of the mean value of 0.09 m/ns for the velocity fluctuations. The type of heterogeneity seen in Fig. 1(a) can be regarded as a realistic first-order abstraction of many surficial alluvial environments (e.g., [61]–[63]), and hence serves as a pertinent and challenging test case for the waveform inversion algorithm. To generate synthetic crosshole georadar data corresponding to Fig. 1(a) in each test case, we used the same FDTD solution of Maxwell’s equations in 2-D Cartesian coordinates that is employed in the waveform inversion procedure. We considered 41 transmitter antenna locations and 41 receiver antenna locations spaced every 0.5 m down the left and right edges of the 10-m-by-20-m modeling domain. Each of the FDTD simulations was run at less than 70% of the stability limit as defined by the Courant criterion. The spatial discretization interval was set to be 0.02 m, which ensured at least 15 grid points per minimum wavelength and thus allowed us minimize numerical artifacts related to grid dispersion in the resulting data [35], [37]. Please

note that the geometrical setup and dimensions of this model as well as the number of transmitters and receivers can, for all practically relevant intents and purposes, be regarded as being representative of the vast majority of experimental setups used for crosshole georadar imaging. This far-reaching and effective standardization of the experimental design for crosshole georadar and seismic data acquisition has been largely guided by a combination of methodological considerations of optimal resolution as well as by practical considerations of operational and computational efficiency [16], [17], [43]. In this context, it is also important to note that the experimental setup for the acquisition of crosshole georadar is bistatic and thus consists of one transmitting and one receiving antenna. As a result, a transmitter gather is formed by keeping the transmitting antenna at a given location in one borehole while covering all of the recording positions by moving the receiver antenna along the other borehole. Fig. 1(b) shows the tomographic image that we obtained from the ray-based inversion of the first-arrival traveltimes corresponding to Fig. 1(a). The latter were obtained by applying a semi-automatic picking procedure to the synthetic data. The tomogram in Fig. 1(b) served as the initial model for both our source wavelet estimation and waveform inversion procedures. Note the inherent smoothness and associated limitations in resolution of the ray tomogram compared to the true model in Fig. 1(a). It is these limitations that we aim to overcome through the use of waveform inversion. We first consider the case of weak variability in the true source wavelet as a function of transmitter position. Fig. 2(a) shows the 41 source wavelets that were randomly distributed over the transmitter locations and used to generate the synthetic data for this example. They are Ricker wavelets having different amplitudes and slightly different frequency content as characterized by a standard deviation of 4.5% with regard to the mean dominant frequency of 110 MHz. In Fig. 2(b), we plot the source wavelets that were estimated from the corresponding synthetic data. The single average estimated wavelet obtained using the data from all transmitter gathers is shown in red, whereas the multiple estimated source wavelets considering each transmitter gather separately are shown in blue. Overall, there is good agreement between the true and multiple estimated source wavelets, and the single average wavelet is seen to provide a good representative pulse. Note, however, that all of the estimated source wavelets are of slightly lower amplitude than the true wavelets. This is because the ray-based tomogram that we used to approximate the Earth’s impulse response is too smooth to completely account for all of the effects of scattering [Fig. 1(b)]. That is, there exists scattering attenuation in the crosshole data that is not predicted when using the approximate Earth impulse response based on the ray tomogram, and thus the estimated wavelets must be reduced in amplitude to match the true data. Fig. 3(a)–(c) shows the waveform tomographic reconstructions of dielectric permittivity that were obtained from the synthetic data using the true source wavelets, the single estimated wavelet based on all transmitter gathers, and the multiple estimated wavelets based on each transmitter gather separately. Fig. 3(d)–(f) shows the corresponding residual tomograms,

BELINA et al.: INFLUENCE OF SOURCE WAVELET VARIABILITY

4615

Fig. 2. Weak variability in the true source pulse with transmitter location. (a) True source wavelets randomly assigned to the 41 transmitter locations. (b) Source wavelet estimated using the data from all transmitter gathers together (red) versus multiple source wavelets estimated using the data from each transmitter gather separately (blue). Amplitudes are normalized by the overall maximum amplitude of the true and estimated wavelets.

Fig. 3. Weak variability in the true source pulse with transmitter location. Top row: Waveform permittivity tomograms obtained using (a) the true source wavelets, (b) a single estimated source wavelet based on all transmitter gathers, and (c) multiple estimated source wavelets based on each transmitter gather separately. Bottom row: Percentage deviation between the tomograms in (a)–(c) and the true dielectric permittivity model in Fig. 1(a).

which we define as the percentage deviation of the waveform tomograms from the true model in Fig. 1(a). Looking at these results, we see that all three waveform inversions allow for a greatly superior reconstruction of the true permittivity model

compared to the ray-based traveltime inversion result shown in Fig. 1(b). This clearly illustrates the significant improvements in resolution that can be obtained through waveform inversion and provides further motivation for work in this domain.

4616

IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 50, NO. 11, NOVEMBER 2012

Fig. 4. Moderate variability in the true source pulse with transmitter location. (a) True source wavelets randomly assigned to the 41 transmitter locations. (b) Source wavelet estimated using the data from all transmitter gathers together (red) versus multiple source wavelets estimated using the data from each transmitter gather separately (blue). Amplitudes are normalized by the overall maximum amplitude of the true and estimated wavelets.

Looking at Fig. 3(b) and (c), we see that a number of small artifacts appear in the waveform reconstructions near the left model edge where the transmitter borehole is located. These artifacts are characteristic of the inversion scheme attempting to compensate for inaccuracies in the forward model [7], here caused by differences between the true and estimated wavelets. That is, the algorithm has placed heterogeneity near the transmitter locations in order to compensate for errors in the assumed source wavelets. The artifacts in Fig. 3 are visibly worse for the case of a single, or average, estimated wavelet [Fig. 3(b)] compared to using multiple wavelets estimated for each transmitter location separately [Fig. 3(c)], because the former cannot address the fact that the true source pulse varied with the transmitter antenna location. Nevertheless, the artifacts in Fig. 3(b) are minor, and we can conclude that, in the case of weak variability of source wavelet, a single estimated wavelet is adequate to obtain an excellent waveform-based tomographic reconstruction. Next, we consider the case of moderate variability in the true source wavelet with transmitter position and perform the same analysis as above. Fig. 4(a) shows the true source pulses considered in this example, which again were distributed randomly among the 41 transmitter locations and used for the synthetic modeling. We obtained these source wavelets by applying our wavelet estimation procedure to field data collected at the Boise Hydrogeophysical Research Site near Boise, Idaho, USA [7], [45]. The field data are considered in further detail in Section III-B. The wavelet amplitudes in Fig. 4(a) can be seen to vary significantly, and the differences in frequency content are quantified by a standard deviation of 12% with regard to the average center frequency of 89 MHz. The somewhat “ringy” character of the wavelets is related to the fact the boreholes were filled with water [7], [35]. In Fig. 4(b), we show the single wavelet estimated from the corresponding synthetic data based on all transmitter gathers in red, and the multiple estimated source wavelets considering each transmitter gather separately in blue. Again, the overall agreement between the multiple estimated source wavelets compared to the true source wavelets

is good, and the single estimated wavelet can be considered to be adequately representative. Fig. 5 shows the waveform tomograms obtained using the true source wavelets, the single estimated source wavelet, and the multiple estimated source wavelets, as well as the percent difference between these tomograms and the true permittivity model in Fig. 1(a). Again, all three tomographic reconstructions clearly represent a significant improvement in resolution compared to the ray-based inversion result in Fig. 1(b), and both a single estimated average source wavelet and multiple estimated wavelets for each transmitter gather are able to provide reconstructions that are very close to the best possible image obtained when the true source wavelets are known. Artifacts again appear along the left model edge in Fig. 5(b) and (c) as a result of the inversion scheme trying to compensate for errors in the estimated source wavelets. These artifacts are slightly greater in Fig. 5(b) than in Fig. 3(b) because of the greater amount of source wavelet variability in this example, but they can still be considered as rather benign. For our third test, we consider the case of extremely strong variability in the true source wavelet as a function of transmitter position. To do this, we generate synthetic data using highly variable wavelets characterized by phase differences of up to 180◦ and significantly different maximum amplitudes and center frequencies ranging between 79 and 116 MHz [Fig. 6(a)]. This set of source wavelets represents a mix of the wavelets used in the two previous tests. The corresponding multiple estimated source wavelets based on each transmitter gather again show good agreement with the true source wavelets [Fig. 6(b)]. However, the single average source wavelet estimated using the data from all transmitter gathers now clearly fails to be representative of any of the true wavelets in Fig. 6(a). Consequently, the waveform tomographic reconstruction obtained using this wavelet [Fig. 7(b)] is inferior to that obtained using the true source wavelets [Fig. 7(a)] and that obtained using the estimated wavelets based on each transmitter gather separately [Fig. 7(c)]. In this case, using multiple estimated source wavelets for each transmitter location is critical for obtaining

BELINA et al.: INFLUENCE OF SOURCE WAVELET VARIABILITY

4617

Fig. 5. Moderate variability in the true source pulse with transmitter location. Top row: Waveform permittivity tomograms obtained using (a) the true source wavelets, (b) a single estimated source wavelet based on all transmitter gathers, and (c) multiple estimated source wavelets based on each transmitter gather separately. Bottom row: Percentage deviation between the tomograms in (a)–(c) and the true dielectric permittivity model in Fig. 1(a).

Fig. 6. Strong variability in the true source pulse with transmitter location. (a) True source wavelets randomly assigned to the 41 transmitter locations. (b) Source wavelet estimated using the data from all transmitter gathers together (red) versus multiple source wavelets estimated using the data from each transmitter gather separately (blue). Amplitudes are normalized by the overall maximum amplitude of the true and estimated wavelets.

a high-quality waveform tomographic image. Indeed, Fig. 7(c) is highly similar to Fig. 7(a) and matches very well the true dielectric permittivity model in Fig. 1(a).

For our final test, we assume the true source pulse to be constant for each transmitter location. Fig. 8 then compares this pulse with the corresponding single estimated source wavelet

4618

IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 50, NO. 11, NOVEMBER 2012

Fig. 7. Strong variability in the true source pulse with transmitter location. Top row: Waveform permittivity tomograms obtained using (a) the true source wavelets, (b) a single estimated source wavelet based on all transmitter gathers, and (c) multiple estimated source wavelets based on each transmitter gather separately. Bottom row: Percentage deviation between the tomograms in (a)–(c) and the true dielectric permittivity model in Fig. 1(a).

Fig. 8. Comparison of the constant true source pulse used at all of the transmitter locations (black), single source wavelet estimated using the data from all transmitter gathers together (red), and multiple source wavelets estimated using the data from each transmitter gather separately (blue). Amplitudes are normalized by the overall maximum amplitude of the true and estimated wavelets.

based on all transmitter gathers and the multiple estimated source wavelets based on each gather separately. The corresponding tomographic reconstructions and residual tomograms

are shown in Fig. 9. Notice in Fig. 8 that the single estimated wavelet based on all transmitter gathers matches the true source pulse quite well, whereas the multiple estimated wavelets differ rather notably. This is due to the fact that far fewer traces are available to infer a separate wavelet for each transmitter position, which reduces the robustness of the deconvolution-based estimate when the source wavelet does not actually vary with location [7]. As a result, in this case, the waveform tomographic reconstruction obtained using the multiple estimated wavelets [Fig. 9(c)] contains more artifacts than the reconstruction based on a single estimated wavelet [Fig. 9(b)]. This clearly illustrates the caveats and potential pitfalls of the procedure explored in this study. Specifically, our results indicate that estimating multiple source wavelets based on fewer data is only truly warranted when there is convincing evidence for pronounced variability of the source wavelet as a function of the transmitter position. This is arguably a key finding of this sensitivity study and an important result per se. Please note that that the persistent blue feature in all of the residual plots [Figs. 3(d)–(f), 5(d)–(f), 7(d)–(f), and 9(d)–(f)] at around 1-m depth is related to the fact that, in this region, we have very poor constraints on the subsurface velocity structure because of poor coverage by the propagating georadar energy. In other words, the recorded georadar data are not very sensitive to features in this region. The particular mentioned feature

BELINA et al.: INFLUENCE OF SOURCE WAVELET VARIABILITY

4619

Fig. 9. No variability in the true source pulse with transmitter location. Top row: Waveform permittivity tomograms obtained using (a) the true source wavelets, (b) a single estimated source wavelet based on all transmitter gathers, and (c) multiple estimated source wavelets based on each transmitter gather separately. Bottom row: Percentage deviation between the tomograms in (a)–(c) and the true dielectric permittivity model in Fig. 1(a).

represents a small low-velocity zone in the true subsurface model [Fig. 1(a)] that is not well resolved in the various reconstructions. The model cell values in this region tend to remain at those of the ray-based starting model [Fig. 1(b)]. The poor resolution due to limited coverage is likely exacerbated by the higher velocity zone located above, through which most of the georadar energy will propagate. Conversely, the previously discussed artifacts in the vicinity of the transmitter locations in Figs. 3, 5, 7, and 9 are inherent to the backpropagation of the residual wavefield used for the model adjustment and the nonamenability of this inversion approach to standard regularization procedures [55]. In an attempt to quantify the findings of this sensitivity study, we have evaluated the root-mean-square (RMS) error of the tomographically reconstructed distributions of the dielectric permittivity [Figs. 3(a)–(c), 5(a)–(c), 7(a)–(c), and 9(a)–(c)] with regard to the underlying true model [Fig. 1(a)]. The corresponding results are shown in Table I. Unfortunately, these RMS values are strongly biased by the aforementioned artifacts in the vicinity of the transmitter locations. This in turn results in a generally increasing RMS with decreasing variability of

TABLE I RMS D IFFERENCE B ETWEEN THE U NDERLYING T RUE M ODEL OF THE D IELECTRIC P ERMITTIVITY S TRUCTURE [F IG. 1(a)] AND W ITH R EGARD TO THE C ORRESPONDING T OMOGRAPHIC R ECONSTRUCTIONS FOR W EAK [F IG. 3(a)–(c)], M ODERATE [F IG. 5(a)–(c)], S TRONG [F IG. 7(a)–(c)], AND NO VARIABILITY OF THE S OURCE WAVELET [F IG . 9(a)–(c)]

the source wavelet, which is clearly counterintuitive to the qualitative visual perception of the corresponding tomographic reconstructions [Figs. 3(a)–(c), 5(a)–(c), 7(a)–(c), and 9(a)–(c)] and their corresponding percentage deviations [Figs. 3(d)–(f), 5(d)–(f), 7(d)–(f), and 9(d)–(f)] with regard to the underlying model [Fig. 1(a)].

4620

IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 50, NO. 11, NOVEMBER 2012

Fig. 10. Schematic diagram of the crosshole georadar experiment conducted at the Boise Hydrogeophysical Research Site. “Trn” and “Rec” denote the transmitter and receiver, respectively; only the first and last antenna locations are shown. The unit boundaries (dashed yellow lines) were determined based on borehole porosity logs. Adapted from Ernst et al. [7].

B. Field Data Example We now investigate the use of single versus multiple estimated source wavelets when performing waveform inversion of a real crosshole georadar data set collected at the Boise Hydrogeophysical Research Site near Boise, Idaho, USA. This site consists of a dense network of boreholes that have been used over the past ∼15 years to conduct a wide range of geological, hydrological, geochemical, and geophysical experiments (e.g., [45], [63]–[67]). The geology at the site is characterized by an approximately 20-m-thick alluvial layer consisting of braidedriver gravels and sands with minimal fractions of silt and clay underlain by a layer of red clay, which serves as the basal aquitard for the corresponding unconfined surficial aquifer. At the time the georadar data were acquired, the groundwater table was at a depth of 2.96 m below measuring point (bmp) or 2.36 m below the land surface. A schematic diagram of the crosshole experiment is shown in Fig. 10. The georadar data were acquired between boreholes C5 and C6, which have a diameter of approximately 10 cm. Both boreholes are slightly deviated and are roughly 8 m apart. The data acquisition was carried out from ∼4 to ∼19.5 m depth bmp, which ensured that the antenna elements were completely below the groundwater table. Over this depth range, the boreholes penetrate three different pebble- and cobble-dominated zones with avearge

Fig. 11. (a) Comparison of the single estimated source wavelet based on all transmitter gathers (red) with the multiple estimated source wavelets based on each transmitter gather separately (blue). (b) The same as (a), except that the estimated source wavelets are based on receiver gathers.

porosities of approximately 21%, 26%, and 23%. This zonation, which is largely based on the interpretation of crosshole georadar data [7], [45], is, with minor differences in detail, consistent with evidence from high-resolution neutron porosity logs, [64], capacitive conductivity logs [67], and analyses of core material [65]. A total of 77 transmitter locations and 40 receiver locations were used for the crosshole georadar survey, with spacings of 0.2 m and 0.4 m down the C6 and C5 boreholes, respectively. Although a 250-MHz RAMAC borehole georadar system was employed, the recorded data have a dominant frequency of approximately 80 MHz because the borehole fluid and surrounding water-saturated sediments make the antenna “electrically longer” than when it is employed in air (e.g., [35], [47]). Processing of the georadar data involved removal of the DC component and application of a 0–250 MHz zero-phase lowpass filter. For the waveform inversions, we only consider data where the absolute transmitter-receiver angle was less than 45◦ from the horizontal in order to avoid potential problems with high-angle measurements [7], [34]. Anisotropy in the georadar velocity was previously determined to be negligible at the site and the probed subsurface region was found to be largely nondispersive [7]. As in our synthetic examples, we only invert for the dielectric permittivity and assume that the magnetic

BELINA et al.: INFLUENCE OF SOURCE WAVELET VARIABILITY

4621

Fig. 12. (a) Ray-based traveltime inversion tomogram that was used as an initial dielectric permittivity model for the source wavelet estimation and waveform inversions for the field data from Boise Hydrogeophysical Research Site. This is compared with the waveform permittivity tomograms obtained using: (b) a single estimated source wavelet based on all transmitter gathers; (c) multiple estimated source wavelets based on each transmitter gather separately; (d) a single estimated source wavelet based on all receiver gathers; (e) multiple estimated source wavelets based on each receiver gather separately.

permeability and electrical conductivity are known and equal to constant values of µ0 and 2 mS/m, respectively. Evidence from extensive earlier work clearly suggests that the almost universal assumption of an essentially nonmagnetic subsurface (e.g., [39]–[43]) is justified for this site [45], [J.H. Bradford, personal communication, 2012]. Indeed, magnetic susceptibility values in SI from corresponding logs run in boreholes C5 and C6 range approximately between 0.001 and 0.005 (W. Barrash, personal communication, 2012). This in turn implies that the maximum error of the assumption µ = µ0 is of the order of 0.5% or less, which for all practical intents and purposes is negligible. Moreover, a previous sensitivity analysis has shown that, for the physical model considered in this study, variations in the electrical conductivity only result in corresponding variations

of the absolute amplitude of the estimated source wavelet, but do not affect its phase [26]. Given that we are working with normalized source wavelets, the assumption of a constant conductivity, corresponding approximately to the average value of the probed volume [45], is therefore adequate. As in our synthetic examples, we invert the field data considering a single average source wavelet estimated using the data from all transmitter gathers together, and multiple estimated source wavelets based on the data from each transmitter gather separately. In addition, we also consider reciprocity [68] in the field example in order to investigate the influence of variations in the recorded georadar signal due to differences in receiver antenna coupling with position. In doing so, we invert the data using a single source wavelet estimated from all of the receiver

4622

IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 50, NO. 11, NOVEMBER 2012

gathers together and multiple estimated source wavelets based on each receiver gather separately. For the latter two cases, the receiver positions were used as the transmitter positions during the inversion procedure. For each inversion, we compare the resulting waveform tomograms with neutron porosity logs acquired along boreholes C5 and C6 [45] by converting the dielectric permittivity values obtained along the boreholes to porosity using the complex refractive index method formula (CRIM) [69]. Fig. 11(a) and (b) show the different source wavelets that were estimated from the field data based on the transmitter and receiver gathers, respectively. We see that estimating a wavelet based on each receiver gather separately leads to more variability in the estimated wavelets than when each transmitter gather is considered separately. Also, note that the estimated source wavelet based on all transmitter gathers together is similar to the estimated source wavelet based on all receiver gathers together and can be considered to be an adequate representative wavelet in both cases. Fig. 12(a) shows the ray-based traveltime tomogram that was used as an initial model for the waveform inversion and source wavelet estimation procedures. Fig. 12(b)–(e) then shows the waveform tomograms that were obtained in each case. Notice that the waveform tomogram obtained using the estimated wavelet based on all transmitter gathers [Fig. 12(b)] is almost identical to its reciprocal counterpart based on all receiver gathers [Fig. 12(d)]. In contrast, the waveform tomogram obtained using the estimated source wavelets based on each transmitter gather separately [Fig. 12(c)] shows slight differences mainly along the transmitter borehole C6, where some high permittivity zones, such as that at ∼15.5 m depth bmp, are underestimated compared to Fig. 12(b) and (d). The same holds true, possibly in an even more pronounced manner, for the waveform tomogram obtained using the estimated source wavelets based on each receiver gather separately [Fig. 12(e)]. In this case, the most prominent example is the high permittivity zone at approximately 6 m depth bmp, which is clearly underestimated compared to Fig. 12(b)–(d). The latter observations are corroborated by Fig. 13, where we compare porosities obtained from the permittivity values along the boreholes in each waveform tomogram in Fig. 12 with the corresponding neutron porosity log data in boreholes C5 and C6. Our explanation for the systematic underestimation of highpermittivity features along the boreholes when using multiple estimated source wavelets is the following: If a different source wavelet is estimated for each transmitter or receiver gather, then strong fluctuations in the dielectric permittivity in the immediate vicinity of the corresponding transmitter or receiver locations tend to be compensated through the source wavelet estimation procedure since the traces of the corresponding transmitter or receiver gather used to estimate the wavelet are primarily affected by this local heterogeneity. In other words, just as we saw how errors in the source wavelet can be compensated in the waveform tomograms by near-borehole heterogeneity, the effects of near-borehole heterogeneity can be compensated through the multiple estimated source wavelets. Conversely, when estimating a single average source wavelet based on all transmitter or receiver gathers together, the corresponding effects are limited to relatively few traces and tend to

Fig. 13. Comparison of neutron porosity logs along boreholes C5 and C6 (black) with porosities along these boreholes derived from the waveform permittivity tomograms shown in Fig. 5(b) (light blue), Fig. 5(c) (green), Fig. 5(d) (red), and Fig. 5(e) (blue).

be averaged out. In this context, it is important to emphasize that the variability observed in the estimated source wavelets based on each transmitter or receiver gather separately may indeed not reflect actual variability in the true source wavelets. Rather these variations might represent stochastic fluctuations resulting from a smaller number of traces used for the estimation, which we explored in the last synthetic test case (Fig. 8), or strong near-borehole heterogeneity unaccounted for in the starting model. The fact that the waveform tomograms obtained using a single estimated wavelet are more consistent with the porosity-log data indicates that the latter may indeed occur in the considered field data set. The stronger degree of variability seen when we estimate a source wavelet for each receiver gather separately [Fig. 11(b)], as opposed to each transmitter gather separately [Fig. 11(a)], probably results from stronger fluctuations in the dielectric permittivity close to borehole C5, which we indeed see in the corresponding porosity log data. C. Computational Aspects The source wavelet estimation and waveform inversion procedures described and tested above are inherently parallel in nature and hence lend themselves to be carried out on a Beowulf-type high-performance computing cluster. The nodes of the cluster we used for this purpose consist of 4 AMD Opteron 6136 processors with 12 cores clocked at 2.2 GHz. The inversion algorithm was parallelized by associating each transmitter gather with a different core. An additional “master” core then coordinates these “slave” cores. The synthetic examples considered above contained 41 transmitter locations, which corresponds to 42 cores (1 master and 41 slaves). The grid spacings used for the forward and inverse problems are 0.02 and 0.06 m, corresponding to 549 × 1050 and 183 × 350 cells, respectively. Convergence was

BELINA et al.: INFLUENCE OF SOURCE WAVELET VARIABILITY

reached after 40 iterations, each of which took approximately 25–30 min. For the field example, the 40 and 77 transmitter gathers were allocated to 41 and 78 cores, respectively. Here, the grid spacings for the forward and inverse problems were 0.02 and 0.14 m, corresponding to 504 × 833 and 72 × 119 cells, respectively. Approximately 45 iterations were needed to reach convergence with each iteration taking approximately 6–7 min. The convergence time and number of iterations evidently depend on the adequacy of the starting model. In this context, it would be interesting to explore the relative benefits of standard ray-based starting models versus more advanced methods, such as, for example, tomographic inversions based on the Born approximation [70], [71]. Additionally, it might be interesting to consider a deconvolution-based Born-type approach for the simultaneous estimation of the initial source wavelet and starting model of spatial distribution of the electrical parameters. IV. C ONCLUSION While the presence of spatial and/or temporal fluctuations in the source wavelet in crosshole georadar surveying is widely known and acknowledged, the importance of these fluctuations for the newly emerging practice of waveform inversion of such data has been as-of-yet largely unexplored. In particular, it was not clear whether the common practice of using a single estimated source wavelet averaged over all transmitter positions is valid when source wavelet variability exists, and in this case, whether the quality and resolution of the waveform tomographic reconstructions could be improved by estimating a separate source wavelet for each transmitter position. To address these questions, we have tested and compared these two different source wavelet estimation strategies through a series of synthetic tests and application to a pertinent field data set. The results of our synthetic tests indicated that the estimation of separate source wavelets for each transmitter position is only beneficial and warranted in the presence of very-strong-toextreme variability in the true source wavelets. In the presence of moderate-to-strong source wavelet variability, the standard approach of estimating a single average wavelet is entirely adequate. Indeed, in the latter case, which is likely to be representative of most situations of practical importance, the estimation and use of multiple source wavelets may even result in tomographic images of inferior quality compared to those based on a single wavelet assumption. This observation can be attributed to the fact that multiple source wavelets estimated for each transmitter position separately are subject to significantly greater stochastic uncertainty in the course of the estimation procedure. Our synthetic results and their interpretation were found to be fully consistent with findings based on the field data set. Our tests on the field data also indicated that the use of multiple source wavelets is likely to be negatively affected by the fact that the wavelet estimation procedure may be dominated by pronounced heterogeneity in the immediate vicinity of the considered transmitter location. Finally, it is important to note that due to the far-reaching analogies between georadar and visco-acoustic wave phenomena, the results and conclusions of this study can be expected to be also fully valid for the waveform-based inversion of crosshole seismic data.

4623

ACKNOWLEDGMENT The authors thank H. Maurer and the support staff at ETH Zurich for facilitating work at the local high-performance computing center, the Center for the Geophysical Investigation of Shallow Subsurface at Boise State University for granting the permission to work with the crosshole georadar and neutron porosity log data collected at the Boise Hydrogeophysical Research Site, and two anonymous referees for thoughtful and constructive comments, which helped to improve the quality of this manuscript. The authors are also particularly grateful to W. Barrash of Boise State University for his careful proofreading of the revised version of this manuscript and for providing access to unpublished magnetic susceptibility logs for boreholes C5 and C6 at the Boise Hydrogeophysical Research Site. R EFERENCES [1] F. Soldovieri, A. Brancaccio, G. Prisco, and G. Leone, “A kirchhoff-based shape reconstruction algorithm for the multimonostatic configuration: The realistic case of buried pipes,” IEEE Trans. Geosci. Remote Sens., vol. 46, no. 10, pp. 3031–3038, Oct. 2008. [2] E. Pettinelli, A. Di Matteo, E. Mattei, L. Crocco, F. Soldovieri, J. D. Redman, and A. P. Annan, “GPR response from buried pipes: Measurement on field site and tomographic reconstructions,” IEEE Trans. Geosci. Remote Sens., vol. 47, no. 8, pp. 2639–2645, Aug. 2009. [3] L. Lo Monte, D. Erricolo, F. Soldovieri, and M. C. Wicks, “Radio frequency tomography for tunnel detection,” IEEE Trans. Geosci. Remote Sens., vol. 48, no. 3, pp. 1128–1137, Mar. 2010. [4] J. Minet, S. Lambot, E. C. Slob, and M. Vanclooster, “Soil surface water content estimation by full-waveform GPR signal inversion in the presence of thin layers,” IEEE Trans. Geosci. Remote Sens., vol. 48, no. 3, pp. 1138–1150, Mar. 2010. [5] F. Soldovieri, O. Lopera, and S. Lambot, “Combination of advanced inversion techniques for an accurate target localization via GPR for demining applications,” IEEE Trans. Geosci. Remote Sens., vol. 49, no. 1, pp. 451– 461, Jan. 2011. [6] J. R. Ernst, H. Maurer, A. G. Green, and K. Holliger, “Full-waveform inversion of crosshole radar databased on 2-D finite-difference time-domain solutions of Maxwell’s equations,” IEEE Trans. Geosci. Remote Sens., vol. 45, no. 9, pp. 2807–2828, Sep. 2007. [7] J. R. Ernst, H. Maurer, A. G. Green, and K. Holliger, “Application of a new 2D time-domain full-waveform inversion scheme to crosshole radar data,” Geophysics, vol. 72, no. 5, pp. J53–J64, Aug. 2007. [8] F. A. Belina, J. R. Ernst, and K. Holliger, “Inversion of crosshole seismic data in heterogeneous environments: Comparison of waveform and raybased approaches,” J. Appl. Geophys., vol. 68, no. 1, pp. 85–94, May 2009. [9] A. Tarantola, “Inversion of seismic reflection data in the acoustic approximation,” Geophysics, vol. 49, no. 8, pp. 1259–1266, Aug. 1984. [10] A. Tarantola, The Seismic Reflection Inverse Problem, in Inverse Problems of Acoustic and Elastic Waves. Philadelphia, PA: SIAM, 1984. [11] A. Tarantola, “A strategy for nonlinear elastic inversion of seismicreflection data,” Geophysics, vol. 51, no. 10, pp. 1893–1903, Oct. 1986. [12] A. Tarantola, Inverse Problem Theory and Methods for Model Parameter Estimation. Philadelphia, PA: SIAM, Dec. 2005. [13] O. Gauthier, J. Virieux, and A. Tarantola, “Two-dimensional nonlinear inversion of seismic wave-froms—Numerical results,” Geophysics, vol. 51, no. 7, pp. 1387–1403, Jul. 1986. [14] P. Mora, “Nonlinear two-dimensional elastic inversion of multi-offset seismic data,” Geophysics, vol. 52, no. 9, pp. 1211–1228, Sep. 1987. [15] A. Pica, J. P. Diet, and A. Tarantola, “Nonlinear inversion of seismicreflection data in a laterally invariant medium,” Geophysics, vol. 55, no. 3, pp. 284–292, Mar. 1990. [16] R. G. Pratt and M. H. Worthington, “Inverse-theory applied to multisource cross-hole tomography, Part 1: Accoustic wave-equation method,” Geophys. Prospect., vol. 38, no. 3, pp. 287–310, Apr. 1990. [17] R. G. Pratt, “Inverse-theory applied to multisource cross-hole tomography, Part 2: Elastic wave-equation method,” Geophys. Prospect., vol. 38, no. 3, pp. 311–329, Apr. 1990. [18] R. G. Pratt, “Seismic waveform inversion in the frequency domain, Part 1: Theory and verification in a physical scale model,” Geophysics, vol. 64, no. 3, pp. 888–901, May/Jun. 1999. [19] R. G. Pratt and R. M. Shipp, “Seismic waveform inversion in the frequency domain, Part 2: Fault delineation in sediments us-

4624

[20] [21] [22] [23] [24]

[25] [26]

[27] [28] [29] [30] [31] [32]

[33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44] [45]

IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 50, NO. 11, NOVEMBER 2012

ing crosshole data,” Geophysics, vol. 64, no. 3, pp. 902–914, May/ Jun. 1999. B. Zhou and S. A. Greenhalgh, “Crosshole seismic inversion with normalized full-waveform amplitude data,” Geophysics, vol. 68, no. 4, pp. 1320– 1330, Jul./Aug. 2003. T. Watanabe, K. T. Nihei, S. Nakagawa, and L. R. Myer, “Viscoacoustic wave from inversion of transmission data for velocity and attenuation,” J. Acoust. Soc. Amer., vol. 115, no. 6, pp. 3059–3067, Jun. 2004. S. Kuroda, M. Takeuchi, and H. J. Kim, “A full-waveform inversion algorithm for interpreting crosshole radar data,” in Proc. Exp. Abst. SEG 75th Annu. Meeting, Soc. Expl. Geophys, Houston, TX, 2007, pp. 1176–1179. S. Kuroda, M. Takeuchi, and H. J. Kim, “Full-waveform inversion algorithm for interpreting crosshole radar data: A theoretical approach,” Geosci. J., vol. 11, no. 3, pp. 211–217, Sep. 2007. G. A. Meles, J. Van der Kruk, S. A. Greenhalgh, J. R. Ernst, H. Maurer, and A. G. Green, “A new vector waveform inversion algorithm for simultaneous updating of conductivity and permittivity parameters from combination crosshole/borehole-to-surface GPR data,” IEEE Trans. Geosci. Remote Sens., vol. 48, no. 9, pp. 3391–3407, Sep. 2010. C. Zhou, G. T. Schuster, S. Hassanzadeh, and J.-M. Harris, “Elastic wave equation traveltime and waveform inversion of crosswell data,” Geophysics, vol. 62, no. 3, pp. 853–868, May/Jun. 1997. F. A. Belina, J. Irving, J. R. Ernst, and K. Holliger, “Analysis of an iterative deconvolution approach for estimating the source wavelet during waveform inversion of crosshole georadar data,” J. Appl. Geophys., vol. 78, pp. 20–30, Mar. 2012. R. G. Pratt and N. R. Goulty, “Combining wave-equation imaging with traveltime tomography to form high-resolution images from crosshole data,” Geophysics, vol. 56, no. 2, pp. 208–224, Feb. 1991. D. R. Velis and T. J. Ulrych, “Simulated annealing wavelet estimation via fourth-order cumulant matching,” Geophysics, vol. 61, no. 6, pp. 1939– 1948, Nov./Dec. 1996. L. Wang and J. M. Mendel, “Adaptive minimum prediction-error deconvolution and source wavelet estimation using Hopfield neural networks,” Geophysics, vol. 57, no. 5, pp. 670–679, May 1992. Q. Cheng, R. Chen, and T. Li, “Simultaneous wavelet estimation and deconvolution of reflection seismic signals,” IEEE Trans. Geosci. Remote Sens., vol. 34, no. 2, pp. 377–384, Mar. 1996. A. B. Weglein and B. G. Secrest, “Wavelet estimation for a multidimensional acoustic or elastic earth,” Geophysics, vol. 55, no. 7, pp. 902–913, Jul. 1990. F. A. Belina, J. Irving, J. R. Ernst, and K. Holliger, “Evaluation of the reconstruction limits of a frequency-independent crosshole georadar waveform inversion scheme in the presence of dispersion,” J. Appl. Geophys., vol. 78, pp. 9–19, 2012. B. Lampe and K. Holliger, “Effects of fractal fluctuations in topographic relief, permittivity and conductivity on ground-penetrating radar antenna radiation,” Geophysics, vol. 68, no. 6, pp. 1934–1944, Nov./Dec. 2003. J. Irving and R. Knight, “Effect of antennas on velocity estimates obtained from crosshole GPR data,” Geophysics, vol. 70, no. 5, pp. K39–K42, 2005. J. R. Ernst, H. Maurer, A. G. Green, and K. Holliger, “Realistic FDTD modeling of borehole georadar antenna radiation: Methodology and application,” Near Surf. Geophys., vol. 4, pp. 19–30, 2006. K. S. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propag., vol. AP-14, no. 3, pp. 302–307, May 1966. T. Bergmann, J. O. A. Robertsson, and K. Holliger, “Numerical properties of staggered finite-difference solutions of Maxwell’s equations for groundpenetrating radar modelling,” Geophys. Res. Lett., vol. 23, pp. 45–48, 1996. M. Sato and R. Thierbach, “Analysis of borehole radar in cross-hole mode,” IEEE Trans. Geosci. Remote Sens., vol. 29, no. 6, pp. 899–904, Nov. 1991. G. Turner and A. F. Siggins, “Constant Q attenuation of subsurface radar pulses,” Geophysics, vol. 59, no. 8, pp. 1192–1200, Aug. 1994. F. Hollender and S. Tillard, “Modeling ground-penetrating radar wave propagation and reflection with the Jonscher parametrization,” Geophysics, vol. 63, no. 6, pp. 1933–1942, Nov./Dec. 1998. R. J. Knight, “Ground penetrating radar for environmental applications,” Annu. Rev. Earth Planet Sci., vol. 29, no. 1, pp. 229–255, 2001. A. Neal, “Ground-penetrating radar and its use in sedimentology: Principles, problems and progress,” Earth Sci. Rev., vol. 66, no. 3/4, pp. 261– 330, Aug. 2004. A. P. Annan, “GPR methods for hydrogeological studies,” in Hydrogeophysics, Y. Rubin and S. S. Hubbard, Eds. New York: Springer-Verlag, 2005, ch. 7, p. 532. A.Taflove and S. C. Hagness, Computational Electrodynamics: The FiniteDifference Time-Domain Method. Norwood, MA: Artech House, 2000. J. Tronicke, K. Holliger, W. Barrash, and M. D. Knoll, “Multivariate analysis of cross-hole georadar velocity and attenuation tomograms for

[46] [47] [48] [49] [50] [51] [52] [53] [54] [55] [56] [57] [58] [59]

[60] [61] [62]

[63] [64] [65]

[66] [67]

[68] [69]

[70] [71]

aquifer zonation,” Water Resources Res., vol. 40, p. W01 519, 2004. doi:01510.01029/02003WR002031. K. Holliger and T. Bergmann, “Numerical modelling of borehole georadar data,” Geophysics, vol. 67, no. 4, pp. 1249–1257, Jul./Aug. 2002. B. Lampe and K. Holliger, “Resistively loaded antennas for groundpenetrating radar: A modeling approach,” Geophysics, vol. 70, no. 3, pp. K23–K32, May 2005. J. P. Berenger, “A perfectly matched layer for the absorption of electromagnetic-waves,” J. Comput. Phys., vol. 114, no. 2, pp. 185–200, Oct. 1994. J. Fang and Z. Wu, “Generalized perfectly matched layer for the absorbing of propagating and evanescent waves in lossless and lossy media,” IEEE Trans. Microw. Theory Tech., vol. 44, no. 12, pp. 2216–2222, Dec. 1996. E. Polak and G. Ribière, “Note on convergence of conjugate direction methods,” Revue Francaise d’Informatique de Recherche Operationnelle, vol. 3, pp. 35–43, 1969. P. Hashemzadeh, A. Fhager, and M. Persson, “Experimental investigation of an optimization approach to microwave tomography,” Electromagn. Biol. Med., vol. 25, no. 1, pp. 1–12, 2006. A. Fhager and M. Persson, “Using a priori data to improve the reconstruction of small objects in microwave tomography,” IEEE Trans. Microw. Theory Tech., vol. 55, no. 11, pp. 2454–2462, Nov. 2007. F. Costen, J. P. Berenger, and A. K. Brown, “Comparison of FTTD hard source with FDTD soft source and accuracy assessment in Debye media,” IEEE Trans. Antennas Propag., vol. 57, no. 7, pp. 2014–2022, Jul. 2009. W. Menke, Geophysikal Data Analysis: Discrete Inverse Theory. New York: Academic, 1984. A. Fichtner, Full Seismic Waveform Modelling and Inversion. Berlin, Germany: Springer-Verlag, 2011. A. T. de Hoop, Handbook of Radiation and Scattering of Waves: Acoustic Waves in Fluids, Elastic Waves in Solids, Electromagnetic Waves. San Diego, CA: Academic, 1995. J. A. Goff and T. H. Jordan, “Stochastic modeling of seafloor morphology—Inversion of sea beam data for 2nd-order statistics,” J. Geophys. Res.—Solid Earth Planets, vol. 93, pp. 13 589–13 608, 1988. H. H. Hardy and R. A. Beyer, Fractals in Reservoir Engineering. Singapore: World Scientific, 1994. K. Holliger and J. A. Goff, “A generalised model for the 1/f-scaling nature of seismic velocity fluctuations,” in Heterogeneity in the Crust and Upper Mantle—Nature, Scaling, and Seismic Properties. Norwell, MA: Kluwer, 2003, pp. 131–154. M. Kelkar and G. Perez, Applied Geostatistics for Reservoir Characterization. Richardson, TX: Soc. Petroleum Eng., 2002. L. W. Gelhar, Stochastic Subsurface Hydrology. Englewood Cliffs, NJ: Prentice-Hall, 1993. J. Tronicke and K. Holliger, “Quantitative integration of hyrogeophysical data: Conditional geostatistical simulation for characterizing heterogeneous alluvial aquifers,” Geophysics, vol. 70, no. 3, pp. H1–H10, May 2005. B. Dafflon, J. Irving, and K. Holliger, “Simulated-annealing-based conditional simulation for the local-scale characterization of heterogeneous aquifers,” J. Appl. Geophys., vol. 68, pp. 60–70, 2009. W. Barrash and T. Clemo, “Hierarchical geostatistics and multifacies systems: Boise hydrogeophysical reserach site, Boise, Idaho,” Water Resources Res., vol. 38, no. 10, p. 1196, 2002. doi:10.1029/2002WR001436. W. Barrash and E. C. Reboulet, “Singificance of porosity for stratigraphy and textural composition in subsurface, coarse fluvial deposits: Boise hydrogeophysical reserach site,” Geol. Soc. Amer. Bull., vol. 116, no. 9/10, pp. 1059–1073, Sep./Oct. 2004. W. P. Clement and W. Barrash, “Crosshole radar tomography in a fluvial aquifer near Boise, Idaho,” J. Environ. Eng. Geophys., vol. 11, no. 3, pp. 171–184, Sep. 2006. C. J. Mwenifumbo, W. Barrash, and M. D. Knoll, “Capacitive conductivity logging and electrical stratigraphy in a high-resistivity aquifer, Boise hydrogeophysical research site,” Geophysics, vol. 74, no. 3, pp. E125– E133, Apr. 2009. L. Knopoff and A. F. Gangi, “Seismic reciprocity,” Geophysics, vol. 24, no. 4, pp. 681–691, Oct. 1959. K. R. Roth, R. Schulin, H. Fluhler, and W. Attinger, “Calibration of time domain reflectometry for water content measurement using a composite dielectric approach,” Water Resources Res., vol. 26, no. 10, pp. 2267– 2273, 1990. A. Abubakar, P. M. van den Berg, and S. Y. Semenov, “A robust iterative method for Born inversion,” IEEE Trans. Geosci. Remote Sens., vol. 42, pp. 342–354, Feb. 2004. R. Persico and F. Soldovieri, “Effects of background removal in linear inverse scattering,” IEEE Trans. Geosci. Remote Sens., vol. 46, no. 4, pp. 1104–1114, Apr. 2008.

BELINA et al.: INFLUENCE OF SOURCE WAVELET VARIABILITY

Florian A. Belina received the M.Sc. degree in geophysics from the Swiss Federal Institute of Technology (ETH) in Zurich, Switzerland, in 2006, and the Ph.D. degree in the same field from the University of Lausanne, Lausanne, Switzerland, in 2011. His research interests focus on the full-waveform inversion of crosshole georadar data. A particular emphasis of his recent work has been on the estimation of the source wavelet and its impact on the tomographic imaging process.

James Irving holds a Senior Lecturer position in the Institute of Geophysics at the University of Lausanne, Lausanne, Switzerland. He received the M.Sc. degree in geophysics from the University of British Columbia, Vancouver, in 2000, and the Ph.D. degree in geophysics from Stanford University, Stanford, CA, in 2006. From 2007 to 2009, he worked as a Postdoctoral Scientist at the University of Lausanne. After spending one year as an Assistant Professor in the School of Engineering at the University of Guelph in Canada, he returned to Switzerland to take up his current faculty position in February 2011. Since 2009, Dr. Irving has served as an Associate Editor for the journal Geophysics. He is currently also serving as the President of the Near-Surface Geophysics Section of the Society of Exploration Geophysicists. In April 2011, he was awarded the Environmental and Engineering Geophysical Society and Geonics Limited Early Career Award, which acknowledges academic excellence and research contributions in the field of near-surface geophysics.

4625

Jacques R. Ernst (S’07–A’07)received the M.Sc. and Ph.D. degrees in applied and environmental geophysics from the Swiss Federal Institute of Technology (ETH) in Zurich, Switzerland, in 2002 and 2007, respectively. Following the Ph.D. degree, he was with UBS AG Zurich. Since April 2010, he has been with EGL, a major energy trading company. His research interests include simulation of realistic ground-penetrating radar (GPR) antennas and full-waveform modeling and inversion of crosshole-borehole GPR data. Dr. Ernst is a member of the Society of Exploration Geophysicists.

Klaus Holliger received the M.Sc. and Ph.D. degrees in geophysics from the Swiss Federal Institute of Technology (ETH) Zurich, Zurich, Switzerland, in 1987 and 1990, respectively, and a postgraduate degree in economics (2000) from the University of London, London, U.K. After an extended postdoc at Rice University in Houston, TX, he joined ETH’s newly founded applied and environmental geophysics group in 1994. He has also worked for shorter periods of time at the United States Geological Survey, Imperial College, and the University of Cambridge. In 2005 he moved to a Chaired Professorship at the University of Lausanne, Lausanne, Switzerland. He finished a term as a Vice-Dean of research and is now the Director of the university’s geophysics institute. Since 2008, he has also been serving as the Editor-in-Chief of the Journal of Applied Geophysics. He has broad scientific interests and has worked in a variety of fundamental and applied research domains. To date, he has coauthored more than a hundred peer-refereed publications. Currently, his main research interest is the emerging and inherently interdisciplinary field of hydrogeophysics.