Underground Cavity Detection Through Spectral Distortion of a GPR Signal

We present a novel method to detect underground cavities using cross-borehole ground-penetrating radar (GPR). We model the propagation of a GPR signal across a cavity and find that the spectrum of the signal will be distorted in a specific low-pass manner. Comparing the spectra of received signals with the predicted spectrum using a hypothesis-testing approach allows us to detect cavities with high probability. The proposed method is validated on both simulated and real measurements. Our approach generally remains effective even when conventional methods fail; our method can also be combined with conventional methods in order to more successfully distinguish between cavities and geologic clutter.


I. INTRODUCTION AND RELATED WORK
W E consider the problem of using cross-borehole ground-penetrating radar (GPR) to detect underground cavities. GPR has been applied to this problem since at least the 1980s [1], and methods exist to detect cavities and other anomalies using both zero-offset profiling (ZOP) and multiple-offset gathers (MOG). Most cavity-detection methods are based around traveltime (e.g. [1]- [6]); since an electromagnetic (EM) wave will propagate more quickly in air than in soil, when a signal arrives at the receiver borehole earlier it may indicate that it propagated through a cavity. Imaging using full-waveform inversions (FWI) [7]- [11] is also increasingly used for cavity-detection applications. Other useful features include received energy [1], [12], [13] and time-frequency methods [14]- [16].
An interesting feature for cross-borehole GPR (not applied specifically to the cavity-detection problem) is offered by [17]. They would have liked to use the received energy as a feature for tomographic inversion. Unfortunately, they found that for a variety of reasons, direct measurements of the received energy are often unacceptably noisy (we have also found this to be true in our own measurements). To this end, they measure the centroid frequency of the received signal and use this quantity as a proxy for the attenuation of the signal.
We model the propagation of an EM wave through a cavity and show that this will tend to have a low-pass effect on the spectrum of the signal. This spectral distortion can be used as a feature for cavity detection. Indeed, we show that cavities can sometimes be detected using only the centroid frequency of [17]. However, having a model of what the spectrum of the received signal should look like both in the case in which the signal passed through a cavity and in the case in which it did The authors are with Tel-Aviv University, School of Electrical Engineering. not pass through a cavity, we can adopt a hypothesis-testing approach. As was found by [18], adopting a hypothesis-testing approach dramatically improves the efficacy of our method.
In this paper we restrict our treatment to cavities which are significantly extended in the direction perpendicular to the line connecting the two boreholes and which have a width of around a meter or two. These are the types of cavities which have been considered in [1]- [3], [13], [15], [16], [18], [19], among many others. A loose three-view drawing of such a cavity can be found as supplemental material to this paper. Our treatment in this paper may not be valid for cavities which are significantly smaller or larger than this -for small fractures, we expect that diffraction effects, which we neglect, might dominate, while for large caverns, the dominant effect might not be the propagation through the cavity upon which we base our model, but rather the absence of propagation in soil. A treatment of these types of cavities is beyond the scope of this paper.
In the following section we model the propagation of an EM wave across a cavity. We subsequently present our method of cavity detection and validate it using both simulated and real GPR measurements. Finally, we conclude and offer directions for future research.

