Rao–Blackwell dimension reduction applied to hazardous source parameter estimation

Rao–Blackwell dimension reduction applied to hazardous source parameter estimation

Signal Processing 132 (2017) 177–182 Contents lists available at ScienceDirect Signal Processing journal homepage: www.elsevier.com/locate/sigpro c...

617KB Sizes 0 Downloads 32 Views

Signal Processing 132 (2017) 177–182

Contents lists available at ScienceDirect

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

crossmark

Rao–Blackwell dimension reduction applied to hazardous source parameter estimation ⁎

Branko Ristica, , Ajith Gunatilakab, Yan Wangc a b c

RMIT University, School of Engineering, Melbourne, VIC3000, Australia Defence Science and Technology Group, Land Division, Melbourne, Australia Tsinghua University, Institute of Public Safety Research, Department of Engineering Physics, Beijing, China

A R T I C L E I N F O

A BS T RAC T

Keywords: Bayesian estimation Chemical–biological–radiological defence Monte Carlo estimation Plume dispersion

Parameter estimation of a source of chemical, biological or radiological emissions is a problem of great importance for public safety. The key parameters of interest are the source intensity and its location. This paper applies the concept of Rao–Blackwell dimension reduction to solve the posterior probability distribution function of source intensity, conditioned on source location, analytically. The paper is cast in the context of a source of a hazardous release of particles or gas and its turbulent transport through the medium. Numerical results, obtained by simulations and using an experimental dataset, demonstrate the statistical efficiency of the proposed method.

1. Introduction The threat of a hazardous chemical–biological–radiological (CBR) attack, either in a form of a release of toxic biochemical aerosols into the atmosphere or an improvised nuclear device, has been well documented [1,2]. For the sake of public safety, it is very important to rapidly detect and localise the CBR source so that the mitigation actions can be carried out promptly. The problem of CBR source localisation has been studied for quite some time. The standard solutions are based on optimisation techniques, such as nonlinear least squares [3,4] and maximum likelihood estimation [5]. These approaches can fail due to local minima or poor convergence and in addition provide only point estimates without uncertainty attached to it. The alternatives are the Bayesian techniques [6–12], which estimate the posterior probability density function (PDF) of a source, thereby providing an uncertainty measure to any point estimate derived from it. Among the Bayesian techniques, Markov chain Monte Carlo (MCMC) is the dominant method for estimation of the posterior PDF, applied in [6,8–10]. Other approaches to posterior estimation are numerical integration [7], importance sampling [11], and approximate Bayesian computation [12]. In formulating the problem, one needs to specify the parameter vector, the dispersion model and the measurement model. The parameter vector θ typically includes the source release rate or intensity, Q0 ∈ + , the position of the source r0 ∈ d (d=2 or 3), and possibly environmental parameters (e.g. the mean wind speed and canopy



Corresponding author. E-mail address: [email protected] (B. Ristic).

http://dx.doi.org/10.1016/j.sigpro.2016.10.005 Received 3 March 2016; Received in revised form 22 July 2016; Accepted 10 October 2016 Available online 12 October 2016 0165-1684/ © 2016 Elsevier B.V. All rights reserved.

characteristics). The dispersion or propagation model describes via mathematical equations the mean concentration of an emitted biochemical substance, or the mean radiation field, at a given sensor location, as a function of the parameter vector θ . All dispersion/ propagation models used in practice are nonlinear with respect to the source location r0 , but linear with respect to intensity Q0. This linear relationship is the key feature exploited in the paper. Finally, the measurement model relates the mean concentration to sensor measurements that are affected by stochastic fluctuations. The most adequate model for this purpose has been experimentally found to be the Poisson distribution – either in the context of measuring the count of the radiated photons [13] or the count of dispersed biochemical particles [14]. The source parameter estimation problem is approached in the Bayesian framework. The main novelty of the paper, compared to the earlier Bayesian approaches, is that it applies the concept of Rao– Blackwell dimension reduction [15,16] to the problem. In doing so, we derive the analytic expressions for: (a) the posterior PDF of source intensity Q0 (conditioned on the source location) and (b) the likelihood of the source location parameter (only). Numerical analysis, using both simulated and experimental data, is carried out to demonstrate the gains of the proposed, dimension reduced, estimation method. The presentation of the paper, without loss of generality, will focus on a hazardous biochemical substance release into the atmosphere, using an open-field two-dimensional (d=2) dispersion model.

Signal Processing 132 (2017) 177–182

B. Ristic et al.

dimension of the parameter vector space. Using the chain rule one can write the posterior PDF as follows:

2. Models and problem formulation We adopt a dispersion model of turbulent transport through the medium from [14], used in a number of recent publications [17–21]. Consider a source of particle release into the environment (atmosphere and water), characterised by a constant emission rate Q0 > 0 and located at r0 = (x 0 , y0)⊤ . The particles propagate through the medium with the isotropic diffusivity D, but can also be advected by flow (due to wind or current), whose mean direction and average speed V are known. Let us adopt the convention that the mean flow (wind) direction coincides with the direction of the x-axis. The average lifetime of a particle (before being absorbed) is τ. A spherical sensor of small size a at a location with coordinates r = (x, y)⊤, non-coincidental with the source location r0 , will experience a series of encounters with the released particles at the rate [14]:

R θ (r) =

⎡ (x − x ) V ⎤ ⎛ (x − x 0 )2 + ( y − y0)2 Q0 exp ⎢ 0 ⎥ ·K0 ⎜⎜ ⎛λ⎞ ⎣ ⎦ ⎝ 2D λ ln ⎜ ⎟ ⎝a⎠

Suppose for a moment that we can calculate the posterior p (Q0 |r0, z) analytically. Then we need to apply a Monte Carlo estimation method to compute only p (r0 |z), which according to Bayes’ rule is given by:

p (r0 |z) =

Dτ . V 2τ 1+ 4D

⎞ ⎟ ⎟ ⎠ (1)

π (Q0 ) = . (Q0 ; η0 , ϑ0) =

μi

e−μi

ℓ(z|θ ) π (θ ) . ∫ ℓ(z|θ ) π (θ ) dθ

Q0(η0 −1) e−Q0 /ϑ0 . ϑ 0η0 Γ (η0 )

∏ 7 (zi; Q0·ρr0 (ri)) i =1

(7)

(8)

where

ρr0 (ri) = t0 R θ (ri)/ Q0

(9)

is independent of Q0. This is a consequence of the dispersion model being linear1 with respect to Q0. Recall that the conjugate prior of the Poisson distribution is the Gamma distribution [22]. Therefore, we can expect the posterior p (Q0 |r0, z) to be also a Gamma distribution, p (Q0 |r0, z) = . (Q0 ; η, ϑ), with parameters η and ϑ that can be calculated analytically as a function of r0 and z . Note that the fact that the likelihood is a product of Poisson distributions (rather than a single Poisson) does not change the scheme, because one can think of this product as an arbitrary order sequence of updates of the Gamma distributed random variable with Poisson distributed measurements, which results in a sequence of Gamma distributed random variables. Proposition: The parameters η and ϑ of the posterior p (Q0 |r0, z) = . (Q0 ; η, ϑ) can be calculated as follows:

(3)

where μi = R θ (ri) t0 is the mean concentration at ri and i = 1, …, S . All sensors are assumed to be of the same (and known) size a. Because the environmental parameters τ, D and V, are also known, (3) represents the full specification of the likelihood function of θ , given a measurement zi (taken at location ri ). Assuming the sensor measurements, conditioned on θ , are independent, the likelihood function of the measurement vector z = [z1, z2, …, zS ]⊤ can be written as a product S ℓ(z|θ ) = ∏i =1 7 (zi ; t0 R θ (ri)). The parameter estimation problem is formulated in the Bayesian framework. The goal is to compute the posterior PDF p (θ|z), which provides a complete probabilistic description of the information contained in z about θ . To compute the posterior distribution, in addition to ℓ(z|θ ) one needs to specify the prior distribution of the parameter vector π (θ ). Using Bayes' rule

p (θ|z) =

(6)

S

ℓ(z|θ ) =

zi

zi !

.

For suitably chosen hyperparameters η0 and ϑ0, this prior can be diffuse, with the support covering a large span of possible values of Q0. Recall that the likelihood function ℓ(z|θ ) is a product of Poisson distributions, which we can write in a slightly different form as:

(2)

The rate of particle encounters R θ (r) is expressed in (1) as a function of the (unknown) parameter vector θ = [Q0 , r⊤0 ]⊤ . Suppose a network of S spatially distributed sensors is measuring the concentration of emitted particles. The stochastic process of sensor encounters with emitted particles is modelled by a Poisson distribution: the probability that ith sensor at location ri = (xi , yi)⊤ encounters zi ∈ + particles during a time interval t0 is then:

7 (zi ; μi ) =

g (z|r0) π (r0)

∫ g (z|r0) π (r0) d r0

The problem with (6) is that the likelihood function g (z|r0) is also unknown – only ℓ(z|θ ) is known. Hence, in order to apply the Rao– Blackwell dimension reduction, we need not only the analytic expression for p (Q0 |r0, z), but also for g (z|r0). Due to the independence assumption between the source intensity and its location, π (θ ) = π (Q0 ) π (r0). Let us assume the prior π (Q0 ) is a Gamma distribution with shape parameter η0 and scale parameter ϑ0, that is

