HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS		QUICK SEARCH:   [advanced] Author: Keyword(s): Year:  Vol:  Page:

This Article Abstract Figures Only Full Text (PDF) Free Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Hansen, T. M. Articles by Nielsen, L. Search for Related Content PubMed Articles by Hansen, T. M. Articles by Nielsen, L. GeoRef GeoRef Citation Agricola Articles by Hansen, T. M. Articles by Nielsen, L. Related Collections Ground Penetrating Radar, GPR Tomography Inverse Procedures/Parameter Estimation

SPECIAL SECTION: GROUND PENETRATING RADAR IN HYDROGEOPHYSICS

Inferring the Subsurface Structural Covariance Model Using Cross-Borehole Ground Penetrating Radar Tomography

^a Univ. of Copenhagen, Niels Bohr Institute, Juliane Maries Vej 28, DK-2100 Copenhagen Ea, Denmark
^b Univ. of Copenhagen, Geological Institute, Oester Voldgade 10, DK-1350 Copenhagen K, Denmark
^* Corresponding author (tmh{at}gfy.ku.dk).

All rights reserved. No part of this periodical may be reproduced or transmitted in any form or by any means, electronic or mechanical, including photocopying, recording, or any information storage and retrieval system, without permission in writing from the publisher.

Received 29 September 2006.

	ABSTRACT

TOP ABSTRACT INTRODUCTION Theory--Posterior Validation of... Case Study: Inferring the... Real Data Example Conclusions Appendix REFERENCES

 
We address a fundamental problem inherent in least squares based ground penetrating radar tomography problems, and linear inverse Gaussian problems in general: how should the a priori covariance model be chosen? The choice of such a prior covariance model is most often a very subjective task that has major implications on the result of the inversion. We present a method that allows quantification of the likelihood that a given choice of prior covariance model is consistent with the observed tomography data. This is done by comparing statistical properties of samples of the prior and posterior probability density function of the tomographic inverse problem. In essence, if samples of the posterior are unlikely samples of the prior, then such a choice of a priori covariance model is deemed unlikely. This enables one to quantify the consistency of a number of equally probable prior covariance models to data observations. A synthetic data set was used to describe and validate the approach. We determined how a known covariance model could be inferred from a synthetic tomography problem. The methodology was then applied to a nonlinear ground penetrating radar tomography case study. The covariance model deemed most likely was consistent with nearby ground penetrating radar reflection profiles. The method provides useful results even if just a subset as small as 10% of the available data is considered.

Abbreviations: GPR, ground penetrating radar • HPD, highest probability density • McMC, Markov chain Monte Carlo • pdf, probability density function

	INTRODUCTION

TOP ABSTRACT INTRODUCTION Theory--Posterior Validation of... Case Study: Inferring the... Real Data Example Conclusions Appendix REFERENCES

 
Cross-borehole tomography has become a widely popular method to obtain structural information about the near-surface aquifer and reservoir properties (Knight et al., 1997; Rea and Knight, 1998; Hubbard et al., 1999; Tercier et al., 2000). The most widely used method of inverting such tomography data is a parametric approach based on least squares inversion, which implies that a multi-Gaussian mathematical model has been selected to describe the inverse problem at hand. The implicit assumption is that the subsurface can be described by a multi-Gaussian probability density function (pdf) with a given mean, m_p, and covariance, C_Mp. The covariance model describes the spatial correlations of the subsurface. We refer to this pdf as the prior pdf, describing the a priori information. Measurement errors are assumed Gaussian, given by C_d.

In a Bayesian formulation, for the linear case, where data observations, d, are related to the subsurface model parameter space, m, through the linear operator g, such that d = gm, the least squares solution is a Gaussian pdf with mean m_est and covariance C_Mest (Tarantola and Valette, 1982, Eq. [41–42]):

[1]

[2]

where g' is the transpose of g.

We refer to this pdf as the a posteriori pdf. Samples of this posterior pdf can be used to describe and analyze the solution of linear(ized) inverse tomographic problems. Equations [1] and [2] explicitly take into account the data and model covariance matrices. This Bayesian approach to solving least squares based inverse problems is referred to as generalized inversion using the least squares criterion (Tarantola and Valette, 1982) or stochastic regularization (Maurer et al., 1998; Linde et al., 2006).

Inverse problems associated with cross-borehole tomography are typically underdetermined. There exist infinitely many sets of m that will fit the data within their uncertainties. The Bayesian approach described by Eq. [1] and [2] reflect this nonuniqueness in that the solution to the inverse problem is a number of (infinitely many possible) realizations of the posterior pdf. Until recently, it has been unfeasible to generate such realizations of the posterior. Instead, the focus has been on generating one model, and fitting the data observations within their uncertainties, subject to some regularization constraint.

