Comparison of Two Approaches for Detection and Estimation of Radioactive Sources

Abstract

This paper describes and compares two approaches for the problem of determining the number of radioactive point sources that potentially exist in a designated area and estimating the parameters of these sources (their locations and strengths) using a small number of noisy radiological measurements provided by a radiation sensor. Both approaches use the Bayesian inferential methodology but sample the posterior distribution differently: one approach uses importance sampling with progressive correction and the other a reversible-jump Markov chain Monte Carlo sampling. The two approaches also use different measurement models for the radiation data. The first approach assumes a perfect knowledge of the data model and the average background radiation level, whereas the second approach quantifies explicitly the uncertainties in the model specification and in the average background radiation level. The performances of the two approaches are compared using experimental data acquired during a recent radiological field trial.

1. Introduction

There is growing concern in recent years regarding the risk of the smuggling and illicit trafficking of stolen radiological material. This state of affairs has heightened the spectre of radiological terrorism involving the use of an improvised radiological dispersal device (dirty bomb) to disperse radiological material over a large area using the force of conventional explosives (see, Panofsky [1]). In this context, the detection, localisation, and characterization of sources of radioactive material using radiation sensors (or, networks of such sensors) have important implications for safety and public security.

Various researchers have focussed on the problem of recovering the physical location and strength of radioactive sources using measurements obtained from a number of radiological detectors. Howse et al. [2] applied a nonlinear, recursive least-squares algorithm for tracking the positions of a moving radioactive point source. Stephens Jr. and Peurrung [3] conducted a cost-benefit analysis involving a variety of (conflicting) operational requirements that need to be considered in the design and deployment of a network of radiation sensors for use in the detection of a moving radioactive source. Nemzek et al. [4] and Brennan et al. [5] investigated the application of a Bayesian methodology for the detection and estimation of mobile radioactive sources, with a specific focus on examination of the signal-to-noise ratio envelope (and the concomitant detection limits) expected from utilization of a distributed sensor network for detection of a radiological source in a moving vehicle. Gunatilaka et al. [6] proposed the application of a maximum likelihood estimation algorithm, extended Kalman filter, and unscented Kalman filter for the localisation of a single static radiological point source. Morelande et al. [7, 8] applied a Bayesian methodology realized using a particle filtering technique with progressive correction to localize multiple static radiological point sources. Finally, Mendis et al. [9] applied binary and continuous genetic algorithms to solve a maximum likelihood estimation problem for recovery of the location and strengths of an a priori known number of static radiological point sources.

In this paper, we describe two approaches for radiological source localisation for the case when the number of sources is unknown a priori and compare their performance by applying them to an experimental data set collected in a field trial in Puckapunyal Military Area (Victoria, Australia). Both approaches use the Bayesian inferential methodology but sample the posterior distribution using different algorithms. In the first approach, the sampling algorithm uses importance sampling with progressive correction (PC), while the second approach uses the Markov chain Monte Carlo (MCMC) algorithm for the sampling. While the first algorithm uses the minimum description length (MDL) to estimate the unknown number of sources, the second algorithm uses the reversible-jump MCMC approach for this purpose. Finally, the measurement model for the radiation data differs for the two approaches in the following way: the first approach uses a form of the posterior distribution that implicitly assumes that the average background radiation count rate is exactly known and that the error in the model for the radiation data is negligible, whereas the second approach utilizes a form of the posterior distribution that relaxes these assumptions.

The experimental data for evaluating the two approaches were acquired using a Geiger-Müller detector in the presence of multiple radiological point sources of unequal strengths. Background radiation data in the absence of these radiation sources were also acquired to provide an estimate of the average background radiation level.

The organisation of the paper is as follows. The radiological source localisation problem is formulated in Section 2. Section 3 describes the two approaches we use for radiological source estimation. Section 4 describes the radiation field trial and experimental data sets. Some implementation aspects of the two algorithms are discussed in Section 5. Section 6 presents the results of the application of the two approaches for source reconstruction to the experimental data. Conclusions are drawn in Section 7.

2. Problem Formulation

The problem of the detection and estimation of the characteristics (e.g., location, activity) of an unknown number of static radiological point sources from a finite number of “noisy” measurements (data) obtained from radiation sensors is addressed using the Bayesian formalism. This formalism involves the application of Bayes' rule:𝑝(𝚯∣𝐳,𝐼)∝𝑝(𝚯∣𝐼)𝑝(𝐳∣𝚯,𝐼),(2.1) where Θ are the parameters that we are trying to determine and 𝐳 is the data. In (2.1), 𝐼 is the background (contextual) information available in the problem, “∣” denotes “conditional upon”, and the various factors that appear in this relation have the following interpretation: 𝑝(Θ∣𝐼) is the prior probability for the parameters Θ, 𝑝(𝐳∣Θ,𝐼) is the likelihood function which specifies the probability that we observe data 𝐳 when Θ is known exactly, and 𝑝(Θ∣𝐳,𝐼) is the posterior probability for the parameters Θ in light of the new information introduced through the acquired data 𝐳. In essence, Bayes' rule informs us how to update our knowledge about the parameters Θ encoded in 𝑝(Θ∣𝐼) following the acquisition of data 𝐳. To apply the Bayesian formalism to the radiological source reconstruction problem, we need to assign appropriate functional forms for the prior probability and the likelihood function, which in turn determines the posterior probability.

2.1. Assignment of Prior Probability

Here, we assume that a number 𝑟 of point sources of gamma radiation are present in a flat open area that is free of any obstructions. The number of sources, 𝑟, is unknown. If 𝑟>0, the 𝑖th point source of gamma radiation (where 𝑖=1,2,…,𝑟) is parameterised by (i)its location 𝐱_𝑠,𝑖≡(𝑥_𝑠,𝑖,𝑦_𝑠,𝑖) in a Cartesian coordinate system with (𝑥_𝑠,𝑖,𝑦_𝑠,𝑖)∈𝒜, where 𝒜⊂ℝ² is some large region that is assumed to contain all the sources, (ii)its intensity rate (source strength), 𝑄_𝑠,𝑖.

The parameter vector of source 𝑖 is thus given by𝜽_𝑖=𝑥_𝑠,𝑖𝑦_𝑠,𝑖𝑄_𝑠,𝑖^⊺,(2.2) where ⊺ denotes the matrix transpose. The source parameter vectors are collected into a stacked vector 𝜽(𝑟)=(^{𝜽⊺1⋯𝜽⊺𝑟})^⊺, and let Θ≡(𝑟𝜽(𝑟)^⊺)^⊺.

If we assume the logical independence of the various source parameters (or equivalently, the components of Θ), the prior probability for the parameters factorizes as follows:𝑝(𝚯∣𝐼)=𝑝(𝑟∣𝐼)^𝑟_𝑖=1𝑝𝑄_𝑠,𝑖∣𝐼𝑝𝐱_𝑠,𝑖∣𝐼.(2.3) Now, we need to assign explicit functional forms for each of the component prior probabilities in (2.3).

The prior distribution of 𝑟, 𝑝(𝑟∣𝐼), is assigned as a binomial distribution:𝑝(𝑟∣𝐼)=𝑟^max−𝑟_min!𝑟−𝑟_min!𝑟_max−𝑟!𝑝^{∗(𝑟−𝑟min)}1−𝑝^∗^{𝑟max−𝑟},(2.4) (𝑟=𝑟_min,𝑟_min+1,…,𝑟_max). In (2.4), 𝑝^∗∈(0,1) is the binomial rate, and, 𝑟_min and 𝑟_max are, respectively, the minimum and maximum number of sources.

Two alternative functional forms are used for the prior distribution of 𝑄_𝑠,𝑖 (𝑖=1,2,…,𝑟). The first of the two alternatives assigns the prior for the source strength 𝑄_𝑠,𝑖 to be a gamma distribution:𝑝𝑄_𝑠,𝑖∣𝐼=𝑄^{𝜅−1𝑠,𝑖}exp−𝑄_𝑠,𝑖/𝜓Γ(𝜅)𝜓^𝜅,𝑖=1,2,…,𝑟,(2.5) where 𝜅 and 𝜓 are the shape and scale parameters, respectively, and Γ(𝑥) is the gamma function. For the second of the two alternatives, the prior distribution for 𝑄_𝑠,𝑖 is assigned a Bernoulli-uniform mixture model:𝑝𝑄_𝑠,𝑖∣𝐼=(1−𝛾)𝛿𝑄_𝑠,𝑖+𝛾𝕀^(0,𝑄max)𝑄_𝑠,𝑖𝑄_max,𝑖=1,2,…,𝑟,(2.6) where 𝛾 is the probability that the source is active (Pr{𝑄_𝑠,𝑖>0}=𝛾), 𝛿(𝑥) is the Dirac delta function, and 𝑄_max is an a priori upper bound on the expected source intensity rate. In (2.6), 𝕀_𝐴(𝑥) denotes the indicator (characteristic) function for set 𝐴, with 𝕀_𝐴(𝑥)=1 if 𝑥∈𝐴 and 𝕀_𝐴(𝑥)=0 if 𝑥∉𝐴.

Finally, the prior distribution for the source location 𝐱_𝑠,𝑖 is assigned to be uniformly distributed over the region 𝒜⊂ℝ² that is assumed a priori to contain the source:𝑝𝐱_𝑠,𝑖∣𝐼=𝕀^𝒜𝐱_𝑠,𝑖meas(𝒜),𝑖=1,2,…,𝑟,(2.7) where meas(𝒜) is the area of the region 𝒜.

2.2. Assignment of Likelihood Function

To prescribe a functional form for the likelihood function, we need to relate the hypotheses of interest about the unknown source(s)¹ as encoded in Θ to the available radiation data 𝐳 measured by detectors emplaced in the vicinity of the source(s). Towards this objective, we need to formulate a measurement model for the radiation data 𝐳.

The radiation counts from nuclear decay obey Poisson statistics (see, Tsoulfanidis [10]), implying that the probability that a gamma radiation detector registers 𝑧∈ℕ∪{0} counts (ℕ∪{0} being the set of natural numbers including zero) in an exposure time of 𝜏 seconds, from a source that emits on average 𝜇 counts per second, is𝒫(𝑧;𝜆)=𝜆^𝑧𝑧!exp(−𝜆),(2.8) where 𝜆=𝜇𝜏 is the parameter (both the mean and the variance) of the Poisson distribution.

