Detecting Landmines
Subsurface imaging by space-time deconvolution for a GPR array.
The UN states that around 80-110 million landmines are laid
throughout the world. Landmines can remain active for decades, so
they continue to be a threat even after the conflict has long been
re-solved. Since mines can not distinguish between soldier and
civili-ans, they cause a high number of civilian casualties. Annually
thou-sands of people are killed or maimed by landmines, so the reason
for the need for clearing them is quite obvious. Various methods
exist to detect the landmines; each having his own advantages and
drawbacks.
Author: Niels van Tol, MSc
Ground Penetrating Radar
My thesis work treats the signal process-ing part of data acquired by an ultra-wideband Ground Penetrating Radar as described in the previous edition of the Maxwell. I will briefly discuss the basics of the particular GPR that I’ve used. The system is based on a time-domain ra-dar that sends out a very short pulse emit-ted by the transmitter antenna locaemit-ted at the top. Because the pulse is in the order of nanoseconds the bandwidth is in the order of gigahertz, hence the term
ultra-wide band. The high bandwidth leads to a resolution in the order of centimetres, making it suitable to detect the smaller landmines as well. The receiving array consists of 13 receiving loops all placed under the transmitter antenna at 7 cm in-tervals, totalling to a
width of 84 cm. This means that the array can scan a 2D area by only having to move the array over one di-mension. The set-up is shown in Figure 1.
Need for
sub-surface imaging
U n f o r t u n a t e l y there’s a fundamen-tal problem with the data acquired by any time-based GPR. Tosimplify the situation a 2D scenario is used. In the case of our radar we could say that we would only take data from one loop. Now, let’s take a look at Fig-ure 2. The antenna is located at the top and there’s an object located beneath it, which is denoted in red. It’s important to
realise that the unit of the y-axes is time, where increasing time is located lower on the axes. Imagine that the radar is mov-ing from left to right over the scan line. At position 1 the radar sends out the pulse and the reflection made by the “red” ob-ject is recorded by the receiving loop. The reflection in this case is simplified as a dot which is detected after a certain amount of time. This time depends among others on the distance between the radar and the object. Now the radar moves to position 2 (directly above the object) and the pulse is sent out again. This time the distance between the radar and the object is small-er leading to a reflection which is earlismall-er in time. When we extend this procedure to all over the scan-line this results in a hyperbola over the scan-line which has it’s apex at the antenna position where the distance between antenna and object is the smallest. So basically, a very small
object creates a large hyperbola in the data. A real example of this is shown in Figure 3. Here a small sphere was placed under the antenna and the antenna was scanned from left to right. The top two bars are the so-called antenna crosstalk; because the transmitter is above the Figure 1: Array of 13 receiving antenna loops
elling; the synthetic B-scan is formed by interpreting the entire system as a cas-cade of linear responses. In this case the system is the combination of the influ-ences of the radar, the air, the soil and the target itself. A schematic of this model is shown in Figure 4, where you can clearly see all the components that make up the reflection of a certain target. An example of a point-spread function (PSF) that has its scatterer at the center of the scan-line is shown in Figure 5.
Inverse Wiener filter and quality
criteria
So we need to solve Λz0 out of Equation
1. A very straightforward way would be receiving array the pulse first passes by
the receiving loop before reaching the tar-get. The hyperbola can be clearly seen in the image.
This distribution of energy does not give good localisation and resolution of the targets and therefore needs to be reversed. For that we use an imaging process called migration.
Problem description of my thesis
Several methods exist to do this migra-tion, but they are very computationally intensive and use only a limited number of characteristics of the radar and soil. One of these methods is called the dif-fraction stack algorithm. This algorithm is based on the travel times from trans-mitter to object to receiver and basically searches for all possible hyperbolas in the data. It gives good migration results but is very computationally intensive.A new method was introduced by [Scheers] that did not have these draw-backs. My job was to develop this method for our GPR and to extend it to a 3D array version. The diffraction stack algorithm was also implemented which served as a benchmark for the to be developed imag-ing process.
TyPES OF SCANS
When the GPR’s position is fixed and a reflection is recorded this has only time-information. This is a 1D-scan, also known as an A-scan. If we move an an-tenna pair over the scan-line with small intervals the recorded data can be com-bined to form a 2D-scan, also known as a B-scan. In our case, when we move the array over the scan-line, the scanned area is 2D; that together with the depth-information in the time-domain leads to a 3D-scan, also known as a C-scan.
Basic concept of migration by
deconvolution
Let’s say that the reflection of any object can be modeled by the combined reflec-tion of a set of point scatterers. So that objects in the ground are modeled by a certain number of point scatterers in a certain place.
If the combining of these reflections is linear this is a convolution in space and can be described by
b(y, t) = w(y, z0, t)⊗y,tΛz0(y, t)
where b(y, t) is the recorded 2D-scan,
w(y, z0, t) is the point-spread function
(point-target response) and Λz0 is the
scattering matrix (the objects). The re-corded 2D-scan (for an example, look at Figure 3) is also called a B-scan.
Λz0 is of course what we’re searching for.
It’s the collection of point scatterers that would result in a particular B-scan. Basi-cally it’s our processed dataset.
Now we’ll discuss how the point-spread
function w(y, z0, t) is determined. This
is the 2D-result we would expect to get back from the radar if there was only a very small point target located beneath the radar array. The point-spread func-tion is determined by using forward
mod-Figure 3: Original B-scan (signal level) Figure 5: Point-spread function [signal level]
0 0.15 0.3 0.45 0.6 y [m] 0 2.5 5 7.5 10 t [ns] 0 -5 -10 -15 -20 0 0.15 0.3 0.45 0.6 y [m] 0 2.5 5 7.5 10 t [ns] 0 -5 -10 -15 -20
to deconvolve the point-spread function
w(y, z0, t) out of the B-scan, but doing
this in the space-time-domain is very computationally intensive. So to perform this deconvolution we use an Inverse Wiener filter in the frequency-wavenum-ber domain which is described by
ˆΛ(ky, ω) =
B(ky, ω)· W∗(ky, ω)
W (ky, ω)· W (ky, ω)∗+ β
where B(kx, ky, ω) is the Fourier
trans-formed acquired data-set, W (kx, ky, ω)
the Fourier transformed point-spread
function and β the regularisation
param-eter (1/SNR).
So after this filtering operation, all that
remains is to transform the ˆΛ back to
the time-domain, right? Not quite. This
filter needs the parameter β in order to
function. Mathematically speaking, this parameter is the inverse of the Signal-to-Noise Ratio of the raw input data (the original B-scan). But since the signal is unknown, this is difficult to determine beforehand.
To overcome this issue two quality crite-ria were introduced by [Savelyev] to con-trol the deconvolution result.
The first one is the Energy ratio between the processed result and the original B-scan and can be described by
γ =
ˆΛ(y, t)
b(y, t)
The second one, the Error, is the measure of difference between the original signal and deconvolution result convolved with the point-spread function. Theoretically this convolution should result in exactly the recorded B-scan. The error can be de-scribed by
=
b(y, t) − ˆΛ(y, t) ⊗ w(y, t)
b(y, t) +ˆΛ(y, t) ⊗ w(y, t) Ideally the error would be zero, so the natural approach is to vary the regulari-sation parameter to minimise the error. The maximum Energy Ratio functions as a constraint for the parameter, since not only the Error matters but the quality of the resulting image as well. This maxi-mum Energy Ratio depends on the sce-nario (2D- or 3D-data), but should always be lower than 100% since that would mean that the processed result contains more energy than the original dataset. The regularisation parameter is varied it-eratively until the actual Energy Ratio is close to the maximum Energy Ratio Figure 4: Scheme
Figure 7: Deconvolution result of sub-surface landmine [dB] Figure 6: Diffraction stack result of sub-surface landmine [dB]
Q
0 0.25 0.5 0.75 1 y [m] 0 3 6 9 12 depth [cm] 0 -2.5 -5 -7.5 -10 0 0.25 0.5 0.75 1 y [m] 0 3 6 9 12 depth [cm] 0 -2.5 -5 -7.5 -10constraint. The Error is calculated after-wards to give an indication on the quality of the result.
Results
The B-scan that’s shown in Figure 3 is the result of a single sphere with a 2 cm diameter sphere placed at 41 cm from the transmitting antenna. In Figure 6 you can see the migrated result using the diffrac-tion stack algorithm. Figure 7 shows the result of the new algorithm. It’s clear that the deconvolution algorithm gives a more condensed image at a higher resolution. The diffraction stack needs 55 seconds to process while the deconvolution opera-tion needs less than 5 seconds!
Extending to a 3D array method
To determine the performance of the algo-rithm a 3D dataset is used that contained two identical targets at around 80 cm on the scan-line at the centre and the side of the array. All figures are now 2D-projec-tions of 3D datasets. So now the y-axes represents the scan-line and the x-axes represents the dimension over the width of the array (aperture dimension) and the time/depth-information is omitted. To extend this algorithm to a 3D version is relatively straightforward. The raw da-taset is now 3D, e.g. C-scan. So the PSFand the filter need to be extended to 3D as well. The challenge of making a 3D ar-ray version is two-fold. First there is the course resolution over the width of the array (remember that there are only 13 loops to cover 84 cm). Second the finite width of the radar all together (you can make the scan-line as long as you want, but the size of the radar is fixed).
To “increase” the resolution linear inter-polation was applied to the PSFs and raw C-scans over the array dimension. The result when taking these steps is shown in Figure 8 and the diffraction stack result in Figure 9. Only the cen-tre target is visible in the deconvolution result, you can still see a little image at the right side but it is much weaker than the centre target. Both targets are clear-ly visible in the diffraction stack result, but still the side target is weaker. This is because the finite width of the array: a target located at the side will have a part of its hyperbola cut off because the array just doesn’t detect it. On the other hand, when a target is located precisely at the centre of the array a much larger part of the hyperbola will be detected; resulting in a higher total energy for that particular target. This explains why the right target is a bit weaker in the diffraction stack
result, but not why it almost disappears with the deconvolution result!
The use of multiple PSFs
Well remember that the PSF was formed with one scatterer in the middle in the 2D-case. The same was done in the 3D-case. One scatterer was placed at the cen-tre of the scan-line and the cencen-tre of the array. This means that the filter is basi-cally looking for hyperbolae that are lo-cated in the centre, so centre targets are biased and will have stronger results than targets at the side.
To resolve this a method was developed that used multiple PSF’s; each having their scatterer at a different location. Using more PSF’s also results in more processed datasets, which have to be combined in order to receive one imaged result. With the help of some trial and error I found that the optimal method consisted of 13 PSF’s, where each PSF had their scatterer located under a different loop. The result is shown in Figure 10. Both targets can now clearly be seen, but naturally (as with the diffraction stack result) the side target is still perceived weaker due to the fact that there’s less hyperbola energy of a side target in the original C-scan.
Figure 8: ASWEP deconvolution result [dB] Figure 9: ASWEP diffraction stack result [dB]
-40 -20 0 20 40 x [cm] 100 75 50 25 0 scan-line [cm] 0 -5 -10 -15 -20 -40 -20 0 20 40 x [cm] 100 75 50 25 0 scan-line [cm] 0 -5 -10 -15 -20
But just because we have this same prob-lem with another algorithm, doesn’t make it right. Fact of the matter is that two identical targets should give two identical responses.
Weighting the result
To counter this a weighting scheme was developed that multiplied the signal lev-els of the deconvolution results from each separate PSF with a linear weighting fac-tor before combining the results. This weighting is of course relative and the deconvolution result obtained with the one PSF that has its scatterer at the cen-tre will have a weighting factor of 1. But how do we determine with which factor to multiply?
It is necessary to normalise the PSF’s be-fore filtering; a.o. in order to make sure that the quality criteria make sense, but this does remove the energy
informa-tion. You can imagine that a PSF with its scatterer at the side has less total energy than the PSF with its scatterer at the cen-tre. So the normalisation factors of the 13 different PSF’s are used to come up with a weighting factor for each decon-volution result. The factors are relative, which makes the lowest factor 1. The distribution is shown in Figure 11. The 13 weighting factors each correspond to a PSF which has their scatterer at a differ-ent position over the width of the array. Combining the 13 results after weighting leads to the final result shown in Figure 12.
Conclusions
The complete deconvolution algorithm needed less than 8 minutes to do the same as the diffraction stack algorithm did in over 22 hours. The shapes of the targets in the diffraction stack result are
nicer and the deconvolution result also has more artefacts. But the new algo-rithm does show both targets at compa-rable strengths. The good image quality together with the very big cut in compu-tational effort make the newly developed method very interesting.
Interested?
You want to do a challenging thesis work with good guidance? You do not want to create a report that disappears in a cabi-net? You want to work with a GPR that is designed to detect landmines? Anywhere from antenna design to signal processing? Contact Dr. Alexander Yarovoy
a.yarovoy@irctr.tudelft.nl
or T.G. Savelyev
t.g.savelyev@irctr.tudelft.nl
at IRCTR.
References
[Savelyev] T. G. Savelyev, L. van Kempen, and H. Sahli. Ground Penetrating Radar, 2nd edition, chap-ter “Deconvolution techniques”, pages 298–310. D. J. Daniels ed., IEE Radar, Sonar, Navigation and Avion-ics Series 15, 2004.
[Scheers] B. Scheers and M. Acheroy. Ground
Penetrating Radar, 2nd edition, chapter “Migration
technique based on deconvolution”, pages 283–293. D. J. Daniels ed., IEE Radar, Sonar, Navigation and Avionics Series 15, 2004.
A
Figure 10: ASWEP result with median filtering [dB]
Figure 11: Weights for each deconvolution result related to position of scatterer
Figure 12: ASWEP combined deconvolution result [dB]
-40 -20 0 20 40 x [cm] 100 75 50 25 0 scan-line [cm] 0 -5 -10 -15 -20 -40 -20 0 20 40 x [cm] 100 75 50 25 0 scan-line [cm] 0 -5 -10 -15 -20