II. MODELING PROPAGATION ACROSS A CAVITY
Consider a plane wave propagating between two boreholes as shown in Figure 1. The wave propagates a distance s 1 through soil, across a cavity of width a 1 containing only free space, and a further distance s 2 through soil before arriving at the receiver borehole. At each interface a portion of the wave is transmitted and a portion is reflected. It is well known [20] that a plane wave with initial magnitude E 0 , after propagating a distance h through a medium, is given by where is the propagation constant, α is the attenuation constant, and β is the phase constant. Proceeding similarly to [21], we write equations for the waves shown in Figure 1, where the wave is measured at Fig. 1. A vertically-polarized (parallel to the interfaces) EM plane wave propagates across a cavity between two boreholes. The wave is reflected at every interface. the tail of the arrow (the wave w 0 is ignored as it does not affect the field measured at the receiver borehole). Letting E rc denote the wave at the receiver borehole after propagation through a cavity, k = ω c be the wavenumber of the wave in free space, t sa and t as be the coefficients of transmission from soil to air and from air to soil respectively, and r as be the coefficient of reflection from air to soil, and suppressing the e jωt time dependence as is customary, we write and Solving and making explicit the e jωt time dependence, we find that if we let v 0 = A 0 e jωt , and letting d = s 1 + s 2 + a 1 be the distance between the boreholes, then By linearity, it is clear that if our transmitted signal is a sum of various A n e jωnt (with A n possibly complex), then the wave at the Rx borehole will be given by a sum of solutions to (5). It will be convenient to rewrite (5) in terms of the wave which would have been received in the absence of a cavity. Denoting this wave by E ra = A 0 e jωt e −γd , we can rewrite (5) as Note that (6) is based on the following two assumptions: • The wave can be treated as a plane wave (the far-field assumption). • The effect of diffraction around the cavity can be neglected.
We attempted to validate (6) using GPRmax [22], [23] by simulating measurements in simple homogenous soil. When we prevented diffraction around the cavity (by extending the cavity in the appropriate directions until it reached the perfectly-matched layers bounding the simulation), we found the measurements shown in Figure 2, showing that (6) is indeed a reasonable model of the magnitude spectrum. However, Fig. 2. The magnitude spectra of signals which did and did not pass through a cavity, as well as the magnitude spectrum predicted for a signal which passes through a cavity. All magnitude spectra are normalized to have unit energy, and diffraction around the cavity was prevented. The cavity had a width of one meter, and the source was a derivative of the Gaussian waveform centered at 30 MHz. The background medium was homogeneous soil with relative permittivity of 10 and conductivity 0.002 S/m -the parameters given by [24] for "dry, sandy coastal land." (6) utterly fails to model both the phase spectrum and the received energy. We attribute this to the failure of the planewave assumption. Additionally, when we allow diffraction around the cavity, even small changes in depth (relative to the top or bottom of the cavity) can dramatically distort the received magnitude spectrum. While these distortions are difficult to model analytically 1 , if we can sample multiple depths at which the signal passes through the cavity we can hope that these distortions will average out. A version of Figure 2 based on real measurements appears as supplemental material for this paper.
Having seen that we are unable to model the phase spectrum but must suffice with the magnitude spectrum, we can make two additional simplifications: • While the Fresnel coefficients t sa , t as and r as are complex functions of frequency, in practice we can approximate them as purely real constants over the relevant frequency range. • We approximate e −γa ≈ e −j ωa v const -that is, we assume that the wave propagates through the distance a of soil at the speed v const regardless of frequency and suffers only negligible attenuation while doing so. The results of applying these simplifications are visually indistinguishable from those presented in Figure 2.