Probably the most used regularization techniques, used with cross-borehole tomography throughout the vadose zone, is based on smoothing (i.e., Occam's inversion as given by Constable et al. [1987]), which seeks the simplest model that will fit the data, and damping (Marquardt, 1970), which seeks a model with similarity to a specific (prior) model. The smoothing properties of Occam's inversion can be described a priori through C_Mp using Eq. [1] and [2] above (Gouveia and Scales, 1997). Damping can likewise be specified a priori by altering C_d appropriately. In stochastic regularization, a specific choice of covariance model is selected that imposes some smoothness on the model parameters, and m_est (Eq. [1]) is then selected as the solution (see, e.g., Linde et al., 2006).

The main advantage using the generalized least squares approach as given in Eq. [1] and [2] is that the prior chosen C_Mp relates directly to the underground geophysical and geological parameters, whereas the use of smoothing or damping provides algorithm-dependant regularization that is not directly linked to the subsurface geology and geophysics.

There are two major reasons why smoothed and damped regularized inversion sometimes is preferred to the general least squares inversion of Eq. [1] and [2]. First, for large three-dimensional problems, it may not be computationally feasible to calculate C_Mp^–1 (Maurer et al., 1998; Linde et al., 2006). Second, C_Mp may be poorly known in many cases. The problem of computing C_Mp^–1 can be circumvented. Linde et al. (2006) suggested to use circulant matrices to efficiently calculate C_Mp^–1. A different approach is to solve the least squares inversion problem sequentially (point by point), as suggested by Hansen et al. (2006). This involves solving many, but small (and hence computationally fast), linear inverse systems. The second problem of finding a good prior choice of C_Mp was the focus of this study.

An additional problem using least squares linear inversion is that, in practice, m_est is most often used as "the" solution to the inverse tomography problem (Phillips and Fehler, 1991; Cassiani et al., 1998; Day-Lewis and Lane, 2004; Linde et al., 2006). While m_est is a realization of the linear inverse problem, as it is the center of the a posteriori pdf, it is also the so-called minimum variance realization of the a posteriori pdf. Of all possible realizations of the posterior pdf, it is the one sample with the least spatial variation. In this sense, it may not be a very realistic model of the subsurface. In practice, the model with maximum posterior likelihood may, in fact, be of little practical interest, as discussed by Tarantola (2005, p. 54).

This smoothing effect of least squares estimation imposes a serious problem dealing with models of the vadose zone. For example, co-kriging can be used to estimate hydrologic parameters (such as moisture content) using m_est as secondary data. Thus, a linear relation between the hydrologic parameter of interest and the velocity model, m_est, is assumed (Cassiani et al., 1998). Therefore the smoothness of the velocity model will be propagated into the estimate of the hydrologic parameter, which will also be too smooth. In effect, extreme values are disregarded, which may have a great impact on subsequent flow modeling. Day-Lewis and Lane (2004) and Day-Lewis et al. (2005) concluded that the smoothness effect observed when applying regularization techniques has a strong negative effect on the quality of the correlation obtained between hydrologic properties and the tomographic least squares estimate.

The observed smoothing effect is an effect of using only m_est as a solution. Gloaguen et al. (2005) showed how one realization consistent with data observations and the a priori chosen covariance model can be obtained. Hansen et al. (2006) described how actual samples of the posterior pdf can be generated, enabling stochastic inversion, as opposed to the deterministic inversion. Using this method, samples of the posterior pdf will show the correct (unsmoothed) spatial variability as given by the a priori chosen covariance model.

Whether one is interested in performing least squares based deterministic or stochastic inversion, the choice of prior covariance model will have a great effect on the inversion result. In addition, it is also a choice that cannot be avoided: any use of least squares based linear inversion requires one to specify such a prior covariance model, whether it be directly or through the use of damping or smoothing. Therefore it is critical that the used covariance model be chosen with the greatest care. This problem of selecting the prior covariance model where only tomography data are available is what we address here.

Inferring the Subsurface Covariance Model
Several researchers have acknowledged the importance of using a correct choice of covariance model, and some work has been done to estimate such a covariance matrix. Knight et al. (1997), Rea and Knight (1998), and Tercier et al. (2000) computed experimental covariance models directly from two- and three-dimensional ground penetrating radar (GPR) reflection profiles. Hubbard et al. (1999) inferred the spatial correlation length by investigating the power spectra of well log data and linear average data in the form of a smooth minimum variance estimate of a tomographic inverse problem. The longer wavelengths were obtained from the smooth inverted tomography data and the shorter wavelengths were obtained from the well log data. Oldenborger et al. (2003) showed that the horizontal correlation length obtained from semivariogram analysis of GPR reflection data can be used to directly infer the horizontal correlation structure of hydraulic conductivity.

All the approaches considered above assumed that information additional to tomographic data are available in the form of either well log data or GPR reflection data. Gloaguen et al. (2005) suggested a method to infer the subsurface covariance that relies only on observed travel time delays in a cross-borehole setup. This method was implemented by Giroux et al. (2007). They suggested comparing an experimental data covariance to a theoretical data covariance obtained from a linear data kernel and an assumed model covariance model. When the theoretical and experimental data covariances matched to a certain level, the corresponding model covariance model was selected for subsequent inversion.

Here we present a new approach to infer the subsurface covariance model, based only on available travel time data as observed from cross-borehole GPR tomography. By comparing samples of the a posteriori pdf to samples of the a priori pdf, a number of equally probable a priori models can be distinguished in terms of how consistent they are with respect to data observations. We show the method applied to linear least squares based cross-borehole GPR stochastic inversion and discuss the differences between our approach and the method suggested by Gloaguen et al. (2005). The method is viable for use with any (non)linear inverse problem for which samples can be generated from both the prior and posterior pdf.

	Theory—Posterior Validation of Prior Information

TOP ABSTRACT INTRODUCTION Theory--Posterior Validation of... Case Study: Inferring the... Real Data Example Conclusions Appendix REFERENCES

 
Consider some data observations, d, which are related to some model parameter space, m, through the (possibly nonlinear) relation, g. The forward problem is the problem of estimating d from g and m:

[3]

The inverse problem consists of determining properties of the model parameter space given measurements, d, the (often physically based) relation, g, and some prior information. A general probabilistic solution to the inverse problem is given by Tarantola (2005):

[4]

where k is a normalization factor, (d,m) represents information on the physical relation between d and m, (d,m) represents the a priori pdf, and (d,m) represents the a posteriori pdf. We refer to the solution of the inverse problem as the a posteriori pdf (d,m). Sometimes the choice of prior information is not straightforward, and several equiprobable choices of prior model may exist. We suggest validating the choice of a priori model by comparing the a posteriori pdf to the a priori pdf. Specifically, if samples of the a posteriori pdf are probable (above some probability threshold) samples of the chosen a priori pdf, we say that the choice of prior information is consistent with the observed data. If, on the other hand, a sample of the posterior is a very unlikely sample of the a priori pdf, we say that the choice of a priori information is inconsistent with data observations.

Equation [4] is completely general and applies to any inverse problem that can be formulated in such a probabilistic way. Here we will consider the special group of Gaussian linear inverse problems, where g is a linear operator such that d = gm. Both (d,m) and (d,m) can be described by a multi-Gaussian pdf; (d,m) is also multi-Gaussian and can be calculated from data observations with assumed Gaussian errors. A variant of this formulation of an inverse problem is typically used in GPR tomography as discussed above. The linear and Gaussian assumption leads to the simple and beautiful mathematics of Eq. [1] and [2] that completely describes the Gaussian a posteriori pdf, that is, the solution to a linear inverse Gaussian problem. Adopting the linear Gaussian approach, one must a priori choose a Gaussian pdf describing the a priori information. In GPR tomography, this amounts to selecting a covariance model that describes the spatial variation in the subsurface velocity field.

We show how to build a probabilistic model of the ergodic variations of the mean and the covariance of a prior choice of Gaussian pdf. Then we compute the likelihood that a sample from the posterior pdf is a realization of the prior distribution. In this way, we show how to effectively validate the consistency of choices of Gaussian prior models with data observations.

Validation of the Prior Gaussian Probability Density Function
It may seem trivial to consider how consistent a specific choice of prior covariance model is with data observations. For Gaussian inverse problems, by definition, all samples of the a posteriori pdf are possible samples of the a priori pdf. It is only when a uniform prior distribution is used, however, that all values are equally probable samples of the a priori distribution.

Consider a simple one-dimensional Gaussian inverse problem. Both the prior and posterior are given by one-dimensional Gaussian pdfs. From here on, N(µ,C) refers to a Gaussian distribution with mean µ and covariance C. Two geologists provide their best choice of a priori distribution of a single model parameter. The first geologist, G1, chooses a prior Gaussian pdf as ₁ = N(4000,300), and the other geologist, G2, chooses ₁ = N(6000,300). The change in shape between the prior and posterior is caused by the addition of data observations into the inverse problem. Then 1000 samples of the posterior are generated for each of these two inverse problems, where the prior information differs and data observations and the transfer functions are identical. The corresponding posterior pdfs are described by the one-dimensional Gaussian pdfs ₁ = N(4100,100) and ₂ = N(4300,100). Figure 1 gives a graphical representation of the prior and posterior pdfs for the two cases.

View larger version (8K): [in this window] [in a new window]   FIG. 1. Prior probability density function (pdf) given by geologist no. 1 (light yellow) and geologist no. 2 (light gray). Posterior pdf using the prior information given by geologist no. 1 (dark yellow) and geologist no. 2 (dark gray).

[5]

[6]

Thus for G1, we see that samples of the posterior pdf, ₁, are also highly probable samples of the prior pdf, ₁. For G2, we can see that samples of the posterior pdf, ₂, are extremely unlikely samples of the prior pdf, ₂. Having computed the posterior of both inverse problems (using the same data observations, but different a priori information), we can conclude that G1's initial choice of prior model is much more consistent with data observations than the choice of prior model of G2. We can, in effect, falsify the prior model choice of G2 as a valid choice of prior model. The method we propose to quantify how well a specific choice of prior covariance model is consistent with data observations is but a generalization of this one-dimensional Gaussian example.

Sequential simulation can be used to generate actual samples of both the prior and posterior pdf of Gaussian linear inverse problems such as, for example, GPR cross-borehole tomography (Hansen et al., 2006; Hansen and Mosegaard, 2007). For a large number of samples from the a priori pdf, a Gaussian probability model can be built that describes the expected variation of both the mean and variance using a specific choice of a priori model. From this probability model, it is trivial to calculate the likelihood that a specific realization is a sample of the a priori pdf. Then we generate a number of samples from the a posteriori pdf. From these samples, we calculate the likelihood that each posterior realization is a sample from the a priori distribution. This enables us to calculate the average probability that a sample of the a posteriori pdf is a sample of the a priori pdf.

Building a Probabilistic Model of Ergodic Variability of the Prior Probability Density Function
Samples from the Gaussian a priori and a posteriori pdfs are equivalent to samples of a multivariate Gaussian random function, Z(x), for which the first- and second-order moments are the mean and covariance, N(m_p,C_Mp) (see Goovaerts [1997] for an introduction to the multivariate Gaussian model). Thus, examining the variability of samples from the a priori (and posterior) pdf is equivalent to examining the variability of the samples of the corresponding Gaussian random functions. The reason we introduce random functions is to apply theory from geostatistics, which is largely based on the concept of random functions.

Due to the finite size of the model parameter space, the mean and covariance of samples of a Gaussian random function fluctuates as ergodic variations around the a priori chosen values. The random function is said to be ergodic in the mean if the mean of a realization of the random field tends to m_p as the size of the model parameter space tends to infinity, and to be ergodic in the covariance if the covariance of a realization of the random field tends to C_Mp as the size of the model parameter space tends to infinity (Chiles and Delfiner, 1999).

A measure of the ergodicity can be calculated theoretically (Ortiz and Deutsch, 2002; Emery, 2004; Marchant and Lark, 2004), but as the observed ergodicity is always also a product of the algorithm used, we will make use of an experimentally assessed measure of ergodicity. As an example of ergodic fluctuations, we consider the ergodicity in the mean and covariance as observed from drawing 100 samples from a Gaussian random function (i.e., a choice of prior model) with a mean of 0.13 m/ns and a spherical anisotropic covariance with ranges 10, 20, and 40 m, and a variance of 2 x 10^–4 m²/ns². For bounded semivariogram models, the variance is also referred to as the sill value (see Goovaerts [1997] for a detailed description of choices of covariance models). From here on, we will refer to the semivariogram interchangeable with the covariance as the semivariogram and the covariance is equivalent following (h) = C(0) – C(h) for stationary random functions, which is what we consider here.

Figure 2 shows three independent samples of an unconditional simulation in a 100 by 100 grid with cell size 1 m and an isotropic range of 10, 20, and 40 m. The distribution of the mean of the 100 realizations is shown in Fig. 3 , and the variation of the semivariogram of the realizations is shown in Fig. 4 . The nature of the ergodic fluctuations is heavily dependent on the size of the model grid and the choice of semivariogram model. Specifically, the amplitude (range) of variability around the a priori chosen semivariogram model increases as the range or model size fraction increases.

View larger version (100K): [in this window] [in a new window]   FIG. 2. Three unconditional realizations using a range of 10, 20, and 40 m in a model of size 100 by 100 m, using a grid cell size of 1 m.

View larger version (11K): [in this window] [in a new window]   FIG. 3. Ergodic variation of the mean computed from 100 unconditional realizations using a range of 10, 20, and 40 m. The red curve is the best-fitting normal distribution to the observed data. The vertical line indicates the a priori choice of mean velocity (0.13 m/ns).

View larger version (17K): [in this window] [in a new window]   FIG. 4. Ergodic variation of the semivariance computed from 100 unconditional realizations using a range of [10,20,40], with experimental semivariogram (black lines) and theoretical semivariogram (blue line). Mean and 95% confidence intervals for the variance of the experimental semivariogram are the red solid and dashed lines, respectively.

Consider again the 100 realizations considered in Fig. 2 to 4. As expected, the ergodic variability of the mean, Fig. 3, can be reasonably approximated by a Gaussian distribution. The red line on Fig. 3 corresponds to the Gaussian pdf, N(µ_m*,_m*), found to match best the observed variability in the mean. The experimental ergodicity in the semivariance is described by the multivariate Gaussian pdf, N(*,*), obtained from the observed ergodic variation in the semivariance from 100 realizations (Fig. 4). The mean semivariance for each of N separation distances k is given by * = [*¹...*^N]. The covariance matrix, *, describing the covariance between the semivariance at different lags, is computed from the semivariance values for each separation distance, _k, as _kk' = Cov(_k,_k'). The diagonal of * is but the observed variance of the semivariance for each separation bin, as illustrated by the 95% confidence interval in Fig. 4.