The measurements of the radiation field are made using a Geiger-Müller (GM) counter. Let 𝑧_𝑗∈ℕ∪{0} (𝑗=1,…,𝑚) be a measured count from the GM counter, taken at (receptor) location (𝑥_𝑙,𝑗,𝑦_𝑙,𝑗)∈ℝ². The following assumptions are made: (1) the GM counter has a uniform directional response; (2) the radiation sources are point sources; (3) there is negligible attenuation of gamma radiation due to air; (4) the radiation measurements are independently distributed, and, (5) the exposure time (or the sampling interval) 𝜏 is constant for all measurements.

2.2.1. Measurement Model: Known Background Rate

The likelihood function is simply the joint density of the measurement vector 𝐳≡(^{𝑧1⋯𝑧𝑚})^⊺ conditional on the parameter vector 𝜽 and the knowledge that 𝑟 sources are present. With this identification, the likelihood function can then be written as a product of Poisson distributions (Martin and Harbison [11]; Tsoulfanidis [10]) as follows:𝑝𝐳∣𝜇_𝑏,𝚯,𝐼≡𝑝𝐳∣𝜇_𝑏,𝑟,𝜽,𝐼≡𝑙𝐳∣𝜇_𝑏,𝑟,𝜽=^𝑚_𝑗=1𝒫𝑧_𝑗;𝜆_𝑗(𝚯),(2.9) where 𝜆_𝑗(Θ) is the mean radiation count for the 𝑗th sensor location:𝜆_𝑗(𝚯)=𝜇_𝑏+^𝑟_𝑖=1𝑄_𝑠,𝑖𝑑^2𝑗𝑖𝜏,(2.10) with𝑑_𝑗𝑖=𝑥_𝑙,𝑗−𝑥_𝑠,𝑖²+𝑦_𝑙,𝑗−𝑦_𝑠,𝑖²^1/2(2.11) being the Euclidean distance between the 𝑖th source and the 𝑗th sensor.

The constant 𝜇_𝑏 in (2.10) is the mean background count rate (namely, the average count rate due to the background radiation only which includes the contribution from cosmic and terrestrial radiation). In the measurement model (and concomitant form of the likelihood function given by (2.9) and (2.10)), the mean background count rate 𝜇_𝑏 and the number of radiological point sources 𝑟 are assumed to be known a priori. In consequence, 𝜇_𝑏 and 𝑟 have been added as quantities to the right of the vertical bar in the likelihood function of (2.9) to indicate explicitly that these quantities are known exactly.

The mean signal count rate 𝜇^𝑗𝑠 at the 𝑗th sensor due to the radioactive point sources only is modelled in (2.10) as𝜇^𝑗𝑠=^𝑟_𝑖=1𝑄_𝑠,𝑖𝑥_𝑙,𝑗−𝑥_𝑠,𝑖²+𝑦_𝑙,𝑗−𝑦_𝑠,𝑖².(2.12) In (2.10), it is implicitly assumed that the model for the mean signal count rate in (2.12) is perfect (no model error) so 𝜇^𝑗𝑠=𝜇^𝑗𝑠.

2.2.2. Measurement Model: Unknown Background Rate

The key assumptions in the formulation of the measurement model of (2.9) and (2.10) are that (1) the mean background count rate 𝜇_𝑏 is known exactly, and (2) the model error associated with the model for the mean signal count rate 𝜇_𝑠 of (2.12) is negligible. Now, we will relax these critical assumptions and approach the (perhaps) more realistic case when we must incorporate both the model error and the effects of an uncertain mean background count rate in the measurement model.

Towards this objective, the following model is assumed for the (true though unknown) mean signal count rate 𝜇^𝑗𝑠 (𝑗=1,2,…,𝑚). The mean signal count rate can be described by the radiation data 𝐳 and by the quantitative model (cf. (2.12)), so 𝜇^𝑗𝑠=𝜇^𝑗𝑠+𝑒^𝑗(1)=𝜇^𝑗𝑠+𝑒^𝑗(2),𝑗=1,2,…,𝑚,(2.13) where 𝜇^𝑗𝑠 is an estimate for the mean signal count rate at the 𝑗th sensor obtained from the data 𝐳 (more specifically, 𝑧_𝑗), 𝜇^𝑗𝑠 is the mean signal count rate at the 𝑗th sensor obtained from the model given in (2.12), 𝑒^𝑗(1) is the uncertainty in the estimate 𝜇^𝑗𝑠 of the mean signal count rate, and 𝑒^𝑗(2) represents the model error incurred by using (2.12) to predict the mean signal count rate. The model errors consist of a number of contributions which include (1) the effects of absorption of the gamma radiation by the air and (2) the effects of backscattering from the ground surface which may (possibly) influence the inverse-square law dependence in the model of (2.12).

The measurement model of (2.13) can be rearranged to give𝜇^𝑗𝑠−𝜇^𝑗𝑠=𝑒^𝑗(2)−𝑒^𝑗(1)≡𝜖_𝑗,𝑗=1,2,…,𝑚,(2.14) where 𝜖_𝑗 is the composite error. If the expectation value ⟨𝜖_𝑗⟩ of the composite error is zero and the variance ⟨𝜖^2𝑗⟩ is assumed to be given by 𝜎^2𝑗, then application of the principle of maximum entropy (see, Jaynes [12]) provides the following explicit form for the likelihood function (assuming that at least one radiological source is present, so 𝑟>0):𝑝(𝐳∣𝚯,𝐼)=1∏^𝑚𝑗=1√2𝜋𝜎_𝑗exp⎧⎪⎨⎪⎩−12^𝑚_𝑗=1𝜇^𝑗𝑠−𝜇^𝑗𝑠𝜎_𝑗²⎫⎪⎬⎪⎭.(2.15) It is noted that the Gaussian likelihood (which is effectively the probability distribution of the “noise” 𝜖_𝑗) obtained from application of the maximum entropy principle is simply the most conservative noise distribution that is consistent with the given facts (namely, the first two moments of the “noise” are known). In the Gaussian likelihood function of (2.15), the variance of the composite error 𝜖_𝑗 is obtained as𝜎^2𝑗=𝜎^21,𝑗+𝜎^22,𝑗,𝑗=1,2,…,𝑚,(2.16) where 𝜎^21,𝑗 is the variance of the uncertainty 𝑒^𝑗(1) in the estimate 𝜇^𝑗𝑠 of the mean signal count rate (obtained from the data 𝐳) and 𝜎^22,𝑗 is the variance of the model error 𝑒^𝑗(2). In this paper, the model error standard deviation (or square root of the variance) is specified as 𝜎_2,𝑗=(𝜇^𝑗𝑠)^1/2.

To complete the specification for the likelihood function in (2.15), we need to provide the estimate 𝜇^𝑗𝑠 along with the variance 𝜎^21,𝑗 in this estimate. To this purpose, we assume that we have available a measurement of background radiation (obtained in the absence of any point radiation sources) given to us as 𝑧^𝑏 counts obtained over an exposure time 𝜏_𝑏. Now, suppose that a radiation sensor at a fixed receptor location obtains a measurement of 𝑧 counts in an exposure time 𝜏. We will now show how to use this information to obtain an estimate for the mean signal count rate 𝜇_𝑠, as well as to characterize the uncertainty in this estimate in terms of the standard deviation 𝜎₁.

To accomplish this, we utilize a Bayesian solution for the extraction of weak signals in strong background interference described by Loredo [13]. We will briefly outline the application of this methodology to the problem in this paper. To begin, we can apply Bayes' rule to compute the joint probability for the mean signal count rate 𝜇_𝑠 and the mean background count rate 𝜇_𝑏 conditioned on the measurement 𝑧. This gives𝑝𝜇_𝑠,𝜇_𝑏∣𝑧,𝐼∝𝑝𝜇_𝑠∣𝜇_𝑏,𝐼𝑝𝜇_𝑏∣𝐼𝑝𝑧∣𝜇_𝑠,𝜇_𝑏,𝐼.(2.17) It is stressed that the contextual information 𝐼 in (2.17) includes the contextual information 𝐼_𝑏 concerning the measurement of the background radiation and the contextual information 𝐼_𝑠 concerning the measurement of the radiation field which may include the presence of one or more radiation sources. Hence, the contextual information 𝐼 that is available in the problem can be decomposed as 𝐼=𝐼_𝑏,𝐼_𝑠.

We will assign a least informative prior probability for 𝜇_𝑠 (the so-called Jeffreys prior; see Jaynes [12]) which has the following form (conditioned on a known mean background rate 𝜇_𝑏):𝑝𝜇_𝑠∣𝜇_𝑏,𝐼=1𝜇_𝑠+𝜇_𝑏.(2.18) The prior probability for 𝜇_𝑏 is chosen to be the following informative prior (using the information that is available, namely, that a measurement of the background radiation yielded 𝑧^𝑏 counts in an exposure time 𝜏_𝑏):𝑝𝜇_𝑏∣𝐼≡𝑝𝜇_𝑏∣𝑧^𝑏,𝐼_𝑏,𝐼_𝑠=𝑝𝜇_𝑏∣𝑧^𝑏,𝐼_𝑏∝𝑝𝜇_𝑏∣𝐼_𝑏𝑝𝑧^𝑏∣𝜇_𝑏,𝐼_𝑏,(2.19) since 𝐼_𝑠 does not have any bearing on the determination of 𝜇_𝑏. The prior 𝑝(𝜇_𝑏∣𝐼_𝑏) is assigned the nonsinformative Jeffreys prior so 𝑝(𝜇_𝑏∣𝐼_𝑏)=1/𝜇_𝑏 and the likelihood function 𝑝(𝑧^𝑏∣𝜇_𝑏,𝐼_𝑏)=(𝜇_𝑏𝜏_𝑏)^𝐳𝑏exp(−𝜇_𝑏𝜏_𝑏)/𝑧^𝑏!, so𝑝𝜇_𝑏∣𝐼=𝑝𝜇_𝑏∣𝑧^𝑏,𝐼_𝑏∝𝜏_𝑏𝜇_𝑏𝜏_𝑏^𝑧𝑏−1exp−𝜇_𝑏𝜏_𝑏.(2.20)