where K0 is the modified Bessel function of order zero and

λ=

(5)

p (θ|z) = p (Q0 |r0, z) p (r0 |z)

S

η = η0 +

∑ zi ,

(10)

i =1

(4)

ϑ=

ϑ0 S

1 + ϑ0 ∑i =1 ρr0 (ri)

Quantities of interest related to θ (e.g., the posterior mean and variance) can be computed from p (θ|z). Note that the prior π (θ ) is typically non-Gaussian: the source position is often restricted to polygon regions, while Q0 is strictly positive. Optimal Bayesian estimation is generally impossible because the posterior PDF cannot be found in closed-form. This is certainly the case for the signal model described above.

. (11)

Proof. The proof is based on two properties of the Gamma distribution: (i) If X ∼ . (η, θ ) then for any constant c > 0 , cX ∼ . (η, cθ ) [23]. (ii) If . (μ; η, θ ) is a prior distribution and n is a sample from the Poisson distributed likelihood function with parameter μ, then the posterior is [22] . (μ; η + n, θ /(1 + θ )). Consider for simplicity only the first measurement, z1, collected at position r1. The likelihood function ℓ(z1 |θ ) is Poisson distributed with mean μ1 = Q0 ρr0 (r1), where Q0 ∼ . (Q0 ; η0 , ϑ0). According to Properties (i) and (ii),

3. Rao–Blackwell dimension reduction This section presents the key result: if the source intensity and its location are independent, the posterior PDF of source intensity, conditioned on the source location, can be calculated analytically. Consequently, Monte Carlo estimation can be applied to a reduced

1 All commonly used dispersion/propagation models are linear with respect to Q0. While in the paper we adopt the model from [14], the method is universal.

178

Signal Processing 132 (2017) 177–182

B. Ristic et al.

⎛ ρr0 (r1)ϑ0 ⎞ ⎟⎟ . μ1 ∼ . ⎜⎜μ; η0 + z1, 1 + ρr0 (r1)ϑ0 ⎠ ⎝

previous stage as an importance distribution. A target distribution which can be used in this context at iteration = = 1, …, H of ISPC is:

If we again use property (i), the posterior distribution of Q0 (after processing z1 collected at r1) becomes:

p= (r0 |z) ∝ [g (z|r0)]Γ= π (r0)

with γj ∈ (0, 1] and ΓH = 1. Because Γ= is an where Γ= = increasing function of ℏ, upper bounded by one, the intermediate likelihood is broader than the true likelihood, particularly in the early stages. Thus, the sequence of target distributions in (17) gradually introduces the correction imposed by the measurement z on the prior π (r0). To derive any benefits from ISPC, it is required after each stage to remove the lowly weighted members of the sample and diversify the remaining ones. Lowly weighted members are removed by resampling [28], while diversification of the remaining samples is performed by Markov transitions whose stationary distribution is the target distribution p= (r0 |z). The outcome is a diverse sample located in the region of the parameter space where the intermediate likelihood has non-zero values. The choice of correction factors γ1, …, γH as well as the number of iterations H are design parameters. We choose γ= parameters adaptively in such a manner that the effective sample size at each iteration is greater than a certain threshold. The pseudo-code of the Rao–Blackwellised (RB) ISPC algorithm, given in Algorithm 1, is best explained step-by-step. The initial source location sample at ISPC stage = = 1, that is {r(0n,1) }1≤ n ≤ N , is drawn from the prior in line 2; N is the sample size. The likelihood of this initial draw is computed in line 3. The while-loop between the lines 5 and 24 carries out the progressive correction of the source location samples. The correction factors γ= are computed adaptively, by ensuring the effective sample size, calculated in line 9, is greater than a threshold χ. The lines 16–19 implement the resampling step, followed by diversification via regularisation of samples [24] in line 20. This step simply jitters each source position sample with ΩDh ϵn , where ϵn is a twodimensional draw from the standard normal distribution, Dh is a square root of the sample covariance matrix of {r(0n, = +1) }1≤ n ≤ N , and Ω is the optimal bandwidth of the Gaussian kernel, see [28, Eq. (3.52)]. The last stage of RB-ISPC, = + 1, is renamed H in line 25. This is followed by the computation of the parameters of the posterior p (Q0 |r 0n, H ), z) in lines 26 and 27. The output of the algorithm is the sample at the final stage of RB-ISPC (see line 28). For the sake of error performance analysis, we also calculate the point estimates as the N 1 l0 = 1 ∑N η(n, H ) ·ϑ(n, H ) . sample means: ^r0 = (xl0 , yl0)⊤ = N ∑n =1 r(0n, H ) ; Q n =1 N