Likelihood that One Realization Is a Sample of the Prior
As the ergodic variation in both the mean and the semivariance can be described by Gaussian statistics, we can compute the likelihood that any given model, m*, is a realization of a specific multivariate Gaussian pdf (describing the observed ergodicity). In other words, we can compute the likelihood that any realization of the posterior pdf is a realization of the prior pdf.

The ergodicity in the mean is described by a one-dimensional Gaussian pdf, with mean µ_m* and variance _m*, and therefore the likelihood that a given realization is a sample of a Gaussian pdf with observed mean m_exp is given by

[7]

In a similar manner, the likelihood that any given sample, m*, with a specific experimental semivariogram _exp = [₁, ..., _N], is a sample of the prior pdf in terms of reproducing the ergodic variability of the semivarigran given by N(*, *) can be found using Tarantola (2005):

[8]

Using Eq. [8], we compute the likelihood that any of the 100 realizations, generated using a range of 10, 20, and 40 m, is a sample of a prior distribution with the range of [5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75] m. The computed logarithmic likelihood is plotted in Fig. 5 . There is a clear maximum around the a-priori-selected correlation length. Of the 100 realizations, 100, 100, and 97 find a maximum likelihood corresponding to the a priori chosen range of 10, 20, and 40 m, respectively. As expected, the ergodic fluctuations are larger when the fraction between the width of the model parameter space considered and the range becomes smaller. Therefore it will be more difficult to infer the covariance model when the subsurface range is relatively long, which is reflected by the slightly worse inference of a correlation of 40 m.