The likelihood function in (2.17), associated with the measurement of 𝑧 counts in exposure interval 𝜏 coming from a radiation field with mean count rate (𝜇_𝑠+𝜇_𝑏), is described by the following Poisson distribution:𝑝𝑧∣𝜇_𝑠,𝜇_𝑏,𝐼=𝜏^𝑧𝜇_𝑠+𝜇_𝑏^𝑧𝑧!exp−𝜇^𝑠+𝜇_𝑏𝜏.(2.21) Now, if we substitute (2.18), (2.20), and (2.21) in (2.17), the joint probability for the mean signal and background count rates is given by𝑝𝜇_𝑠,𝜇_𝑏∣𝑧,𝐼∝𝜇_𝑠+𝜇_𝑏^𝑧−1𝜇^{𝑧𝑏𝑏−1}exp−𝜇_𝑠𝜏exp−𝜇_𝑏𝜏_𝑏+𝜏.(2.22) The marginal probability distribution for the mean signal count rate can be derived by marginalizing with respect to the (unknown) mean background count rate:𝑝𝜇_𝑠∣𝑧,𝐼=^∞0𝑝𝜇_𝑠,𝜇_𝑏∣𝑧,𝐼𝑑𝜇_𝑏,(2.23) which on using (2.22) and performing the integration finally gives (on normalization of the probability density function)𝑝𝜇_𝑠∣𝑧,𝐼=^𝑧_𝑖=1Φ_𝑖𝜏𝜇_𝑠𝜏^𝑖−1(𝑖−1)!exp−𝜇^𝑠𝜏,(2.24) whereΦ_𝑖=1+𝜏^𝑏/𝜏^𝑖𝑧+𝑧^𝑏−𝑖−1!/(𝑧−𝑖)!∑^𝑧𝑗=11+𝜏_𝑏/𝜏^𝑗𝑧+𝑧^𝑏−𝑗−1!/(𝑧−𝑗)!.(2.25) Note that (2.24) and (2.25) provide an explicit expression for the posterior probability of the mean signal count rate 𝜇_𝑠 that is independent of the (unknown) mean background count rate 𝜇_𝑏. More specifically, the precise value of 𝜇_𝑏 is not known, so it is treated as a nuisance parameter here and removed by integration (marginalization).

Given the posterior probability of 𝜇_𝑠 summarized by (2.24) and (2.25), a best estimate for the mean signal count rate can be obtained as the posterior mean:𝜇_𝑠≡⟨𝜇_𝑠⟩=^∞0𝜇_𝑠𝑝𝜇_𝑠∣𝑧,𝐼𝑑𝜇_𝑠=1𝜏^𝑧_𝑖=1𝑖Φ_𝑖.(2.26) The second moment of 𝜇_𝑠 about zero can be evaluated as𝜇^2𝑠=^∞0𝜇^2𝑠𝑝𝜇_𝑠∣𝑧,𝐼𝑑𝜇_𝑠=1𝜏^2𝑧_𝑖=1𝑖(𝑖+1)Φ_𝑖,(2.27) from which the uncertainty 𝜎₁ in the estimate for 𝜇_𝑠 given in (2.26) can be determined as𝜎₁=𝜇^2𝑠−𝜇_𝑠²^1/2.(2.28) Equations (2.26), (2.27) and (2.28) when used in conjunction with (2.16) and (2.15) fully determine the form of the likelihood function (assuming implicitly that 𝑟>0). This functional form for the likelihood function accounts explicitly for the model error and the sampling error (the latter of which is compounded by an uncertain mean background count rate).

For the case 𝑟=0, the likelihood function needs to be determined as follows (again assuming that our state of knowledge consists of a background radiation measurement yielding 𝑧^𝑏 counts in exposure time 𝜏_𝑏):𝑝𝐳∣𝑧^𝑏,𝐼=^𝑚_𝑗=1^∞0𝒫𝑧_𝑗∣𝜇_𝑏𝜏𝑝𝜇_𝑏∣𝑧^𝑏,𝐼_𝑏𝑑𝜇_𝑏,(2.29) which on insertion of the normalized form of (2.20) for 𝑝(𝜇_𝑏∣𝑧^𝑏,𝐼_𝑏) and evaluation of the integral gives𝑝𝐳∣𝑧^𝑏,𝐼=^𝑚_𝑗=1𝜏^𝑧𝑏_𝑏𝜏^𝑧𝑗𝜏+𝜏_𝑏^{𝑧𝑗+𝑧𝑏}𝑧_𝑗+𝑧^𝑏−1!𝑧_𝑗!𝑧^𝑏−1!.(2.30) In (2.29) and (2.30), the radiation measurements 𝑧_𝑗 are assumed to be made using the same exposure time 𝜏 (without any loss in generality).

Finally, the likelihood function of (2.15) valid for 𝑟>0 and the likelihood function of (2.30) valid for 𝑟=0 can be combined to give the following general likelihood function for the measurement model assumed in this section:𝑝𝐳∣𝑧^𝑏,𝚯,𝐼=𝕀^{𝑟>0}(𝑟)∏^𝑚𝑗=1√2𝜋𝜎_𝑗exp⎧⎪⎨⎪⎩−12^𝑚_𝑗=1𝜇^𝑗𝑠−𝜇^𝑗𝑠𝜎_𝑗²+𝕀_{𝑟=0}⎫⎪⎬⎪⎭(𝑟)^𝑚_𝑗=1𝜏^𝑧𝑏_𝑏𝜏^𝑧𝑗𝜏+𝜏_𝑏^{𝑧𝑗+𝑧𝑏}𝑧_𝑗+𝑧^𝑏−1!𝑧_𝑗!𝑧^𝑏−1!,(2.31) where the assumption of the availability of a background radiation count 𝑧^𝑏 (obtained over an exposure time 𝜏_𝑏) has been included in the conditioning of the likelihood function (and of the posterior distribution based upon this likelihood function).

3. Computational Framework

This section describes the computational procedures that were used for extracting the source parameter estimates for source reconstruction. To this purpose, two different methodologies for sampling from two different forms of the posterior distribution for the source parameters are described.

3.1. Importance Sampling Using Progressive Correction

Importance sampling using progressive correction is a computational methodology for sampling from a posterior distribution. In this paper, this sampling procedure is applied to a form of the posterior distribution for the source parameters in which the number of sources 𝑟 and the mean background count rate 𝜇_𝑏 are assumed to be known a priori (see (5.1)). To this purpose, the information contained in the prior distribution for the source parameters 𝑝(𝜽(𝑟)∣𝑟,𝐼) and in the measurements 𝐳 assimilated through the likelihood function 𝑝(𝐳∣𝜇_𝑏,𝑟,𝜽(𝑟),𝐼)≡𝑙(𝐳∣𝜇_𝑏,𝑟,𝜽(𝑟)) is combined to obtain the posterior distribution for 𝜽(𝑟) (assuming that the number of sources 𝑟 and the background count rate 𝜇_𝑏 are known):𝑝𝜽(𝑟)∣𝜇_𝑏,𝑟,𝐳,𝐼∝𝑙𝐳∣𝜇_𝑏,𝑟,𝜽(𝑟)𝑝(𝜽(𝑟)∣𝑟,𝐼).(3.1) The minimum mean-squared error (MMSE) estimate of 𝜽 is then the posterior expectation,𝜽(𝑟)≡⟨𝜽(𝑟)∣𝜇_𝑏,𝑟,𝐳⟩=𝜽(𝑟)𝑝𝜽(𝑟)∣𝜇_𝑏,𝑟,𝐳,𝐼𝑑𝜽(𝑟).(3.2)

Because the posterior PDF 𝑝(𝜽(𝑟)∣𝜇_𝑏,𝑟,𝐳,𝐼) and, hence, the posterior expectation cannot be found exactly for the measurement model of Section 2, an approximation of the integral in (3.2) is computed via importance sampling. This involves drawing 𝑁 samples of the parameter vector from an importance density and approximating the integral by a weighted sum of the samples:𝜽(𝑟)≈^𝑁_𝑛=1𝑤^𝑛𝜽^𝑆,𝑛(𝑟),(3.3) where {𝜽^𝑆,𝑛(𝑟),𝑛=1,…,𝑁} is a “particle” set consisting of samples drawn from the importance density and {𝑤^𝑛,𝑛=1,…,𝑁} is the set of importance weights. The sum of these weights is equal to one (namely, ∑^𝑁𝑛=1𝑤^𝑛=1).

Because the prior distribution will often be more diffuse than the likelihood, many samples drawn from the prior may fall outside the region of parameter space favoured by the likelihood function. Therefore, a straightforward application of (3.1) may result in poor source parameter estimates. The progressive correction approach overcomes this problem by using a multistage procedure in which samples are obtained from a series of posterior distributions which become progressively closer to the true posterior distribution. This is achieved by adopting an approximate likelihood at each stage which is somewhat more diffuse than the true likelihood.

Let 𝑆 denote the number of stages and 𝑝_𝑠(𝜽(𝑟)∣𝜇_𝑏,𝑟,𝐳,𝐼), 𝑠=1,…,𝑆 denote the posterior PDF at stage 𝑠. A series of 𝑝_𝑠(𝜽(𝑟)∣𝜇_𝑏,𝑟,𝐳,𝐼) can be constructed by setting𝑝_𝑠𝜽(𝑟)∣𝜇_𝑏,𝑟,𝐳,𝐼∝𝑙_𝑠𝐳∣𝜇_𝑏,𝑟,𝜽(𝑟)𝑝(𝜽(𝑟)∣𝑟,𝐼),(3.4) where the intermediate likelihood at stage 𝑠=1,…,𝑆 is 𝑙_𝑠(𝐳∣𝜇_𝑏,𝑟,𝜽(𝑟))=𝑙(𝐳∣𝜇_𝑏,𝑟,𝜽(𝑟))^𝐺𝑠 with 𝐺_𝑠=∑^𝑠ℓ=1𝛾_ℓ, such that 𝛾_𝑠∈[0,1), for 𝑠=1,…,𝑆 and 𝐺_𝑆=𝛼≤1. Assume that a random sample {𝑤^𝑛,𝜽^{𝑠−1,𝑛}(𝑟)} from 𝑝_𝑠−1(𝜽(𝑟)∣𝜇_𝑏,𝑟,𝐳,𝐼) is available, and it is desired to produce a sample from 𝑝_𝑠(𝜽(𝑟)∣𝜇_𝑏,𝑟,𝐳,𝐼). The steps of the PC algorithm are given in Algorithm 1. Note that for 𝑠=1, the sample {𝑤^𝑛,𝜽^{𝑠−1,𝑛}(𝑟)} is drawn directly from the prior 𝑝(𝜽(𝑟)∣𝑟,𝐼) with 𝑤^𝑛=1/𝑁.

