Published online 25 February 2008
Published in Vadose Zone J 7:263-271 (2008)
DOI: 10.2136/vzj2007.0008
© 2008 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
SPECIAL SECTION: GROUND PENETRATING RADAR IN HYDROGEOPHYSICS
Accounting for Correlated Data Errors during Inversion of Cross-Borehole Ground Penetrating Radar Data
Knud S. Cordua,
Majken C. Looms and
Lars Nielsen*
Dep. of Geography and Geology, Univ. of Copenhagen, Øster Voldgade 10, DK-1350 Copenhagen K, Denmark
* Corresponding author (ln{at}geol.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 14 January 2007.
 |
ABSTRACT
|
|---|
Least-squares tomographic inversion of cross-borehole ground penetrating radar (GPR) data was used to estimate the radar wave velocity distribution in the upper 
10 m of the vadose zone at a field site in northern Zealand, Denmark. The radar wave velocities were transformed to values of water saturation and formed a basis for hydrologic studies of flow characteristics of the vadose zone. Cross-borehole GPR data are likely to be contaminated by correlated data errors that may give rise to significant artifacts in the inverse estimate of the velocity distribution if they are not properly accounted for during the inversion process. We analyzed two sources of correlated data errors (unknown cavities and small-scale clay-enriched zones close to the borehole walls), which we assumed to play a significant role in the cross-borehole GPR data sets collected at our field site. The correlated errors may be accounted for by specification of data error covariance matrices, which are included as a priori knowledge in the inverse operator used in the tomographic algorithm. The study indicates that proper accounting for correlated data errors significantly suppresses the effects of these data errors and results in trustworthy inverse estimates of the radar wave velocities between the boreholes. Suppression of the influence of the correlated data errors is a necessity for obtaining realistic hydrologic models of the vadose zone. Using popular inverse methods that disregard the data error correlation properties can lead to severe artifacts in the inverse model estimate and, therefore, can lead to incorrect hydrologic interpretations. Our findings are based on both synthetic tests and a real data example.
Abbreviations: GPR, ground penetrating radar
|
INTRODUCTION
|
|---|
C ross-borehole ground penetrating radar (GPR) methods are increasingly used in integrated geophysical–hydrogeologic studies of the uppermost parts (upper 20 m) of the subsurface. The tomographic GPR images are subsequently used for estimating physical parameter values of the geologic layers between the boreholes (Binley et al., 2001, 2002; Alumbaugh et al., 2002; Cassiani and Binley, 2005; Winship et al., 2006). The primary information, which is most often extracted from the radar wave signals, is the set of first arrival travel times from the transmitter to the receiver antenna positions. From these travel times, radar wave velocities can be estimated by the use of tomographic inversion algorithms. The radar wave velocity field is transformed to values of water saturation of the vadose zone using empirical relations (e.g., Topp et al., 1980). Furthermore, models of radar wave attenuation may be estimated from the cross-borehole GPR data if not only the travel times but also the waveforms of the radar wave signals are included in the inversion (e.g., Peterson, 2001).
Tomographic schemes, originally developed for inversion of seismic data, have been adopted and modified to fit cross-borehole GPR problems (e.g., McMechan, 1983; Hole, 1992; Spakman and Bijward, 2001). The forward algorithms typically rely on some kind of ray approximation of the radar waves. The simplest schemes assume straight ray paths between the source and receiver positions. Straight-ray algorithms give reliable results if the velocity fluctuations of the subsurface are moderate and the distances between sources and receivers are relatively small. If strong velocity fluctuations are expected or the distances between sources and receivers are large, algorithms that take bending of the rays into account will produce more reliable results. Regardless of the choice of ray approximation, the tomographic algorithm often relies on a least-squares inversion scheme.
Least-squares inversion algorithms require input from the user regarding a priori model covariance and data error covariance. Studies have been directed toward choosing proper a priori models and making sound estimates of the covariance of the model parameters, because the inverse estimates obtained during inversion strongly depend on the a priori model properties (e.g., Hansen et al., 2008). Likewise, the assumptions made regarding the data errors and their correlation have significant influence on the inverse estimate (e.g., Tarantola, 1987). Ideally, sources of significant correlated data errors should be corrected before inversion (Peterson, 2001). This is often a difficult task because the exact location of the data errors in the data sets may be difficult or practically impossible to assess, especially if the data sets are large. It is often possible, however, to make sound estimates of the standard deviation and correlation properties of the data errors. These data error properties can be described in the data error covariance matrix of the least-squares inverse operator, and thereby the data errors may be accounted for during inversion (Jackson, 1979; Tarantola, 1987; Nielsen and Jacobsen, 1996).
We used cross-borehole GPR investigations for studying water saturation and flow characteristics of the vadose zone in a sandy environment in northern Zealand, Denmark (Fig. 1
). In our study area, we expect unknown cavities in the borehole walls and small-scale low-velocity anomalies close to the source and receiver positions to give rise to strongly correlated data errors. These types of data errors are often referred to as static errors and they may have significant influence on the quality of cross-borehole tomographic inverse estimates (e.g., Squires et al., 1992; Maurer and Green, 1997; Vasco et al., 1997). In this study, we investigated the influence of such correlated data errors on our GPR data, and we assessed the standard deviation and correlation properties of the errors. We outline below how data error covariance matrices, which are included in the inverse operator, may be specified for the expected correlated data errors. In a synthetic study, we demonstrate how the data errors propagate into the inverse estimate and affect model resolution. Further, we applied the method to a field data example from our study area. We show that the assumptions made about data errors and their correlation during inversion are crucial for the quality of the resulting inverse estimate.
View larger version (26K):
[in this window]
[in a new window]
|
FIG. 1. (a) Location map of the Danish area. The black circle indicates the location of the field site in northern Zealand. (b) Layout of boreholes for ground penetrating radar (GPR) recording and electrical resistivity tomography (ERT) measurements at the field site (see also Looms et al., 2008). Tests made in the field have shown that the electrodes placed in the boreholes for the ERT measurement do not significantly affect the GPR signal recorded in the other boreholes.
|
|
Static-like errors in cross-borehole seismic or GPR data sets are traditionally treated as source- or receiver-dependent extra model parameters, which are estimated from the travel time data during the inversion process (Maurer and Green, 1997; Vasco et al., 1997). We discuss the main differences between our approach and the traditional approaches used for suppressing static-like errors.
|
Site and Data Acquisition Description
|
|---|
Our field site is located on the peninsula of Arrenæs, northern Zealand, Denmark (Fig. 1). Here, the vadose zone is approximately 25 m thick. The sediments of this zone are loosely packed glacial deposits, which consist mainly of sandy or gravely material. Minor intercalations of clay-rich sediments exist in the area. These sand- and gravel-dominated glacial deposit sediments have an average radar wave velocity of 0.13 to 0.14 m/ns (Looms et al., 2008). The boreholes used for GPR investigations have been placed in a regular geometry (Fig. 1b). The boreholes were drilled with a 0.10-m-wide auger drill stem, and 0.05-m-wide polyvinyl chloride (PVC) tubes for GPR antenna access were placed in these boreholes shortly after drilling. The antennae, which were provided by Sensors & Software Inc., Mississauga, ON, Canada, have a nominal frequency of 100 MHz, and are 0.026 m in diameter. The boreholes were drilled to a depth of 14 m. When the PVC tubes were lowered into the boreholes, it was discovered that the lowermost 1 to 2 m of the boreholes were filled with sediments. Thus, a significant amount of the loosely packed sediments must have fallen from the borehole walls to the bottom of the borehole shortly after drilling. The void between the borehole walls and the PVC tubes was backfilled with the material brought to the surface during the drilling process. It is very likely, however, that some significant cavities are still present close to the borehole walls.
We used a ray geometry consisting of a total of 702 rays for our cross-borehole tomographic GPR setup in this study (Fig. 2
). The spacing between neighboring transmitter and receiver positions in the two boreholes is 1.0 and 0.25 m, respectively. All rays dipping >45° were omitted from the data set. The two-dimensional model section between the two boreholes was parameterized using a regular grid of cells, where each cell has constant slowness (inverse velocity). The grid cell length is 0.50 m. For each transmitter–receiver pair, we picked the first arrival travel time to be able to estimate the radar wave velocity by tomographic inversion. In general, the signal/noise ratio of the recorded waveforms was high, and the picking uncertainties were very low and mainly determined by the temporal sampling interval of the data set.
View larger version (68K):
[in this window]
[in a new window]
|
FIG. 2. Ray geometry. A total of 702 data points were used. Rays dipping >45° were omitted from the study. Transmitter positions were 1 m apart, whereas receiver positions were 0.25 m apart.
|
|
|
Correlated Data Errors and Covariance Specification
|
|---|
Several sources of data error exist in cross-borehole GPR tomographic problems. Some examples of important error sources are (see Peterson [2001] for an exhaustive analysis of different error types that may influence cross-borehole GPR data sets):- (i) unknown cavities in the borehole walls close to the receiver or transmitter antenna,
- (ii) unknown small-scale anomalies close to the antenna positions,
- (iii) incorrect positioning of the transmitter or receiver antennae during data acquisition,
- (iv) time jumps or time drift due to problems with transmitter calibration,
- (v) incorrect information about the geometry of the two boreholes,
- (vi) mispicking of the travel time of the first arrival,
- (vii) electrical background noise from the surroundings or the GPR instrumentation.
Data error sources (i) through (v) will result in correlated data errors. Only error sources (vi) and (vii) may be expected to result in uncorrelated data errors. The unknown cavities and small-scale anomalies close to the borehole walls are treated as sources of data error because they are believed to be smaller than the grid cell size, and, therefore, they cannot be captured by the model grid. Furthermore, they are positioned in areas where the spatial resolution provided by the cross-borehole data is low due to a limited number of crossing ray paths. The influence of the different error types will vary from data set to data set, since they strongly depend on the local geology, instrumentation, and background noise level.
In the following, we discuss why we expect errors sources (i) and (ii) to significantly contaminate our cross-borehole GPR data sets, and we show how these error types may be described through data error covariance specification. As described above, observations made during drilling of the boreholes strongly indicated that material had fallen from the borehole walls to the bottom of the boreholes shortly after the drilling process had ended, and small-scale clay-rich anomalies are known from other investigations to exist in the study area. The cavities in the borehole walls will generate correlated errors for data collected using receiver or transmitter positions close to the cavities. In the vadose zone, the cavities will primarily be filled with air, which has a radar wave velocity of 0.3 m/ns, i.e., approximately twice the average velocity of the sediment that is studied here. If the radius of the unknown cavities is on the order of 0.05 to 0.25 m, then the unknown cavities may cause all travel times collected close to the cavities to be 0.5 to 1.5 ns early. Furthermore, we know from wells drilled in the area that minor intercalations of clay-rich material may be found in the overall sandy–gravelly sediments of the study area (well information from the area is available in the Jupiter borehole database at www.geus.dk/jupiter/index-dk.htm [in Danish, verified 13 Sept. 2007]). The clay-enriched material is expected to have a radar wave velocity of about 0.06 m/ns, which is about half the velocity of the overall sedimentary material under study. Thus, even small (0.05–0.25 m in diameter), unknown anomalies enriched in clay may cause all the data collected using a specific receiver or transmitter position to be about 0.5 to 2 ns late if these anomalies are located close to the receiver or transmitter antenna position. These anomalies are too small to be captured by the model parameterization. Therefore, we choose to regard the influence of such small-scale anomalies as errors of the data set.
A priori information about data errors and data error correlation is specified in the data error covariance matrix, CD, which is part of the least-squares inverse estimate of the model parameters (e.g., Tarantola, 1987):
 | [1] |
where mprior contains the a priori values of the slowness distribution of the grid of cells describing the two-dimensional section between the two boreholes, CM is the model covariance matrix, dobs is the observed travel time data set, and G is the linear mapping matrix expressing the relation between model parameters and observed data. Matrix G contains the distances that the different rays have traveled in the individual cells.
The data error covariance matrix is a square and invertible matrix consisting of N by N elements, where N is the total number of picked travel times. For more details regarding the characteristics of data error covariance matrices, we refer to Tarantola (1987). The choice of CD has significant impact on the inverse estimate. The diagonal elements of CD contain the variance of the individual data errors, whereas the off-diagonal elements describe the correlation between the individual data errors. In the special case where all data errors are assumed to be uncorrelated and have the same standard deviation, the data error covariance matrix takes on the familiar form
 | [2] |
where I is the identity matrix, and 
is the standard deviation of the data errors.
Error types I and II are spatially correlated, and the spatial correlation properties must be specified in the off-diagonal elements of the data error covariance matrix. Clearly, the unknown cavities and low-velocity anomalies strongly affect data collected with a transmitter or receiver position close to the location of a cavity or a low-velocity anomaly. We specify separate data error covariance matrices related to the transmitter and receiver positions (CTr and CRe, respectively). For simplicity, we refer to the corresponding correlated data errors as the transmitter static errors and the receiver static errors, respectively, in the following. The term static is borrowed from the reflection seismic literature, where this term is used for describing seismic travel time data errors, which are (almost) the same for all seismic data recorded at a specific location. We choose a covariance function to describe the different elements of CTr and CRe (e.g., Walden and Hosken, 1985). For CTr we have
 | [3] |
where s(i,j) is the distance between the ith and the jth transmitter position, and L is the spatial correlation length. Matrix CRe is specified in a similar way.
The covariance matrix for uncorrelated data errors related to mispicking of the first arrival may be set up using Eq. [2]. As mentioned above, the picking uncertainty is low for our data set, which offers a high signal/noise ratio for the recorded waveforms. In the synthetic study below, where we focused on the effects of expected correlated error types, we chose a low standard deviation of 0.1 ns for the uncorrelated data errors.
|
Simulation of Correlated Data Errors
|
|---|
In a synthetic test study, we used a spatial correlation length, L, of 5 m and a standard deviation, , of the correlated data errors of 0.5 ns for setting up CTr and CRe using Eq. [2]. Matrices CTr and CRe are illustrated in Fig. 3a and 3c
, respectively. For simplicity, only 200 measurements were used to create these figures. Note the blocky appearance of CTr. Each of the white blocks following the diagonal of CTr describe the perfect correlation of the data errors of all the travel times coming from a specific transmitter position. The gray tones of the other blocks show the amount of correlation of the data errors originating from different transmitter points. Matrix CRe also has a blocky appearance but shows different correlation properties within the individual blocks due to the way the travel-time data set is organized: The transmitter antenna is lowered to a certain position in the borehole. From that position, radar signals are transmitted to all the different receiver antenna locations in the other borehole. Then the transmitter is lowered 1.00 m to the next transmitter position, and signals are again transmitted to all receiver locations in the other borehole. The receiver positions are only 0.25 m apart. This procedure is continued until all transmitter locations in both boreholes have been visited (Fig. 2). The data error covariance matrix related to the uncorrelated data errors (standard deviation of 0.1 ns) is specified using Eq. [2] and illustrated in Fig. 3e.
View larger version (47K):
[in this window]
[in a new window]
|
FIG. 3. Data error covariance matrices and corresponding simulated data noise: (a–b) spatially correlated transmitter static errors; (c–d) spatially correlated receiver static errors; (e–f) uncorrelated errors related to mispicking of first arrivals. In (a–d), the spatial correlation length, L, is 5 m. The data error covariance matrices were set up using Eq. [3], and the different types of noise were simulated as outlined by Jacobsen (1993). The total data set consists of 702 travel times (Fig. 2). For illustration purposes, the data errors are shown for only 200 data points.
|
|
Based on the specified data error covariance matrices, we may simulate the corresponding data noise following the scheme outlined by Jacobsen (1993). Figure 3 shows simulations of transmitter static errors, receiver static errors, and the uncorrelated data errors. The strongly correlated nature of the static errors is evident. The transmitter static errors have a step-like appearance. In contrast, neighboring receiver static errors are typically different; however, the receiver static errors are repeated in a cyclic manner through the data set. These appearances of the static error types are a consequence of the data sampling scheme outlined above.
|
Inversion of Synthetic Data with Correlated Errors
|
|---|
In a synthetic test, we assessed the influence of the correlated errors on the least-squares tomographic inverse estimate of the slowness distribution between a set of boreholes placed 5 m apart. Using the ray geometry of Fig. 2 and an anomaly test pattern with a known slowness perturbation in the central part of the model grid (Fig. 4d
, 4h, and 4l), we generated a noise-free first arrival travel-time data set. We then set up the data error covariance matrices for the transmitter and receiver static errors as well as the uncorrelated mispicking errors using Fig. 3. Subsequently, we added the three data error covariance matrices and created the total error budget, which we then added to the noise-free synthetic travel-time data set to obtain a synthetic data set including correlated data errors.
View larger version (56K):
[in this window]
[in a new window]
|
FIG. 4. Synthetic test study: (a–c), (e–g), and (i–k) inversion results; (d, h, and l) true slowness test pattern that should ideally be fully recovered by inversion. In (a–c), data error correlation was not taken into account during the inversion process. In (e–g) and (i–k), the inversion was performed assuming that the data errors were correlated; different assumptions about the correlation properties and the standard deviations of the correlated data errors were made. Only in (j) was the inversion performed using the correct data error covariance matrix of the total error budget.
|
|
The synthetic data set contaminated by the correlated transmitter and receiver static errors as well as the uncorrelated mispicking errors was inverted using Eq. [1]. We used an a priori model in which all cells were assumed to have a slowness equal to the background slowness of the true anomaly test pattern. We performed nine different inversions that differed only in the choice of a priori data error covariance matrix.
First, we performed the inversion disregarding the two correlated data error types, i.e., we used an a priori data error covariance matrix that was a diagonal matrix containing only the variance of the uncorrelated data error (Fig. 3e). The inversion result is shown in Fig. 4a. The test anomaly in the central part of the model grid is indicated in the inverse estimate, although it is smeared and its shape is distorted (compare Fig. 4a and 4d). The correlated data errors have a significant effect on the inverse estimate, and they give rise to pronounced artifacts in the estimated velocity model, especially, but not only, close to the boreholes. Such effects of static-like errors are expected and have also been noted for seismic cross-borehole investigations (e.g., Squires et al., 1992). When severe errors are assumed to be present in a given data set used for tomographic inversion, the standard procedure is often to simply increase the standard deviation of the a priori data uncertainty above the level of the expected picking uncertainty to damp the least-squares solution. Therefore, we again performed the inversion just accounting for uncorrelated data errors, but this time we increased the standard deviation of the a priori data uncertainties from 0.1 to 0.7 ns to account for the total amplitude of the composite error budget (Fig. 4b). The artifacts related to the correlated data errors are damped but still evident. The estimation of the central test anomaly has suffered from additional smoothing. The approach where we only account for uncorrelated data errors can partly suppress the unfortunate effects of the correlated data errors if the a priori standard deviation is increased to 1.4 ns, which is twice the standard deviation of the total error budget of the synthetic data set (Fig. 4c). In this case, however, the resolution of the central test anomaly is dramatically reduced due to smearing caused by the strong damping of the least squares solution.
We then performed the inversion, assuming that a part of the data errors are correlated, and we tested how the inversion scheme performs on the noisy synthetic data for different choices of the correlation length, L, and the standard deviation of the correlated errors. In all the remaining synthetic examples, we correctly assumed that the standard deviation of the uncorrelated mispicking errors was 0.1 ns. First, we correctly assumed that L was 5 m, but we incorrectly assumed that the standard deviations of the transmitter statics and the receiver statics were only 0.1 ns (only one-fifth of the actual standard deviation of the synthetic transmitter and receiver statics). The inverse estimate is shown in Fig. 4e. Although not fully recovered and slightly distorted, the test anomaly in the central part of the model grid is better recovered than in any of the three cases in which data error correlation was disregarded. Significant undesired artifacts are, however, present even in the central parts of the model grid. The inverse estimate of Fig. 4f was found assuming that L = 5 m and that the standard deviation of both the transmitter and receiver static errors was 0.25 ns (half that of the actual synthetic static errors). In this case, the resolution of the central test anomaly is significantly improved and the amplitude of the artifacts is reduced compared with the previous inversion result (compare Fig. 4f to Fig. 4e). The main artifacts are now observed along the top and bottom of the model grid where the resolution of the data set is poor due to the limited ray coverage in these areas (cf. Fig. 2). For the estimate of Fig. 4g, we again made the correct assumption that L is 5 m, but we overestimated the standard deviations of the correlated static errors by a factor of 4 (2 ns). The central test anomaly is well resolved, and the main artifacts observed in the inverse estimate are limited to the areas where the ray coverage was poor. For the inverse estimates in Fig. 4i to 4k, we made the correct assumptions regarding the standard deviation of the different error types. These three estimates only differ in the choice of L. From Fig. 4i to 4k, L increases in steps: 0, 5, and 10 m. All three inversion results provide good resolution of the central test anomaly, although the estimate shown in Fig. 4k (L = 10 m) exhibits slightly more smearing than the estimates of Fig. 4i and 4j; however, the differences are minor. Significant artifacts of the inverse estimates are, again, limited to the regions of the model with the poorest ray coverage (top and bottom).
Taking the correlation of the data errors into account resulted in significantly better recovery of the central test anomaly than the popular inversion approach in which correlation of the data errors is disregarded. The synthetic test examples indicate that the approach of taking data correlation into account is robust. Even if the standard deviation of the correlated data errors is underestimated by a factor of 2 or overestimated by a factor of 4, the inversion scheme provides good resolution of the central test anomaly and the influence of artifacts is limited. Over- or underestimating the value of the correlation length, L, by 100% results in essentially the same inversion results as the case where the correct choice of L is made.
|
Tomographic Inversion of Field Data
|
|---|
Tomographic data sets were obtained using the GPR boreholes at the Arrenæs field site (Fig. 1b). We inverted and interpreted first arrival travel-time data acquired between Boreholes 1 and 3 using the source–receiver geometry illustrated in Fig. 2. Looms et al. (2007) described several data sets obtained at the Arrenæs field site and discussed the hydrogeologic implications of the geophysical interpretations. We inverted the field data, assuming different standard deviations and correlation properties of the data errors (Fig. 5
). The inverse estimates that we found and discuss below only differ in the choices we made regarding a priori data error correlation and standard deviation.
View larger version (37K):
[in this window]
[in a new window]
|
FIG. 5. Results from inversion of the data set acquired between GPR Boreholes 1 and 3 at the Arrenæs field site (Fig. 1b). Borehole 1 is to the left and Borehole 3 to the right. In (a–c), the data errors were assumed to be uncorrelated. In (d–h), the data errors were assumed to be correlated; different assumptions about the correlation properties and the standard deviations of the correlated data errors were made.
|
|
First, we inverted the tomographic data using the popular approach where data error correlation is not taken into account (Fig. 5a–5c). We used an a priori data error standard deviation of 0.8 ns to obtain the result of Fig. 5a. Overall, the velocity model of Fig. 5a is characterized by relatively high GPR wave velocities above 7.5-m depth and lower velocities below this depth. In the uppermost 7.5 m, a high-velocity anomaly is present in the central part of the model grid. Marked high-velocity anomalies are found at 2- to 3-m depth near both borehole walls. Below 3- to 4-m depth, low-velocity anomalies are found along both boreholes. In the area of lower velocities below 7.5-m depth, a pronounced low-velocity anomaly is found just below 8-m depth about 1 m to the right of the leftmost borehole. Thus, marked two-dimensional variations are observed in the velocity field. The high-velocity areas could be interpreted as sandy or gravelly sediments with low water saturation, whereas the low-velocity areas could be interpreted as areas with an increased content of fine-grained material or increased water content (e.g., Reynolds, 1997). The strong lateral variations in the velocity field are inconsistent, however, with results from interpretation of GPR reflection data sections, which show that, overall, the geologic layering of the subsurface is horizontal to subhorizontal, with average dips of 6° and maximum dips of 10° (see Hansen et al. [2008] for illustration of the GPR reflection data). We suspect that the observed strong lateral velocity variations may be caused by unknown cavities in the borehole walls or local, small-scale, low-velocity anomalies near the boreholes, as discussed above. We sought inverse estimates where these unrealistic anomalies are suppressed. The inverse estimate illustrated in Fig. 5b was found assuming uncorrelated data errors with a standard deviation of 1.9 ns. The two high-velocity anomalies near the boreholes are now partly suppressed and the central high-velocity anomaly above 7-m depth is smeared. Relatively strong lateral variations are still present in the model, however, although all velocity structures appear more smoothly smeared than the model of Fig. 5a. Increasing the standard deviation of the assumed uncorrelated data errors to 3.8 ns resulted in the inverse image of the velocity field shown in Fig. 5c. Here, all local anomalies are significantly suppressed, and, in general, the velocity field is smooth, with more limited lateral velocity variations.
To suppress the unrealistic velocity anomalies in the inverse estimates while keeping a high degree of model resolution with limited smearing in the central parts of the model grid, we inverted the travel-time data while accounting for correlated, static-like errors through a priori data error covariance specification as outlined in the synthetic studies above. Inversion results obtained assuming different characteristics of the static-like data errors are shown in Fig. 5d to 5h. In these five inversions, we also assumed uncorrelated data errors with a standard deviation of 0.8 ns. Thus, we assumed all three data error types illustrated in Fig. 3 to be present in the field data. For the three inverse estimates of Fig. 5d to 5f, we used L = 5 m, whereas we made different choices regarding the standard deviation of the correlated data errors. In Fig. 5d, we assumed a standard deviation of 0.5 ns for the two correlated data error types. Although partly suppressed, the unrealistic anomalies near the boreholes and the high-velocity anomaly in the central part of the model grid are still evident in the tomographic image. Increasing the a priori standard deviation of the assumed correlated data errors to 1.2 ns (Fig. 5e) and 2.0 ns (Fig. 5f) resulted in velocity models showing less overall lateral variation. Note that in Fig. 5f the anomalies near the boreholes have been effectively suppressed. Furthermore, the high velocities above 7.5-m depth are not centered around the middle of the model grid. Instead, we observe a more laterally continuous high-velocity body terminating at about 7.5-m depth, consistent with the observations made in the reflection GPR data. Note that the low-velocity anomaly just below 8-m depth seems to be well constrained by independent crossing ray paths and has not suffered smearing although the standard deviation of the total error budget has been significantly increased. This anomaly could represent a small anomaly enriched in fine-grained material (clay and silt). Setting L to 0 m while keeping the standard deviation of the correlated data errors at 2.0 ns resulted in the inverse estimate of Fig. 5g, whereas increasing L to 10 m (still keeping the correlated data error standard deviation at 2.0 ns) resulted in the velocity distribution shown in Fig. 5h. The results of Fig. 5g and 5h show much resemblance to the results of Fig. 5e and 5f, although slightly more lateral variation is observed in the velocity field of Fig. 5g. Thus, in this real data example also, our approach seems to be robust: a broad selection of a priori standard deviations and correlation lengths of the correlated data errors resulted in inverse estimates in which the general features are almost identical.
|
Discussion
|
|---|
Several sources of data error exist in cross-borehole GPR studies, and some error types are strongly correlated. Based on studies related to our field site in northern Zealand, Denmark, we investigated the influence of selected sources of correlated errors (unknown irregularities of the borehole walls and small-scale anomalies close to the borehole walls), which are believed to have a significant effect on our cross-borehole GPR data. We found that these correlated error types had a significant negative effect on the inverse estimate of the radar wave velocity distribution if not properly accounted for by careful a priori data error covariance specification. They introduced pronounced positive and negative velocity anomalies not only close to the boreholes but also in other parts of the model grid if a standard least-squares operator, which does not take data error correlation into account, was used. The synthetic studies show that if we account for data error correlation, significantly improved inverse estimates may be obtained even if the correlation properties are only approximately known. The inversion approach used to suppress the effects of the correlated data errors is robust: a wide range of a priori standard deviations and correlation lengths of the correlated data errors resulted in inverse estimates in which the recovered test anomaly patterns are similar. In the field data example, strong isolated high-velocity anomalies are present in the inverse estimates if data error correlation is not accounted for. These high-velocity anomalies could well represent cavities in the borehole walls, especially since they occur in the uppermost 2 to 3 m of the boreholes, where we expect the sediments to be most loosely packed due to lack of compaction. The undesired high-velocity anomalies can only be suppressed by the popular approach where error correlation is disregarded if very high a priori data uncertainties are assumed, in which case the entire model is significantly smoothed and the velocity structures appear to be poorly resolved. Furthermore, the popular approach results in relatively strong lateral variations in the main velocity anomalies of the inverse estimate. This observation is in conflict with observations made in reflection GPR data sections from the field site (Hansen et al., 2008). Consistent with our findings, Squires et al. (1992) showed that static-like errors introduce unrealistic, strong lateral variations in the inverse estimates of seismic wave velocities in a seismic cross-borehole study. Accounting for data error correlation suppresses the effects of the anomalies close to the borehole walls and reduces the unrealistically strong lateral variation of the main velocity anomalies, consistent with observations made in the reflection GPR data sections from the field site. As for the synthetic test examples, the inversion scheme accounting for data error correlation appears to be robust.
Other inversion schemes have been used for suppression of static-like errors. A common feature of these schemes is that the sources of the correlated data errors are included as extra model parameters that have to be estimated during the inversion process. For example, Crosson (1976) used a source-relocation technique to account for incorrect information about the position of earthquake sources in large-scale seismic investigations. Maurer and Green (1997) used a similar technique to suppress the systematic data errors resulting from incorrect knowledge about the source and receiver positions in a high-resolution cross-well seismic tomographic experiment. Vasco et al. (1997) used an inversion scheme where source and receiver static terms were estimated as extra model parameters in a cross-borehole tomographic GPR study. Although GPR tomographic problems normally involve more data points than model parameters, the inverse problem cannot be solved without using some kind of regularization technique. Vasco et al. (1997) minimized the roughness of the velocity model during inversion using a finite-difference approximation of the spatial derivatives of the velocity field, and they did not allow for extremely large-amplitude velocity variations and unrealistically large source and receiver static terms. Thus, also in the approach of Vasco et al. (1997), prior assumptions were made regarding the amplitude of the static-like errors. In our study, we chose to account for the expected data errors through specification of data error covariance matrices, which are included in the inverse operator. With this approach, we do not need to increase the number of model parameters that must be estimated from the same number of travel times, and we can specify geologically reasonable standard deviations and correlation properties of the data errors. Furthermore, our approach is not limited to accounting for static-like features in the data sets. The scheme is flexible and any error type for which a covariance matrix may be specified can be accounted for using this approach.
We have focused on the correlated data errors generated by unknown cavities or small-scale anomalies near the boreholes, which we assume may play an important role at our field site. As discussed above, several other sources of correlated data errors occur. We have not identified any systematic data errors that may unequivocally be linked to erroneous positioning of the source and receiver antennae. The effects generated by such errors may, however, at least partly be accounted for by the covariance matrices that we set up in this study. Before, during, and after each cross-borehole survey, we performed transmitter and receiver calibration (adjustment of time zero). We did not note any drift or time jumps, which could cause significant miscalibration. Incorrect knowledge about borehole geometry will also produce significant systematic errors (Peterson, 2001). We have measured the distance and topographical differences between the boreholes using a high-precision theodolite instrument, resulting in uncertainties of <0.001 m. In this study, we assumed that our uncertainties linked to borehole geometry were insignificant.
The error types that we have studied in detail here may not have a significant impact on cross-borehole GPR data sets collected in other environments. For example, unknown cavities around the borehole walls may be less likely to arise in more compact sediments such as limestone or sandstone. Moreover, the sources of error should be carefully analyzed for each separate data set.
Peterson (2001) discussed a number of systematic (or correlated) errors likely to contaminate cross-borehole GPR data sets and outlined how they may be corrected. Correction of known errors should, of course, always be made before inversion. Sometimes purely static errors can be identified by visual inspection of travel times or travel-time residuals (e.g., Squires et al., 1992). Manual inspection and correction of all travel-time picks may be not be a feasible approach, however, especially for large data sets containing 700 travel-time picks or more, and the risk of missing significant errors during the manual correction process is large. We have demonstrated that if the variance and correlation properties of correlated data errors can be approximately estimated, then the effect on the inverse estimate of these correlated data errors may be greatly damped through careful specification of a priori data error covariance matrices.
In the synthetic test examples shown in this study, it is evident that the inversion results obtained accounting for data error correlation provide better resolution of the test anomaly pattern. In general, tomographic algorithms are known to produce severe smearing of anomalies (e.g., Hole, 1992). In particular, in areas with uneven ray coverage (few rays in one direction and many rays in a different direction) standard inversion algorithms, which disregard spatial correlation, tend to cause smearing of anomalies in the direction of the many rays. Nielsen and Jacobsen (1996) demonstrated that taking spatial correlation into account suppresses the smearing effects and increases model resolution.
Most program packages used for cross-borehole GPR tomographic inversion do not offer the possibility of taking correlated data errors into account by flexible specification of the data error correlation properties in the a priori data error covariance matrix, which is part of the least-squares inverse operator. Here we have demonstrated how correlated data errors can be taken into account through flexible data error covariance specification. The approach is simple in the way that it only requires basic reasoning regarding data errors likely to occur in a given data set. As soon as the interpreter has got an estimate of the standard deviation of the data errors and their spatial (or temporal) correlation properties, the data error covariance matrix can be set up. The incorporation of the data error covariance matrix in the least squares operator is straightforward (Eq. [3]; Tarantola, 1987).
|
Conclusions
|
|---|
Cross-borehole GPR data used for studying water saturation and hydrologic properties of the vadose zone are likely to be contaminated by correlated data errors. Some of the correlated data errors may be difficult or practically unfeasible to correct before inversion of the GPR data. Standard software for interpretation of cross-borehole GPR data do not normally offer the possibility of taking correlated data errors into account during the tomographic least-squares inversion process. Disregarding the correlation properties during inversion may lead to severe artifacts in the inverse model estimate and, therefore, lead to incorrect results regarding water saturation and hydrologic properties.
At our field site in northern Zealand, Denmark, we expected that unknown cavities and unknown small-scale low-velocity anomalies existed at or close to the walls of the boreholes used for the GPR data collection. These two sources of error may significantly affect the recorded first arrival travel times of the cross-borehole GPR data sets, and the resulting travel-time errors have a strongly correlated nature. Through analysis of the correlation properties of these error types, we have specified data error covariance matrices that describe the main characteristics of the errors. In synthetic inversion examples, we have shown that disregarding these correlated error types during inversion leads to significant contamination of the inverse estimate by undesired anomalies close to the borehole walls. The influence of the correlated data errors was effectively damped and the resolution properties of the inverse estimate were significantly improved if data error correlation was accounted for during the inversion process. The proposed method is robust in the sense that even significant overestimation of the data error standard deviation or incorrect assumptions about the spatial correlation length of the data errors results in inverse estimates that are similar to the ideal case where the exact correct standard deviation and spatial correlation length are known.
The proposed inversion scheme was applied to a tomographic GPR data set acquired at the field site at Arrenæs. The inverse estimates so obtained show more geologically plausible velocity structures and significant suppression of spurious high-amplitude anomalies near the borehole walls than inversion results calculated using standard inversion schemes disregarding data error correlation.
|
ACKNOWLEDGMENTS
|
|---|
We wish to thank the Water Department, Section for Water Quality, at Copenhagen Energy for all their help and for the use of their field location at Arrenæs infiltration plant. The cross-borehole instrumentation was purchased at Sensors & Software Inc., Mississauga, ON, Canada.
|
REFERENCES
|
|---|
- Alumbaugh, D., P.Y. Chang, L. Paprocki, J.R. Brainard, R.J. Glass, and C.A. Rautman. 2002. Estimating moisture contents in the vadose zone using cross-borehole ground penetrating radar: A study of accuracy and repeatability. Water Resour. Res. 38(12):1309, doi:10.1029/2001WR000754.[CrossRef]
- Binley, A., G. Cassiani, R. Middleton, and P. Winship. 2002. Vadose zone flow model parameterisation using cross-borehole radar and resistivity imaging. J. Hydrol. 267:147–159.[CrossRef]
- Binley, A., P. Winship, R. Middleton, M. Pokar, and J. West. 2001. High-resolution characterization of vadose zone dynamics using cross-borehole radar. Water Resour. Res. 37:2639–2652.[CrossRef]
- Cassiani, G., and A. Binley. 2005. Modeling unsaturated flow in a layered formation under quasi-steady state conditions using geophysical data constraints. Adv. Water Resour. 28:467–477.[Medline]
- Crosson, R.S. 1976. Crustal structure modeling of earthquake data: 1. Simultaneous least squares estimation of hypocenter and velocity parameters. J. Geophys. Res. 81:3036–3046.
- Hansen, T.M., M.C. Looms, and L. Nielsen. 2008. Inferring the subsurface structural covariance model using cross-borehole ground penetrating radar tomography. Vadose Zone J. 7:249–262 (this issue).[Abstract/Free Full Text]
- Hole, J.A. 1992. Nonlinear high-resolution three-dimensional seismic travel time tomography. J. Geophys. Res. 97:6553–6562.[CrossRef]
- Jackson, D.D. 1979. The use of a priori data to resolve non-uniqueness in linear inversion. Geophys. J. R. Astron. Soc. 57:137–157.
- Jacobsen, B.H. 1993. Practical methods of a priori covariance specification for geophysical inversion. p. 1–10. In K. Mosegaard (ed.) Proc. Interdisciplinary Inversion Workshop 2, Copenhagen. 19 May 1993. Niels Bohr Inst., Univ. of Copenhagen, Copenhagen.
- Looms, M.C., K.H. Jensen, A. Binley, and L. Nielsen. 2008. Monitoring unsaturated flow and transport using cross-borehole geophysical methods. Vadose Zone J. 7:238–248 (this issue).[Abstract/Free Full Text]
- Maurer, H., and A.G. Green. 1997. Potential coordinate mislocations in crosshole tomography: Results from the Grimsel test site, Switzerland. Geophysics 62:1696–1709.[CrossRef][Web of Science]
- McMechan, G.A. 1983. Seismic tomography in boreholes. Geophys. J. R. Astron. Soc. 74:601–612.
- Nielsen, L., and B.H. Jacobsen. 1996. Resolution and error propagation analysis for tomographic data with correlated errors. In B.H. Jacobsen et al. (ed.) Inverse methods: Interdisciplinary elements, methodology, computation, and applications. Lecture Notes in Earth Sci. 63. Springer-Verlag, Berlin.
- Peterson, J.E. 2001. Pre-inversion corrections and analysis of radar tomographic data. J. Environ. Eng. Geophys. 6:1–18.
- Reynolds, J.M. 1997. An introduction to applied and environmental geophysics. John Wiley & Sons, New York.
- Spakman, W., and H. Bijward. 2001. Optimization of cell parameterization for tomographic inverse problems. Pure Appl. Geophys. 158:1401–1423.[CrossRef]
- Squires, L.J., S.N. Blakeslee, and P.L. Stoffa. 1992. The effects of statics on tomographic velocity reconstructions. Geophysics 57:353–362.[CrossRef]
- Tarantola, A. 1987. Inverse problem theory: Methods for data fitting and model parameter estimation. Elsevier Sci., New York.
- Topp, G.C., J.L. Davis, and A.P. Annan. 1980. Electromagnetic determination of soil water content: Measurements in coaxial transmission lines. Water Resour. Res. 16:574–582.[CrossRef]
- Vasco, D.W., J.E. Peterson, and K.H. Lee. 1997. Ground-penetrating radar velocity tomography in heterogeneous and anisotropic media. Geophysics 62:1758–1773.[CrossRef][Web of Science]
- Walden, A.T., and J.W.J. Hosken. 1985. An investigation of the spectral properties of reflection coefficients. Geophys. Prospect. 33:400–435.[CrossRef]
- Winship, P., A. Binley, and D. Gomez. 2006. Flow and transport in the unsaturated Sherwood sandstone: Characterization using cross-borehole geophysical methods. In R.D. Barker and J.H. Tellam (ed.) Fluid flow and solute movement in sandstones: The onshore UK Permo-Triassic red bed sequence. Geol. Soc. Spec. Publ. 263. Geol. Soc., London.
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]
|
|
|
TOP
ABSTRACT
INTRODUCTION