View larger version (21K): [in this window] [in a new window]   FIG. 5. Likelihood of 100 samples of a Gaussian random function with a specific range of 10, 20, and 40 m being a realization of a Gaussian random function with range from 5 to 75 m.

Likelihood That a Sample from the a Posteriori Probability Density Function Is a Sample of the a Priori Probability Density Function
In reality, inverse problems are never completely resolved, and therefore a number of models, all realizations of the a posteriori pdf, must be considered. Therefore our aim is to calculate the likelihood that a sample (a set of realizations) of a Gaussian pdf (specifically, the a posteriori pdf) can be seen as a sample of another Gaussian pdf (specifically, the a priori pdf). We do this by calculating the average likelihood, L_av, that each of N realizations from the a posteriori pdf, with experimental covariance _exp = [_exp1, ..., _expN], is a realization from the prior pdf described by the Gaussian pdf N(*,*) as

[9]

In the case of no data observations, the prior and the posterior are the same, and * will be constant and independent of the choice of covariance model, L_av(*) = L_av, meaning that all choices of a priori models are equally probable.

Using Eq. [9], the consistency of a choice of prior covariance model to data observations can now be quantified. We will make use of this in the following case study to quantify the consistency of a number of choices of prior covariance models to travel time delay observations, considering both a synthetic and a real data set.

Equation [9] is appropriate to use in the case where the covariance model is isotropic (when the correlation range is independent of direction). The appendix contains the mathematical basis for calculating a likelihood measure based on ergodicity observed using anisotropic covariance models, and the joint use of ergodicity in the mean and the covariance.

Comparison with Other Techniques
Of the methods known to infer the subsurface covariance model before performing inversion, the method developed by Gloaguen et al. (2005) and Giroux et al. (2007) also relies solely on travel time delay data. The method is based on the comparison of an experimental and synthetic data covariance matrix. The researchers did not, however, in these studies, document how well the proposed method works. It has not been tested on synthetic data, and the covariance models eventually used in the case studies of both Gloaguen et al. (2005) and Giroux et al. (2007) were not presented. The method suggested by Gloaguen et al. (2005) relies on calculating the experimental data covariance function from the smooth m_est model. The smoothness of m_est entails that high-frequency spatial variation is absent, and, therefore, the covariance model of the least squares estimate cannot be used to capture models with significant roughness. Due to the smoothing effect, we expect that using the method of Gloaguen et al. (2005) to infer the subsurface covariance model will result in underestimating the sill and overestimating the range of anomalies of the subsurface.

We suggest modifying the approach of Gloaguen et al. (2005) such that the experimental data covariance matrix is calculated from a (large) number of realizations from the posterior pdf. This will ensure that the data from which the experimental data covariance matrix is calculated possess the spatial features requested by the prior choice of covariance model, if allowed by the observed data.

The approach of Gloaguen et al. (2005) is viable only for Gaussian inverse problems, whereas the method we propose here is viable for any Bayesian inverse problem where samples can be generated from both the a priori and the a posteriori pdf.

	Case Study: Inferring the Subsurface Covariance Model Using Cross-Borehole GPR Tomography

TOP ABSTRACT INTRODUCTION Theory--Posterior Validation of... Case Study: Inferring the... Real Data Example Conclusions Appendix REFERENCES

 
As shown, the proposed methodology can be used to successfully infer the subsurface covariance model of a model of the underground in the special case where the subsurface is completely resolved. We use the methodology to infer the subsurface covariance function in case no direct samples of the model parameter space are available, but where linear average measurements of the subsurface are available in the form of recorded travel time data from a cross-borehole GPR tomography experiment. To validate the approach, we will initially consider a synthetic data set making use of the actual geometry used to collect the true data. Subsequently, we will apply the method to the real observed data.

Data Geometry
A cross-borehole GPR tomography experiment is currently being performed in Arrenæs (55.982002° lat., 12.079760° long.), northern Sjælland, Denmark. The aim of this integrated geophysical–hydrogeologic project is to estimate the radar wave velocity distribution between a number of boreholes. Subsequently, the moisture content of the sedimentary structures can be estimated using an empirical relationship between radar wave velocity and moisture content (Looms et al., 2007). The collected moisture content variations are subsequently used to estimate the subsurface infiltration rate, which is important to the rate of degeneration of nutrients in the unsaturated zone.