(1) 𝑠=0, 𝐺0=0, 𝜇𝑏=𝜇∗𝑏, 𝑟=𝑟∗; (2) for 𝑛=1,…,𝑁 do (3)   Draw 𝜽0,𝑛(𝑟)∼𝑝(𝜽(𝑟)∣𝑟,𝐼); (4)   Compute 𝑙(𝐳∣𝜇𝑏,𝑟,𝜽0,𝑛(𝑟)); (5) end for (6) while 𝐺𝑠<𝛼 and 𝑠<𝑆 do (7)   𝑠←𝑠+1; (8)   Select 𝛾𝑠; (9)   𝐺𝑠=𝐺𝑠−1+𝛾𝑠; (10)  for 𝑛=1,…,𝑁 do (11)    Weights: 𝑤𝑠,𝑛=𝑙(𝐳∣𝜇𝑏,𝑟,𝜽𝑠−1,𝑛(𝑟))𝛾𝑠/∑𝑁𝑗=1𝑙(𝐳∣𝜇𝑏,𝑟,𝜽𝑠−1,𝑛(𝑟))𝛾𝑠; (12)  end for (13)  Resample {𝜽𝑠−1,𝑛(𝑟)}𝑁𝑛=1 with weights {𝑤𝑠,𝑛}𝑁𝑛=1 to give {𝜽𝑛∗}𝑁𝑛=1 with uniform weights (14)  for 𝑛=1,…,𝑁 do (15)    𝜽𝑠,𝑛(𝑟)=𝜽𝑛∗+𝝐𝑠,𝑛, where 𝝐∼𝑔𝑠; (16)  end for (17) end while

Algorithm 1 

Importance sampling with progressive correction algorithm.

The performance of the procedure depends somewhat on the number of steps 𝑆 and the expansion factors 𝛾₁,…,𝛾_𝑆. The expansion factors 𝛾_𝑠 were selected in line 8 of Algorithm 1 using an adaptive scheme proposed in Musso et al. [14]. In this adaptive scheme, the expansion factors are selected after each step rather than being selected a priori. Line 13 in Algorithm 1 performs a resampling on the samples {𝑤^𝑠,𝑛,𝜽^{𝑠−1,𝑛}(𝑟)}^𝑁𝑛=1 to give samples {1/𝑁,𝜽^𝑛∗}^𝑁𝑛=1 that have uniform weights. This resampling is undertaken using the systematic resampling algorithm (see, Kitagawa [15]), which is summarized in the form of a pseudocode given in Table  3.2 of Ristic et al. [16]. Line 15 in Algorithm 1 performs regularisation (jittering) of samples in order to avoid their duplication. Here, 𝑔_𝑠 is a Gaussian regularisation kernel (see, Musso et al. [14]).

The case 𝛼=1 in the progressive correction algorithm should be used only when the likelihood function model in (2.9) is perfectly correct and its parameters precisely known. As we discussed in Section 2.2.2, this is unrealistic in practice: there are many potential errors in the model, such as the background radiation level, air attenuation, sensor location, and propagation loss. In order to make the PC algorithm robust against the modelling errors, it is recommended to adopt 𝛼<1 (this change affects the “while” loop in line 6 of Algorithm 1). In this way, the measurement likelihood is effectively approximated by a fuzzy membership function, but the Bayesian estimation framework holds as the weights in Step  9 of Algorithm 1 are normalised. This is similar to the treatment of the ambiguously generated unambiguous measurements in ([17], Chapter 7). Parameter 𝛼 is a tuning parameter, which needs to be determined by field trials and calibration.

The PC estimate of (3.2) is approximated as (after the resampling step of line 13 and the regularization step of line 15 in Algorithm 1)𝜽(𝑟)=1𝑁^𝑁_𝑛=1𝜽^𝑆,𝑛(𝑟).(3.5)

3.1.1. Source Number Estimation

The importance sampling algorithm using progressive correction described above assumes that the number of sources is known a priori (see line 1 of Algorithm 1) which initializes 𝑟=𝑟_∗, where 𝑟_∗ is the number of sources. However, the number of sources will not be known in advance in practice, requiring this number also to be estimated from the data.

For this purpose, we use the minimum description length (MDL) criterion which chooses the value of 𝑟 that minimises (Kay [18] page 225]); Rissanen [19]):𝜒_𝑟=−log𝑙𝐳∣𝜇_𝑏,𝑟,𝜽_ML(𝑟)+𝑛^𝑟2log𝑚,(3.6) where 𝜽_ML(𝑟) is the maximum likelihood estimate of 𝜽(𝑟) assuming 𝑟 sources are present and 𝑛_𝑟≡3𝑟 is the dimension of 𝜽(𝑟) under the hypothesis 𝑟∈ℳ. Here, ℳ={0,…,𝑟_max} denotes a set of candidate source numbers up to some maximum 𝑟_max. While MDL is derived based on the maximum likelihood estimate 𝜽_ML(𝑟) (see (3.6)), we will use the MDL in conjunction with the PC by simply replacing 𝜽_ML(𝑟) with the PC estimate 𝜽(𝑟) in (3.6), so𝜒_𝑟=−log𝑙𝐳∣𝜇_𝑏,𝑟,𝜽(𝑟)+𝑛^𝑟2log𝑚.(3.7) Finally, the estimate of the number of sources 𝑟 is then determined as follows:𝑟̂=argmin_𝑟∈ℳ𝜒_𝑟.(3.8)

To estimate the number of sources using the MDL criterion, we first run the importance sampling using progressive correction for all values of 𝑟=1,…,𝑟_max. For each of these values of 𝑟, the estimated value for 𝜽(𝑟) in (3.5) is used to compute the corresponding value for 𝜒_𝑟 from (3.7). The case 𝑟=0 (namely, the hypothesis that there are no sources present) is treated separately. There is no need to run the importance sampling using PC in this case, as (3.7) reduces to𝜒₀=−log𝑙𝐳∣𝜇_𝑏(3.9) when 𝑟=0.

3.2. Reversible-Jump Markov Chain Monte Carlo Algorithm

In Section 3.1, we described an importance sampling algorithm for sampling from the posterior distribution based on a sequence of proposal densities (involving, as such, a sequence of modified forms of a likelihood function that has been raised to powers 𝛾_𝑗<1, 𝑗=1,2,…,𝑆). In its application here, calculations with the importance sampling algorithm using progressive correction were made separately for each model structure (namely, for each value of 𝑟), and comparisons of the different models were realized using the MDL criterion. In this section, we describe an alternative for sampling from the posterior distribution that involves application of an MCMC algorithm that allows sampling over both model and parameter spaces simultaneously (namely, 𝑟 which indexes the various model structures is treated simply as another parameter, so that it is estimated jointly with the other model parameters 𝜽(𝑟)).

Towards this objective, we apply a reversible-jump MCMC (RJMCMC) algorithm. The formalization of RJMCMC algorithms for dealing with variable dimension models has been described in the seminal work of Green [20]. The application of the RJMCMC algorithm to the problem of inverse dispersion has been described previously by Yee [21, 22]. The algorithm involves constructing a Markov chain whose stationary distribution is the posterior distribution 𝑝(Θ∣𝑧^𝑏,𝐳,𝐼) of the source parameters Θ≡(𝑟𝜽(𝑟)^⊺)^⊺∈ℝ^1+3𝑟. To this purpose, let {Θ^(𝑡)}≡{(𝑟^(𝑡)𝜽(𝑟(𝑡))^⊺)^⊺} (𝑡=0,1,2,…) be the state vector of a Markov chain that is so constructed that its stationary distribution coincides with 𝑝(Θ∣𝑧^𝑏,𝐳,𝐼) (see (5.2) for the explicit form of the posterior distribution that is sampled from using the RJMCMC algorithm in this paper). The construction of the Markov chain uses an RJMCMC algorithm in which the MCMC moves over the model and parameter spaces are separated into two categories: (1) propagation moves which do not change the dimensionality of the parameter space; and, (2) transdimensional jump moves which change the source distribution by ±1 radiological point source and, in so doing, changes the parameter space dimension by ±3.

In the first category of moves (which are dimension-conserving), the parameter vector 𝜽(𝑟) (𝑟 fixed) is partitioned as follows: 𝜽(𝑟)=(𝜽^1⊺𝜽^2⊺)^⊺ where 𝜽^1⊺≡(𝑄_𝑠,1𝑄_𝑠,2⋯𝑄_𝑠,𝑟)∈ℝ^𝑟 and 𝜽^2⊺≡(^{𝐱⊺𝑠,1𝐱⊺𝑠,2⋯𝐱⊺𝑠,𝑟})∈ℝ^2𝑟. The parameters collected together in 𝜽¹ are linearly related to the radiation “measurements” (cf. (2.12)). For these parameters, we apply a Gibbs sampler for the update (propagation) move involving drawing samples directly from the univariate full conditional distribution for 𝑄_𝑠,𝑖 (which happens to be a truncated Bernoulli-Gaussian distribution in this case). The remaining parameters in 𝜽² are related nonlinearly to the radiological measurements (cf. (2.12)). Owing to the fact that the univariate conditional posterior distribution for these parameters cannot be determined analytically, we apply a Metropolis-Hastings (M-H) sampler to update the parameters in 𝜽². Towards this purpose, we update 𝜽^2𝑘∈𝜽² (𝑘=1,2,…,2𝑟) by generating a new value 𝜽^𝑘2 from a proposal transition kernel that is chosen here to be a Gaussian mixture, with each component of the mixture having a mean 𝜃^2𝑘 and different variances 𝛽^2𝑘. The variances of the components in this Gaussian mixture are chosen to cover several orders of magnitude from 0.001 to 0.1 times the “length” of the domain of definition for the parameter (assigned using the prior distribution for the parameter).

In the second category of moves (which are dimension-changing), the source distribution is modified by ±1 radiological point source by (1) a birth move that adds a single radiological point source and (2) a death move that removes a single radiological point source to/from the current source distribution. If a birth move is selected, the “coordinates” of the new radiological point source are generated by drawing a random sample from the prior distribution for each coordinate. Alternatively, if a death (reverse) move is selected, then a radiological point source in the current source distribution is randomly selected and removed. Following the recommendation of Green [20], the probabilities for the birth (𝑝_𝑏) and death (𝑝_𝑑) moves are specified as follows:𝑝_𝑏(𝑟)=12min1,𝑝(𝑝𝑟(+𝑟∣1𝐼∣)𝐼),𝑟=𝑟^min𝑝,…,𝑟_max−1,_𝑑(𝑟+1)=12min1,𝑝(𝑝𝑟(+𝑟∣1𝐼∣)𝐼),𝑟=𝑟^min,…,𝑟_max−1,(3.10) ensuring that the probability of a jump move (either birth or death) 𝑝_𝑏+𝑝_𝑑 lies in [0.5,1] at each iteration. We note that for 𝑟=𝑟_min, 𝑝_𝑑(𝑟)=0 and for 𝑟=𝑟_max, 𝑝_𝑏(𝑟)=0.

