Information driven search for point sources of gamma radiation

ARTICLE IN PRESS Signal Processing 90 (2010) 1225–1239

Contents lists available at ScienceDirect

Signal Processing journal homepage: www.elsevier.com/locate/sigpro

Information driven search for point sources of gamma radiation Branko Ristic a,, Mark Morelande b, Ajith Gunatilaka c a

ISRD, Defence Science and Technology Organisation, 506 Lorimer Street, Fishermans Bend, VIC 3207, Australia Melbourne Systems Laboratory, Department of Electrical and Electronic Engineering, The University of Melbourne, VIC 3010, Australia c HPPD, Defence Science and Technology Organisation, 506 Lorimer Street, Fishermans Bend, VIC 3207, Australia b

a r t i c l e in fo

abstract

Article history: Received 10 March 2009 Received in revised form 3 July 2009 Accepted 9 October 2009 Available online 22 October 2009

The problem is to estimate the number of radioactive point sources in a specified area and to estimate their parameters (locations and magnitudes), using measurements ¨ collected by a low-cost Geiger–Muller counter. The measurements are Poisson distributed with the mean proportional to the radiation field intensity. The radiation field represents a superposition of background radiation and the source contributions subjected to the inverse distance squared attenuation. The solution is based on an information gain driven search which comprises a sequential Bayesian estimator coupled with a sensor/observer control unit. The control unit directs the observer(s) to move to new locations and acquire measurements that maximise the information gain in the Re´nyi divergence sense. The performance of the proposed information driven search, including a comparison with a unform search along a predefined path, is studied by simulations. A successful application of the proposed technique to experimental datasets, recently collected in the field trials, verifies the measurement model and the theoretical considerations. Crown Copyright & 2009 Published by Elsevier B.V. All rights reserved.

Keywords: Nuclear search Radioactive sources Sequential Monte Carlo estimation Sensor management Particle filter

1. Introduction Of growing concern for national security of many nations are numerous incidents of lost or stolen radioactive materials [1,2]. The ability to rapidly detect and localise illegal radioactive materials is an important prerequisite for the emergency responders to react quickly. Consequently, the topic of detection, localisation and search for radioactive material has created a lot of interest in the scientific community recently. Howse et al. [3] considered localisation of a moving radioactive point source based on the least squares algorithm. Hjerpe et al. [4] proposed to detect a gamma radiation source using statistical tests (Mann–Whitney and Kolmogorov–Smirnov). Nicol et al. [5]

 Corresponding author.

E-mail addresses: [email protected] (B. Ristic), [email protected] (M. Morelande), [email protected] (A. Gunatilaka).

developed several geometry inspired solutions to detect vehicles carrying radioactive sources over the border crossings. Another geometry-based solution to localisation of a single radioactive source was proposed in Rao et al. [6]. Both Nemzek et al. [7] and Brennan et al. [8] explore the effectiveness of a large sensor network deployed along a roadway for detection of nuclear materials in moving vehicles. A somewhat similar problem was also considered in [9]. Distributed detection of static radioactive sources using a small number of radiation counters is proposed in [10]. The Cramer–Rao lower bound and a statistically/ computationally efficient algorithm for detection and localisation of multiple static radioactive sources are presented in [11]. The subject of search for radiological sources has been less explored. The main characteristic of search is the ability to control the observer(s) in order to detect/localise sources. The conventional approach to search is based on a predefined path that scans the area in a uniform manner (e.g. parallel sweep search) [12]. Ziock and Goldstein [13]

0165-1684/$ - see front matter Crown Copyright & 2009 Published by Elsevier B.V. All rights reserved. doi:10.1016/j.sigpro.2009.10.006

ARTICLE IN PRESS 1226

B. Ristic et al. / Signal Processing 90 (2010) 1225–1239

study the design of parallel sweeps for a constant velocity observer in order to detect a radioactive source of known strength. Klimenko at al. [14] also consider a nuclear search strategy along a predefined path, but they adaptively change the exposure time, using the ideas from the sequential testing theory. The modern approach to search is information-driven. The observer is controlled sequentially, based on the current posterior probability density function (pdf) and the anticipated information content of future measurements, to select the best sensing action. Information-driven observer control, or sensor management, is widespread in many fields of engineering, such as bearings-only motion control [15,16], target tracking [17–19], classification [20], computer vision [21], machine learning [22], robotics [23], etc. In the context of search for nuclear material, information-driven observer control has been applied in [24,25]. Cortez et al. [24] divide the search area into a regular grid of cells. Using a simplified measurement model, where the counts observed in a particular cell depend only on the radiation level in that cell, they propose a search strategy where: (a) the observer stays in one cell until the variance of the measurement falls below a certain threshold; (b) the observer moves to a neighbouring cell selected based on its ability to reduce the uncertainty in the estimated radiation field. Ristic and Gunatilaka [25] consider a single-sensor single-source set-up where the observer control is based on the Fisher information gain of future measurements. This paper proposes an information driven search algorithm for an unknown number of point sources of gamma radiation, carried out by a team of coordinated multiple mobile observers. The characteristics of sources, such as their type (e.g. uranium, cesium, etc.) and strength of the radioactive material, are unknown. Each observer is ¨ equipped with a low-cost Geiger–Muller (GM) counter, a precise self-localisation instrument (e.g. a differential GPS receiver) and a two-way communication link with a data fusion and control (DFC) unit. The DFC unit sequentially performs estimation and, based on its current estimates, controls the observers to take actions. An action vector specifies to an observer its new location (i.e. where to go) and its radiation exposure time (i.e. how long to stay there) in order to acquire the next measurement. Once the observer completes its action, it sends back to the DFC processor the readings from its GM counter. The DFC responds with a new action vector and the cycle continues. The search area is assumed to be open, flat and polygon shaped. The DFC processor is assumed to know the boundaries of the search area as well as the average background radiation. The measurements are Poisson distributed with the mean proportional to the radiation field intensity. The radiation field represents a superposition of background radiation and the source contributions subjected to the inverse distance squared law. The estimation of the number of sources and the parameters of sources are carried out in a sequential Bayesian framework implemented via a particle filter in a manner similar to the multi-target track-before-detect [26,27].

This framework for hybrid estimation is particularly suitable to situations where one deals with raw (unprocessed) intensity measurements. The observer control is performed by a sensor management policy which maximises the information gain in the Re´nyi divergence sense. The contribution of the paper is thus the development of a radiological search algorithm, which combines a multi-target track-before-detect particle filter with information-driven observer control. The resulting algorithm, which uses multiple mobile sensors to perform a coordinated search for an unknown number of sources, is a significant improvement on existing radiological search algorithms which operate with a single sensor and use a simplified measurement model [24] or assume only a single source [25]. The proposed algorithm has been tested and verified both by simulations and using experimental datasets. Although the paper is devoted to a niche application of radiological search, the proposed framework of multisource parameter estimation with information driven data collection may be applicable in a variety of signal processing applications, provided the appropriate measurement model is used. The paper is organised as follows. A formal problem description and the measurement model are presented in Section 2. The sequential Bayesian joint estimator of the number of sources and the source parameters is described in Section 3. The observer control is studied in Section 4, with important implementation details described in Section 5. The simulation and experimental results are presented in Sections 6 and 7, respectively. Conclusions and future research directions are summarised in Section 8. 2. Problem description and the measurement model The radiation counts from nuclear decay obey Poisson statistics [14,28]: the probability that a gamma radiation detector registers z 2 N counts (N is the set of natural numbers including zero) in t seconds, from the source that emits on average m counts per second is Pðz; lÞ ¼

lz z!

el ;

ð1Þ