A transmitter is lowered into one of two boreholes. With the transmitter position fixed, a receiver is lowered into the second borehole. For every 0.25 m the receiver is lowered, the transmitter emits an electromagnetic pulse, and the propagating signal is recorded at the receiver. Once the receiver has reached the base of the borehole, the transmitter is lowered 1 m, and the recording process is repeated. This continues until the transmitter reaches the base of the borehole, at which point the antennas are interchanged. In this way, GPR data sections are obtained. For each ray, the first arrival time, along with the ray geometry, gives a measure of the mean velocity along a ray path. The observed mean velocity is a linear average of the subsurface. Sequential Gaussian simulation can now be used to draw samples from the a posteriori pdf of the linear inverse tomography problem (Hansen et al., 2006). We consider a model parameter space of 12 by 5 m using a regular grid of 0.25- by 0.25-m cells.

To illustrate how well the proposed method can be expected to work, we initially treat the GPR problem as linear, for the synthetic case study, assuming that rays traverse the model parameter space along straight lines, as shown in Fig. 6 . As the method is based on linear least squares based inversion, it provides optimal results only when the inverse problem is linear or can be linearized. In the real data case study, we make use of bending rays, obtained through linearization of the nonlinear inverse problem, and discuss the implications of using the straight ray approximation.

View larger version (78K): [in this window] [in a new window]   FIG. 6. Geometry of tomography setup. Each line represents a ray. The color of each line indicates the mean velocity observed along the ray.

Inferring the Subsurface Horizontal and Vertical Correlation Lengths
First, we assume the dip is known (i.e., we assume horizontal layers, and β = 0°). We consider four cases of data availability: 10, 30, 70, and all 702 ray data are in turn considered. The 10-, 30-, and 70-ray data were randomly selected from all 702 available data. For each data geometry, Fig. 7 shows the 5, 50, and 95% highest probability density (HPD) regions (Box and Tiao, 1992). These plots are based on approximately 200 sets of evaluated likelihood values, according to Eq. [9] and the use of anisotropy as defined in the appendix. The HPD plot is a convenient way to visualize the two-dimensional likelihood function obtained. Pairs of (h_max,h_min) located outside the 5% HPD contour line are located in an area where only 5% of the total mass of probability is present. In effect, pairs of (h_max,h_min) located in this white area are unlikely to be consistent with data observations, whereas the prior choice pairs of (h_max,h_min) within the 95% HPD area are highly consistent with data observations. Figure 8 shows the corresponding one-dimensional marginal posterior pdfs. From these figures it is evident that, as the number of data increase, the HPD intervals shrink, implying better constraints on h_max and h_min. It proves difficult to obtain a precise estimate of h_max because of the large horizontal range, which is on the order of the width of the subsurface model space across which ray data is sampled. In addition, rays traversing the model parameter space run mostly horizontally, and, therefore, the sensitivity of the observed data is highest with respect to the vertical variability. We conclude that h_max is larger than h_min, and that h_max is probably >3.5 m.

View larger version (21K): [in this window] [in a new window]   FIG. 7. The 5% (light gray), 50% (gray), and 95% (black) highest probability density regions using (a) 10, (b) 30, (c) 70, and (d) 702 rays. The red circle denotes the choice of horizontal and vertical correlation lengths (hmax, hmin) = (6 m, 1.5 m) used to generate the reference model.

View larger version (13K): [in this window] [in a new window]   FIG. 8. One-dimensional marginal distributions of (left) horizontal correlation length (hmax) and (right) vertical correlation length (hmin), using data observations for 10, 30, 70, and 702 rays.

Inferring the Subsurface Correlation Lengths and Dip
Now we consider both the correlation lengths (h_max and h_min) and the dip (β) to be unknown, using the same synthetic data set as considered above. We make use of a Metropolis–Hastings algorithm to sample the three-dimensional likelihood function (Hastings, 1970). The Metropolis–Hastings algorithm is a random walk sampler that ensures that all accepted points are sampled proportional to the likelihood density. Therefore a one-dimensional histogram of accepted (h_max, h_min, β) is simply the one-dimensional marginal distribution of the likelihood function.

Figure 9 shows the one-dimensional marginal a posteriori distribution for h_max, h_min, and β based on 200 sampled likelihood values. The correlation length and dip used to generate the reference model are located in areas of nonzero likelihood. Other set of parameters, however, are deemed much more likely. From the one-dimensional marginal pdfs, the optimal correlation lengths and dip are estimated at (h_max, h_min, β) = (4 m, 1.2 m, 15°). The vertical range is thus found relatively well (1.2 m compared with the reference of 1.5 m), whereas the horizontal range and dip show some variation from the parameters chosen for the reference model. Again, the bias may be due to the ergodic variability: the reference model may be a more probable realization of a covariance model with (h_max, h_min, β) = (4 m,1.2 m, 15°) even though it is generated as a realization of a random function with a covariance model with (h_max, h_min, β) = (6 m,1.5 m, 0°). The effect of ergodicity will be most pronounced in the case where the correlation length is relatively long with respect to the scale of the model parameter space. Therefore it is no surprise that the vertical correlation length is the easiest to infer. This latter example also shows that, in a case where the dip is rather well known, the value of the dip should be fixed and only the correlation lengths inferred. And, in a case where the horizontal correlation length can be inferred from, for example, GPR reflection profiles, this could also be fixed, such that only the vertical correlation length should be inferred. The more parameters to infer, the more uncertain the inferred parameters become.

View larger version (10K): [in this window] [in a new window]   FIG. 9. One-dimensional marginal distributions of (left) horizontal correlation length (hmax), (middle) vertical correlation length (hmin), and (right) dip (β), using data observations for all 702 rays. The dashed line indicates the value used to generate the reference random field.

	Real Data Example

TOP ABSTRACT INTRODUCTION Theory--Posterior Validation of... Case Study: Inferring the... Real Data Example Conclusions Appendix REFERENCES

 
We now consider the actual observed data collected with the same source–receiver geometry as considered in the synthetic case study. The central frequency of the antennae used to obtain the data was 100 MHz.