In practical terms, the RJMCMC algorithm as applied to the problem of source reconstruction can be set up as follows. (1)Specify values for the hyperparameters (𝑟_min,𝑟_max,𝑄_max,𝛾,𝑝^∗) and 𝒜 which define the prior distribution 𝑝(Θ∣𝐼), as well as a value for 𝑡_upper (maximum number of iteration steps to take). (2)Set the iteration counter 𝑡=1 and choose an initial state Θ^(𝑡−1) for the Markov chain by sampling from 𝑝(Θ∣𝐼). (3)Starting from Θ^(𝑡−1), conduct the following sequence of moves to update the state vector to Θ^(𝑡): 𝚯−−−→^{ℳ(𝑡−1)𝑏,𝑑}𝚯−−→_⋆ℳ¹𝚯−−→_ℳ⋆⋆²𝚯^(𝑡),(3.11) where Θ_⋆ and Θ_⋆⋆ denote some intermediate transition states between iterations (𝑡−1) and 𝑡. (4)Change the counter from 𝑡 to 𝑡+1, and return to step 3 until a maximum number of steps (𝑡_upper) has been taken.

In Step  3, ℳ_𝑏,𝑑 denotes a jump move involving either the birth (addition) of a radiological point source at a random locations or the death (removal) of an existing radiological point source with probabilities 𝑝_𝑏 and 𝑝_𝑑, respectively; ℳ₁ denotes the update of the source intensity rates using a Gibbs sampler; and, ℳ₂ denotes the update of the source locations using an M-H sampler.

3.2.1. Simulated Annealing

To improve convergence of the Markov chain, we have implemented a simulated annealing procedure. This procedure is very similar to the progressive correction that is used in conjunction with the importance sampling (as described in Section 3.1). In a nutshell, simulated annealing mimics the quasistationary evolution of a system through a sequence of equilibrium samplings, involving an infinitely slow switching between the initial and final states. More concretely, let us consider an ensemble of 𝑁_∗ different source distributions (system) that has been randomly drawn from a modified posterior having the following form: 𝑝_𝜆(Θ∣𝑧^𝑏,𝐳,𝐼)∝𝑝(Θ∣𝐼)𝑝^𝜆(𝐳∣𝑧^𝑏,Θ,𝐼), where 𝜆∈[0,1] is a “coolness” parameter. These samples will be labelled Θ_𝑘(𝜆) (𝑘=1,2,…,𝑁_∗). When 𝜆=0 (initial state), the likelihood function is switched off and the modified posterior distribution reduces exactly to the prior distribution. On the other hand, when 𝜆=1 (final state) the modified posterior distribution is exactly the posterior that we wish to sample from. In between these two extremes, with 𝜆∈(0,1), the effects of the radiation measurements 𝐳 are introduced gradually through the modified (or “softened”) likelihood function 𝑝^𝜆(𝐳∣𝑧^𝑏,Θ,𝐼).

The parameter 𝜆 can be interpreted as an inverse temperature parameter 𝑇 (𝜆≡1/𝑇), so 𝜆∈[0,1] implies that 𝑇∈[1,∞]. The modified posterior 𝑝_𝜆(Θ∣𝑧^𝑏,𝐳,𝐼) (𝜆∈[0,1)) corresponds to “heating” the posterior distribution to a temperature 𝑇=1/𝜆>1. This heated distribution is flattened relative to the posterior distribution 𝑝(Θ∣𝑧^𝑏,𝐳,𝐼), with the result that moves required to adequately cover the parameter space become more likely. In other words, this facilitates Markov chain mobility in the remote regions of the posterior distribution, ensuring a more reliable relaxation of the chain into the potentially (exponentially) small regions of the hypothesis space where the posterior distribution is large.

We initialize the stochastic sampling scheme at 𝜆=0 (infinite temperature) by randomly drawing 𝑁_∗ samples of source distributions Θ_𝑘(0) (𝑘=1,2,…,𝑁_∗) from the prior distribution 𝑝₀(Θ∣𝑧^𝑏,𝐳,𝐼). Given an ensemble of 𝑁_∗ samples Θ_𝑘(𝜆) that has achieved equilibrium (at temperature 𝑇=1/𝜆) with respect to the modified posterior 𝑝_𝜆(Θ∣𝑧^𝑏,𝐳,𝐼), an ensemble of 𝑁_∗ samples Θ_𝑘(𝜆+𝛿𝜆) that is consistent with 𝑝_{𝜆+𝛿𝜆}(Θ∣𝑧^𝑏,𝐳,𝐼) (at the reduced temperature 𝑇=1/(𝜆+𝛿𝜆), 𝛿𝜆>0) can be obtained by using the weighted resampling method (see, Gamerman and Lopes [23]) applied to Θ_𝑘(𝜆) (𝑘=1,2,…,𝑁_∗). To this purpose, each sample Θ_𝑘(𝜆) is assigned a normalized importance weight 𝑤_𝑘 as follows: 𝑤_𝑘=𝑝^{𝜆+𝛿𝜆}𝚯_𝑘(𝜆)∣𝑧^𝑏,𝐳,𝐼𝑝_𝜆𝚯_𝑘(𝜆)∣𝑧^𝑏,𝐳,𝐼×^𝑁∗_𝑗=1𝑤_𝑗⁻¹=𝑝^𝛿𝜆𝐳∣𝑧^𝑏,𝚯_𝑘(𝜆),𝐼×^𝑁∗_𝑗=1𝑤_𝑗⁻¹,(3.12) for 𝑘=1,2,…,𝑁_∗. We note that this procedure of computing weights is analogous to line 11 in Algorithm 1 for the importance sampling algorithm using progressive correction.

Next, a resample is drawn from the discrete distribution concentrated at the samples {Θ_𝑘(𝜆)}^𝑁∗_𝑘=1 with weights {𝑤_𝑘}^𝑁∗_𝑘=1. This resampling involves mapping the random discrete measure defined by {𝑤_𝑘,Θ_𝑘(𝜆)}^𝑁∗_𝑘=1 into the new random discrete measure defined by {1/𝑁_∗,Θ_𝑘^∗}^𝑁∗_𝑘∗=1 having uniform weights (namely, Pr{Θ_𝑘^∗=Θ_𝑘(𝜆)}=𝑤_𝑘 with the result that Θ_𝑘^∗∼𝑝_{𝜆+𝛿𝜆}(Θ∣𝑧^𝑏,𝐳,𝐼) where ~ denotes “is distributed as”). The resampling conducted here is analogous to line 13 in Algorithm 1 for the importance sampling algorithm using progressive correction. The resampling used to produce the weighted replica ensemble acts as a selection mechanism in the sense that the number of replicas (copies) of a particular sample Θ_𝑘(𝜆) is increased for those samples that are assigned a large weight 𝑤_𝑘 (which, in turn, enhances the number of relevant configurations in the ensemble). Note that the weight determined according to (3.12) provides a measure as to whether the state Θ_𝑘(𝜆) contains new information about areas of the “energy” landscape that are favourable (reflected in larger values of the likelihood function).

There are many ways to implement this resampling (with replacement) from {Θ_𝑘(𝜆)}^𝑁∗_𝑘=1, but in this paper we use the systematic resampling scheme described by Kitagawa [15] (namely, the same resampling scheme that is used in the PC algorithm described in Section 3.1). This scheme requires 𝑂(𝑁_∗) time to execute and minimizes the Monte Carlo variation in the resample. For completeness, we briefly describe this resampling scheme as applied to the current problem. A vector (𝑛₁,𝑛₂,…,𝑛_𝑁∗) of copies of the samples Θ_𝑘(𝜆) (𝑘=1,2,…,𝑁_∗) is obtained by computing a vector (𝜌₁,𝜌₂,…,𝜌_𝑁∗) of the cumulative sums of 𝑁_∗×(𝑤₁,𝑤₂,…,𝑤_𝑁∗), generating a random draw from the uniform distribution 𝑢∼𝒰([0,1]), and determining 𝑛_𝑘 (𝑘=1,2,…,𝑁_∗) from𝑛_𝑘=𝑛⌊𝜌_𝑘+𝑢⌋−⌊𝜌_𝑘−1+𝑢⌋,𝑘=2,3,…,𝑁_∗−1,₁=⌊𝜌₁+𝑢⌋,𝑛_𝑁∗=𝑁_∗−𝜌_𝑁∗−1+𝑢,(3.13) where ⌊⋅⌋ denotes the “integer part of”. After resampling, the weights for each member of the resample are reset to 𝑤_𝑘=1/𝑁_∗ (namely, an equal weight is assigned to each member of the resample). After each resampling step to give an ensemble of 𝑁_∗ samples Θ_𝑘(𝜆+𝛿𝜆) that is in equilibrium at temperature 𝑇=1/(𝜆+𝛿𝜆) (approximately or better) with respect to 𝑝_{𝜆+𝛿𝜆}(Θ∣𝑧^𝑏,𝐳,𝐼), we apply 𝑁_𝑟 (typically 10) Markov chain transitions 𝑇_{𝜆+𝛿𝜆} to each member of the ensemble obtained from the resampling operation. These Markov chain transitions utilize the update scheme given by (3.11) and leave 𝑝_{𝜆+𝛿𝜆}(Θ∣𝑧^𝑏,𝐳,𝐼) invariant.

It is necessary to specify the sequence of 𝛿𝜆's used to move the ensemble of samples {Θ_𝑘(𝜆)}^𝑁∗_𝑘=1 from 𝜆=0 to 𝜆=1. This sequence defines the annealing schedule for 𝜆∈[0,1]. Taking a hint from the physical process of annealing, it is useful to allow the system to cool slowly. This allows the system to transition through a sequence of quasiequilibrium states and, in so doing, minimizes the probability of the system being trapped into a metastable state along the path of annealing. To mimic this process for the case of simulated annealing, at each step of the annealing schedule (starting from 𝜆=0 where the stochastic sampling scheme is exploring the prior) the step size 𝛿𝜆 is chosen so that1𝑁_∗𝑁^∗_𝑘=1||𝑤_𝑘−1||≈𝑅,(3.14) where 𝑤_𝑘 are the weights computed in accordance to (3.12) and 𝑅 is a small constant (𝑅≪1). The constant 𝑅 is chosen to provide a limiter on the maximum permissible size of 𝛿𝜆 in each step of the annealing schedule. Note that unlike the PC algorithm, the sequence of 𝜆_𝑖∈[0,1] (𝑖=0,1,2,…,𝐾) with 𝜆_𝑖<𝜆_𝑗 for 𝑖<𝑗 and 𝜆₀=0, 𝜆_𝐾=1 that define the annealing schedule need not sum to 𝛼 (namely, ∑^𝐾𝑖=0𝜆_𝑖≠𝛼 in general).

