ON PARTICLE FILTERS FOR LANDMINE DETECTION USING IMPULSE GROUND PENETRATING RADAR

  ON PARTICLE FILTERS FOR LANDMINE DETECTION USING IMPULSE GROUND.pdf (Size: 1.05 MB / Downloads: 0)
ABSTRACT
In this paper, we present an online stochastic approach for landmine
detection based on ground penetrating radar (GPR) signals
using sequential Monte Carlo (SMC) methods. Since the existence
of true landmines is unknown and random, we propose to
use the reversible jump Markov chain Monte Carlo (RJMCMC)
in association with the SMC methods to jointly detect and localize
landmines in the light of observations. Computer simulations
on real GPR measurements demonstate the robust and consistent
performance of the proposed method.
I. INTRODUCTION
Owing to the good penetration, depth resolution and excellent detection
of metallic and nonmetallic objects, ground penetrating
radar (GPR) has become an emerging technique for landmine detection
[1–3]. A GPR system receives returned electromagnetic
signal from the ground by which landmines can be located, if
present. In reality, the signals originating from various types of
ground surfaces, like soil or clay, are nearly indistinguishable from
those of the genuine landmines. Thus robust and intelligent approaches
for the problem are needed.
Typical landmine detection approaches are based on background
removal and the corresponding techniques include adaptive
background subtraction [4], background modelling using a timevarying
linear prediction [5] and its improvement [6]. Zoubir et
al. [7] have proposed to combine a Kalman filter (KF) for state
estimation and a detection method based on statistic testing, but
this approach heavily relies on the appropriate selection of parameters
in order to perform properly. Tang et al. have suggested
in [8] using sequential Monte Carlo (SMC) methods, also known
as particle filtering (PF), [9–11] for landmine detection application,
where the ground bounce signals are estimated and removed
prior to localizing landmines.
In this paper, three major contributions are presented. Even
though the data model adopted in this paper primarily follows that
in [7], a simplified version is proposed with a smaller number of
variables without sacrificing the overall localization performance.
In a given scan of surface, not only is the number of objects unknown,
but it is also varying. The second contribution is to employ
the reversible jump Markov chain Monte Carlo (RJMCMC) [12] to
perform a soft rather than a hard detection of possible landmines,
followed by state estimation using PF. The last one is the evalua-
¤The work described in this paper was supported by a grant from the
Research Grants Council of the Hong Kong Special Administrative Region,
China (Project No. CityU 119605)
tion and demonstration of the superior and consistent performance
of the proposed method using real GPR measurements [13].
The rest of the paper is organised as follows. In Section II, we
provide the state-space model and the prior distribution and likelihood
functions. In Section III we present the development of the
SMC method with the RJMCMC for landmine detection application.
Computer simulations and evaluations on real GPR data are
included in Section IV, and conclusions are given in Section V.
II. STATE-SPACE MODEL
The radargram or B-scan of a GPR scan usually consists of
K distance measurements along the x-axis on every time index
n 2 f1; :::;Ng. Let Ál;k = [®T
l;k; ¯T
l;k]T be a combined
state vector containing the background and target vector
signals, ®l;k 2 RM£1 and ¯l;k 2 RM£1, where L =
bN=Mc is the number of strips, M is the strip size or successive
times samples at a given distance index k, and ()l;k =
[()Ml+1;k; ()Ml+2;k; :::; ()Ml+M;k]T . In addition, we define a
stochastic existence variable sl;k 2 f0; 1g, indicating whether or
not within the current M samples a landmine is present. Then a
general state evolution model can be expressed as follows
Ál;k = f(Ál;k¡1) + Bsl;kvl;k; (1)
where f(¢) may be a linear or nonlinear dynamical function,
Bsl;k =
·
IM 0M
0M sl;k £ IM
¸
with IM and 0M being an
M £ M identity and zero matrices. The vector vl;k 2 RM£1
consists of zero-mean, white Gaussian random variables with covariance
matrix §l
v =
·
¾2
v;0;lIM 0M
0M ¾2
v;1;lIM
¸
for strip l.
When comparing the above model in (1) with that in [7], one
may realize that there are two main differences. Firstly, the background
signal ®l;k is always present and varies in every scan k, regardless
whether a target exists (sl;k = 1) or not (sl;k = 0). Secondly,
the bias term suggested in the model in [7] has been combined
with the target signal ¯l;k. The motivations of our changes
are to simplfiy the dynamical model by reducing one unknown
variable, and hence to improve the overall performance of target
detection and localization.
For the random variable sl;k, we model it by the stochastic
relationship sl;k = sl;k¡1 + ²s [14], where ²s is an i.i.d. random
variable such that its prior distribution function is given by
p(sl;kjsl;k¡1) =
8><
>:
Pr(²s = ¡1) = pd;
Pr(²s = 0) = 1 ¡ pb ¡ pd;
Pr(²s = 1) = pb;
(2)
978-1-4244-2241-8/08/$25.00 ©2008 IEEE 225
where pb; pd 2 f0; 1g are respectively the probabilities of incrementing
and decrementing the number of targets such that pb = 0
if sl;k¡1 = 1 and pd = 0 if sl;k¡1 = 0.
From this point onwards, our parameters of interest are denoted
by µl;k = fÁl;k; sl;kg for strip l, whose joint prior distribution
function becomes
p(µl;kjµl;k¡1) =p(Ál;k; sl;kjÁl;k¡1; sl;k¡1);
=p(Ál;kjÁl;k¡1; sl;k) £ p(sl;kjsl;k¡1);
(3)
where p(Ál;kjÁl;k¡1; sl;k¡1) and p(sl;kjsl;k¡1) are the prior distribution
functions of Ál;k and sl;k, respectively.
It is further assumed that the states at different strips are statistically
independent as considered in [7] such that the joint prior
distribution function of µk becomes
p(µkjµk¡1) =
LY¡1
l=0
p(µl;kjµl;k¡1); (4)
where µk = [µ0;k; :::; µL¡1;k]T and p(µl;kjµl;k¡1) is the prior
distribution function in (3).
For the observation model, we have yl;k = g(µl;k) + ²l;k,
where the function g(¢) may be linear or nonlinear, and ²l;k is a
zero-mean, white Gaussian random variable vector with covariance
matrix §²l = ¾2
²lIM. As a result, the likelihood of the
observation yl;k due to µl;k can be expressed as p(yl;kjµl;k) =
N(yl;kjg(µl;k);§²l ).
Denoting yk = [y0;k; ::::; yL¡1;k]T and µk = [µ0;k; :::; µL¡1;k]T ,
we would like to estimate µl;k = fÁl;k; sl;kg 8l sequentially upon
the receipt of yk.
III. SEQUENTIAL MONTE CARLO METHODS
In the context of online parameter estimation using SMC methods
or PFs, we are interested in approximating the posterior distribution
¼(µkjy1:k) using Ns number of Monte Carlo samples fµ(i)
k g
for i = f1; :::;Nsg with their associated importance weights
[9–11, 15–17]. These samples are propagated and updated using
a sequential version of importance sampling as new measurements
become available. Hence statistical inferences, such as expectation,
maximum a posteriori (MAP) estimates, and minimum mean
square error (MMSE), can be computed from these samples.
III.1. Sequential Importance Sampling
Since we have assumed that the states from different strips are statistically
independent, we may separately rather than jointly compute
the particles and their associated weights for every individual
strip. In other words, we will have L independent PFs, each
of which estimates the posterior distribution function of µl;k of a
particular strip l. From a large set of Ns particles fµ(i)
l;k¡1gNs
i=1
with their associated importance weights fw(i)
l;k¡1gNs
i=1, we approximate
the posterior distribution function ¼(µl;k¡1jyl;1:k¡1)
as ¼(µl;k¡1jyl;1:k¡1) ¼
PNs
i=1 w(i)
l;k¡1 ±(µl;k¡1¡µ(i)
l;k¡1); where
±(¢) is the Dirac delta function and every particle µ(i)
l;k is generated
from µ(i)
l;k » q(µl;kjµ(i)
l;k¡1; yl;1:k), for i = f1; :::;Nsg. Among
many practical PFs, we choose to use the bootstrap PF [9] and
assign q(µl;kjµ(i)
l;k¡1; yl;1:k) = p(µl;kjµ(i)
l;k¡1) as in (3). The associated
importance weights w(i)
l;k can be recursively updated as
follows
w(i)
l;k / w(i)
l;k¡1 £
p(yl;kjµ(i)
l;k)p(µ(i)
l;kjµ(i)
l;k¡1)
q(µ(i)
l;kjµ(i)
l;k¡1; yl;1:k)
; (5)
with
PNs
i=1 w(i)
l;k = 1. It follows that the new set of particles
fµ(i)
l;kgNs
i=1 with the associated importance weights fw(i)
l;kgNs
i=1 is
then approximately distributed according to ¼(µl;kjyl;1:k).
Having the set of particles fµ(i)
l;k;w(i)
l;kg for i 2 f1; :::Nsg,
we may compute the residual energies for all strips as ´l;k =
PNs
i=1 w(i)
l;k £ VAR[yl;k ¡ g(Á(i)
l;k; 0)], where g(Á(i)
l;k; 0) is an estimate
of observation yl;k based solely on the background signals
and VAR[¢] is the variance operator. That is, if the observation
yl;k contains background signal only, the residual energy
´l;k is expected to be a small value. On the other hand, if the
observation yl;k contains target signal, using the background signal
alone to estimate yl;k is going to yield a much larger value
of ´l;k. We also denote the aggregate residual energy signal by
´k =
PL
l=1 ´l;k, where the residual energy signals are summed
over all strips l along the x-axis. This quantity will be used later
when receiver operating characteristic (ROC) curves are prepared
for performance evaluation in Section IV.
As the particle filter operates through time, only a few particles
contribute significant importance weights in (5), leading to
the well-known problem of degeneracy [11, 16]. To avoid this,
one needs to resample the particles according to their importance
weights. That is, those particles with more significant weights
will be selected more frequently than those with less significant
weights. More detailed discussion of degeneracy and resampling
can be found in [11].
III.2. Reversible Jump Markov chain Monte Carlo
Unlike some existing methods, like [7], rely on the computation
of test statistics in an ad hoc fashion and compare it with some
pre-determined threshold to deterministically decide if a landmine
is present, we propose to devise the RJMCMC, a variation of
Metropolis-Hastings (MH) algorithm [18,19], to explore all possible
model spaces and softly determine if a landmine is present in a
given strip of the radargram.
In our application for a given strip l at every k and particle
we randomly generate a new candidate µ?
l;k = fµ?
l;k; s?
l;kg
from a distribution function d(¢j¢), which is conditional on µ(i)
l;k =
fµ(i)
l;k; s(i)
l;kg. The candidate will be accepted with probability
» = minf1; rg, where
r =
¼(µ?
l;kjy) d(µ(i)
l;kjµ?
l;k)
¼(µ(i)
l;kjy) d(µ?
l;kjµ(i)
l;k)
£ J; (6)
with J being the Jacobian of the transformation from µ(i)
l;k to µ?
l;k.
In effect, the proposed RJMCMC method jumps between parameter
subspaces, thus visiting all relevant models. In particular three
different moves are randomly selected to enable the exploration of
the parameter subspace. For a given strip l at distance index k,
1. For each particle i = f1; :::;Nsg:
² Sample u » U[0;1].
² If (u < pb), then “birth move”.
226
Distance (k)
Time (n)
10 20 30 40 50 60 70 80 90
20
40
60
80
100
120
140
Fig. 1. The contour plot of the radargram from channel 46.
² Else, if (u < pd + pb), then “death move”.
² Else, update all parameters where sl;k = sl;k¡1.
2. k à k + 1, goto step 1.
In short, in the proposed method all Ns particles assume different
models according to the RJMCMC at a given scan k, and a
histogram of the available models can be constructed from fs(i)
l;kg
and a detection can be made if necessary.
III.2.1. Birth, Death and Update Move
When the birth move is selected, it is assumed that a landmine is
present in the current region, i.e., sl;k = 1 given sl;k¡1 = 0.
A candidate particle fµ?
l;k; s?
l;k = 1g is generated according to
(3) and is accepted with probability »birth = minf1; rbirthg,
where rbirth = exp
½
¡ 1
2
µ
e?T
§¡1
² e?¡e
0T§¡1
² e
0T
¶¾
, where
e? = yl;k ¡ g(µ?
l;k) and e0 = yl;k ¡ g(µ
0
l;k), and µ
0
l;k »
p(µl;kjµl;k¡1; s0
l;k = 0; sl;k¡1 = 0). When the death move is
selected, similar procedures in the birth move are taken in which
a candidate particle µ?
l;k with s(i)
l;k = 0 and another particle µ
0
l;k
with s(i)
l;k = 1 are generated. The candidate particle will be accepted
with probability »death = minf1; 1
rbirth
g. The particle
µ(i)
l;k = fµ(i)
l;k; s(i)
l;kg, where fµ(i)
l;k; s(i)
l;kg = fµ?
l;k; s?
l;kg if the candidate
is accepted. Otherwise, fµ(i)
l;k; s(i)
l;kg = fµ
0
l;k; s0
l;kg. The new
particle µ(i)
l;k will then be used to update the importance weights
w(i)
l;k. If, however, the update move is selected, we simply generate
the particle µ(i)
l;k » q(µl;kjµ(i)
l;k¡1; yl;1:k) with s(i)
l;k = s(i)
l;k¡1,
which will then be used to compute the importance weights w(i)
l;k.
IV. SIMULATION RESULTS
In this section the performance of the proposed algorithm on a set
of real GPR measurements [13] is evaluated. Furthermore, the results
from the proposed method will be compared with those from
Distance (k)
Strip (l)
10 20 30 40 50 60 70 80 90
1
2
3
4
5
6
7
8
9
10
Fig. 2. The contour plot of localization results of the radargram in
Fig. 1 from the proposed method averaged over 100 independent
trials.
the Zoubir’s method in [7]. Here we specify the form of the dynamical
function f(¢) and observation function g(¢). In particular,
we choose to use the same linear model as in [7] in order
to carry out a fair performance comparison. Accordingly, the dynamical
function is f(Ál;k¡1) = I2MÁl;k¡1 and the observation
function becomes gj(µl;k) = Hj =
£
IM j £ IM
¤
Ál;k,
for j = f0; 1g. Note that for other forms of functions, say nonlinearity,
chosen for f(¢) and g(¢) would not drastically affect
the performance of the proposed SMC algorithms as they are developed
to tackle problems that are nonlinear and non-Gaussian
[9–11, 15–17].
In the data set [13] there are a total of 91 channels, and in
this evaluation we consider one sub-set of channels from 42 ¡ 49,
where there are two landmine objects centered at the range of x-
positions: 28¡29 cm and 72¡75 cm, respectively. Fig. 1 exhibits
the radargram with K = 91 distance measurements and N = 150
collected in one of the channels, where the first 50 rows of data
have been taken away as the signals in this block correspond to
the responses to the ground surface which are inappropriate for
localization purposes.
The proposed method is investigated in this experiment for localizing
the objects withM = 10, Ns = 500, and pd = pb = 0:2.
A total of 100 independent trials are conducted on these channels.
According to the contour plots in Fig. 2, the proposed method
clearly localizes the objects when compared with the true radagrams
in Fig. 1.
To quantitatively evaluate the localization performance on the
real measurements, the ROC curves are constructed for the results
from the proposed method. Prior to constructing ROC curves for
the evaluation on both methods, we need to first define the assumed
width of a genuine landmine. In the absence of this piece of information,
we assume the width of a landmine is 7 distance indicies,
i.e., 3 indicies left and right to its assumed center. For comparison
purposes the Zoubir’s approach is also studied with these sets of
measurements. From Fig. 3, it is seen that the proposed method
outperforms the Zoubir’s approach when real measurements are
considered.
In summary according to the ROC curves the proposed ap-
227
0 0.2 0.4 0.6 0.8 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
False alarm probability
Detection probability
Zoubir’s method
Proposed method
Fig. 3. The ROC curves obtained from the proposed and the
Zoubir’s methods on the data sets over 100 independent trials.
proach is able to consistently localize the landmine objects from
real GPR data using the PF and RJMCMC method, and is superior
to the Zoubir’s method. Nevertheless, it requires higher computational
load when compared with that of the Zoubir’s and other
competing non-PF based methods.
V. CONCLUSION
A stochastic online landmine detection method for ground penetrating
radar data using the sequential Monte Carlo approach with
reversible jump Markov chain Monte Carlo (RJMCMC) has been
presented. The proposed method in association with the simplied
Zoubir’s data model takes the advantage of the RJMCMC to explore
different model spaces and to expend the extensive computation
only on the most possible model space. The benefit is obvious
that no hard or predetermined thresholds are needed to decide
which model should be used given the data. Computer simulations
demonstrate that the proposed approach is able to successfully localize
the landmine objects from real GPR data, and it outperforms
the Zoubir’s method at the expense of higher computational load.