A Priori Information
Some information about the subsurface is available in the form of GPR reflection profiles, well logs, and some actual samples, all recorded or obtained close to the borehole locations used for this case study. Thus, we know that the subsurface primarily consists of gravel and sand, and we have some knowledge about the spatial connectivity of the sedimentary layers.

Horizontal and Vertical Variability
A GPR reflection profile, Fig. 10 , collected close and parallel to the case study area, with a vertical resolution of about 0.3 to 0.5 m, gives indications about the expected spatial variability. An automatic gain control operator with an operator length of 50 ns was applied to the traces of the GPR reflection section to correct for spherical spreading, and frequencies above 160 MHz have been suppressed by the application of a high-cut frequency filter. The GPR reflection section was migrated using a constant root-mean-square velocity of 0.11 m/ns. Note that the average velocity of 0.13 m/ns (see below) of the zone sampled by the cross-borehole GPR data is higher than the migration velocity, because the cross-borehole data do not sample the uppermost 2 m of the subsurface, which have a relatively lower GPR wave velocity. The GPR profile shows a layered subsurface with clear evidence of layers dipping about 6.5° toward the northeast. From visual inspection of the records, we estimate the lateral continuity of the individual layers to vary between 4 and 20 m, and the layer thickness typically is interpreted to be on the order of 1 to 5 m.

View larger version (42K): [in this window] [in a new window]   FIG. 10. Ground penetrating radar reflection profile recorded close and parallel to the cross-borehole profile considered. The green area indicates the approximate location of the model parameter space considered.

Data Uncertainty
The first arrival times of the GPR data section have been picked with an accuracy of about 0.4 ns, which was chosen as the level of measurement uncertainty. For some GPR traces, the uncertainty is even higher, due to a refracted wave traveling along the soil–air interface. Fortunately, these refracted signals are strongly correlated and can be easily identified in the GPR records (Peterson, 2001). The GPR data recorded in the topmost 1.5 m of the boreholes were therefore excluded in the inversion. For all subsequent simulation, these measurement errors were incorporated as independent data uncertainties by adding the measurement variance to the diagonal of C_d in Eq. [1] and [2].

GPR Tomography as a Nonlinear Inverse Problem
Previously, we dealt with GPR tomography as a purely linear inverse problem, i.e., making use of straight rays. In reality, GPR tomography is a weakly nonlinear inverse problem in the sense that the nonlinear problem can be linearized using, for example, iterative techniques (see Tarantola, 2005). Once linearized, VISIM can be used to generate samples of the posterior pdf of the nonlinear inverse problem. This will, of course, only work properly to the extent that the GPR tomography problem can be properly linearized. In the remainder of this study (except where indicated), all results have been based on linearization of the inverse tomography problem. Specifically, for each considered covariance model, the inverse problem was linearized before samples were generated from the posterior pdf. Figure 11 indicates the degree of nonlinearity by comparing a linear ray kernel to a linearized sensitivity kernel obtained using a covariance model with h_max = 10 m and h_min = 3 m. The ray kernel corresponds to g in Eq. [1].

View larger version (71K): [in this window] [in a new window]   FIG. 11. Sensitivity kernel using a (a) linear and (b) nonlinear (linearized) approach, in the case considering a covariance model with horizontal and vertical correlation lengths and dip (hmax, hmin, β) = (10 m, 3 m, 6.5°).

View larger version (75K): [in this window] [in a new window]   FIG. 12. Eight realizations from the prior probability density function (pdf) using horizontal and vertical correlation lengths of (a) (hmax, hmin) = (2 m, 1 m), (b) (hmax, hmin) = (10 m, 3 m), and (c) (hmax, hmin) = (18 m, 6 m).

View larger version (76K): [in this window] [in a new window]   FIG. 13. Eight realizations from the posterior probability density function (pdf) using horizontal and vertical correlation lengths of (a) (hmax, hmin) = (2 m, 1 m), (b) (hmax, hmin) = (10 m, 3 m), and (c) (hmax, hmin) = (18 m, 6 m).

Inferring the Subsurface Horizontal and Vertical Correlation Lengths
We performed a similar analysis as in the synthetic case, except that we also compared the effect of assuming straight vs. bending rays. Figure 14a shows the HPD region obtained in the case where a pure linear approach was taken, i.e., where the straight-ray approximation was used. Figure 14b show the corresponding HPD region for the linearized case. In both cases, all 702 ray observations were used. The difference between Fig. 14a and 14b is not very pronounced. In both cases, a maximum likelihood for the horizontal range is found between 5 and 10 m. The vertical correlation length has maximum likelihood in the range between 1.5 and 3.0 m, using the linear and linearized approaches, respectively. Figure 15 shows the one-dimensional marginal pdfs corresponding to the HPDs of Fig. 14. Table 2 summarizes the found maximum likelihood model as picked from the two-dimensional HPD plots (Fig. 14) and the one-dimensional marginal pdfs (Fig. 14). The most prominent effect of using a linear sensitivity kernel seems to be that there is no apparent upper limit for the horizontal correlation length. One can only infer that the horizontal correlation length is above 3 to 5 m. The relatively small effect of using a wrong linear sensitivity kernel is probably due to the fact that, for the present case study, the straight-ray approximation is in fact relatively good. Little to no water is expected to be present in the vadose zone. If water had been present, much larger velocity contrast would be present, resulting in more ray bending than is apparent for this case study (see Fig. 11). This is in good agreement with Kowalsky et al. (2005), who concluded that the straight-ray approximation is good for heterogeneous velocity fields, but bad in case a water plume (implying locally large velocity contrast) is present.

View larger version (14K): [in this window] [in a new window]   FIG. 14. The 20% (white), 40%, 60%, and 80% (black) highest probability density (HPD) intervals for the likelihood function using (a) 702 rays and a linear sensitivity kernel, (b) 702 rays and nonlinear sensitivity kernels, and (c) 70 rays and nonlinear sensitivity kernels, as a function of the horizontal (hmax) and vertical (hmin) correlation lengths. The HPD plots are based on approximately 300 sample points.

View larger version (17K): [in this window] [in a new window]   FIG. 15. One-dimensional (1D) marginal probability density function distributions of the likelihood of the horizontal (hmax) and vertical (hmin) correlation lengths using (a) 702 rays and a linear sensitivity kernel, (b) 702 rays and nonlinear sensitivity kernels, and (c) 70 rays and nonlinear sensitivity kernels.