III. THE CAVITY-DETECTION METHOD
Let us assume that there exist either zero or one cavities and treat the problem as a hypothesis-testing problem. Let Y d (ω) be the magnitude of the spectrum of the normalized (with respect to energy) signal received at depth d and G(ω) be the magnitude of the spectrum received in the absence of a cavity. At each depth in a ZOP measurement there are two possibilities: The reason for using proportions rather than equalities is that, as mentioned above, our model is capable of predicting the received magnitude spectrum only up to a multiplicative constant. Additionally, we found, as did [17], that in practice our measurements of the received energy are unacceptably noisy. Having made the decision to work only with the proportions and having made the simplifying assumption that the Fresnel coefficients are independent of frequency, we can ignore the t sa t as factor in H 1 .
To use these two hypotheses for Neyman-Pearson-style hypothesis testing, it is necessary to first specify a distribution for the noise n(ω). As we work with magnitude spectra, the most natural choice for the distribution of n(ω) is probably the Rayleigh distribution. We have found, however, that that distribution will tend to lead to poor results. This is most probably because the noise is mostly not additive white Gaussian noise (AWGN) caused by e.g. thermal (Johnson-Nyquist) noise in the receiver and added to each measurement. Rather, the "noise" is mostly caused by the failure of the assumptions upon which we built our model and by the inhomogeneities of the soil. The Central Limit Theorem justifies the use of the normal distribution when the noise is due to the sum of these many different sources of error. Additionally, we note that the normal distribution is "well-behaved" in the sense that, no matter the value of the variance, the probability of an observation decreases monotonically with the distance of the observation from the mean. Use of the Rayleigh distribution, in contrast, might plausibly lead to the "failure mode" where the estimate of the scale parameter σ is too far from the correct value; an observation might then be penalized for having too little noise. If the estimate of the scale parameter σ is correct, then it is of course correct to penalize the observation for having too little noise, just as it is correct to penalize the observation for having too much noise. Still, if the performance of the algorithm is not robust with respect to the estimate of the scale parameter σ, that might partially explain why it yields poor results.
This being the case, we make the following four additional assumptions: • The noise n(ω) is normally distributed. • The width of the cavity is at least approximately known.
• If a cavity is present, then the wave impinges upon the cavity at approximately normal incidence (or more generally, the angle of incidence is at least approximately known). • The noise at each frequency is independent (the naïve assumption).
The last assumption does not hold in practice. Ideally we would compute likelihoods using a multivariate normal distribution, but since we work with a large number of frequencies we would then have to invert large matrices, rendering our method intractable. This is reminiscent of the problem in machine learning where a classifier is expected to make a decision based on a large number of features, yet estimating the covariance of the features or using it to compute probabilities is computationally intractable. One solution is to assume that all features are independent. This assumption, known as the naïve assumption, has been used in various machine learning algorithms, most notably the Naive Bayes algorithms. These algorithms tend to correctly (and tractably) classify data, though they are also known for being overconfident -since they may consider each of a large number of correlated features as an independent confirmation of a certain classification, these classifiers often return probabilities which are very close to zero or very close to one. While the probabilities returned by these classifiers are often next to useless, the classifications themselves tend to be correct. We should therefore expect that as we make use of our own naïve assumption, the likelihood ratios we find will be similarly overconfident. However, we hope that we too will find large values of the likelihood ratio only at depths at which a cavity does indeed exist.
We are almost ready to perform Neyman-Pearson hypothesis testing on the measurements received at each depth. All that remains is to estimate r as and G(ω). We will find r as by maximizing over all plausible values, and we will find G(ω) by a process inspired by cell-averaging Constant False Alarm Rate (CA-CFAR) radar.
We first slide a window over all depths in our profile. At each step, we will take the mean measurement of the window (at each ω) as the measurement of the center of the window (we hope that this step will allow us to "average out" deviations from our model caused by diffraction effects). When the window is centered at depth d, this value is Y d (ω), and we must determine its likelihood under H 0 and under H 1 . To do this, we will ignore all depths within a certain distance guardSize of the window (similar to the manner in which CA-CFAR radar commonly ignore the guard cells adjacent to the cell under test), and we will take the mean value (at each ω) of the measurements greater than a distance guardSize from the window to be G(ω). We will refer to these depths as the training depths. As discussed above, we do not wish our method to be influenced by the total energy of the received signals, so we normalize both Y d (ω) and G(ω). We similarly find the variance, for each ω, of G(ω) from the training depths. At this point, using our assumptions that the noise at each frequency is distributed normally and that the noise at each frequency is independent, we can simply compute the likelihood (or in practice, the log-likelihood) of measuring Y d (ω) under H 0 . We can similarly compute the likelihood (or log-likelihood) of measuring Y d (ω) under H 1 for any value of r as , so we take the likelihood under H 1 to be the maximum over all plausible 2 r as of the likelihood of H 1 . At this point, we are able to compute the likelihood 2 It is only necessary to compute the maximum over geologically plausible values of ras. In this paper, we have computed the maximum over all values of ras computed using the Fresnel equations when we take the refractive index n of soil to be a real number between one and five. All greater values of n were deemed to be geologically implausible. This article has been accepted for publication in IEEE Transactions on Geoscience and Remote Sensing. This is the author's version which has not been fully edited and content may change prior to final publication. Citation information: DOI 10.1109/TGRS.2022.3185802 ratio (in practice, the difference of the log-likelihoods) of the two hypotheses for every depth in the profile. A pseudocode description of this method appears in the appendix.
We must consider only a finite number of frequencies ω. By default, we use all frequencies sampled by the FFT of the entire received signal. In certain circumstances, it might plausibly be advantageous to consider only a subset of the possible frequencies. In the limiting case of considering only a single frequency, the naïve assumption holds trivially; if only a small set of widely-spaced frequencies are considered it might hold approximately 3 . Of course, in these cases we would be ignoring a large amount of potentially-useful information that resides in the portions of the spectrum not considered.
In this method we have attempted to estimate the "background signal" by calculating the mean from among the training depths. The justification for this is that if we assume that at each of the training depths, we measured the true background signal contaminated by multivariate Gaussian noise, then the maximum likelihood estimator (MLE) for the background signal is indeed the mean signal calculated over the background depths. This model is simplistic; in some cases where this model is less correct it might be possible for a more sophisticated method to more accurately estimate the background signal. In related research where we have used similar CFAR-style techniques to estimate the background signal we have sometimes found that we can improve results by taking the background signal to be the largest singular vector computed from the singular-value decomposition (SVD) of all background measurements. In this case the SVD method seems to interact poorly with the complexity of real soil.