⎛ ⎞ ϑ0 ⎟⎟ , Q0 |r0, z1 ∼ . ⎜⎜Q0 ; η0 + z1, 1 + ρ ( r )ϑ ⎝ r0 1 0 ⎠ which corresponds to (10) and (11) with S=1. By repeating the sequence of update steps using the remaining measurements z2, …, zS , it is easy to verify that (10) and (11) hold for any S.□ The analytic expression for the likelihood g (z|r0), which features in (6), is derived next. Note first that the posterior PDF p (Q0 |r0, z), according to Bayes's rule, is:

p (Q0 |r0, z) =

ℓ(z|Q0 , r0) π (Q0 ) g (z|r0)

(12)

where

g (z|r0) =

∫ ℓ(z|Q0, r0) π (Q0) dQ0

(13)

Based on the material presented above, the analytic expressions for p (Q0 |r0, z), ℓ(z|Q0 , r0) and π (Q0 ) are already available. Hence we can write from (12), using (8), (10) and (11):

g (z|r0) =

π (Q0 )ℓ(z|Q0 , r0) p (Q0 |r0, z)

(14)

S

g (z|r0) =

. (Q0 ; η0 , ϑ0) ∏i =1 7 (zi ; Q0 ·ρr0 (ri)) . (Q0 ; η, ϑ)

(15)

where η and ϑ were given by (10) and (11), respectively. Upon substitution of the expressions for Gamma and Poisson distributions in (15), one can show that Q0 cancels out, resulting in S

S ∑i =1 zi Γ (η0 + ∑i =1 zi )

g (z|r0) = ϑ 0

Γ (η0 )

S

∏ i =1

[ρr0 (ri)]zi zi !

(17)

= ∑ j =1 γj

. (16)

Now that we have derived the analytic expressions for p (Q0 |r0, z) and g (z|r0), we can apply a Monte Carlo method of choice to estimate the posterior p (r0 |z).

Algorithm 1. RB-ISPC algorithm for CBR source estimation.

4. Rao–Blackwellised importance sampling estimator

1: Input: π (r0), η0, ϑ0, z , N

The posterior PDF p (r0 |z) of (6) needs to be computed numerically. Various Monte Carlo methods are applicable and we adopt a technique, similar to the one reported in [11], referred to as importance sampling with progressive correction (ISPC) [24]. The essence of the ISPC method, explained in detail in [11], has been proposed under different names, e.g. tempering [25, p. 540], annealed importance sampling [26], and bridging densities [27]. A desirable property of an importance density to be used for importance sampling is that it produces weights with a very small variance. Ideally, the importance weights should be even, which can be achieved only if the importance density equals the (unknown) posterior. Importance sampling with the prior PDF is useful only if there is a significant overlap between the likelihood and the prior PDF. In this case the posterior is reasonably similar to the prior and the importance weights are not highly variable. The idea behind ISPC is to mimic this situation in a sequential procedure by constructing a series of target distributions from which to draw samples. The first target distribution should be similar to the prior while the final target distribution is the posterior. The series of target distributions should be such that consecutive target distributions in the series do not differ too greatly. This enables an importance sample with reasonably even weights to be obtained from the target distribution at a particular stage using the target distribution at the

2: Draw r(0n,1) ∼ π (r0) for n = 1, …, N tion

▹ Prior PDF of source loca-

3: Calculate g (z|r(0n,1) ), Eq. (16), for n = 1, …, N 4: Set γ1 = 1, Γ0 = 0 , = = 0 5: while Γ= < 1 do 6: = = = + 1 ∼(n, = ) = [g (z|r(n, = ) )]γ= for n = 1, …, N 7: Calculate weights w 0 8: 9: 10: 11:

leff = [∑N (w (n, = ) )2 ]−1 Compute N k n =1 leff < χ do while N γ= = γ= /2

12:

∼(n, = ) = [g (z|r(n, = ) )]γ= for n = 1, …, N Calculate weights w 0

13:

−1 ∼(n, = )·⎛⎜∑N w ∼( j, = ) ⎞⎟ , for n = 1, …, N Normalisation: w(n, = ) = w j =1 ⎠ ⎝

14:

leff = [∑N (w (n, = ) )2 ]−1 Compute N k n =1

15: 16: 179

−1 ∼(n, = )·⎛⎜∑N w ∼( j, = ) ⎞⎟ , for n = 1, …, N Normalisation: w(n, = ) = w j =1 ⎠ ⎝