When 𝜆=1, the annealing phase is complete (corresponding, as such, to the burn-in phase of the algorithm) and probabilistic exploration of the hypothesis space proceeds using the RJMCMC algorithm described at the end of Section 3.2. In the application of this algorithm following simulated annealing, Step  2 is modified as follows: instead of sampling Θ⁽⁰⁾ from the prior distribution 𝑝(Θ∣𝐼), these initial states are simply chosen to coincide with the samples from the ensemble of source distributions {Θ_𝑘(𝜆=1)}^𝑁∗_𝑘=1 obtained at the end of the annealing phase.

4. Experimental Trial

A radiological field trial was conducted on a large, flat, and open area without obstacles within a military site in Puckapunyal, Victoria, Australia. An AN/PDR-77 radiation survey meter equipped with a gamma probe containing two Geiger-Müller tubes to cover both low and high ranges of dose rates was used to collect the radiation data. When connected to the AN/PDR-77 radiation survey meter, this gamma probe is capable of measuring gamma radiation dose rates from background to 9.99 Sv h⁻¹ without saturating [24], and it has a fairly flat response (∼±10%) from 0.1 to above 1 MeV (see [24]; page 7). Measured dose rate data were recorded in 𝜇Sv h⁻¹. We converted these dose rate data into raw count measurements 𝑧_𝑗∈ℕ∪{0} by multiplication by a conversion factor (which was 0.21⁻¹ for our probe) to convert to count rate data and subsequently by multiplication by the exposure time to convert to raw count data (measured over the exposure time).

Three radiation sources were used in the field trial: two cesium sources (137Cs) and one cobalt (60Co) source (see Table 1). While the strengths of the radiation sources are typically characterised in terms of their activity in GBq (cf. Table 1), for simplicity we will characterise the strength (intensity rate) of a radiation point source by the expected dose rate at a distance of one meter from the source. The strengths of our three radiation sources calibrated in this manner were 𝑄_𝑠=1912, 392.28, and 98.07 𝜇Sv m²h⁻¹, for Sources 1, 2, and 3, respectively. We will use these values as the “true” source strengths in our source estimations.

Table 1 

Radiation sources used in the field trial.

The sources were mounted at the same height above the ground as the gamma probe. To ensure that radiation sources appear as isotropic in the horizontal plane, they were placed in a vertical configuration so that the handling rods were pointing up. Radiation dose measurements were collected on a grid of surveyed points that were carefully measured and marked beforehand in the local Cartesian coordinate system on the asphalt surface of the airfield. The data were acquired when the trolley-mounted gamma probe was positioned over individual grid points. During data collection at each grid point, the gamma probe remained stationary until sixty measurements were acquired. The exposure time for each radiation dose measurement (effectively the sampling interval) was kept constant at about Δ=0.8 s. Multiple measurements were collected at each location so that the candidate algorithms could be tested using several different batches of experimental data, or different total exposure intervals can be used in the source reconstruction by choosing the appropriate number of data points at each grid point (e.g., to test the source reconstruction algorithms for raw counts obtained over an exposure interval 𝜏, we can simply select 𝑛_∗ data points at each grid position such that 𝜏≈𝑛_∗Δ).

Data sets were collected with one, two, and three radiation sources emplaced, respectively. These data sets are referred to as test sets 1, 2 and 3. The sources used when collecting these data sets and their locations in the local Cartesian coordinate system are listed in Table 2.

Table 2 

The locations of radiation sources in the field trial.

An areal picture of the Puckapunyal airfield site where the field trial was conducted is shown in Figure 1. The three red stars show the locations where the three radiation sources were emplaced. The white cross symbols indicate the grid points where data were collected. The 𝑥- and 𝑦-axis markings on the plot show the distances along these directions in meters (as well as the local Cartesian coordinate system that is used here to refer to locations of sources and grid positions).

Figure 1 

Aerial image of the Puckapunyal airfield site where the field trial was conducted. The red stars located at coordinates (11,10) m, (3,50) m, and (41,5) m of the local Cartesian coordinate system indicate where the three sources were emplaced for Test 3. The numbered grid points marked with crosses are the measurement points.

In order to estimate the background radiation level in the field trial, measurements were collected at a number of grid positions in the absence of radiation sources. At each of these grid positions, the detector was held stationary until 120 measurements were acquired. The background radiation counts were found to be spatially homogeneous and verify a Poisson distribution with mean background count rate of 𝜇_𝑏≈0.9 counts s⁻¹.

5. Implementational Aspects

5.1. Importance Sampling Algorithm Using PC

The importance sampling algorithm using progressive correction was applied to the experimental field trial data for source reconstruction. Towards this purpose, this algorithm was used to sample from the following posterior distribution of the source parameters:𝑝𝜽∣𝜇_𝑏,𝑟,𝐳,𝐼∝∝𝑝(𝜽∣𝑟,𝐼)𝑝𝐳∣𝜇_𝑏,𝑟,𝜽,𝐼_𝑚_𝑗=1×𝒫𝑧_𝑗;𝜆_𝑗(𝚯)_𝑟_𝑖=1𝑄^{𝜅−1𝑠,𝑖}exp−𝑄_𝑠,𝑖/𝜓Γ(𝜅)𝜓^𝜅×𝕀^𝒜𝐱_𝑠,𝑖meas(𝒜),(5.1) which combines the prior distribution for source strengths given by (2.5) and the prior distribution for source location given by (2.7) with the likelihood function given by (2.9). The quantities 𝜆_𝑗(Θ) are determined in accordance withs (2.10). In this form of the posterior distribution, it is assumed that the mean background count rate 𝜇_𝑏 and the number of sources 𝑟 are known a priori.

The importance sampling algorithm using progressive correction is applied to the experimental field trial data using the following values for the hyperparameters that define the prior distribution in (5.1). The gamma distribution that defines the prior for the source strengths has a shape parameter 𝜅=1.5 and a scale parameter 𝜓=8000 𝜇Sv m² h⁻¹. This choice of the hyperparameters for the gamma distribution provides a prior for the source strengths that is broad enough to cover all likely source strength values. The prior for the source location is a uniform distribution over a specified rectangular-shaped area 𝒜, typically containing all measurement points. For the current application, the area 𝒜 is selected as the following rectangular region: 𝒜≡{(−50,200)⊗(−20,150)}, where ⊗ denotes Cartesian product.

The PC algorithm is initialized by drawing a random sample in 3𝑟-dimensional parameter space from 𝑝(𝜽∣𝑟,𝐼). The number of samples in the PC algorithm is set to 𝑁=2500𝑟.

The expansion factors and the actual number of stages of the PC algorithm are data dependent. We have tuned the PC algorithm so that on average the number of stages is directly proportional to the number of sources. Parameter 𝛼 was set to 1.

5.2. Reversible-Jump MCMC Algorithm

The reversible-jump MCMC algorithm (applied in conjunction with simulated annealing) was used to infer an unknown number of radiological point sources from the experimental field trial data. This involves drawing samples of radiological point source distributions from the following posterior distribution for the source parameters Θ:𝑝𝚯∣𝑧^𝑏,𝐳,𝐼∝∝𝑝(𝚯𝕀∣𝐼)𝑝𝐳∣𝑧^𝑏,𝚯,𝐼_{𝑟>0}(𝑟)∏^𝑚𝑗=1√2𝜋𝜎_𝑗exp⎧⎪⎨⎪⎩−12^𝑚_𝑗=1𝜇^𝑗𝑠−𝜇^𝑗𝑠𝜎_𝑗²×𝑟_max⎫⎪⎬⎪⎭−𝑟_min!𝑟−𝑟_min!𝑟_max−𝑟!𝑝^{∗(𝑟−𝑟min)}1−𝑝^∗^{𝑟max−𝑟}×^𝑟_𝑖=1(1−𝛾)𝛿𝑄_𝑠,𝑖+𝛾𝕀^(0,𝑄max)𝑄_𝑠,𝑖𝑄_max×𝕀^𝒜𝐱_𝑠,𝑖+m𝕀eas_{𝑟=0}(𝒜)(𝑟)^𝑚_𝑗=1𝜏^𝑧𝑏_𝑏𝜏^𝑧𝑗𝜏+𝜏_𝑏^{𝑧𝑗+𝑧𝑏}𝑧_𝑗+𝑧^𝑏−1!×1−𝑝_∗𝑧_𝑗!𝑧^𝑏−1!^𝑟max,(5.2)

which combines the likelihood function given by (2.31) with the prior for the number of sources given by (2.4), the prior for the source strength given by (2.6), and the prior for the source location given by (2.7). For the source reconstruction, the hyperparameters that define the prior distribution in (5.2) are assigned the following values: 𝑟_min=0, 𝑟_max=4, 𝑝^∗=0.25, 𝛾=0.25, 𝑄_max=75,000 𝜇Sv m² h⁻¹, and 𝒜={(−50,200)⊗(−20,150)} which is used to define the prior bounds for the location (𝑥_𝑠,𝑦_𝑠) of any radiological point source. Note that the prior bounds for the source location are exactly those assigned for the application of importance sampling using PC.

The stochastic sampling algorithm was executed with 𝑁_∗=50 members of an ensemble of distributions of radiological point sources (originally, randomly drawn from the prior distribution 𝑝(Θ∣𝐼)). During the burn-in (or annealing) phase of the algorithm, the control strategy used to determine the annealing schedule (see (3.14)) was applied with 𝑅=0.03. After termination of the annealing at 𝜆=1, a further 𝑡_upper=1000 iterations of the RJMCMC procedure were applied to each of the 𝑁_∗=50 members of the ensemble of distributions of radiological point sources obtained at the end of the annealing phase. This gave 50,000 samples of distributions of radiological point sources drawn from the posterior distribution 𝑝(Θ∣𝑧^b,𝐳,𝐼) given in (5.2) (corresponding, as such, to the probabilistic exploration phase of the algorithm).