IV. EVALUATION AND RESULTS
We validated the performance of our method on several simulations, as well as on a series of real measurements. Using GPRmax, we first simulated cross-borehole GPR measurements across homogeneous soil with a relative permittivity of 10 and conductivity of 0.002 S/m, the parameters given by [24] for "Dry, sandy, coastal land." A cavity (containing only free space) with a width of one meter, significantly extended in length (perpendicular to the line connecting the two boreholes), was placed between the depths of 3.85 and 5.5 meters below the surface. The transmitter and receiver boreholes were separated by six meters. ZOP measurements were simulated every 0.25 meters up to a depth of 10 meters. The source waveform was a derivative of the Gaussian waveform with a center frequency of 30 MHz. As seen in Figure 3, our method successfully detects the cavity. It is important to note, however, that perfectly homogeneous soil is unrealistically simple.
We further simulated a set of three ZOP profiles over soils characterized by strata. These simulations attempted to capture the complexity of real-world soils. Each stratum had dielectric properties that conformed to the fractal Peplinski model of 3 Alternatively, in this case it may be tractable to dismiss the naïve assumption and assume the noise obeys a multivariate Gaussian distribution, with the covariance matrix estimated from the training depths in a manner similar to that which we used to estimate the variance at each frequency. A more complete examination of this potentially powerful variation of our method is left for future research. Fig. 3. The results of applying our method to simulated homogeneous soil. The cavity exists between the depths of 3.85 and 5.5 meters below the surface. When applying our method, both the size of the window and the parameter guardSize were four depths (one meter). Fig. 4. The results of applying our method to simulated stratified soil. The cavity exists at the depths of the major peaks of the likelihood ratios given in blue in the first two scenarios. In scenario (c), the cavity exists at the depth of the smaller peak at around 4.3 meters below the surface. When applying our method, the size of the window was three depths (0.72 meters) and the value of the guardSize parameter was 12 depths (2.88 meters). The width of the cavity in all cases was assumed to be 1.5 meters. [25]; the properties in each stratum were chosen randomly using a Brownian process based only on the properties of the previous stratum. The probability distributions for each property were normal and uniform distributions whose parameters were based upon the measurements in Table 1 in [26] and Table 4 in [27]. Each scenario was simulated both with and without a cavity. A complete description of these simulated stratified soils is available as supplemental material for this paper. The results of applying our method to these simulations are shown in Figure 4. While these realistic stratified soils are clearly more challenging, our method was mostly successful in detecting the cavities.
In each of the first two scenarios, the cavity exists at the depth of the large peak in the likelihood ratio. On scenario (c), which is shown in Figure 5, our method fails; the cavity exists at the depth of the peak at around 4.3 meters below the  Figure 4. The simulations were identical other than the cavity which was simulated in the measurements shown on the right. The signals have been cropped to show the most relevant portions. The strange behavior of the signal at the shallowest and deepest depths is most likely a result of numerical artifacts and should have been ignored by our method due to its CFAR-style nature.
surface. The large peak at around 7.2 meters is a false positive; the stratum at that depth seems to act as a low-pass filter such that both our method and the CFDS method of [17] find an anomaly there, characteristic of a cavity, which is stronger than the anomaly caused by the actual cavity. The GPR signal also propagates anomalously quickly in that stratum, as it does at the depth of the actual cavity. Interestingly, more energy arrives at the receiver at the depth of the actual cavity; in general less energy arrives when a cavity is present as much of the energy is reflected away from the receiver at the soilair and air-soil interfaces. In short, the anomalous stratum at 7.2 meters appears cavity-like to every method (though at least traveltime, unlike the other methods, is slightly more anomalous at the depth of the actual cavity).
Finally, to show that our method can indeed detect real cavities, we carried out a series of real measurements at a test site characterized by sedimentary rock. A single cavity with a cross-section of around one square meter intersected two lines of boreholes. Eight ZOP profiles were taken between boreholes which did not bracket the cavity, while a further three profiles were taken between boreholes which did bracket the cavity. In each case, there were approximately seven meters between the boreholes and measurements were taken every 0.25 meters. While the exact signal transmitted is unknown, the energy of the received signal is centered at around 30-40 MHz (the received centroid frequencies are shown in Figure 8). The results of applying our method to these measurements are shown in Figure 6. Our method clearly succeeds in detecting the cavity. Note, however, that changing the parameters used in our method may result in a likelihood ratio at the anomaly at the deepest depths which can exceed the likelihood ratio at the true cavity. This anomaly is discussed further in the next section.
Our method compares favorably with standard cavitydetection methods based on traveltime. While the cavity from the homogeneous simulation in Figure 3 can be easily detected by the early arrival of the signal (not shown), the other examples are more complicated. As we have seen, the measurements in Figure 5 show a stratum in which the signal exhibits low Fig. 6. The results of applying our method to real ZOP measurements, taken from among different pairs of boreholes at the same test site. Blue measurements were taken across pairs of boreholes which bracketed the same cavity at around 19 meters below the surface, while red measurements were taken between pairs of boreholes which did not bracket any cavities. When applying our method, the size of the window was six depths (1.5 meters) and the parameter guardSize was one depth (0.25 meters). The cavity was assumed to have a width of 1.5 meters. Fig. 7. Raw signals from two of the real ZOP profiles used in Figure 6. Unlike the signals shown in Figure 5, here the measurements at each depth were normalized to have unit energy. The signals have dimensions of Volts.
traveltime. The measurements from Figure 6 exhibit similar strata, as shown in Figure 7. Here there are two strata through which the signal exhibits low traveltime; one lies at around 13 meters below the surface, and the other lies at around 19 meters below the surface. Detection is complicated by the fact that the cavity lies in or adjacent to one of these strata. Only by comparison between adjacent profiles can the cavity be detected using traveltime alone, and this can only be done because there are only 7 meters between the two boreholes. As the distance between the boreholes grows, the difficulty of distinguishing the effect of the cavity from the effect of the stratum on traveltime increases; 7 meters is approximately the largest separation between the boreholes at which it is possible to do so. Fig. 8. The centroid frequency of the same ZOP profiles used in Figure 6. It is clear that this feature, though powerful in other contexts [17], cannot be used alone to solve our cavity detection problem.