where l ¼ mt is the parameter (both the mean and the variance) of the Poisson distribution. Let A be a flat open-field polygon-shaped search area, whose map is known to the DFC processor. Assume that an unknown number of sources rZ0 are present in A. The case r ¼ 0 corresponds to the absence of any radiation sources, while for r40, each point source of gamma radiation s (where s ¼ 1; 2; . . . ; r) is parameterised by

 Its Cartesian coordinates ðxs ; ys Þ 2 A.  Its equivalent strength, fs , which is a single parameter that takes into account the activity of the radioactive source, the value of gamma energy per disintegration and the scaling/conversion factors involved [29]. The parameter vector of source s is thus given by xs ¼ ½xs ys fs > ;

ð2Þ

ARTICLE IN PRESS B. Ristic et al. / Signal Processing 90 (2010) 1225–1239

where > denotes the matrix transpose. The source parameter vectors are collected into a stacked vector > > Xr ¼ ½x> 1 ; . . . ; xr  . The measurements of radiation field are made using identical low-cost omnidirectional GM counters1 mounted on mobile observers. Let zik 2 N, be a measured count from a GM counter, taken at time tk by observer i i i 2 f1; . . . ; Og located at ðxk ; zk Þ 2 A. Given the number of point sources r and the parameter vector Xr , the measurement likelihood function can be expressed as [28,29] i

pðzik jr; Xr Þ ¼ Pðzik ; lk ðXr ÞÞ;

is the mean radiation count at location where i i ðxk ; zk Þ given by ! r X fs bi dis i k lik ðXr Þ ¼ mb þ e tk ; ð4Þ is 2 s¼1 ðdk Þ with qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi i i ðxk  xs Þ2 þ ðzk  yi Þ2

ð5Þ

being the distance between the source s and the i th sensor at time tk . An explanation of the terms in (4) follows:

 tik is the radiation exposure interval of sensor i at time 



way that excessively small observer-to-source distances are not allowed. By controlling the motion and the exposure intervals of observers which collect the radiation counts, the goal is to determine in the most efficient manner (the shortest time and the smallest number of measurements) the number of radioactive point sources r, their locations and strengths. The measurement model is given by (4) with airattenuation factors assumed to be bi  0. The estimation part of the search will be described in the next section. The decision logic for measurement locations and exposure intervals will be presented in Section 4.

ð3Þ

lik ðXr Þ

dis k ¼

tk (meaning that the count measurement zik is collected over the interval ½tk  tik ; tk ). bi in (4) is the air absorption coefficient which depends on the energy spectrum of gamma radiation (and thus the type of the radiation source). For cesium and cobalt, for instance, this coefficient equals 9:34  103 m1 and 6:86  103 m1 , respectively [30]. Since we use a low-cost GM counter which does not provide the spectral information and cannot identify the type of the radiation source, we will make an approximation i here. In order to remove the dependence of lk ðXr Þ on the type of radioactive source, in our study we restrict ourselves to relatively small search areas (with typical In this way we distances dis k of a few tens of meters). is can adopt the approximation ebi dk  1 in (4), which appears reasonable for cesium and cobalt. mb represents the average count-rate due to the background radiation (from cosmic rays and from naturally occurring radioactive materials in earth). This number is assumed to be known and constant over the search area (a realistic assumption with A being a small open field area).

Note that according to (4), if the distance between an observer and a source is zero, the measured count would be theoretically infinite. In practice this is not the case since a GM counter cannot measure radiation counts above a certain level (due to the saturation effect). We avoid this anomaly by controlling the observer(s) in such a

3. Sequential Bayesian estimation The Bayesian approach to parameter estimation provides a complete probabilistic description of the parameters of interest. This is important in sensor management since it is desired to select suitable actions in the presence of uncertainty. In fact, the Bayesian representation of parameter uncertainty, the posterior pdf, provides the information state upon which decisions are made in a partially observed Markov decision process (POMDP) [31] (the framework in which the search algorithm is formulated in Section 4). A further useful property of Bayesian estimation is that it can easily be performed in a sequential manner. This is important when sensing actions must be continually re-evaluated based on new data. 3.1. Conceptual solution In order to estimate the number of sources present in the search area A, parameter r is considered as a discretevalued random variable which can take values r 2 f0; 1; . . . ; rmax g. Joint estimation of the number of sources and their parameters (locations and strengths) is then treated as a hybrid (continuous–discrete) sequential Bayesian estimation problem [32, Chapter 11, 33]. Measurements of radiation field are used directly, as intensity measurements, rather than being subjected to a detection threshold. This will result in an estimation framework which effectively corresponds to multi-target track-before detect [27], where the concepts of probability of detection and false alarm on the measurement level do not exist. In order to simplify notation in this section we will drop the superscript i referring to the observer identification. Let the joint posterior density of the number of sources r and the parameter vector Xr , after processing k  1 measurements, be denoted by pðr; Xr jz1:k1 Þ, where z1:k1 ¼ fz1 ; . . . ; zk1 g is the accumulated set of count measurements. Given a new k th measurement zk , the goal is to construct the updated posterior density pðr; Xr jz1:k Þ using the measurement likelihood model described in the previous section. The update step is done via pðr; Xr jz1:k Þ ¼

1

1227

Low-cost GM counters are unable to analyse the energy spectrum of gamma radiation.

pðzk jr; Xr Þpðr; Xr jz1:k1 Þ ; pðzk jz1:k1 Þ

ð6Þ

with the likelihood function pðzk jr; Xr Þ introduced in (3).

ARTICLE IN PRESS 1228

B. Ristic et al. / Signal Processing 90 (2010) 1225–1239

Once the posterior pdf pðr; Xr jz1:k Þ is known, the posterior probability Pk ðrÞ ¼ Prfrjz1:k g

ð7Þ

that r sources are present in A is computed as the marginal: Z Pk ðrÞ ¼ pðr; Xr jz1:k Þ dXr : ð8Þ From this, the maximum a posteriori (MAP) or the expected a posteriori (EAP) estimates of the number of sources [34], after processing k measurements, can be computed as MAP ¼ arg max Pk ðrÞ; r^ k

ð9Þ

0rrrrmax

EAP r^ k ¼

rmax X

rPk ðrÞ;

ð10Þ

r¼0

respectively. Similarly, the posterior pdfs of individual source parameter vectors xs can be computed as the marginals of pdf pðs; Xs jz1:k Þ for s ¼ 0; . . . ; r^ k. 3.2. Particle filter The described hybrid state estimator is implemented in the form of a particle filter (PF) [32]. In this framework the posterior pdf is represented by a random sample fwnk ; Xnk ðrkn ÞgN n¼1 as pðr; Xr jz1:k Þ 

N X

wnk dðXr



Xnk ðrkn ÞÞ;

ð11Þ

n¼1

where N is the size of the sample (number of particles), wnk is the normalised importance weight of particle n ¼ 1; . . . ; N and d is the Dirac function. Each particle Xnk ðrkn Þ thus characterises both the number of sources rkn 2 f0; . . . ; rmax g and the stacked parameter vector ½ðxn1;k Þ> ; . . . ; ðxnrn ;k Þ> > . k

3.2.1. Prior pdf Before any measurement is received (at k ¼ 0), the particles are initialised using the prior pdf. The prior probability for the existence of r sources in A is adopted as rmax  r þ 1 P0 ðrÞ ¼ Prfrjz0 ¼ |g ¼ Prmax r¼0 ðr þ 1Þ

ð12Þ

for r ¼ 0; . . . ; rmax. In this way, initially, the majority of particles r0n P0 ðrÞ anticipate that there is no source present in A. For particles characterised by r0n 40, the prior pdf of the components of the source parameter vector in the location–strength ðx; y; fÞ space are adopted as follows:

 The source location coordinates ðxns;0 ; yns;0 Þ, s ¼ 1; . . . ; r0n ,