These intervals of high probability obtained as h_max = [5 ...10] m and h_min = [1.5 ... 3] m, as well as the model of maximum likelihood of (h_max,h_min) = (8.5, 2.4) (see Table 2), are thus highly consistent with the information found from the nearby GPR reflection profiles of Fig. 10, which suggests a horizontal correlation length of about 4 to 20 m and a vertical correlation length of about 1 to 5 m.

Inferring the Subsurface Correlation Lengths and Dip
While the average dip of subsurface layers is assumed relatively well known from the nearby GPR reflection profile, we will, for illustration purposes, consider the dip unknown in addition to the variability in correlation length as given above. We make use of the Metropolis–Hastings algorithm to obtain the one-dimensional marginal pdfs of the correlation lengths and the dip as shown in Fig. 16 . As can be seen, the choice of β = 6.5° can be considered a probable parameter of an anisotropic subsurface covariance function, although the most probable set of parameters for the subsurface covariance model is (h_max, h_min, β) = (11 m, 2.5 m, 9°). This result is quite similar to the results obtained when fixing the dip. The variability of the one-dimensional marginal distribution of h_max and h_min of Fig. 16 is remarkably similar to that of Fig. 15, obtained for a fixed dip.

View larger version (13K): [in this window] [in a new window]   FIG. 16. One-dimensional (1D) marginal probability density function (marg. pdf) distributions of the likelihood of the horizontal (hmax) and vertical (hmin) correlation lengths and the dip using 702 rays, based on 175 samples found using the Metropolis–Hastings algorithm.

	Conclusions

TOP ABSTRACT INTRODUCTION Theory--Posterior Validation of... Case Study: Inferring the... Real Data Example Conclusions Appendix REFERENCES

 
We have presented a method to validate the choice of the prior model in inverse problems. Specifically, we have shown how the prior choice of subsurface structural model (i.e., the prior covariance model) associated with least squares based GPR cross-borehole tomography can be validated directly by travel time data observations.

We suggest quantifying the likelihood that a specific choice of prior model is consistent with data by calculating the probability that a sample of the a posteriori pdf is a sample of the a priori pdf. This approach is applicable to any inverse problem where samples can be generated from both the prior and posterior pdfs. In the specific case of linear inverse Gaussian problems, where both the a priori and a posteriori pdfs are multi-Gaussian, we built a model of observed ergodic variability of the mean and covariance of realizations of the a priori pdf, related to a specific choice of prior covariance model, C_Mp. Then we calculated the likelihood that a sample of the a posteriori pdf is a sample of the a priori pdf in terms of reproducing the C_Mp within the model of ergodic variability. This enables one to quantify how consistent a choice of prior covariance model is with data observations. The methodology is relatively easy to apply, although CPU demanding for larger problems; however, using just a subset of the available data results in useful estimates of the subsurface covariance model.

Through a cross-borehole GPR tomography case study, we investigated how well the method works on both synthetic and real data. For the synthetic case, we have shown that the subsurface correlation lengths can be inferred, even when the range is comparable to the width of the model parameter space considered. For the real case, the correlation lengths inferred were highly consistent with information obtained from a nearby reflection GPR profile. For both cases, we found that inferring the dip of the continuity was less trivial, while we found the inferred distribution of the dip to be consistent with the reference model and reflection GPR profile, respectively. We therefore believe that we have presented a novel approach to locate choices of prior covariance models that are consistent with data observations.

The main use of the proposed method with GPR cross-borehole tomography is in conjunction with stochastic simulation, where a number of realizations of the subsurface can be generated with realistic small- and large-scale variability, consistent with observed travel time data. Such realizations can then be transformed into hydrologic parameters with realistic spatial variability.

	Appendix