V. DISCUSSION
It might appear that (6) predicts that propagation through a cavity will have an effect on the magnitude spectrum of the wave that is cyclical in frequency. However, the soil itself will only pass low frequencies. This property is reflected in the fact that α given in (2) will grow as ω increases.
In practice, the only frequencies which are passed by the soil are those for which propagation through a cavity has a low-pass effect. It would then appear that rather than using our method, cavities may be simply detected by searching for depths at which the received spectrum has an anomalously large amount of power in the lower frequencies relative to the higher frequencies.
Indeed, in some cases we can detect a cavity simply by computing the centroid frequency as defined by [17] at each depth in a ZOP profile. Note that [17] define the centroid frequency as f s = ∞ 0 f P(f )df ∞ 0 P(f )df where P(f ) is the magnitude of frequency f . The centroid frequencies at each depth of the profiles used in Figure 6 are shown in Figure 8. It can be seen that the profiles where the wave passed through the cavity do have anomalously low centroid frequencies at around 19 meters below the surface, which is the approximate depth of the cavity. However, as the centroid frequency is often lower when the signal did not pass through a cavity, it is clear that we cannot use the centroid frequency alone to detect cavities.
As was noted in the preceding section, and as can be seen from Figure 7, our real measurements were taken at a test site in which there exists a stratum at around 13 meters below the surface such that the signal propagates through the stratum anomalously quickly. Yet as can be seen from Figure 6, our method, contrary to traveltime-based methods, does not find an anomaly at this depth (with the exception of a very minor anomaly in a single profile). We would indeed expect our method to generally be robust to these anomalies. Our method is at the greatest risk of returning a false positive when the low-pass effect of propagation through soil is anomalously strong (as it may be confused with the low-pass effect of propagation through a cavity); conversely, our method is at the least risk of returning a false positive when the low-pass effect of propagation through soil is anomalously weak. Since we see from (2) that the phase velocity v phase = ω β increases as decreases, and that the low-pass effect of propagation through soil, whose strength varies with α, increases as increases, we would expect not to see strata in which the wave both propagates anomalously quickly and simultaneously experiences an anomalously strong low-pass effect. This means we would expect our method to be robust to strata which would tend to cause a false-positive in conventional traveltimebased methods, and, conversely, conventional traveltime-based methods would tend to be robust to strata which would tend to cause our method to return a false positive.
One should be careful, however, not to overstate the preceding argument. In particular, we neglected to consider that and σ are themselves (weak) functions of frequency, and that the strata are not necessarily perfectly homogeneous 4 . Thus we see that scenario (c) in Figure 4 does indeed contain a stratum in which the signal propagates anomalously quickly and in which our method also returns a false positive. Additionally, we see from Figure 7 that the signal propagates anomalously quickly at the depth of the actual cavity, and we also see from Figure 6 that even when a cavity does not exist, the likelihood ratio is higher at that depth than it is for other depths (though it is still much smaller when a cavity does not exist than when it does exist).
There is an interesting anomaly at the deepest depths of our real measurements. We see from Figure 7 that the signal exhibits a much greater traveltime at these depths, while we see from Figures 6 and 8 that the soil at these depths has some sort of low-pass effect on the spectra of the signals. Additionally, an anomalously low amount of energy is received at these depths (not shown), though as mentioned we suspect that the measurement of the signal's energy is noisy 5 . It is possible that these effects are simply the result of the properties of the soil at these depths. We conjecture, however, that these effects are the result of borehole drift. While [19] write that "It is well known that about 5 meters are deviated per every 100 meters in digging vertical boreholes," we have seen that the actual deviation can differ from this figure by more than an order of magnitude in either direction. Furthermore, we have seen that the deviation is sometimes gradual and sometimes sudden (presumably caused by a refraction-like effect when drilling through different strata). It is therefore plausible that a sudden deviation in our boreholes occurs at these deepest depths; it is furthermore easy to see that random perturbations 4 Note that there is one profile in Figure 6 which exhibits an anomalous likelihood ratio at around 13 meters, and similarly exhibits an anomalously low centroid frequency at the same depth in Figure 8. This indicates that there is a localized anomaly at that depth which does not affect other profiles at that depth, leading us to conclude that indeed the stratum is not homogenous. 5 A geologist familiar with the lithology of the test site claimed that our measurements of the received energy were implausible. This accords with [17], who found in similar situations that measurements of the received energy were inaccurate. Alternatively, it is possible that our GPR accurately measured the received energy and it was only due to borehole drift that the measurements of the received energy were deemed implausible.
in the positions of our boreholes would tend to increase the distance between them.
If the distance between the boreholes were suddenly to increase, we would expect to see an increase in traveltime, a decrease in the amount of received energy, and a stronger low-pass effect on the spectrum of the signal; while the latter two effects might cause us to believe that a cavity lies at those depths, the first effect would disabuse us of that notion. Conversely, if the boreholes happened to draw closer together, then we would expect to see a smaller traveltime, which would cause us to suspect the existence of a cavity, along with a smaller low-pass effect (and greater received energy) which would allow us to conclude that a cavity does not in fact exist. This latter case of converging boreholes is illustrated in supplemental material to this paper.
Recall that our method performs Neyman-Pearson-style hypothesis testing at each depth with the two hypotheses being either the existence of a cavity in soil otherwise identical to the soil at the training depths, or alternatively that the soil is simply identical to that at the training depths. We could perhaps evaluate, as a third hypothesis, the possibility that the soil is identical to that at the training depths but that the distance between the boreholes has changed. However, this seems unlikely to be worth the additional complexity. We conclude that it is best to use both conventional traveltimebased methods and techniques such as the method presented in this paper, as while both techniques may be susceptible to false positives, it is unlikely that the same phenomena would simultaneously cause false positives in both techniques.