end while for n = 1, …, N do

▹ Resampling

Signal Processing 132 (2017) 177–182

B. Ristic et al.

40

Table 1 RMS errors in source parameter estimation as a function of the sample size N: Comparison between the RB-ISPC and ISPC, using the same threshold χ = 5. Two sets of hyperparameters of the RB-ISPC applied: (a) η0 = 3, ϑ 0 = 5.2 and (b) η0 = 4 , ϑ 0 = 4.5. The theoretical lower bounds (based on the Cramer–Rao analysis) are: σ xl0 = 2.10 , σ yl0 = 1.42 , σ Ql0 = 1.49 .

12

30 10

20 10

xl0

8

y

0

l0 Q

yl0

RB-ISPC

ISPC

RB-ISPC

ISPC

RB-ISPC

ISPC

4.24 2.90 2.37 2.16 2.09

12.63 9.47 6.57 4.43 3.07

2.10 1.82 1.70 1.64 1.62

8.58 6.09 4.04 2.94 2.14

1.49 1.49 1.50 1.49 1.48

5.58 4.55 3.46 2.85 2.25

6

-10 -20

(a) N=25 N=50 N=100 N=200 N=400

4

-30 2

-40 -50 -50

-40

-30

-20

-10

0

x

10

20

30

xl0

40

Fig. 1. A top-down view of the simulation setup for a source placed at x 0 = −14 , y0 = 2

(b) N=25 N=50 N=100 N=200 N=400

with Q0 = 10 : the mean concentration rate R θ (r) shown as a gray-scale intensity plot; the locations of sensors marked by red circles, whose radius is proportional to the corresponding concentration measurement. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

RB-ISPC

ISPC

RB-ISPC

ISPC

RB-ISPC

ISPC

4.26 2.99 2.43 2.19 2.16

11.92 8.09 5.84 3.89 2.79

2.07 1.84 1.75 1.64 1.64

7.73 5.49 3.65 2.70 2.05

1.49 1.49 1.50 1.50 1.50

5.61 4.38 3.46 2.71 2.21

Select index j n ∈ {1, …, N} with probability wk(n, = )

17: 18: 19:

r(0n, = +1) = r(0jn, = ) end for

20:

Regularisation of source position particles {r(0n, = +1) }1≤ n ≤ N

Table 2 RMS errors in source parameter estimation using the RB-ISPC algorithm, as function of the sample size N: the influence of threshold χ. Two sets of hyperparameters of the RBISPC applied: (a) η0 = 3, ϑ 0 = 5.2 and (b) η0 = 4 , ϑ 0 = 4.5. The theoretical lower bounds (based on the Cramer–Rao analysis) are: σ xl0 = 2.10 , σ yl0 = 1.42 , σ Ql0 = 1.49 .

21: Calculate g (z|r(0n, = +1) ), Eq. (16), for n = 1, …, N 22: Γ= = Γ= −1 + γ= ; 23: γ= +1 = 1 − Γ= 24: end while 25: H = = + 1 S

26: η(n, H ) = η0 + ∑i =1 zi , for n = 1, …, N 27: ϑ(n, H )

xl0

(a) N=25 N=50 N=100 N=200 N=400

▹ Eq. (10)

⎞−1 ⎛ S = ϑ0 ⎜1 + ϑ0 ∑i =1 ρr(n, H ) (ri) ⎟ , for n = 1, …, N 0 ⎠ ⎝



Eq. (11) 28: Output: {(r(0n, H ), η(n, H ) , ϑ(n, H ) )}1≤ n ≤ N

(b) N=25 N=50 N=100 N=200 N=400

5.1. Simulations Fig. 1 shows the simulation setup, with the source of intensity2 Q0 = 10 placed at x 0 = −14 , y0 = 2 . The corresponding mean concentration rate (the rate of particle encounters) R θ (r) is displayed in Fig. 1 as a gray-scale intensity plot (darker areas indicate higher concentration rate). The dispersion model parameters are τ = 250 , D=1, V=0.4, the sensor size is a=1 and the integration time is adopted as t0 = 1. The locations of S=63 sensors used for collecting measurements are indicated in Fig. 1 as the centres of red circles. The sizes of circles correspond to one particular set of concentration measurements (the larger the circle, the higher the value of zi, i = 1, …, S ). Table 1 compares the error performance of the proposed RB-ISPC algorithm against the full-parameter ISPC algorithm (applied to the entire three-dimensional space of the parameter vector θ ). For completeness, the pseudo-code of this algorithm is given in the Supplementary material. Each entry in Table 1 represents an RMS error of the estimated source parameter, obtained from 5000 Monte

l0 Q

yl0

χ=5