The samples of radiological point source distributions drawn from the posterior distribution of (5.2) can be used to make various inferences on the values of the source parameters. For example, the number of sources can be obtained from the maximum a posteriori (MAP) estimate as follows:𝑟̂=argmax_𝑟𝑝𝑟∣𝑧^𝑏,𝐳,𝐼,(5.3) where the (marginal) posterior probability for the number of sources is estimated from the samples as follows:𝑝̂𝑟∣𝑧^𝑏,𝐳,𝐼=1𝑁^#𝑡∶𝑟^(𝑡)=𝑟=1𝑁^𝑁′_𝑡=1𝕀_{𝑟}𝑟^(𝑡).(5.4) Here, 𝑁^ is the number of samples Θ⁽¹⁾,Θ⁽²⁾,… drawn from a sampled realization of the Markov chain. Similarly, using the Markov chain samples drawn from the posterior distribution, the posterior means of the various source parameters 𝜃^𝑖(𝑟)̂ (various components of 𝜽(𝑟)̂; see (2.2)) can be estimated by𝜃^𝑖(𝑟)̂=∑^{𝑁′𝑡=1}𝜃^{𝑟𝑖(𝑡)}_(𝑡)𝕀_{𝑟_}̂𝑟^(𝑡)∑^𝑁′_𝑡=1𝕀_{𝑟_}̂𝑟^(𝑡).(5.5) Similar estimates can be constructed for the posterior standard deviation of 𝜃^𝑖(𝑟)̂ which can be used as a measure of uncertainty in the recovery of this quantity. Alternatively, a 𝑝% credible (or highest posterior distribution (HPD)) interval that contains the source parameter 𝜃^𝑖(𝑟)̂ with 𝑝% probability, with the upper and lower bounds specified such that the probability density within the interval is everywhere larger than outside it, can be used as a measure of the uncertainty in the determination of this quantity.

In addition to these statistical quantities, the RJMCMC algorithm applied in conjunction with simulated annealing permits the normalization constant (or evidence) 𝑍 for the posterior distribution of (5.2) to be determined². This can be accomplished using an approach known as thermodynamic integration (see, von der Linden [25]).³ In this approach, 𝑍 can be evaluated as follows:log𝑍=¹⁰log𝑝𝐳∣𝑧^𝑏,𝚯,𝐼_𝜆d𝜆,(5.6) where ⟨⋅⟩ denotes the expectation operation taken with respect to the modified posterior 𝑝_𝜆(Θ∣𝑧^𝑏,𝐳,𝐼). This quantity can be computed during the simulated annealing phase of the stochastic simulation procedure. More specifically, the integral in (5.6) is discretized using the following Riemann sum:log𝑍≈^𝐾_𝑖=1log𝑝𝐳∣𝑧^𝑏,Θ,𝐼_𝜆𝑖𝛿𝜆_𝑖,(5.7) where 𝛿𝜆_𝑖≡(𝜆_𝑖−𝜆_𝑖−1) is determined according to (3.14); and, the annealing schedule so defined provides the set of quadrature points (𝜆₀=0,𝜆₁,…,𝜆_𝐾−1,𝜆_𝐾=1).

The evidence 𝑍 can be used to determine the relative entropy (or, information gain) corresponding to the assimilation of the radiation measurements 𝐳, which involves updating our state of knowledge about the source embodied in the prior distribution to that embodied in the posterior distribution. To this end, the relative entropy 𝐻 can be determined as follows (see, Cover and Thomas [26]): 𝐻≡log𝑝𝚯∣𝑧^𝑏,𝐳,𝐼𝑝(𝚯∣𝐼)𝑝𝚯∣𝑧^𝑏=log𝑝𝐳∣𝑧_𝑏,𝐳,𝐼𝑑𝚯,𝚯,𝐼−log(𝑍).(5.8) The relative entropy (of the posterior with respect to the prior) is a measure of the information gain provided by the receipt of the radiation measurements 𝐳. Geometrically, this can be visualized as the logarithm of the volumetric factor by which the prior has been compressed to give the posterior in the hypothesis space (the greater this compression, the greater the information gain embodied in 𝐳).

6. Results

To test the usefulness of the two approaches for source reconstruction, we have applied them to test data sets 1, 2, and 3 corresponding to one-, two-, and three-source examples, respectively. However, owing to space limitations, we will show only the results of application of the two approaches for test case 3. This test is arguably the most difficult case, but the results from the test are representative also for the other two test cases.

The results for test case 3 will be shown for two different applications: namely, the radiation measurements at the grid positions (see Figure 1) were obtained over exposure intervals of 𝜏=16 and 48 s. In the application of importance sampling using PC, the mean background count rate 𝜇_𝑏 was assumed to be exactly known (𝜇_𝑏=0.9 counts s⁻¹ was used in the algorithm). For the application of RJMCMC used in conjunction with simulated annealing, the first 24 s (𝜏_𝑏=24 s) of a measurement of the background counts 𝑧^𝑏 at receptor location (𝑥_𝑙,𝑦_𝑙)=(−40,−15) m (or, at grid position 1 as shown in Figure 1) was used in the source recovery.

Firstly, we show the results of the two applications of the importance sampling algorithm using PC for test case 3. For each of these two applications, the algorithm (summarized in Algorithm 1) was executed with fixed values of 𝑟=𝑟_∗ with 𝑟_∗=1, 2, and 3. Essentially, the assumption in the algorithm is that the number of sources is exactly known a priori, with 𝑟=𝑟_∗. At the completion of this process, the MDL was computed according to (3.7) for 𝑟∈{0,1,2,3} and a “best” estimate for the number of sources was found using (3.8). In both applications of the importance sampling algorithm using PC here (namely, for 𝜏=16 and 48 s), the minimum value of 𝜒_𝑟 was found to correspond to 𝑟̂=3, an estimate that coincides with the correct number of sources for test case 3.

The results of the source recovery for the two applications of the importance sampling algorithm are summarized in Table 3 for 𝑟=3 (assumed correct number of sources). Here, the posterior mean, posterior standard deviation, and lower and upper bounds for the 95% HPD interval of the source parameters (𝑥_𝑠,𝑦_𝑠,𝑄_𝑠) are shown for 𝜏=16 s (upper portion of Table 3) and for 𝜏=48 s (lower portion of Table 3). In addition, histograms of the source parameters (𝑥_𝑠,𝑦_𝑠,𝑄_𝑠) for the three sources are exhibited in Figures 2 and 3, respectively, for 𝜏=16 and 48 s.

Table 3 

The posterior mean, posterior standard deviation, and lower and upper bounds of the 95% HPD interval of the parameters 𝑥_𝑠,𝑖 (m), 𝑦_𝑠,𝑖 (m), and 𝑄_𝑠,𝑖 (𝜇Sv m² ℎ⁻¹) for 𝑖=1, 2, and 3 calculated from samples of distributions of radiological point sources drawn from the posterior distribution of 𝜃(𝑟)(𝑟=3) using importance sampling with progressive correction. The upper and lower parts of the table show the results of the source recovery using radiological measurements with exposure times of 𝜏=16 and 48 s, respectively.

(a) Source 1

(b) Source 2

(c) Source 3

Figure 2 

Inference of the parameters of the three radiological point sources obtained from samples drawn from the posterior distribution 𝑝(𝜽∣𝜇𝑏,𝑟,𝐳,𝐼) for 𝜏=16 s and 𝑟=3 using the importance sampling algorithm with PC correction. (a, b, c) Histograms for the three parameters, namely 𝑥𝑠-coordinate of source, 𝑦𝑠-coordinate of source, and source intensity rate 𝑄𝑠 that characterize sources 1, 2, and 3. In each frame, the solid vertical line indicates the true value of the parameter and the dashed vertical line corresponds to the best estimate of the parameter obtained as the posterior mean of the marginal posterior distribution for the parameter.

(a) Source 1

(b) Source 2

(c) Source 3

Figure 3 

Inference of the parameters of the three radiological point sources obtained from samples drawn from the posterior distribution 𝑝(𝜽∣𝜇𝑏,𝑟,𝐳,𝐼) for 𝜏=48 s and 𝑟=3 using the importance sampling algorithm with PC correction. (a, b, c) Histograms for the three parameters, namely, 𝑥𝑠-coordinate of source, 𝑦𝑠-coordinate of source, and source intensity rate 𝑄𝑠 that characterize sources 1, 2, and 3. In each frame, the solid vertical line indicates the true value of the parameter and the dashed vertical line corresponds to the best estimate of the parameter obtained as the posterior mean of the marginal posterior distribution for the parameter.

For the application with 𝜏=16 s, the algorithm successfully recovered the true parameters for all three sources to within the stated errors (e.g., the true parameters were all contained within a three-standard deviation interval about the best estimates of the parameters, and/or the true parameters were contained within the 95% HPD interval). However, this is not the case for the application of the algorithm with 𝜏=48 s. More specifically, in this case it is seen that neither the three-standard deviation interval about the best estimates (based on the posterior mean) nor the 95% HPD interval for (𝑥_𝑠,𝑦_𝑠,𝑄_𝑠) for all three sources enclose the true values for these parameters. The reason for this is that the likelihood function in (5.1) that is used in the importance sampling algorithm with PC does not incorporate model error (namely, the model for the radiation measurements given by (2.12) is assumed to be exact). In consequence, for 𝜏=48 s, the sampling error is smaller than the model error, so the latter becomes the dominant component in the error specification. The neglect of the model error in this case leads to underestimates for the precision at which some of the source parameters are recovered by the importance sampling algorithm.

Next, we show the results for source reconstruction for test case 3 using the RJMCMC algorithm (applied in conjunction with simulated annealing) for the two applications with 𝜏=16 and 48 s. Figure 4(a) displays a trace plot for the number of radiological point sources 𝑟 in source distribution models sampled from the posterior distribution of Θ (for 𝜏=48 s). The comparable result for 𝜏=16 s is very similar (and hence, not shown). The trace plot here corresponds to the first 500 iterations of the RJMCMC sampler obtained during the probabilistic exploration phase (post-burn-in iterations following the completion of the annealing phase). Note that death moves for source distribution models, involving transitions from 𝑟=3 to 2 (or from 𝑟=2 to 1 and from 𝑟=1 to 0) do not occur. Furthermore, birth moves involving transitions from 𝑟=2 to 3 do not occur, either. On the other hand, dimension-changing moves (births or deaths) involving transitions from 𝑟=3 to 4 (and vice versa) are seen to occur. The residence time for occupation of a state with 𝑟=3 is seen to be significantly longer than that for a state with 𝑟=4. This is confirmed in Figure 4(b) which exhibits the probability distribution 𝑝(𝑟)≡𝑝(𝑟∣𝑧^𝑏,𝐳,𝐼) for the number of sources (estimated as a histogram from the 50,000 post-burn-in samples in accordance with (5.4)).

(a)

(b)

Figure 4 

Trace plot (a) of the number of radiological point sources 𝑟 in the first 500 samples drawn from 𝑝(Θ∣𝑧𝑏,𝐳,𝐼) during the probabilistic exploration phase using RJMCMC (after the completion of the annealing phase) and the corresponding posterior distribution (b) for the number of sources, 𝑝(𝑟)≡𝑝(𝑟∣𝑧𝑏,𝐳,𝐼).