VI. CONCLUSION
We modeled the effect on the spectrum of a GPR signal of passing through a cavity and found that it will in general be a specific low-pass effect. Using this model, we developed a method to detect underground cavities using GPR ZOP measurements. Our method compares favourably to existing methods. Moreover, we conjecture that while both our method and existing methods may be susceptible to false positives, the phenomena which would tend to cause a false positive in one method would have the opposite effect on the other method. This implies that combining our method with traditional methods would lead to a more robust cavity-detection solution.
In this paper we used the spectral distortion implied by (6) as a feature to detect cavities. While it might seem more natural to have compared the entire signal predicted by (6) with the received signal, so far our attempts to do so have met with only limited success. This may be due to the inability of (6) to predict the phase of the received signal. We leave this as a topic for future research.
We have seen that our use of the naïve assumption destroys the probabilistic significance of the likelihood ratios our method returns. It would clearly be useful to have our method return accurate likelihoods for both hypotheses. We expect finding a more sophisticated replacement for the naïve assumption to be a fruitful topic for future research.
Using other features derived from (6) the way we have used the spectral distortion is an additional topic for future research.
Preliminary research shows that the group dispersion implied by (6) is a strong candidate for such a feature. Significantly, the group dispersion experienced by a GPR signal appears to be relatively unaffected by small deviations in the distance between the boreholes.
Finally, in this paper we have considered only detection of the cavity using ZOP measurements. Using the feature considered here to localize the cavity horizontally or to perform a tomographic inversion remains an interesting problem. Most tomographic inversion algorithms are based around integrating some property of the medium, such as the slowness or opacity, over the path of the ray. The feature we have considered here, however, is influenced by the relative position of the two walls of the cavity; as it is in some sense a nonlocal property, using it as a feature for tomographic inversion is nontrivial. We expect this to be a fruitful direction for future research.

APPENDIX A PSEUDOCODE DESCRIPTION OF OUR METHOD
A pseudocode description of the method described in section III is shown as Algorithm 1.