χ = 10

χ=5

χ = 10

χ=5

χ = 10

4.24 2.90 2.37 2.16 2.09

2.85 2.22 2.10 2.09 2.10

2.10 1.82 1.70 1.64 1.62

1.77 1.61 1.56 1.55 2.51

1.49 1.49 1.50 1.49 1.48

1.49 1.49 1.50 1.49 1.49

xl0

5. Numerical results

2

l0 Q

yl0

l0 Q

yl0

χ=5

χ = 10

χ=5

χ = 10

χ=5

χ = 10

4.26 2.99 2.43 2.19 2.16

2.63 2.20 2.14 2.08 2.09

2.07 1.84 1.75 1.64 1.64

1.78 1.62 1.61 1.55 1.53

1.49 1.49 1.50 1.50 1.50

1.49 1.48 1.50 1.50 1.50

Carlo runs. The number of samples N is varied from 25 to 400. Other algorithm parameters (identical for RB-ISPC and ISPC) were: the prior PDF π (r0) is the square area shown in Fig. 1, i.e. π (x 0 ) = π ( y0) = < [−50, 50]; and threshold χ = 5 (line 10 of Algorithm 1). The hyperparameters of the RB-ISPC are selected as: (a) η0 = 3, ϑ0 = 5.2 and (b) η0 = 4 , ϑ0 = 4.5. The source intensity samples are drawn in ISPC from the prior π (Q0 ) with these hyperparameters. The theoretical lower bounds of source parameter estimation errors have been computed using the Mathematica software (the source code is given in the Supplementary material). These bounds are important for the assessment of the in silico results of Table 1. The values of the bounds for both cases of the hyperparameters, i.e. (a) η0 = 3, ϑ0 = 5.2 and (b) η0 = 4 , ϑ0 = 4.5, are equal, because the Gamma distributions with these sets of parameters have equal variances. The resulting

All physical quantities in Section 5.1 are in arbitrary units (a.u).

180

Signal Processing 132 (2017) 177–182

B. Ristic et al.

500

1.8

Prior PDF Posterior PDF

1.6

400

1.4

y [mm]

300

1.2

200

1

100

0.8 0.6

0

0.4

−100

0.2

−200

0

0

2

4

−300

6

8

10

12

14

16

18

20

Q0 [1/s]

−400 −500 −500

0

Fig. 4. RB-ISPC applied to the experimental dataset: the prior PDF π (Q0 ) plotted in black; the estimated (final) posterior PDF of Q0 plotted in green. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

500

x [mm]

prominent for smaller values of N. (iv) There is a small difference between the results in Table 1(a) and (b), obtained using two different sets of hyperparameters. (v) The source intensity RMS error, obtained using the RB-ISPC, approximately equals the theoretical low bound, independent of the number of samples N. Next we analyse the influence of the threshold χ on the estimation error performance. Table 2 presents the RMS errors of the RB-ISPC obtained with χ = 5 and χ = 10 , using the same setup and the parameters as in Table 1. From Table 2 one can observe somewhat better results (smaller RMS errors) for χ = 10 , than for χ = 5. The price of this improvement is an increased number of stages H (because Γ= grows with smaller increments) and consequently a higher computa-

Fig. 2. A top-down view of the experimental setup. Location of the source is marked by a red asterisk at (−373.5, 0). Sensor locations indicated by blue circles, whose radius is proportional (on the log-scale) to the corresponding concentration measurement. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

bounds are: σ xl0 = 2.10 , σ yl0 = 1.42 , σ Ql0 = 1.49. Several observations can be made from Table 1: (i) The RB-ISPC estimates are consistently more accurate than the ISPC estimates, thus signifying the benefit of Rao–Blackwellisation. (ii) The RMS errors of RB-ISPC approach the theoretical bounds as N is increased. (iii) The improvement achieved by Rao–Blackwell dimension reduction is more

0.5

y [m]

y [m]

0.5

0

-0.5

-1

-0.5

0

0

-0.5

0.5

-1

x [m]

0.5

0

0.5

0.5

y [m]

y [m]

0

x [m]

0.5

0

-0.5

-0.5

-1

-0.5

0

-0.5

0.5

x [m] Fig. 3. RB-ISPC applied to the experimental dataset: the

sample {r(0n, = ) }1≤ n ≤ N

0

-1

-0.5

x [m] shown by green dots at iteration: (a) = = 0 (a draw from the prior PDF), (b) at = = 2 , (c) at = = 3 and (d) at

= = 5 (final). The true source position marked with a blue asterisk. The red circles indicate the positions of sensors. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.) 181

Signal Processing 132 (2017) 177–182

B. Ristic et al.