need to be chosen so that this prior is sufficiently vague to cover a large span of possible strength levels.

are uniformly distributed over the search region A, i.e. ðxns;0 ; yns;0 ÞUðAÞ. As A is polygon-shaped, the initial random placement of particles is carried out by the acceptance–rejection method [35]. The prior source strength is gamma distributed with shape and scale parameters k and c, respectively: fns;0 Gðk; cÞ, s ¼ 1; . . . ; r0n . The parameters k and c

The parameter vector components of the particles characterised by r0n ¼ 0 are not defined. The importance weights are initially uniform, i.e. wn0 ¼ 1=N, n ¼ 1; . . . ; N. 3.2.2. Update step After receiving a measurement zk , ðk ¼ 1; 2; . . .Þ, the update step is carried out according to (6). Recall that in the sequential importance sampling framework, the posterior pdf pðr; Xr jz1:k Þ is approximated by the empirical distribution fwnk ; Xnk ðrkn ÞgN n¼1 as in (11). Assuming that the importance density equals the prior (the prior in this case is the posterior at time k  1, pðr; Xr jz1:k1 Þ), it can be shown that the importance weights are computed as [32] wnk ¼ Cwnk1 pðzk jXnk ðrkn ÞÞ;

ð13Þ

where C is a normalisation constant, wnk1 is the normalised weight of particle n at time k  1 and pðzk jXnk ðrkn ÞÞ is the likelihood of particle n given by (3). In our implementation of the PF, the particles are resampled at every step k. The role of resampling is to eliminate the samples with low importance weights and to replicate the samples with high importance weights [32]. Due to resampling wnk1 ¼ 1=N and (13) is further simplified to wnk ¼ C1 pðzk jXnk ðRnk ÞÞ;

ð14Þ

with C1 ¼ C=N being a normalisation constant. The weights in (14) are used to resample the particles at time k. The straightforward use of (14), followed by resampling, can cause the particle filter to collapse due to particle degeneracy for two reasons. First, the posterior at k  1 can be expected to be much more diffuse than the likelihood. As a result, many (if not all) samples will be drawn in undesirable parts of the parameter space and consequently may be given zero weights. Second, since we consider a static parameter estimation case, the variety of particles will decrease with k because of the accumulation of the repetition in the sample and all particles may ultimately concentrate on a single point of the state-space. In order to overcome the described deficiencies of a straightforward implementation of the update step, we apply the progressive correction (PC) technique with particle regularisation [36]. In order to explain the PC technique, let us rewrite (6) in a shortened form as p ¼ ‘p0 =c, where ‘ is the likelihood function, p and p0 are the posterior and prior pdfs, respectively, and c is the normalisation constant in the denominator. Suppose we decompose the likelihood as ð‘1 ; ‘2 ; . . . ; ‘S Þ

such that ‘ ¼

S Y

‘j :

ð15Þ

j¼1

Observe that the update p ¼ ‘p0 =c is equivalent to the sequence of update steps:

p ¼ ‘S ð‘S1 ð. . . ‘2 ð‘1 p0 =cÞÞÞ:

ð16Þ

ARTICLE IN PRESS B. Ristic et al. / Signal Processing 90 (2010) 1225–1239

Let us adopt for j ¼ 1; . . . ; S that ‘j ¼ ‘gj , with power gj such that S X 0rgj o1 and gj ¼ 1:

Resampling (lines 14 and 26) is based on the systematic resampling algorithm [38], described by pseudo-code in [32, Table 3.2]. Particle regularisation is necessary in order to avoid particle degeneracy. Lines 15–17 describe the regularisation of a subset of particles in the source parameters space ðx; y; fÞ. It involves jittering the samples by a small random vector drawn from a Gaussian kernel Ks (line 16). The kernel bandwidth depends on the sample covariance matrix, for details see [36]. In line 27 of the algorithm we perform regularisation in the source number space r by forcing (in a probabilistic manner) some particles to change their value of rkn as follows. A particle with rkn ¼ rormax is forced, with probability Pb 51 (birth), to add one more source to its concatenated parameter n becomes rkn þ 1 with probability Pb ), vector (hence rkþ1 with the added parameter vector drawn from the prior pdf in the ðx; y; fÞ space. Similarly, a particle with rkn ¼ r40 is forced, with probability Pd 51 (death), to remove one source (chosen randomly) from its concatenated parameter n becomes rkn  1 with probability Pd ). Eq. vector (hence rkþ1 (14) is implemented in line 24 of the algorithm. The particle filter outputs sequentially, at each time index k, probabilities Pk ðrÞ of (8), r ¼ 1; . . . ; rmax , approximated as: Pk ðrÞ  Mr =N. The MAP estimate is then computed using (9) and the EAP estimate (10) is approximated as

j¼1

The update step of the PF will be carried out in S stages. The posterior at stage s ¼ 1; 2; . . . ; S is given by Q ps ¼ ðp0 =cÞ sj¼1 ‘j , with the posterior at the final stage pS ¼ p. As a result, in the early stages (for small values of s), the importance weights are computed using likelihood Qs j¼1 ‘j which is broader than the true likelihood ‘. In the later stages the likelihood approximation sharpens and ultimately becomes the true likelihood. The unnormalised weights are computed at stage s ¼ 1; . . . ; S as wnk ðsÞ ¼ ‘s ¼ ½pðzk jXnk ðrkn ÞÞgs :

1229

ð17Þ

After computing the weights at stage s, the particles are resampled and regularised. The main idea of progressive correction is to help more particles survive the resampling process by using broader likelihood at early stages and gradually moving the samples in the region of the parameter space which is supported by the likelihood function. The PC can be seen as a stochastic version of the particle flow with log-homotopy, recently proposed in [37]. A description of the update step of the particle filter is given in Table 1. The progressive correction is applied separately to each subset of particles characterised by the same value of rkn 40. The formation of these subsets is described in lines 3–8; Mr denotes the cardinality of a subset of particles with r sources ðr ¼ 0; . . . ; rmax Þ. The selection of coefficients gs in line 9 can be done in many ways; we have used the adaptive scheme suggested in [36].

EAP r^ k 

N 1X rn : N n¼1 k

ð18Þ

When the search is complete (the search termination to be discussed in Section 5), and r^ k 40, the PF outputs the

Table 1 The update step of the particle filter at time k ¼ 1; 2; . . .. 1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15: 16: 17: 18: 19: 20: 21: 22: 23: 24: 25: 26: 27: 28: 29:

n function Update StepðfXnk1 ðrk1 ÞgN n¼1 ; zk Þ for r ¼ 1; . . . ; rmax do m¼0 for n ¼ 1; . . . ; N do n ¼ r then if rk1

x Formation of a subset of particles with r sources

n m ¼ m þ 1; IðmÞ ¼ n; Ym;0 ¼ Xnk1 ðrk1 Þ; k1 end if end for S s gs¼1 ,

Select coefficients fg for s ¼ 1; . . . ; S do for m ¼ 1; . . . ; Mr do

s. t. gs 2 ð0; 1 and

PS

s¼1

m;0 Mr x A subset of particles fYk1 gm¼1 is formed

gs ¼ 1

Þgs ; Compute the weights: wm;s ¼ C½pðzk jYm;s1 k1 end for m;s1 Mr m;s Mr Resampling: ½fYm gm¼1  k1 gm¼1  ¼ RESAMPLE½fYk1 ; w for m ¼ 1; . . . ; Mr do m;s m þ em;s Draw em;s Ks and set Yk1 ¼ Yk1 end for end for for m ¼ 1; . . . ; Mr do

XIðmÞ ¼ Ym;S ; rkIðmÞ ¼ r; k k1 end for end for for n ¼ 1; . . . ; N do Compute the weights: wnk ¼ C½pðzk jXnk ðrkn ÞÞ; end for n n n N n N Resampling: ½fXn k ðrk Þgn¼1  ¼ RESAMPLE½fXk ðrk Þ; wk gn¼1  Regularisation in the source number space r n N return fXn k ðrk Þgn¼1 end function

x Progressive Correction in S stages x C such that weights sum to one

x Regularisation in source paramater space ðx; y; I Þ x Ks is a Gaussian kernel

x C such that weights sum to one

ARTICLE IN PRESS 1230

B. Ristic et al. / Signal Processing 90 (2010) 1225–1239

estimate of the concatenated parameter vector computed as the sample mean: N X ^ k ðr^ k Þ ¼ 1 Xn ðr n Þ dK ½rkn ; r^ k ; X Mr^ k n¼1 k k

ð19Þ

where dK ½i; j is the Kroneker delta. The sample covariance matrix is also computed to measure the spread of particles, as one of the search termination criteria. 4. Observer control In this section the problem of observer control is formulated as a POMDP [31]. The elements of a POMDP include an information state, or posterior pdf, a control vector belonging to a set of admissible controls and a reward function. Let uk ðik Þ 2 U k ðik Þ denote the control applied at time tk to observer ik 2 f1; . . . ; Og where U k ðiÞ is the set of admissible control vectors at time tk for observer i. The notation is extended to include the applied control vectors so that pðXr ; rjz1:k ; u0:k1 ði1:k1 ÞÞ is the posterior pdf after applying the control vector sequence u0 ði0 Þ; . . . ; uk1 ðik1 Þ and observing the measurement sequence z1:k . After processing the first k measurements the control vector to be applied at time tk to observer ik is selected as uk ðik Þ ¼ arg max E½Dðv; pðXr ; rjz1:k ; u0:k1 ði0:k1 ÞÞÞ; v2U k ðik Þ

ð20Þ

where D is the real-valued reward function. The reward will typically depend on the future measurement, acquired after the control is applied. Since it is desired to maximise over control vectors v 2 U k ðik Þ in (20) without actually applying them (without acquiring all corresponding future measurements), the expectation E of the reward is taken with respect to the prior measurement pdf. An optimal sensor management strategy would have to look an infinite number of steps ahead, rather than one step ahead (myopic) as in (20). In the absence of obstacles in the search area, however, myopic or one step ahead sensor management is a reasonable tradeoff between the computational requirements and the performance [19]. Approximation computation of the posterior pdf using a particle filter was described in Section 3. In this section we describe the control vectors and reward function. 4.1. Control vectors We assume a centralised fusion architecture, where an observer, upon acquiring a count measurement, sends a message (measurement value, own location and ID) to the

τ1k Observer 1

time 2 τk+2

2 τk+1

Observer 2

time tk

tk+1

tk+2

tk+3

Fig. 1. Asynchronous processing of events.

DFC centre for processing. The observers and the DFC are assumed to be connected by a two-way radio link. Thus after processing the message from the observer, the DFC responds with the next course of action to be taken, represented by the control vector uk ðik Þ for observer ik at time tk . In this study we will neglect the transmission delays of the radio link, the CPU processing time of the DFC and the time it takes the observer to travel from one measurement point to the next. Thus an event when a message is sent to the DFC and the DFC replied to it, is for us instantaneous. Due to the different exposure times of observers, these events will happen asynchronously. An example of the timing axes for two observers is shown in Fig. 1. Note that it is sometimes possible that one observer (which operates with a short exposure time) is more active than the other and communicates more often with the DFC. In our application the control vector tells the observer: (1) where to go next and (2) how long to stay there, in order to acquire the next measurement. The threedimensional control vector is then: i

uk ðik Þ ¼ ½xkk0 0

where k 4k,

zikk0

tikk0 ;

i i ðxkk0 ; zkk0 Þ

ð21Þ

is the next measurement location and

tikk0 is the next exposure time for observer ik. Note that we

use index k0 , rather than k þ 1, to denote the next measurement from observer ik, because of the timing issues described in Fig. 1: the next ðk þ 1Þ st measurement in the sequence may arrive from a different observer. 4.2. Information gain

Information driven search involves using the anticipated information gain of future measurements as a reward function to select the best sensing action. Various measures of information gain have been considered for this task, such as Fisher information, entropy, Kullback–Leibler (KL) divergence, etc. [39]. We used Fisher information in our previous paper [25]. However, since the Fisher information is not defined before the source is detected (and its parameter vector estimated), during the pre-detection phase it was not clear what to do and the search was carried out along a predefined path. The use of the Fisher information in general is unsatisfactory for hybrid estimation problems. A much better solution in this case is to use an information gain formula which involves the entire probability density function. In this paper we adopt for this task the Re´nyi (or alpha) divergence between the current and the future posterior densities. This approach will enable us to control the observers with the maximum information gain throughout the entire search duration, irrespective of the current estimate r^ k of the number of sources [19,39]. The Re´nyi divergence between densities f0 and f1 is defined as [19,39] Z 1 ln f1a ðxÞf01a ðxÞ dx; ð22Þ Da ðf1 jjf0 Þ ¼ a1 where the integral is of the dimension of the state space and aZ0 is a parameter which determines how much we emphasize the tails of two distributions in the metric. In the special cases of a-1 and a ¼ 0:5, the Re´nyi divergence

ARTICLE IN PRESS B. Ristic et al. / Signal Processing 90 (2010) 1225–1239

becomes the Kullback–Leibler divergence and the Hellinger affinity [40], respectively [39]. The reward function used in (20) is taken to be the Re´nyi divergence between the posterior density at time k, f0 ¼ pðr; Xr jz1:k ; u0:k1 ði0:k1 ÞÞ, and the future posterior density, f1 ¼ pðr; Xr jz1:k ; zikk0 ; u0:k ði0:k ÞÞ, computed using the new measurement zikk0 , (with k0 4k) obtained after observer ik has been controlled to take action uk ðik Þ. The time interval tk0  tk ¼ tikk0 is the third component of the control vector uikk , Eq. (21). To reduce the complexity of the notation, in the following the past measurements z1:k and controls u0:k1 ði0:k1 Þ will be omitted. Thus, pðr; Xr Þ denotes the posterior density at time tk and pðr; Xr jzikk0 ; uk ðik ÞÞ denotes the future posterior density. Then from (22) we have the reward function [41]: Dðuk ðik Þ; pðXr ; rÞÞ ¼

1

a1

ln

rmax Z X

pðr; Xr jzikk0 ; uk ðik ÞÞa pðr; Xr Þ1a dx:

r¼0

ð23Þ Expression (23) can be computed relatively simply using the random sample representation of the posterior pdf furnished by the algorithm of Section 3. After using Bayes rule and substituting the particle approximation of the posterior pdf (after resampling and regularisation, hence with equal weights), the Re´nyi divergence between the current and future posterior densities can be approximated by [41] " # g ðzik0 juk ðik ÞÞ 1 ; ð24Þ ln a i k Dðuk ðik Þ; pðXr ; rÞÞ  a  1 g ðz k0 juk ðik ÞÞa 1

k

where, for z 2 N,

ga ðzjuk ðik ÞÞ ¼ 1=N

N X

pðzjXnk ðrkn ÞÞa :

ð25Þ

n¼1

The likelihood of particle n, pðzikk0 jXnk ðrkn ÞÞ was introduced in (3). Note that g1 ðzjuk ðik ÞÞ is the particle approximation to the prior measurement pdf evaluated at z. Recall that we require the expectation of the reward function to select a suitable control vector using (20). This can be approximated as

E½Dðuk ðik Þ; pðXr ; rÞÞ 

  Z2 1 X g ðzjuk ðik ÞÞ : g1 ðzjuk ðik ÞÞ ln a a  1 z¼Z g1 ðzjuk ðik ÞÞa 1

ð26Þ Note that in (26) the theoretically infinite summation over the measurement value z has been reduced to a summation over z ¼ Z1 ; . . . ; Z2. The summation limits Z1 and Z2 are selected such that only values of z for which the summand is significant are considered. 5. Implementation remarks

1231

one-step ahead admissible locations is adopted to be i

i

fðxk þ jr0 cosð‘y0 Þ; zk þ jr0 sinð‘y0 ÞÞ; j ¼ 0; . . . ; Nr  1; ‘ ¼ 0; 1; . . . ; Ny  1g;

ð27Þ

where y0 ¼ 2p=Ny and r0 is a conveniently selected radial step. In this way the observer can stay in its current location ðj ¼ 0Þ, or can move radially in incremental steps. The actual set of admissible locations is a subset of (27): it i i excludes all future locations ðxk0 ; zk0 Þ for which the i i straight-line path from the current location ðxk ; zk Þ steps outside the search region A. In order to reduce the dimension of the control vector optimisation space, for a given new prospective location i i ðxk0 ; zk0 Þ, the exposure time is computed directly. This is done by introducing a maximum ‘‘tolerable’’ value for the i mean radiation count lk ðXr Þ of (4). It can be easily shown [42] that if the Poisson distribution is approximated by an additive Gaussian with the mean and variance equal to lik ðXr Þ, the signal to noise ratio (SNR) in the logarithmic i scale reduces to SNR ¼ 10 log lk ðXr Þ. By fixing the SNR to i SNRmax and substituting lk ðXr Þ with the expression in (4), we can compute tik0 . For this we use the current particle filter estimate of the concatenated vector given in (19). When using the expected Re´nyi divergence as the reward function, the experimental results show that observers are usually directed towards the source. In other words, the Re´nyi divergence acts as an attractive force. While this is in accordance with our intuition, it may create some practical problems in the extreme cases when an observer is too close to a source (recall from Section 2 that the measurement model is not valid for very small observer–source distances). In order to prevent the observers approaching too closely to the radiation sources (and thus avoid prohibitively large radiation counts2), we introduce a repelling force in the reward function. Intuitively, the repelling force should be weak when the radiation is small, and vice versa. We adopt the following definition of the repellant: i

ð28Þ F j ¼ expðj /lkk SÞ; PN ik ik n n where /lk S ¼ ð1=NÞ n¼1 lk ðXk ðrk ÞÞ is the sample mean of the mean radiation count defined in (4); jZ0 in (28) is a user-defined parameter. The Re´nyi divergence is multiplied by the repelling force to give the expected reward for taking action uk ðiÞ as   Z2 Fj X g ðzjuk ðik ÞÞ E½Dðuk ðik Þ; pðXr ; rÞÞ  g1 ðzjuk ðik ÞÞln a a ; a  1 z¼Z g1 ðzjuk ðik ÞÞ 1

ð29Þ which is then used in (20). The search is terminated at time k if any of the following two criteria are satisfied:

 The total search time is greater than a certain interval of time T s ,

First we describe the formation of U k ðiÞ, the admissible action set for observer i at time k. According to (21), the control vector space is continuous-valued and for the sake of implementation on a computer, we will need to discretize it into a finite number of states. If the current i i position of observer i is ðxk ; zk Þ, its provisional set of

k X

ti‘‘ ZT s :

‘¼1

2

This is important when humans act as observers.

ð30Þ

ARTICLE IN PRESS 1232

B. Ristic et al. / Signal Processing 90 (2010) 1225–1239

 For particles characterised by rkn ¼ r^ MAP , the spread of k

short exposures. The red and yellow dots in Fig. 2 indicate the positions of two-source particles, characterised with rkn ¼ 2 (red for one source, yellow for the other source). Fig. 2(a) also contains a few blue dots (not easily seen, because their number is very small) that indicate the locations of particles with rkn ¼ 1. None of the sources has been precisely localised after 45 s, although the EAP estimate r^ k in Fig. 3(c) is close to 2. The algorithm categorically decides that two sources are present after about 70 s. Afterwards, both observers are circulating around respective sources with small exposure times. The parameters used by the search algorithm in simulations are: mb ¼ 1 cnt=s, rmax ¼ 2, N ¼ 2000, Ny ¼ 16, Nr ¼ 3, r0 ¼ 5 m, k ¼ 2, c ¼ 8000, Pb ¼ Pd ¼ 0:003, j ¼ 0:15, SNRmax ¼ 16 dB. The Re´nyi divergence was computed with a ¼ 0:99999 (thus effectively we used the KL divergence). The termination criteria parameters were T s ¼ 350 s and Sp ¼ 1:6 m. The search terminated after 100 s and resulted in both sources being detected and localised: Fig. 2(b) indicates that concentrations of red and yellow particles are formed in vicinity of the true source locations. The final source y^ 1 ¼ 81:3 m, f^ 1 ¼ estimates were: x^ 1 ¼ 25:3 m, 4340 cnt=s for source 1 and x^ 2 ¼ 120:8 m, y^ 2 ¼ 68:3 m, f^ 2 ¼ 8476 cnt=s for source 2.

MAP sources) is particles in the x2y plane (for each of r^ k smaller than a certain threshold Sp .

The irregular shape of the search region A affects both the implementation of the particle filter and the selection of the control vectors. The key problem in placing the particles or selecting the observer action is to determine if a certain point belongs to the interior of A. For this purpose we apply the ray casting algorithm [43] which has been originally developed in computer graphics. To determine if a point ðxp ; yp Þ belongs to the interior of a polygon (possibly with polygonal holes), one casts a horizontal ray emanating from ðxp ; yp Þ to the right. If the number of times this ray intersects the line segments making up the polygon is an even number, then the point is outside the polygon, and vice versa. 6. Simulation results and comparisons 6.1. Simulation set-up and an illustrative run

6.2. Monte Carlo simulations First we study the performance of the search algorithm for different values of parameter a. The simulation set-up and the parameters of the search algorithm are identical to those in Section 6.1, except that a takes values 0.1, 0.5 and 0.99999. The smaller values of a emphasize the tails between the distributions. The results, obtained over 50 Monte Carlo runs, are presented in Table 2. Overall the choice of a makes only a minor impact. Using a ¼ 0:1, the search is the quickest (both the search time and the number of acquired measurements are the smallest), but the least accurate (both in the source position and strength). Using a ¼ 0:5 appears to be the best compromise: parameter estimation is the most accurate, while the search time is only slightly longer than if a ¼ 0:1. The KL divergence ða-1Þ, which is

120

120

100

100

80

80 y [m]

y [m]

The search area A with two radioactive sources (indicated by green asterisks) is shown in Fig. 2. The parameters of source 1 are: x1 ¼ 25 m, y1 ¼ 80 m, f1 ¼ 4000 cnt=s and of source 2 are: x2 ¼ 120 m, y2 ¼ 70 m, f2 ¼ 7000 cnt=s. The search displayed in Fig. 2 involves O ¼ 2 observers entering the search area at locations (30,0) and (45,0) m. Their paths are shown after 45 s and 100 s in Fig. 2(a) and (b), respectively. The colour coding is: pink for observer 1 and cyan for observer 2. The exposure times, measured counts and the current EAP estimates r^ k are shown in Fig. 3(a), (b) and (c), respectively, using the same colour coding. In the simulations the assumption is that the message exchange between the DFC and the observers is instantaneous. Every such message exchange is indicated on the time axis in Fig. 3. Note that the cyan-coloured path is considerably longer in Fig. 2(a). This is due to the much longer exposure intervals of the pink observer in the beginning of the search (see Fig. 3(a)). The strategy of the information driven search appears to be to keep one observer initially fairly static, while the other is traversing the area with

60

60

40

40

20

20

0

0 0

20

40

60 80 100 120 140 x [m]

0

20

40

60 80 x [m]

100 120 140

Fig. 2. The search paths and random samples approximating the posterior pdf after: (a) 45 s and (b) 100 s (when search completed).

ARTICLE IN PRESS

Exposure (s)

B. Ristic et al. / Signal Processing 90 (2010) 1225–1239

1233

20 15 10 5 0

0

10

20

30

40

50

60

70

80

90

100

60

70

80

90

100

60

70

80

90

100

z (counts)

Time [s]

60 40 20 0 0

10

20

30

40

50 Time [s]

Est, rk

2 1.5 1 0.5 0 0

10

20

30

40

50 Time [s]

Fig. 3. The search results (pink for observer 1, cyan for observer 2). (a) Exposure intervals tik . (b) Measured count zk . (c) The EAP estimate of rk . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 2 Search performance for different values of a (50 Monte Carlo runs).

a ¼ 0:1 a ¼ 0:5 a ¼ 0:99999 EAP Average final r^ k RMS error in x (m) RMS error in y (m) RMS error in f (cnt/s) Average search duration (s) Average number of measurements

2

2

2

1.407 0.850 740 87.2 63

0.652 0.772 536 98.0 80

0.858 0.884 1059 97.4 82

Table 3 Search performance: one versus two observers (50 Monte Carlo runs). O¼2 EAP Average final r^ k RMS error in x (m) RMS error in y (m) RMS error in f (cnt/s) Average search duration (s) Average number of measurements

O¼1

2

2

0.652 0.772 536 98 80

0.667 1.268 505 210 112

quite a popular information measure, results in slightly higher estimation errors than the Hellinger affinity ða ¼ 0:5Þ. According to [44], the choice of a ¼ 0:5 stresses the minor differences in the tails between two distributions. While searching for radioactive sources, the future and the current posterior are very similar and the discrimination in tails appears to be the most effective.

Similar results with regards to the choice of a were reported in [45,39]. Next, using the same simulation set-up as before, we compare the search performance of one observer versus two observers. The results obtained from 50 Monte Carlo runs and using a ¼ 0:5, are presented in Table 3. Not surprisingly the most remarkable difference is in the average duration of the search time (see Table 3): it has increased from 98 s (for O ¼ 2) to 210 s (for O ¼ 1). In summary, a single observer uses more than twice of the search time required by a pair of observers. The last set of Monte Carlo simulations compares the performance of the described information driven search with the performance of (1) a uniform-deterministic search along a pre-specified path; (2) a uniform random search. In order to simplify the analysis, we consider the case with a single source, single observer and a rectangular search area. The source of strength f ¼ 3000 cnt=s is placed randomly at each run (but always inside the search area). In every Monte Carlo run, we run the info-driven search first and based on its search duration T and the number of measurements K, we create the corresponding paths for two uniform searches. Both the uniformdeterministic path and the random-walk path cover uniformly the search area, with the number of measurement locations approximately equal to K, and with the exposure time at each measurement location equal to T =K, so that the total search durations are approximately equal. Figs. 4 and 5 show a single run of the info-search and the uniform-deterministic search, respectively

ARTICLE IN PRESS 1234

B. Ristic et al. / Signal Processing 90 (2010) 1225–1239

15

100

10

90

5

80

0

70

0

20

40 60 Time [s]

80

0

20

40 60 Time [s]

80

0

20

40

80

y [m]

60 50

40

40

20

30

0

20 10 0 0

20

40

60

80

100

x [m]

2 1.5 1 0.5 0 60

Time [s] Fig. 4. Info-search. (a) Search area and the final path. (b) Exposure intervals tik . (b) Measured counts zk . (c) EAP estimates r^ k .

8 6

100

4

90

2 0

80 70

0

20

40

60

80

Time [s]

y [m]

60 150

50

100

40

50

30

0

20

0

20

40 60 Time [s]

80

0

20

40 60 Time [s]

80

10 0 0

20

40

60 x [m]

80

100

2 1.5 1 0.5 0

Fig. 5. Uniform search. (a) Search area and the final path. (b) Exposure intervals tik . (b) Measured counts zk . (c) EAP estimates r^ k .

ARTICLE IN PRESS B. Ristic et al. / Signal Processing 90 (2010) 1225–1239

1235

Table 4 Search performance: info-search versus uniform search (50 Monte Carlo runs, single source).

EAP Average final r^ k RMS error in x (m) RMS error in y (m) RMS error in f (cnt/s) Average search duration (s) Average number of measurements

Info-driven search

Uniform deterministic search

Uniform random search

1

1

1

0.75 0.72 165 89.9 20.7

3.55 2.18 595 90.1 20.5

2.72 2.81 521 90.2 20.6

160 140

6

13

5

12

42

120

65

72

64

71

63

70

41

40

100

39 4

Y (m)

80

11 27 49 55 20 3138 54 26 48 19 37 53

60

25 47 52 18 3036 51 24 46 50 17 35

40 20 0 −20

3

10

23 45 16 2934

9 2

22 15

62 69 61 60 68 59 58 67 57 56 66

44 33

8 1

7

−40

−20

21 43 14 2832

0

20

76 80

84

75 79

83

74 78

82

73 77

81

87

91

94

97

93

96

88

92

95

100

120

90 86

40

60 X (m)

89 85

80

140

160

Fig. 6. A top view of the data collection set-up.

(random walk not shown). For a given search duration, the info-search collects more informative measurements and consequently results in more accurate estimates (note the concentration of particles around the true source location marked by green asterisks in Figs. 4 and 5). This observation is confirmed by Monte Carlo averaged results presented in Table 4: the information search is superior to both uniform searches. If we were to require the uniform search to achieve the same level of accuracy as the info-search, the duration of the uniform search would be much longer (due to a longer search path and/or a longer exposure interval). 7. Experimental results A radiological field trial was conducted in Puckapunyal Military Area (Victoria, Australia) on a large, flat, and open area without any obstacles. The satellite image of the area

is shown in Fig. 6. Datasets were collected using the DSTO-developed low cost advanced radiological survey (LCAARS) system in the presence of one and two radiological point sources. An additional dataset was recorded in the absence of any sources, so that the average background radiation can be estimated. The experimental measurements were collected without using the proposed search algorithm, along a predefined path following the yellow circles at enumerated grid points indicated by white crosses, see Fig. 6. The true source locations are indicated by green stars. The data were acquired when the trolley-mounted GM counter was positioned over individual grid points. During data collection at any grid point, the GM counter was held stationarily until approximately sixty measurements, with a constant exposure time of D  1 s, were acquired. Two sets of data were used for the verification of the proposed info-search algorithm. Dataset 1 refers to the

ARTICLE IN PRESS 1236

B. Ristic et al. / Signal Processing 90 (2010) 1225–1239

100 80

y

60 40 20 0 −20 −40 −20

0

20

40

60

80 100 120 140 160

Exp.time [s]

x 20 15 10 5

Estimate rk

z (counts)

0 0

5

10

15

20

25 30 Time [s]

35

40

45

50

0

5

10

15

20

25 30 Time [s]

35

40

45

50

0

5

10

15

20

35

40

45

50

120 100 80 60 40 20 0

1.2 1 0.8 0.6 0.4 0.2 0 25

30

Time [s] Fig. 7. Dataset 1 (single source) final results (single run). (a) Area with the final path (pink line) and the estimate (red dots are particles, green asterisk is the true location). (b) Exposure intervals. (c) Measured counts. (d) EAP estimate of the number of sources. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

set-up with a single 137 Cs source (referred to as source 1) of approximately 26.8 GBq emitting 0.662 MeV, which corresponds to source strength, as defined in this paper, of f1 ¼ 9104:8 cnt=s at a distance of 1 m. The location of this source was x1 ¼ 11 m and y1 ¼ 10 m in the local coordinate system, see Fig. 6. Dataset 2 refers to the set-up with both sources present: source 1 (as above) and source 2 (another 137 Cs source of strength 1868 cnt/s, placed at

x2 ¼ 3 m and y2 ¼ 50 m). More technical details about the experimental set-up are given in [46]. Since the data were acquired along a predefined path with a constant exposure time, in order to be able to apply and verify the proposed info-search algorithm we had to modify the selection of admissible future locations in the control vector. Instead of the procedure described by (27), the admissible locations are the actual

ARTICLE IN PRESS B. Ristic et al. / Signal Processing 90 (2010) 1225–1239

1237

100 80

y [m]

60 40 20 0 −20 −40 −20

0

20

40

60

80 100 120 140 160

x [m]

Exp.time (s)

8 6 4 2 0

0

10

20

30 40 Time [s]

50

60

0

10

20

30 40 Time [s]

50

60

0

10

20

50

60

z (counts)

150 100 50

Est Num Src r

0

2 1.5 1 0.5 0 30

40

Time [s] Fig. 8. Dataset 2 (two sources) final results (single run). (a) Area with the final path (pink line) and the estimate (red and yellow dots are particles, green asterisks are the true location). (b) Exposure intervals; (c) Measured counts. (d) EAP estimate of the number of sources. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

data collection grid points, which fall inside the circle of radius 35 m centered at a current observer location. The controlled exposure time was computed as a multiple d of the quantum of D ¼ 1 s, while a cumulative count measurement is formed by summing up d counts collected at the grid point. In order to avoid a repeated use of the same data, all count measurements that were used during the search are marked as such, so that

they are not re-used if the observer returned to an already visited location. Only a single observer is used on real data. 7.1. Single source results First we used Dataset 1 to verify the algorithm. The parameters of the algorithm were set as follows:

ARTICLE IN PRESS 1238

B. Ristic et al. / Signal Processing 90 (2010) 1225–1239

mb ¼ 0:9 cnt=s (measured independently in the absence of any source), rmax ¼ 1, N ¼ 2000, k ¼ 3, c ¼ 5000, Pb ¼ Pd ¼ 0:003, a ¼ 0:1, j ¼ 0:12, SNRmax ¼ 14:5 dB. The

the higher level of difficulty3 involved in processing Dataset 2.

algorithm terminates when the spread of particles in x2y plane was less than Sp ¼ 1:5 m. Fig. 7(a) shows the search area with the entire path of the observer (pink line), the data collection points (black dots), and the final estimated location of the source (red dots indicate particle locations). The true source location is indicated by the green asterisk. The estimate of the true number of sources r, the measured counts, and the controlled exposure times are shown in Fig. 7(b), (c) and (d), respectively. The observer entered the search area at point ð150; 15Þ m with the exposure time of 2 s. Since the first count measurement at this location was very small (one count only), the observer is controlled to go towards the centre of the search polygon with a large exposure time of about 18 s. The second count measurement, collected at (120, 0), was 40, and the source is at this point of time positively detected (the estimate of r is equal to 1 after 20 s, when the second measurement is processed). In the following four steps, the observer is controlled to go towards the estimated source location, and then to go around it. Finally the source is localised with a positional error of 0.8 m. The total search time was only 52 s (effectively this is the total exposure time, since we ignore the observer travel time). Since the particle filter runs differently every time, we have repeated 20 times the search using Dataset 1 (always with the same entry point). The obtained RMS errors in x and y were 0.06 m and 0.87 m, respectively. The RMS error in f was 3599 cnt/s. The mean search time was 54.8 s. The experimental results thus confirm the validity of the measurement model and the algorithm performance studied by simulations.

8. Conclusions

7.2. Results with two sources

References

Fig. 8 shows a single run of the algorithm using Dataset 2. Note that the measurement locations in Dataset 2 are slightly different from those in Dataset 1. The pink circles indicate the size of a measured count. The parameters of the algorithm were set as follows: mb ¼ 0:9 cnt=s, rmax ¼ 2, N ¼ 3000, k ¼ 1:5, c ¼ 8000, Pb ¼ Pd ¼ 0:003, a ¼ 0:1, j ¼ 0:12, SNRmax ¼ 14 dB. The algorithm terminates when the spread of both red and yellow particles in x  y plane was less than Sp ¼ 2:25 m. The entry point was (120,15) m with the initial exposure time of 2 s. Sources 1 and 2 were localised with positional errors of 2.6 m and 2.8 m, respectively. The total search time was 62 s. As before, we have repeated 20 times the search using Dataset 2 (with the same entry point). The obtained RMS errors in x and y for source 1 were 1.94 m and 2.18 m, respectively. RMS errors in x and y for source 2 were 3.13 m and 1.70 m, respectively. The RMS errors in strength f were 2432 cnt/s and 1074 cnt/s, for sources 1 and 2, respectively. The mean search time was 93.4 s. Note that the estimation errors and the mean search time are somewhat higher than with Dataset 1, reflecting

The paper presented a probabilistic information driven search algorithm which jointly performs the estimation of the number of gamma-radiation point sources and their parameters (locations and strengths). The search area is an open-field polygon-shaped area, where measurements are collected by observers equipped with low-cost ¨ Geiger–Muller counters and precise self-localisation units. By controlling the observers, the collected measurements are more informative in comparison with a predefined path search, thus effectively reducing the search time. The observer control is done sequentially based on the current estimate of the posterior probability density function and by computing the information gain (in the Re´nyi divergence sense) of the candidate future measurements. Successful application of the proposed technique to two sets of experimental data verifies the measurement model and simulation studies. Future work will consider situations where the search area is populated by obstacles. It is anticipated that some prior knowledge of the map of the area, the likely types of radioactive sources and obstacle materials will be required. Multiple-step ahead sensor management is expected to be beneficial in this case.

Acknowledgements The authors would like to thank their DSTO colleagues A. Skvortsov, R. Gailis, A. Eleftherakis, D. Marinaro, P. Harty and I. Mitchell for useful technical discussions and involvement in the field trail.

[1] W.K.H. Panofsky, Nuclear proliferation risks, new and old, Issues in Science and Technology 19 (2003) /http://www.issues.org/19.4/ panofsky.htmlS. [2] IAEA, Trafficking in nuclear and radioactive material in 2005 /http://www.iaea.org/NewsCenter/News/2006/trafficking stats2005.htmlS, August 2006. [3] J.W. Howse, L.O. Ticknor, K.R. Muske, Least squares estimation techniques for position tracking of radioactive sources, Automatica 37 (2001) 1727–1737. [4] T. Hjerpe, R.R. Finck, C. Samuelsson, Statistical data evaluation in mobile gamma spectrometry: an optimisation of on-line search strategies in the scenario of lost point sources, Health Physics 80 (6) (2001) 563–570. [5] D.M. Nicol, R. Tsang, H. Ammerlahn, M.M. Johnson, Sensor fusion algorithms for the detection of nuclear material at border crossings, in: Proceedings of the SPIE, vol. 6201, 2006. [6] N. Rao, M. Shankar, J.C. Chin, D.K.Y. Yau, C.Y.T. Ma, Y. Yang, J.C. Hou, X. Xu, S. Sahni, Localization under random measurements with application to radiation sources, in: Proceedings of the International Conference on Information Fusion, Cologne, Germany, 2008. [7] R.J. Nemzek, J.S. Dreicer, D.C. Torney, T.T. Warnock, Distributed sensor networks for detection of mobile radioactive sources, IEEE Transactions on Nuclear Science 51 (4) (2004) 1693–1700.

3 Recall that the observer is typically tasked to move around a source. With Dataset 2 this was not possible for source 2, because of a very restrictive set of measurement locations, see Fig. 8(a).

ARTICLE IN PRESS B. Ristic et al. / Signal Processing 90 (2010) 1225–1239

[8] S.M. Brennan, A.M. Mielke, D.C. Torney, Radioactive source detection by sensor networks, IEEE Transactions on Nuclear Science 52 (3) (2005) 813–819. [9] D.L. Stephens, A.J. Peurrung, Detection of moving radioactive sources using sensor networks, IEEE Transactions on Nuclear Science 51 (5) (2004) 2273–2278. [10] A. Sundaresan, P.K. Varshney, N.S.V. Rao, Distributed detection of a nuclear radioactive source using fusion of correlated decisions, in: Proceedings of the 10th International Conference on Information Fusion, Que´bec, Canada, July 2007. [11] M. Morelande, B. Ristic, A. Gunatilaka, Detection and parameter estimation of multiple radioactive sources, in: Proceedings of the 10th International Conference on Information Fusion, Que´bec, Canada, July 2007. [12] B.O. Koopman, Search and Scanning, Pergamon Press, Oxford, 1980. [13] K.P. Ziock, W.H. Goldstein, The lost source, varying backgrounds and why bigger may not be better, in: J.I. Trombka, D.P. Spears, P.H. Solomon (Eds.), Unattended Radiation Sensor Systems for Remote Applications, CP632, American Institute of Physics, 2002. [14] A.V. Klimenko, W.C. Priedhorsky, N.W. Hengartner, K.N. Borozin, Efficient strategies for low-level nuclear searches, IEEE Transactions on Nuclear Science 53 (3) (June 2006) 1435–1442. [15] J.M. Passerieux, D. Van Cappel, Optimal obsever maneuver for bearings-only tracking, IEEE Transactions on Aerospace and Electronic Systems 34 (3) (1998) 777–788. [16] J.-P. Le Cadre, S. Laurent-Michel, Optimizing the receiver maneuvers for bearings-only tracking, Automatica 35 (4) (1999) 591–606. [17] S. Blackman, R. Popoli, Design and Analysis of Modern Tracking Systems, Artech House, 1999. [18] A. Doucet, B. Vo, C. Andrieu, M. Davy, Particle filtering for ;multitarget tracking and sensor management, in: Proceedings of the 5th International Conference on Information Fusion,vol. 1, Annapolis, Maryland, USA, 2002, pp. 474–482. [19] C.M. Kreucher, A.O. Hero, K.D. Kastella, M.R. Morelande, An information based approach to sensor management in large dynamic networks, Proceedings of the IEEE 95 (5) (May 2007) 978–999. [20] K. Kastella, Discrimination gain to optimize classification, IEEE Transactions on Systems, Man, and Cybernetics—Part A 27 (1) (1997) 112–116. [21] M.A. Sipe, D. Casasent, Feature space trajectory methods for active computer vision, IEEE Transactions on Pattern Analysis and Machine Intelligence 24 (12) (2002) 1634–1643. [22] J. Denzer, C.M. Brown, Information theoretic sensor data selecetion for active object recognition and state estimation, IEEE Transactions on Pattern Analysis and Machine Intelligence 24 (2) (2002) 145–157. [23] T. Furukawa, F. Bourgault, B. Lavis, H.F. Durrant-Whyte, Recursive Bayesian search-and-tracking using coordinated UAVs for lost targets, in: Proceedings of the IEEE International Conference on Robotics and Automation, Orlando, FL, USA, May 2006, pp. 2521–2526. [24] R.A. Cortez, X. Papageorgiou, H.G. Tanner, A.V. Klimenko, K.N. Borozdin, R. Limia, W.C. Priedhorsy, Smart radiation sensor management, IEEE Robotics & Automation Magazine, September 2008, pp. 85–93. [25] B. Ristic, A. Gunatilaka, Information driven localisation of a radiological point source, Information Fusion 9 (2) (April 2007) 317–326.

1239

[26] D.J. Salmond, H. Birch, A particle filter for track-before-detect, in: Proceedings of the American Control Conference, Arlington, VA, June 2001, pp. 3755–3760. [27] Y. Boers, H. Driessen, Multitarget particle filter track before detect application, IEE Proceedings of the Radar, Sonar and Navigation 151 (6) (December 2004) 351–357. [28] N. Tsoulfanidis, Measurement and Detection of Radiation, Taylor & Francis, Washington, DC, 1995. [29] A. Martin, S.A. Harbison, An Introduction to Radiation Protection, Chapman & Hall, London, 1987. [30] H.E. Johns, Physics of Radiology, fourth ed., Charles C Thomas, 1983. [31] E.K.P. Chong, C. Kreucher, A.O. Hero, POMDP approximation using simulation and heuristics, in: A.O. Hero, D. Castanon, D. Cochran, K. Kastella (Eds.), Foundations and Applications of Sensor Management, Springer, Berlin, 2008 (Chapter 8). [32] B. Ristic, S. Arulampalam, N. Gordon, Beyond the Kalman Filter: Particle Filters for Tracking Applications, Artech House, 2004. [33] J. Czyz, B. Ristic, B. Macq, A particle filter for joint detection and tracking of color objects, Image and Vision Computing 25 (2007) 1271–1281. [34] Y. Bar-Shalom, X.R. Li, T. Kirubarajan, Estimation with Applications to Tracking and Navigation, Wiley, New York, 2001. [35] C.P. Robert, G. Casella, Monte Carlo Statistical Methods, Springer, Berlin, 2004. [36] C. Musso, N. Oudjane, F. LeGland, Improving regularised particle filters, in: A. Doucet, N. deFreitas, N.J. Gordon (Eds.), Sequential Monte Carlo Methods in Practice, Springer, New York, 2001. [37] F. Daum, J. Huang, Particle flow for nonlinear filters with loghomotopy, in: O.E. Drummond (Ed.), Proceedings of the SPIE, vol. 6969, 2008. [38] G. Kitagawa, Monte Carlo filter and smoother for non-Gaussian non-linear state space models, Journal of Computational and Graphical Statistics 5 (1) (1996) 1–25. [39] A.O. Hero, D. Castanon, D. Cochran, K. Kastella, Foundations and Applications of Sensor Management, Springer, Berlin, 2007. [40] D. Pollard, A User’s Guide to Measure Theoretic Probability, Cambridge University Press, Cambridge, 2001. [41] C.M. Kreucher, K.D. Kastella, A.O. Hero, Sensor management using an active sensing approach, Signal Processing 85 (2007) 607–624. [42] A. Gunatilaka, B. Ristic, R. Gailis, On localisation of a radiological point source, in: Proceedings of the Information, Decision and Control (IDC), Adelaide, Australia, February 2007. [43] J.D. Foley, Computer Graphics: Principles and Practice, AddisonWesley, Reading, MA, 1995. [44] A.O. Hero III, B. Ma, O. Michel, J.D. Gorman, Alpha divergence for classification, indexing and retrieval, Technical Report 328, Communications and Signal Processing Lababoratory (CSPL), Department of EECS, University of Michigan, 2001. [45] C.M. Kreucher, An information-based approach to sensor resource allocation, Ph.D. Thesis, The University of Michigan, 2005. [46] A. Gunatilaka, B. Ristic, LCAARS radiological field trial and validation of source localisation algorithms, Technical Report, DSTO, Melbourne, Australia, 2009.