Figure 4(b) shows that the most probable number of sources (MAP estimate) for this test case is 3 (𝑟̂=3), which corresponds to the correct number of sources for this test. More specifically, 𝑟̂=3 is favored with a probability of about 0.7. Interestingly, it is noted that many of the distributions having four radiological point sources (which from Figure 4(b) are favored with probability of about 0.3) have one of the point sources “turned” off (namely, the source strength 𝑄_𝑠 for this source is identically zero). In fact, these distributions of radiological point sources should be considered as being identical to 𝑟=3 sources (and not 𝑟=4 sources owing to the fact that 𝑄_𝑠=0 for one of the sources implies that that source is not present, at least as far as the detection of that source by a radiological sensor is concerned). If this classification for the number of radiological sources is used instead, then 𝑟̂=3 is now favored with a probability greater than about 0.995.

Given the fact that our best estimate for the number of sources is 𝑟̂=3, we can now estimate the source parameters (𝑥_𝑠,𝑦_𝑠,𝑄_𝑠) corresponding to each of these three radiological point sources. Towards this purpose, we extracted all samples of distributions having exactly three sources (including all four-source samples, with one of the sources “turned off”). It is noted that an important issue in the interpretation of the results is the identifiability problem, which arises owing to the fact that the posterior distribution of Θ (cf. (5.2)) is invariant under a relabelling of the identifiers used for each of the radiological point sources in the distribution. In consequence, we relabel the radiological point sources so that sources 1, 2, and 3 correspond to the (arbitrary) labelling of sources used in Table 2. After this re-labelling, Table 4 summarizes the posterior mean and standard deviation, as well as the lower and upper bounds for the 95% HPD interval of the source parameters, for each of the three radiological point sources for radiological measurements with exposure times 𝜏=16 and 48 s.

Table 4 

The posterior mean, posterior standard deviation, and lower and upper bounds of the 95% HPD interval of the parameters 𝑥_𝑠,𝑖 (m), 𝑦_𝑠,𝑖 (m), and 𝑄_𝑠,𝑖 (𝜇Sv m² ℎ⁻¹) for 𝑖=1, 2, and 3 calculated from samples of distributions of radiological point sources drawn from the posterior distribution of Θ using the RJMCMC algorithm. The upper and lower parts of the table show the results of the source recovery using radiological measurements with an exposure times of 𝜏=16 and 48 s, respectively.

From Table 4, it is seen that the estimated values for the locations of the three sources agree well with the true locations of the sources in the reconstruction undertaken for both 𝜏=16 and 48 s. Comparing these estimated values for the source locations with the actual source locations, we see also that the algorithm has adequately recovered the source locations to within the errors as represented by the 95% HPD interval for 𝜏=16 and 48 s. Furthermore, the intensity rates of the sources are estimated quite well for both 𝜏=16 and 48 s, although the accuracy in the recovery of these parameters is not as good as for the source location. Although the 95% HPD interval does not contain the true estimate for 𝑄_𝑠 for sources 1 and 3 for both 𝜏=16 and 48 s, it should be noted that the 99% HPD interval does enclose the true value of 𝑄_𝑠 for these sources.

An examination of Table 4 shows that the recovery of the source parameters for 𝜏=16 and 48 s using the RJMCMC algorithm gives results that are very similar, both in terms of their accuracy and their precision. This is in contrast with the results obtained using the importance sampling algorithm where it was found that there was an underestimate in the precision of recovery of the source parameters for 𝜏=48 s. This difference in performance arises because the importance sampling algorithm was applied to the simpler posterior distribution of (5.1) which does not incorporate the effects of model error, whereas the RJMCMC algorithm was applied to the posterior distribution of (5.2) where the effects of model error are explicitly incorporated.

A comparison of Tables 3 and 4 shows that the source parameters are generally recovered with greater accuracy and better precision using the importance sampling algorithm than those recovered using the RJMCMC algorithm. The reason for this is because the importance sampling algorithm uses a posterior distribution (cf. (5.1)) that is much more informative than the posterior distribution (cf. (5.2)) used by the RJMCMC algorithm. In particular, the number of sources is assumed known a priori in the former, but is an unknown parameter in the latter; the mean background rate is assumed to be known in the former, but is unknown in the latter (where it is assumed that only a background count measurement 𝑧^𝑏 taken over an exposure time 𝜏_𝑏 is available), and the effects of model error on the source recovery are taken into account in the latter, but not in the former.

Histograms of the source parameters (approximating the marginal probability distributions for the parameters), associated with each of the three identified radiological point sources, are exhibited in Figure 5 for 𝜏=48 s. Again, the comparable results for 𝜏=16 s are not shown owing to the fact that they are very similar. A perusal of this figure suggests that the localisation of the three sources in both the x_s- and y_s-directions is generally very good. The recovery of the source intensity rates Q_s is quite good, although both the accuracy and precision in this recovery are poorer than the recovery of the source locations.

(a) Source 1

(b) Source 2

(c) Source 3

Figure 5 

Inference of the parameters of the three radiological point sources obtained from samples drawn from the posterior distribution 𝑝(Θ∣𝑧𝑏,𝐳,𝐼) for τ = 48 s using the RJMCMC algorithm. (a, b, c) Histograms for the three parameters, namely, xs-coordinate of source, 𝑦𝑠-coordinate of source, and source intensity rate 𝑄𝑠 that characterize sources 1, 2, and 3. In each frame, the solid vertical line indicates the true value of the parameter and the dashed vertical line corresponds to the best estimate of the parameter obtained as the posterior mean of the marginal posterior distribution for the parameter.

Finally, the information gain 𝐻 (see (5.8)) obtained from the radiation measurements 𝐳 for 𝜏=16 and 48 s (for test case 3) was found to be 𝐻=53.4 and 60.1 natural units (nits), respectively.⁴ This implies that the information contained in the radiation measurements with exposure times 𝜏=16 and 48 s allowed the “posterior volume” of the hypothesis space (volume of hypothesis space of reasonably large plausibility after receipt of the radiation data) to decrease by a factor of exp(𝐻)≈1.55×10²³ and 1.26×10²⁶, respectively, relative to the “prior volume” of the hypothesis space (volume of the hypothesis space of reasonably large plausibility before the receipt of the radiation data). It is seen that using radiation data measured with exposure time of 𝜏=48 s provided about 7 nits more information gain than that using radiation data measured with an exposure time of 𝜏=16 s.

7. Conclusions

This paper presented two approaches for estimating an unknown number of point radiation sources and their parameters employing radiation measurements collected using a gamma radiation detector. The first approach used importance sampling with PC to sample the posterior distribution, and the second approach used the reversible-jump MCMC technique in conjunction with simulated annealing for this purpose. The two approaches were compared by applying them to experimental data collected in a recent field trial.

It was demonstrated that both approaches perform well on the experimental data: in general, the number of sources is correctly determined and the parameters (e.g., location and intensity rates) that characterize each radiological point source are reliably estimated, along with a determination of the uncertainty (e.g., standard deviation, credible intervals) in the inference of the source parameters. However, it was shown that application of the first approach with 𝛼=1 did not adequately characterize the uncertainty in the source parameter estimates for the radiation measurements obtained with an exposure time 𝜏=48 s. The defect here was not due to the progressive correction algorithm per se, but rather the result of the fact that the simple form of the posterior distribution of 𝜽 used with this algorithm assumed that the model error was negligible. This turned out not to be the case, and application of the RJMCMC with a more sophisticated form of the posterior distribution of Θ that incorporated explicitly the model error was shown to correct this defect.

The results of this paper suggest that Bayesian probability theory is a powerful tool for the formulation of methodologies for source reconstruction. It should be stated that the application of importance sampling using PC to the posterior distribution of (5.1) is fully Bayesian in character (assuming that 𝑟 is known a priori), as is the application of the RJMCMC algorithm to the posterior distribution of (5.2). However, the use of MDL in conjunction with the importance sampling algorithm to determine an estimate of the number of sources is not consistent with a strict Bayesian formalism.

From this perspective, it should be noted that if the Bayes factor is used instead to determine the number of sources (namely, to determine which model of a distribution of radiological point sources to use), then the methodology would be fully consistent with a Bayesian formalism (see, Morelande and Ristic [8]). In fact, as discussed in Jaynes [12], application of the Bayes factor for model selection automatically incorporates Ockham's razor (namely, selecting the simplest model for a distribution of radiological point sources that “best” fits the available radiation data).

Importance sampling using PC may also be carried out over the joint space of 𝑟 and 𝜽 in the same way as the RJMCMC to estimate the number of sources and parameters associated with these sources simultaneously and in a manner that is fully consistent with the Bayesian formalism (see Ristic et al. [27]).

Acknowledgments

The authors would like to thank A. Eleftherakis, D. Marinaro, P. Harty, and I. Mitchell for their involvement in the field trial and A. Skvortsov, R. Gailis, and M. Morelande for useful technical discussions. This work has been partially supported by Chemical Biological Radiological-Nuclear and Explosives Research and Technology Initiative (CRTI) program under project number CRTI-07-0196TD.

Endnotes

1. This hypothesis is of the form “the unknown value of the source parameters is Θ”, so the numerical value of Θ indexes a set of hypotheses of this form.
2. The evidence 𝑍 is the normalization constant in (5.2), required to ensure that the posterior distribution is properly normalized (namely, 𝑍≡∫𝑝(Θ∣𝐼)𝑝(𝐳∣𝑧^𝑏,Θ,𝐼)𝑑Θ, so that 𝑝(Θ∣𝑧^𝑏,𝐳,𝐼)=𝑝(Θ∣𝐼)𝑝(𝐳∣𝑧^𝑏,Θ,𝐼)/𝑍). The value for the evidence 𝑍 should not be too small, as this would imply that the “overlap” of 𝑝(Θ∣𝐼) and 𝑝(𝐳∣𝑧^𝑏,Θ,𝐼) in the hypothesis space (involving hypotheses concerning Θ) is minimal. This would suggest that the prior distribution of Θ and the measurements 𝐳 assimilated through the likelihood function are inconsistent).
3. Actually, the technique of thermodynamic integration dates back to Kirkwood [28] who first used it within the context of the statistical mechanics of fluid mixtures to compute free energy differences.
4. The relative entropy is measured in natural units (nits) rather than the more usual binary units (bits) because a natural logarithm was used in the definition of the relative entropy, rather than a logarithm to the base 2 (see (5.8)).