tional load.

Acknowledgements

5.2. Experimental data

This research is supported in part by the Defence Science and Technology through its Strategic Research Initiative on Trusted Autonomous Systems under Research Agreement `Autonomous Systems for Application in Complex Urban Environments'.

The experimental dataset was collected by COANDA Research & Development Corporation using their large recirculating water channel (10 m long, 1.5 m wide and 0.9 m deep). The floor of the water channel was covered with a metal mesh of height 4 mm to give surface roughness. The source, releasing fluorescein dye at a constant rate from a narrow vertical tube, was placed 4 mm above the bottom of the channel at coordinates x 0 = −373.5 mm , y0 = 0 mm . Concentration data were collected at several downstream positions using a 1D laser induced fluorescence linescan system, at the rate of 300 lines per second, for a total sampling time of 1000 s. The dataset, available in the supplementary material of [12], was extracted from the full recordings averaged over 100 s. It consists of measurements at S=48 sensors (4 rows of 12) at downstream positions, indicated by the circles in Fig. 2. The radius of each circle is proportional (on the log-scale) to the corresponding concentration measurement. The true source location at coordinates (−373.5, 0) is marked by a red asterisk. The height of all sensors was 9.3 mm . In order to apply the dispersion and measurement models described in Section 2, we have converted the raw concentration measurements in the experimental dataset to the nearest integer values. Because the concentration measurements were very small decimal numbers, they were first multiplied by 105. The described l0 (at pre-processing of raw measurements affects only the estimate of Q l0 needs to be divided by 105). the end, Q The results obtained from the proposed RB-ISPC algorithm are shown in Figs. 3 and 4. The following parameters were used: the sample size N=200, integration time t0 = 100 s, sensor size a = 10−3 m , the prior PDFs π (x 0 ) = < [−1.2 m, 0 m], π ( y0) = < [−0.5 m, 0.5 m], and η0 = 2 , ϑ0 = 2.6 ; the RB-ISPC threshold was set to χ = 20 ; the environmental parameters τ = 40 s, D = 3·10−5 m2/s and V=0.01 m/s. The algorithm performed H=4 iterations, with the following values of correction factors: γ1 = 0.015625, γ2 = 0.0615234375, γ3 = 0.230712890625 and γ4 = 0.692138671875. Fig. 3 shows the sample {r(0n, = ) }1≤ n ≤ N at: (a) = = 1 (a draw from the prior, line 2 of Algorithm 1); (b) at = = 2 ; (c) at = = 3 and (d) at the final stage. The true source location is marked with a blue asterisk. Observe how with the growing iteration number = , the members of the sample {r(0n, = ) }1≤ n ≤ N move closer to the true source position. The final ISPC iteration, shown in Fig. 3(d), results in an accurate and precise source position estimate. Fig. 4 plots two PDFs of Q0: the prior PDF, . (Q0 ; η0 , ϑ0), is shown ϑ) is displayed with a with a black line; the posterior PDF, . (Q0 ; ηl , l N N 1 1 green line. Here ηl = N ∑n =1 η(n, H ) and l ϑ = N ∑n =1 ϑ(n, H ) are sample means, where η(n, H ) and ϑ(n, H ) are obtained after the final iteration of the RB-ISPC, see line 28 in Algorithm 1. Because the raw measurements were multiplied by 105, the x-axis in Fig. 4 should be divided by 105.

Appendix A. Supplementary data Supplementary data associated with this article can be found in the online version at http://dx.doi.org/10.1016/j.sigpro.2016.10.005. References [1] R.J. Kendall, S.M. Presley, G.P. Austin, P.N. Smith (Eds.), Advances in Biological and Chemical Terrorism Countermeasures, CRC Press, Boca Raton, 2008. [2] W.K. Panofsky, Nuclear proliferation risks, new and old, Issues Sci. Technol. 19 (4) (2003) 73–74. [3] J. Matthes, L. Gröll, H.B. Keller, Source localization by spatially distributed electronic noses for advection and diffusion, IEEE Trans. Signal Process. 53 (5) (2005) 1711–1719. [4] L.C. Thomson, B. Hirst, G. Gibson, S. Gillespie, P. Jonathan, K.D. Skeldon, M.J. Padgett, An improved algorithm for locating a gas source using inverse methods, Atmos. Environ. 41 (2007) 1128–1134. [5] B. Deb, Iterative estimation of location and trajectory of radioactive sources with a networked system of detectors, IEEE Trans. Nucl. Sci. 60 (2) (2013) 1315–1326. [6] A. Keats, E. Yee, F.-S. Lien, Bayesian inference for source determination with applications to a complex urban environment, Atmos. Environ. 41 (3) (2007) 465–479. [7] R. Humphries, C. Jenkins, R. Leuning, S. Zegelin, D. Griffith, Atmospheric tomography: a Bayesian inversion technique for determining the rate and location of fugitive emissions, Environ. Sci. Technol. 46 (3) (2012) 1739–1746. [8] M. Ortner, A. Nehorai, A. Jeremic, Biochemical transport modeling and Bayesian source estimation in realistic environments, IEEE Trans. Signal Process. 55 (6) (2007) 2520–2532. [9] I. Senocak, N.W. Hengartner, M.B. Short, W.B. Daniel, Stochastic event reconstruction of atmospheric contaminant dispersion using Bayesian inference, Atmos. Environ. 42 (33) (2008) 7718–7727. [10] E. Yee, Inverse dispersion for an unknown number of sources: model selection and uncertainty analysis, ISRN Appl. Math. 2012 (2012) 20 (Article ID 465320). [11] M.R. Morelande, B. Ristic, Radiological source detection and localisation using Bayesian techniques, IEEE Trans. Signal Process. 57 (11) (2009) 4220–4231. [12] B. Ristic, A. Gunatilaka, R. Gailis, A. Skvortsov, Bayesian likelihood-free localisation of a biochemical source using multiple dispersion models, Signal Process. 108 (2015) 13–24. [13] N. Tsoulfanidis, S. Landsberger, Measurement and Detection of Radiation, 4th ed., CRC Press, Boca Raton, 2015. [14] M. Vergassola, E. Villermaux, B.I. Shraiman, Infotaxis as a strategy for searching without gradients, Nature 445 (7126) (2007) 406–409. [15] G. Casella, C.P. Robert, Rao-Blackwellisation of sampling schemes, Biometrica 83 (1) (1996) 81–94. [16] A. Doucet, N. De Freitas, K. Murphy, S. Russell, Rao-Blackwellised particle filtering for dynamic Bayesian networks, in: Proceedings of the Sixteenth Conference on Uncertainty in Artificial Intelligence, 2000, pp. 176–183. [17] J.-B. Masson, M. Bailly-Bachet, M. Vergassola, Chasing information to search in random environments, J. Phys. A: Math. Theor. 42 (2009). [18] A.M. Hein, S.A. McKinley, Sensing and decision-making in random search. Proceedings of the National Academy of Sciences, 109(30), pp.12070-12074. http://dx.doi.org/10.1073/pnas.1202686109. [19] N. Voges, A. Chaffiol, P. Lucas, D. Martinez, Reactive searching and infotaxis in odor source localization, PLoS Comput. Biol. 10 (10) (2014). [20] B. Ristic, A. Skvortsov, A. Gunatilaka, A study of cognitive strategies for an autonomous search, Inf. Fusion 28 (2016) 1–9. [21] H. Hajieghrary, M.A. Hsieh, I.B. Schwartz, Multi-agent search for source localization in a turbulent medium, Phys. Lett. A 380 (20) (2016) 1698–1705. [22] A. Gelman, J.B. Carlin, H.S. Stern, D.B. Rubin, Bayesian Data Analysis, 3rd ed., CRC Press, Boca Raton, 2003. [23] M. Hazewinkel, Gamma-distribution, in: Encyclopedia of Mathematics, Springer, 2001. [24] 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-Verlag, New York, 2001 (Chapter 12). [25] C.P. Robert, G. Casella, Monte Carlo Statistical Methods, 2nd ed., Springer, 2004. [26] R.M. Neal, Annealed importance sampling, Stat. Comput. 11 (2) (2001) 125–139. [27] S. Godsill, T. Clapp, Improvement strategies for Monte Carlo particle filters, in: A. Doucet, N. deFreitas, N.J. Gordon (Eds.), Sequential Monte Carlo Methods in Practice, Springer-Verlag, New York, 2001 (Chapter 7). [28] B. Ristic, S. Arulampalam, N. Gordon, Beyond the Kalman Filter: Particle Filters for Tracking Applications, Artech House, 2004.

6. Conclusions The paper presented an efficient parameter estimation technique for localisation of a chemical, biological or a radiological source. Using Rao–Blackwell dimension reduction, the posterior density of source intensity, conditioned on the source position, is calculated analytically, thereby demanding numerical, Monte Carlo estimation, only in the space of source position. Simulation results demonstrate the statistical benefits of the proposed scheme: the estimates are consistently more accurate, especially for small sample sizes. The practical importance of the proposed scheme is also demonstrated using an experimental dataset characterised by a turbulent flow. In the future, the proposed technique will be integrated with an autonomous robotic search for a source of toxic gas [20].

182