TOP ABSTRACT INTRODUCTION Theory--Posterior Validation of... Case Study: Inferring the... Real Data Example Conclusions Appendix REFERENCES

 
Anisotropy
The theory presented above can easily be extended to anisotropic semivariogram models. Let the experimental semivariance observed in the primary direction of continuity be given _||(h) and in the secondary direction by ₌(h). The joint multi-Gaussian pdf describing the anisotropic ergodicity in the variance is given by N(_2D, _2D), where _2D = (_||,₌) and where _2D = [_{|||| =||}; _=||' _{= =}].

The extension to three-dimensional anisotropy is straightforward, where the multi-Gaussian pdf is given by N(_3D, _3D), where _3D = [_||,₌,], and _3D = [_{|||| ||= ||}; _{=|| = = =}; ₌ ].

Having set up the covariance matrix, , the likelihood that any sample or model is a sample from the multi-Gaussian pdf with a specific choice of spatial covariance model can now be estimated using Eq. [8].

Joint Probability of Both the Mean and the Covariance
For completeness, the joint probability that a realization is a sample of a random function with observed mean m_exp, and covariance _exp, is but a generalization of Eq. [8]:

[10]

where d_exp = [m_exp, _exp], and _d*=[_m*C_(m*,*); C_(*,m*)*], and C_(m*,*)is the observed covariance between the mean value, m, and mean covariance values, .

	ACKNOWLEDGMENTS

 
Thomas Mejer Hansen is funded by the The Danish Research Council for Technology and Production Sciences, Project 57243. Majken C. Looms is funded by the Faculty of Natural Sciences at the University of Copenhagen. We wish to thank the Water Department, Section for Water Quality at Copenhagen Energy for all their help and for use of their field location at the Arrenæs infiltration plant. All computations have been performed using the freely available mGstat (http://mgstat.sf.net/) and VISIM (http://imgp.gfy.ku.dk/visim.php) programs.

	REFERENCES

TOP ABSTRACT INTRODUCTION Theory--Posterior Validation of... Case Study: Inferring the... Real Data Example Conclusions Appendix REFERENCES

 

• Box, G., and G. Tiao. 1992. Bayesian inference in statistical analysis. John Wiley & Sons, New York.
• Cassiani, G., G. Bohm, A. Vesnaver, and R. Nicolich. 1998. A geostatistical framework for incorporating seismic tomography auxiliary data into hydraulic conductivity estimation. J. Hydrol. 206:58–74.[CrossRef]
• Constable, S.C., R.L. Parker, and C.G. Constable. 1987. Occam's inversion: A practical algorithm for generating smooth models from electromagnetic sounding data. Geophysics 52:289–300.[CrossRef][Web of Science]
• Chiles, J.-P., and P. Delfiner. 1999. Geostatistics: Modeling spatial uncertainty. John Wiley & Sons, New York.
• Day-Lewis, F.D., and J.W. Lane. 2004. Assessing the resolution-dependent utility of tomograms for geostatistics. Geophys. Res. Lett. 31(7):L07503, doi:10.1029/2004GL019617.[CrossRef]
• Day-Lewis, F.D., J.W. Lane, and A.M. Binley. 2005. Applying petrophysical models to radar travel time and electrical resistivity tomograms: Resolution-dependent limitations. J. Geophys. Res. 110(B8):B08206, doi:10.1029/2004JB003569.[CrossRef]
• Emery, X. 2004. Testing the correctness of the sequential algorithm for simulating Gaussian random fields. Stoch. Environ. Res. Risk Assess. 18:401–413.[CrossRef]
• Giroux, G., E. Gloaguen, and M. Chouteau. 2007. bh_tomo: A Matlab borehole georadar 2D tomography package. Comput. Geosci. 33:126–137.[CrossRef]
• Gloaguen, E., D. Marcotte, M. Chouteau, and H. Perroud. 2005. Borehole radar velocity inversion using cokriging and cosimulation. J. Appl. Geophys. 57:242–259.[CrossRef]
• Goovaerts, P. 1997. Geostatistics for natural resources evaluation. Oxford Univ. Press, Oxford, UK.
• Gouveia, W.P., and J.A. Scales. 1997. Resolution of seismic waveform inversion: Bayes versus Occam. Inverse Prob. 13:329–349.
• Hansen, T.M., A.G. Journel, A. Tarantola, and K. Mosegaard. 2006. Linear inverse Gaussian theory and geostatistics. Geophysics 71:R101–R111.[CrossRef]
• Hansen, T.M., and K. Mosegaard. 2007. VISIM: Sequential simulation for linear inverse problems. Comput. Geosci. (in press).
• Hastings, W.K. 1970. Monte Carlo sampling methods using Markov chain and their applications. Biometrika 57:17–34.
• Hubbard, S., Y. Rubin, and E. Majer. 1999. Spatial correlation structure estimation using geophysical and hydrogeological data. Water Resour. Res. 35:1809–1825.[CrossRef]
• Knight, R., P. Tercier, and H. Jol. 1997. The role of ground penetrating radar and geostatistics in reservoir description. Leading Edge 16:1576–1583.[Abstract]
• Kowalsky, M.B., S. Finsterle, J. Peterson, S. Hubbard, Y. Rubin, E. Majer, A. Ward, and G. Gee. 2005. Estimation of field-scale soil hydraulic and dielectric parameters through joint inversion of GPR and hydrological data. Water Resour. Res. 41:W11425, doi:10.1029/2005WR004237.[CrossRef]
• Linde, N., A. Binley, A. Tryggvason, L.B. Pedersen, and A. Revil. 2006. Improved hydrogeophysical characterization using joint inversion of cross-hole electrical resistance and ground-penetrating radar traveltime data. Water Resour. Res. 42:W12404, doi:10.1029/2006WR005131.[CrossRef]
• Looms, M.C., K.H. Jensen, A. Binley, and L. Nielsen. 2007. Monitoring unsaturated flow and transport using cross-borehole geophysical methods. Vadose Zone J. 6:[[xxx–xxx]] (this issue).
• Marchant, B.P., and R.M. Lark. 2004. Estimating variogram uncertainty. Math. Geol. 36:867–898.[CrossRef]
• Marquardt, D.W. 1970. Generalized inverses, ridge regression, biased linear estimation, and non-linear estimation. Technometrics 12:591–612.[CrossRef][Web of Science]
• Maurer, H., K. Holliger, and D.E. Boerner. 1998. Stochastic regularization: Smoothness or similarity? Geophys. Res. Lett. 25:2889–2892.[CrossRef]
• Oldenborger, G., R. Schincariol, and L. Mansinha. 2003. Radar determination of the spatial structure of hydraulic conductivity. Ground Water 41:24–32.[CrossRef][Medline]
• Ortiz, J., and C.V. Deutsch. 2002. Calculation of uncertainty in the variogram. Math. Geol. 34:169–183.[CrossRef]
• Pardo-Iguzquiza, E., and P. Dowd. 2001. Variance–covariance matrix of the experimental variogram: Assessing variogram uncertainty. Math. Geol. 33:397–419.[CrossRef]
• Peterson, J.E. 2001. Pre-inversion corrections and analysis of radar tomographic data. J. Environ. Eng. Geophys. 6:1–18.
• Phillips, W.S., and M.C. Fehler. 1991. Traveltime tomography: A comparison of popular methods. Geophysics 56:1639–1649.[CrossRef]
• Rea, J., and R.J. Knight. 1998. Geostatistical analysis of ground-penetrating radar data: A means of describing spatial variation in the subsurface. Water Resour. Res. 34:329–339.[CrossRef]
• Tarantola, A. 2005. Inverse problem theory and methods for model parameter estimation. SIAM, Philadelphia, PA.
• Tarantola, A., and B. Valette. 1982. Generalized nonlinear inverse problems solved using the least squares criterion. Rev. Geophys. Space Phys. 20:219–232.[CrossRef][Web of Science]
• Tercier, P., R. Knight, and H. Jol. 2000. A comparison of the correlation structure in GPR images of deltaic and barrier spit depositional environments. Geophysics 65:1142–1153.[CrossRef]

This article has been cited by other articles:

				S. Lambot, A. Binley, E. Slob, and S. Hubbard Ground Penetrating Radar in Hydrogeophysics Vadose Zone J., February 25, 2008; 7(1): 137 - 139. [Full Text] [PDF]

				K. S. Cordua, M. C. Looms, and L. Nielsen Accounting for Correlated Data Errors during Inversion of Cross-Borehole Ground Penetrating Radar Data Vadose Zone J., February 25, 2008; 7(1): 263 - 271. [Abstract] [Full Text] [PDF]

This Article Abstract Figures Only Full Text (PDF) Free Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Hansen, T. M. Articles by Nielsen, L. Search for Related Content PubMed Articles by Hansen, T. M. Articles by Nielsen, L. GeoRef GeoRef Citation Agricola Articles by Hansen, T. M. Articles by Nielsen, L. Related Collections Ground Penetrating Radar, GPR Tomography Inverse Procedures/Parameter Estimation