International Scholarly Research Notices

International Scholarly Research Notices / 2011 / Article

Research Article | Open Access

Volume 2011 |Article ID 189735 |

Eugene Yee, Ajith Gunatilaka, Branko Ristic, "Comparison of Two Approaches for Detection and Estimation of Radioactive Sources", International Scholarly Research Notices, vol. 2011, Article ID 189735, 30 pages, 2011.

Comparison of Two Approaches for Detection and Estimation of Radioactive Sources

Academic Editor: Y. Dimakopoulos
Received25 May 2011
Accepted17 Jun 2011
Published25 Oct 2011


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 π’œβŠ‚β„2 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 π’œβŠ‚β„2 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)1 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 (π‘₯𝑙,𝑗,𝑦𝑙,𝑗)βˆˆβ„2. 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𝑑𝑗𝑖=ξ‚€ξ€·π‘₯𝑙,π‘—βˆ’π‘₯𝑠,𝑖2+𝑦𝑙,π‘—βˆ’π‘¦π‘ ,𝑖2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+𝑦𝑙,π‘—βˆ’π‘¦π‘ ,𝑖2.(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πœ‹πœŽπ‘—βŽ§βŽͺ⎨βŽͺβŽ©βˆ’1exp2π‘šξ“π‘—=1ξƒ©ξπœ‡π‘—π‘ βˆ’ξ‚πœ‡π‘—π‘ πœŽπ‘—ξƒͺ2⎫βŽͺ⎬βŽͺ⎭.(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 𝜎1.

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π‘ξ€·πœ‡π‘ξ€Έξ€·πœ‡βˆ£πΌ=π‘π‘βˆ£π‘§π‘,πΌπ‘ξ€Έβˆπœπ‘ξ€·πœ‡π‘πœπ‘ξ€Έπ‘§π‘βˆ’1ξ€·expβˆ’πœ‡π‘πœπ‘ξ€Έ.(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 𝜎1 in the estimate for πœ‡π‘  given in (2.26) can be determined as𝜎1=ξ‚€ξ«πœ‡2π‘ ξ¬βˆ’ξ€·ξπœ‡π‘ ξ€Έ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πœ‹πœŽπ‘—βŽ§βŽͺ⎨βŽͺβŽ©βˆ’1exp2π‘šξ“π‘—=1ξƒ©ξπœ‡π‘—π‘ βˆ’ξ‚πœ‡π‘—π‘ πœŽπ‘—ξƒͺ2⎫βŽͺ⎬βŽͺ⎭+𝕀{π‘Ÿ=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

The performance of the procedure depends somewhat on the number of steps 𝑆 and the expansion factors 𝛾1,…,𝛾𝑆. 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πœ’0ξ€·=βˆ’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 𝜽1 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 𝜽2 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 𝜽2. Towards this purpose, we update 𝜽2π‘˜βˆˆπœ½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:𝑝𝑏1(π‘Ÿ)=2ξ‚»min1,𝑝(π‘Ÿ+1∣𝐼)𝑝(π‘Ÿβˆ£πΌ),π‘Ÿ=π‘Ÿmin,…,π‘Ÿmaxπ‘βˆ’1,𝑑1(π‘Ÿ+1)=2ξ‚»min1,𝑝(π‘Ÿβˆ£πΌ)𝑝(π‘Ÿ+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)𝑏,π‘‘βˆ’βˆ’βˆ’β†’πš―β‹†β„³1βˆ’βˆ’β†’πš―β„³β‹†β‹†2βˆ’βˆ’β†’πš―(𝑑),(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; β„³1 denotes the update of the source intensity rates using a Gibbs sampler; and, β„³2 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 𝑝0(Ξ˜βˆ£π‘§π‘,𝐳,𝐼). 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=π‘π›Ώπœ†ξ€·π³βˆ£π‘§π‘,πš―π‘˜ξ€ΈΓ—ξƒ©(πœ†),πΌπ‘βˆ—ξ“π‘—=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 (𝑛1,𝑛2,…,π‘›π‘βˆ—) of copies of the samples Ξ˜π‘˜(πœ†) (π‘˜=1,2,…,π‘βˆ—) is obtained by computing a vector (𝜌1,𝜌2,…,πœŒπ‘βˆ—) of the cumulative sums of π‘βˆ—Γ—(𝑀1,𝑀2,…,π‘€π‘βˆ—), generating a random draw from the uniform distribution π‘’βˆΌπ’°([0,1]), and determining π‘›π‘˜ (π‘˜=1,2,…,π‘βˆ—) fromπ‘›π‘˜=βŒŠπœŒπ‘˜+π‘’βŒ‹βˆ’βŒŠπœŒπ‘˜βˆ’1+π‘’βŒ‹,π‘˜=2,3,…,π‘βˆ—π‘›βˆ’1,1=⌊𝜌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=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 Θ(0) 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βˆ’1 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βˆ’1. We converted these dose rate data into raw count measurements π‘§π‘—βˆˆβ„•βˆͺ{0} by multiplication by a conversion factor (which was 0.21βˆ’1 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 m2hβˆ’1, for Sources 1, 2, and 3, respectively. We will use these values as the β€œtrue” source strengths in our source estimations.

Source Type Activity (MBq)

1 1 3 7 Cs 2 6 Γ— 1 0 3
2 1 3 7 Cs 5 Γ— 1 0 3
3 6 0 Co 0 . 2 Γ— 1 0 3

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.

Test set Source 1 Location Source 2 Location Source 3 Location

1 ( 1 1 , 1 0 )  m β€”β€”
2 ( 1 1 , 1 0 )  m ( 3 , 5 0 )  mβ€”
3 ( 1 1 , 1 0 )  m ( 3 , 5 0 )  m ( 4 1 , 5 )  m

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).

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βˆ’1.

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 m2 hβˆ’1. 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πœ‹πœŽπ‘—βŽ§βŽͺ⎨βŽͺβŽ©βˆ’1exp2π‘šξ“π‘—=1ξƒ©ξπœ‡π‘—π‘ βˆ’ξ‚πœ‡π‘—π‘ πœŽπ‘—ξƒͺ2⎫βŽͺ⎬βŽͺβŽ­Γ—ξ€·π‘Ÿmaxβˆ’π‘Ÿminξ€Έ!ξ€·π‘Ÿβˆ’π‘Ÿminξ€Έ!ξ€·π‘Ÿmaxξ€Έ!π‘βˆ’π‘Ÿβˆ—(π‘Ÿβˆ’π‘Ÿmin)ξ€·1βˆ’π‘βˆ—ξ€Έπ‘Ÿmaxβˆ’π‘ŸΓ—π‘Ÿξ‘π‘–=1𝑄(1βˆ’π›Ύ)𝛿𝑠,𝑖+𝛾𝕀(0,𝑄max)𝑄𝑠,𝑖𝑄maxΓ—π•€π’œξ€·π±π‘ ,𝑖meas(π’œ)+𝕀{π‘Ÿ=0}(π‘Ÿ)π‘šξ‘π‘—=1πœπ‘§π‘π‘πœπ‘§π‘—ξ€·πœ+πœπ‘ξ€Έπ‘§π‘—+𝑧𝑏𝑧𝑗+𝑧𝑏!βˆ’1𝑧𝑗!𝑧𝑏!Γ—ξ€·βˆ’11βˆ’π‘βˆ—ξ€Έπ‘Ÿ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 m2 hβˆ’1, 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 Θ(1),Θ(2),… 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 determined2. This can be accomplished using an approach known as thermodynamic integration (see, von der Linden [25]).3 In this approach, 𝑍 can be evaluated as follows:ξ€œlog𝑍=10𝑝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=0,πœ†1,…,πœ†πΎβˆ’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βˆ’1 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.

Parameter Mean Standard deviation 95% HPD Actual

𝜏 = 1 6  s

  𝑖 = 1
π‘₯ 𝑠 10.7 0.2 ( 1 0 . 3 , 1 1 . 0 ) 11.0
𝑦 𝑠 9.9 0.2 ( 9 . 5 , 1 0 . 4 ) 10.0
𝑄 𝑠 1930 43 ( 1 8 4 8 , 2 0 1 4 ) 1912

  𝑖 = 2
π‘₯ 𝑠 2.1 0.6 ( 1 . 0 , 3 . 2 ) 3.0
𝑦 𝑠 49.8 0.3 ( 4 9 . 3 , 5 0 . 3 ) 50.0
𝑄 𝑠 330 31 ( 2 6 9 , 3 9 2 ) 392

  𝑖 = 3
π‘₯ 𝑠 40.9 0.1 ( 4 0 . 7 , 4 1 . 1 ) 41.0
𝑦 𝑠 5.0 0.3 ( 4 . 5 , 5 . 6 ) 5.0
𝑄 𝑠 92.9 15.2 ( 6 2 . 3 , 1 2 2 . 7 ) 98.1

𝜏 = 4 8  s

  𝑖 = 1
π‘₯ 𝑠 10.9 0.01 ( 1 0 . 8 9 , 1 0 . 9 0 ) 11.0
𝑦 𝑠 9.6 0.01 ( 9 . 6 2 , 9 . 6 3 ) 10.0
𝑄 𝑠 1906 1.3 ( 1 9 0 5 , 1 9 0 7 ) 1912

  𝑖 = 2
π‘₯ 𝑠 2.3 0.01 ( 2 . 3 6 , 2 . 3 8 ) 3.0
𝑦 𝑠 49.6 0.01 ( 4 9 . 6 0 , 4 9 . 6 2 ) 50.0
𝑄 𝑠 336 0.6 ( 3 3 5 . 3 , 3 3 6 . 2 ) 392

  𝑖 = 3
π‘₯ 𝑠 40.8 0.00 ( 4 0 . 8 7 , 4 0 . 8 7 ) 41.0
𝑦 𝑠 5.3 0.01 ( 5 . 2 6 , 5 . 2 8 ) 5.0
𝑄 𝑠 83.6 0.73 ( 8 3 . 0 , 8 4 . 2 ) 98.1

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)).

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.

Parameter Mean Standard deviation 95% HPD Actual

𝜏 = 1 6  s

  𝑖 = 1
π‘₯ 𝑠 8.7 0.8 ( 7 . 2 , 1 0 . 6 ) 11.0
𝑦 𝑠 10.6 0.8 ( 9 . 1 , 1 2 . 1 ) 10.0
𝑄 𝑠 1404 152 ( 1 1 0 9 , 1 7 0 4 ) 1912

  𝑖 = 2
π‘₯ 𝑠 3.0 0.6 ( 3 . 0 , 5 . 0 ) 3.0
𝑦 𝑠 49.8 0.6 ( 4 9 . 0 , 5 0 . 9 ) 50.0
𝑄 𝑠 385 81 ( 2 2 8 , 5 4 8 ) 392

  𝑖 = 3
π‘₯ 𝑠 40.9 0.5 ( 3 9 . 9 , 4 1 . 6 ) 41.0
𝑦 𝑠 5.0 0.9 ( 3 . 5 , 6 . 4 ) 5.0
𝑄 𝑠 194 46 ( 1 0 5 , 2 8 6 ) 98.1

𝜏 = 4 8  s

  𝑖 = 1
π‘₯ 𝑠 9.3 0.9 ( 7 . 6 , 1 1 . 0 ) 11.0
𝑦 𝑠 10.2 0.8 ( 8 . 7 , 1 1 . 9 ) 10.0
𝑄 𝑠 1433 156 ( 1 1 3 0 , 1 7 4 2 ) 1912

  𝑖 = 2
π‘₯ 𝑠 3.8 0.6 ( 3 . 0 , 4 . 9 ) 3.0
𝑦 𝑠 49.8 0.6 ( 4 9 . 0 , 5 0 . 9 ) 50.0
𝑄 𝑠 382 80 ( 2 2 5 , 5 4 3 ) 392

𝑖 = 3
π‘₯ 𝑠 41.1 0.5 ( 4 0 . 0 , 4 1 . 8 ) 41.0
𝑦 𝑠 4.8 0.8 ( 3 . 5 , 6 . 3 ) 5.0
𝑄 𝑠 204 52 ( 1 0 4 , 3 0 9 ) 98.1

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 xs- and ys-directions is generally very good. The recovery of the source intensity rates Qs is quite good, although both the accuracy and precision in this recovery are poorer than the recovery of the source locations.

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.4 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Γ—1023 and 1.26Γ—1026, 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]).


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.


  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)).


  1. W. K. H. Panofsky, β€œNuclear proliferation risks, new and old,” Issues in Science and Technology, vol. 19, no. 4, pp. 73–75, 2003. View at: Google Scholar
  2. J. W. Howse, L. O. Ticknor, and K. Muske, β€œLeast squares estimation techniques for position tracking of radioactive sources,” Automatica, vol. 37, no. 11, pp. 1727–1737, 2001. View at: Publisher Site | Google Scholar | Zentralblatt MATH
  3. D. L. Stephens Jr. and A. J. Peurrung, β€œDetection of moving radioactive sources using sensor networks,” IEEE Transactions on Nuclear Science, vol. 51, no. 5 I, pp. 2273–2278, 2004. View at: Publisher Site | Google Scholar
  4. R. J. Nemzek, J. S. Dreicer, D. C. Torney, and T. T. Warnock, β€œDistributed sensor networks for detection of mobile radioactive sources,” IEEE Transactions on Nuclear Science, vol. 51, no. 4, pp. 1693–1700, 2004. View at: Publisher Site | Google Scholar
  5. S. M. Brennan, A. M. Mielke, and D. C. Torney, β€œRadioactive source detection by sensor networks,” IEEE Transactions on Nuclear Science, vol. 52, no. 3, pp. 813–819, 2005. View at: Publisher Site | Google Scholar
  6. A. Gunatilaka, B. Ristic, and R. Gailis, β€œOn localisation of a radiological point source,” in Proceedings of the Information, Decision and Control, (IDC '07), pp. 236–241, Adelaide, Australia, February 2007. View at: Publisher Site | Google Scholar
  7. M. Morelande, B. Ristic, and A. Gunatilaka, β€œDetection and parameter estimation of multiple radioactive sources,” in Proceedings of the 10th International Conference on Information Fusion, (FUSION '07), Quebec, Canada, July 2007. View at: Publisher Site | Google Scholar
  8. M. Morelande and B. Ristic, β€œRadiological source detection and localisation using Bayesian techniques,” IEEE Transactions on Signal Processing, vol. 57, no. 11, pp. 4220–4231, 2009. View at: Publisher Site | Google Scholar
  9. C. Mendis, A. Gunatilaka, B. Ristic, S. Karunasekera, and A. Skvortsov, β€œExperimental verification of evolutionary estimation algorithms for radioactive source localisation,” in Proceedings of the 5th International Conference on Intelligent Sensors, Sensor Networks and Information Processing, (ISSNIP '09), pp. 151–156, Melbourne, Australia, December 2009. View at: Publisher Site | Google Scholar
  10. N. Tsoulfanidis, Measurement and Detection of Radiation, Taylor & Francis, 1995.
  11. A. Martin and S. A. Harbison, An Introduction to Radiation Protection, Chapman & Hall, 3rd edition, 1986.
  12. E. T. Jaynes, Probability Theory: The Logic of Science, Cambridge University Press, Cambridge, UK, 2003. View at: Publisher Site
  13. T. J. Loredo, β€œFrom Laplace to supernova SN 1987A: Bayesian inference in astrophysics,” in Maximum Entropy and Bayesian Methods, P. F. FougΓ¨re, Ed., pp. 81–142, Kluwer Academic Publishers, 1990. View at: Google Scholar | Zentralblatt MATH
  14. C. Musso, N. Oudjane, and F. LeGland, β€œImproving regularized particle filters,” in Sequential Monte Carlo Methods in Practice, A. Doucet, N. deFreitas, and N. J. Gordon, Eds., Springer, New York, NY, USA, 2001. View at: Google Scholar
  15. G. Kitagawa, β€œMonte Carlo filter and smoother for non-Gaussian nonlinear state space models,” Journal of Computational and Graphical Statistics, vol. 5, no. 1, pp. 1–25, 1996. View at: Publisher Site | Google Scholar
  16. B. Ristic, S. Arulampalam, and N. Gordon, Beyond the Kalman Filter: Particle Filters for Tracking Applications, Artech House, 2004.
  17. R. Mahler, Statistical Multisource Multitarget Information Fusion, Artech House, 2007.
  18. S. M. Kay, Statistical Signal Processing: Detection Theory, Prentice-Hall, 1998.
  19. J. Rissanen, β€œModelling by shortest data description,” Automatica, vol. 14, pp. 465–471, 1978. View at: Publisher Site | Google Scholar
  20. P. J. Green, β€œReversible jump Markov chain Monte Carlo computation and Bayesian model determination,” Biometrika, vol. 82, no. 4, pp. 711–732, 1995. View at: Publisher Site | Google Scholar | Zentralblatt MATH
  21. E. Yee, β€œInverse dispersion of an unknown number of continuous sources,” in Proceedings of the 15th Joint Conference of Air Pollution Meteorology with the A&WMA, p. 17, New Orleans, La, USA, 2008, paper no. 7.1. View at: Google Scholar
  22. E. Yee, β€œTheory for reconstruction of an unknown number of contaminant sources using probabilistic inference,” Boundary-Layer Meteorology, vol. 127, no. 3, pp. 359–394, 2008. View at: Publisher Site | Google Scholar
  23. D. Gamerman and H. F. Lopes, Markov Chain Monte Carlo: Stochastic Simulation for Bayesian Inference, Chapman & Hall, 2nd edition, 2006.
  24. AN/PDR-77 User’s Guide, Canberra Industries, Meriden, Conn, USA, 2000.
  25. W. von der Linden, R. Fischer, and V. Dose, β€œEvidence integrals,” in Maximum Entropy and Bayesian Methods, K. M. Hanson and R. N. Silver, Eds., Kluwer Academic Publishers, 1996. View at: Google Scholar | Zentralblatt MATH
  26. T. M. Cover and J. A. Thomas, Elements of Information Theory, John Wiley & Sons, New York, NY, USA, 1991. View at: Publisher Site
  27. B. Ristic, M. Morelande, and A. Gunatilaka, β€œInformation driven search for point sources of gamma radiation,” Signal Processing, vol. 90, no. 4, pp. 1225–1239, 2010. View at: Publisher Site | Google Scholar | Zentralblatt MATH
  28. J. G. Kirkwood, β€œStatistical mechanics of fluid mixtures,” The Journal of Chemical Physics, vol. 3, no. 5, pp. 300–313, 1935. View at: Publisher Site | Google Scholar | Zentralblatt MATH

Copyright Β© 2011 Her Majesty the Queen in Right of Canada. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

More related articles

 PDF Download Citation Citation
 Download other formatsMore
 Order printed copiesOrder

Related articles

We are committed to sharing findings related to COVID-19 as quickly as possible. We will be providing unlimited waivers of publication charges for accepted research articles as well as case reports and case series related to COVID-19. Review articles are excluded from this waiver policy. Sign up here as a reviewer to help fast-track new submissions.