Localization of a radioactive source in an urban environment using Bayesian Metropolis methods

Localization of a radioactive source in an urban environment using Bayesian Metropolis methods

Nuclear Inst. and Methods in Physics Research, A ( ) – Contents lists available at ScienceDirect Nuclear Inst. and Methods in Physics Research, A ...

4MB Sizes 0 Downloads 26 Views

Nuclear Inst. and Methods in Physics Research, A (

)



Contents lists available at ScienceDirect

Nuclear Inst. and Methods in Physics Research, A journal homepage: www.elsevier.com/locate/nima

Localization of a radioactive source in an urban environment using Bayesian Metropolis methods Jason Hite a ,∗, John Mattingly a , Dan Archer b , Michael Willis b , Andrew Rowe b , Kayleigh Bray b , Jake Carter b , James Ghawaly b a b

Department of Nuclear Engineering, North Carolina State University, Raleigh, NC 27695, United States Oak Ridge National Laboratory, Oak Ridge, TN 37830, United States

ARTICLE

INFO

Keywords: Source localization Bayesian parameter estimation Sensor networks

ABSTRACT We present a method for localizing an unknown source of radiation in an urban environment using a distributed detector network. This method employs statistical parameter estimation techniques, relying on an approximation for the response of a detector to the source based on a simplified model of the underlying transport phenomena, combined with a Metropolis-type sampler that is modified to propagate the effect of fixed epistemic uncertainties in the material macroscopic cross sections of objects in the scene. We apply this technique to data collected during a measurement campaign conducted in a realistic analog for an urban scene using a network of six mobile detectors. Our initial results are able to localize the source to within approximately 8 m over a scene of size 300 m × 200 m in two independent trials with 30 min count times, including a characterization of the uncertainty associated with the poorly known macroscopic cross sections of objects in the scene. In these measurements, the nearest detectors were between 20 m to 30 m from the source, and recorded count rates between approximately 3 and 13 times background. A few detectors had line-of-sight to the source, while the majority were obscured by objects present in the scene. After extending our model to account for the orientation of the detectors and correcting for anomalies in the measurement data we were able to further improve the localization accuracy to approximately 2 m in both trials.

1. Introduction In this article, we investigate determining the location and intensity of an unknown source of radiation in a heterogeneous environment. A source is assumed to exist in the search area and is observed by a network of radiation detectors, which record counts due to both the source and to natural background radiation present in the scene. Using these measurements, we seek to predict the location of the source, as well as to quantify the uncertainty in our estimate. In previous work, we have demonstrated an algorithm based on Markov-Chain Monte Carlo to perform the localization [1,2]. These results were based on a simulated city block, with synthetic measurement data generated using a simplified model for the transport of gamma rays from the source to the detector in a heterogeneous medium. In May 2017, with the assistance of Oak Ridge National Laboratory (ORNL), we performed a series of outdoor measurements intended to serve as a field test of our localization algorithm. These measurements were designed to simulate a real-world search scenario in an urban setting. We present the results of these experiments here, with a focus on the analysis of the

discrepancies observed when moving from synthetic measurement data generated using a simplified transport model to actual measurements taken in the field. 1.1. Prior work Localizing an emissive source based on the measured signal strength of a network of detectors is a topic that has been widely investigated, with a variety of algorithms described in the literature. The simplest method is trilateration, where the signal strength reported by each detector is used to infer a distance from the source to the detector. This creates a set of overlapping circles, each centered on a detector, with radius given by the estimated distance from the source to the corresponding detector. The source is then localized by finding the point of intersection for the set of circles. Variations of this method have been used by the nuclear security community for localizing a radioactive source [3,4] and by other communities for numerous applications, including localization of mobile phones and other wireless devices [5].

∗ Corresponding author. E-mail address: [email protected] (J. Hite).

https://doi.org/10.1016/j.nima.2018.09.032 Received 15 January 2018; Received in revised form 16 August 2018; Accepted 11 September 2018 Available online xxxx 0168-9002/© 2018 Published by Elsevier B.V.

Please cite this article in press as: J. Hite, et al., Localization of a radioactive source in an urban environment using Bayesian Metropolis methods, Nuclear Inst. and Methods in Physics Research, A (2018), https://doi.org/10.1016/j.nima.2018.09.032.

J. Hite et al.

Nuclear Inst. and Methods in Physics Research, A (

Trilateration is essentially a geometric solution to the source localization problem and is sensitive to uncertainties in the signal strength reported by the detectors [6]. This sensitivity is especially important when concerned with radioactive sources, where the count rates measured by the detector can be highly uncertain, particularly for short counting times and/or long standoff distances. An alternative approach is to cast source localization in the framework of statistical parameter estimation. Broadly, this involves the construction of a statistical model describing the detector measurements as random variables, with distributions that depend on the characteristics of the source (e.g., location and activity). The statistical model can also be extended to account for uncertainty and the resulting estimates of the source parameters are typically less sensitive to noise in the measurement data. The classic solution to the statistical parameter estimation problem is the maximum likelihood estimate (MLE), which is defined as the source parameters that maximize the likelihood of observing the detector responses measured in the field. MLE methods are very popular due to their robustness and the reliability of modern numerical optimization algorithms. They have seen wide application, with Ref. [7] demonstrating the use of an MLE for a detector response model dependent on the source activity and the distance from the source to the detector. Recent work has also applied MLE to the complementary problem of estimating the spatial and temporal distribution of background and subsequently using this estimate to detect an anomalous source using a mobile detector network [8]. Bayesian techniques have also seen a growth in applications in recent years, enabled largely by the increasing computational resources available to researchers. Ref. [9] provides an early example, deriving a maximum a posteriori estimator (related to the MLE, but formulated in a Bayesian context, see Section 2.4) for the source location intended for use in real-time tracking. Similarly, Ref. [4] describes the use of Bayesian model selection techniques to distinguish the presence of multiple sources in the scene and localize them via a Monte Carlo estimator. Lastly, there is a separate class of Bayesian algorithms collectively known as particle filtering. This approach can be understood as a type of genetic algorithm, where several ‘‘particles’’ are distributed throughout the search space, with each particle representing a candidate source location. Measured data is used to determine the fitness of each particle by computing its likelihood and the system is evolved to produce a new generation of particles; over several generations, the particles tend to cluster around the true source location. This method has proven especially useful in online monitoring, with Ref. [10] demonstrating the algorithm applied to a vehicle portal monitoring scenario.

2.1. Detector response model Evaluation of 𝐃 (𝐫, 𝐼 ; 𝛴) during the sampling process requires a computational model of the counts recorded by the detector network due to the source, which is a function of the source location and intensity. This model will be evaluated repeatedly during the sampling procedure (hundreds of thousands of times in our application), and in a way that is primarily serial by construction.1 This necessitates the use of a computationally inexpensive model for the underlying transport process; in previous work we have demonstrated the use of a point source/detector transport model, which only computes the uncollided flux arriving at the detector [1]. We assume that the detectors are sufficiently small that at typical standoff distances (tens to hundreds of meters), a gamma ray that is initially emitted within the solid angle of a detector and subsequently scatters is very unlikely to reach the detector, hence the detector response is dominated by the uncollided flux arriving at the detector. Under these assumptions, the model for the response of the 𝑖th detector to the source takes the form of exponential attenuation given in Eq. (2), with parameters defined in Table 1 [13,15], where the total attenuation factor is calculated as a line integral of 𝛴(𝐫 ′ ) along the line drawn from source to detector, denoted 𝐫 → 𝐫𝑖 . Note that this model includes a parametric dependence on 𝛴(𝐫 ′ ), the macroscopic cross section of the scene at position 𝐫 ′ . Clearly the detector response will be influenced by these cross-sections, however, these are treated separately in the sampling algorithm since they are considered to be fixed epistemic uncertainties arising from our imprecise knowledge of the composition and internal dimensions of objects in the scene (see Section 2.3). ( ) 𝐴𝑖 int ′ 𝑑𝑖 (𝐫, 𝐼 ; 𝛴) = 𝐼𝛥𝑡𝑖 ⋅ 𝜖𝑖 ⋅ ⋅ exp − 𝛴(𝐫 ) 𝑑𝑠 (2) ∫𝐫→𝐫𝑖 4𝜋‖𝐫 − 𝐫𝑖 ‖2

We formulate the problem of estimating the source location 𝐫0 and intensity 𝐼0 in the context of Bayesian parameter estimation methods. Assume that we have a network of 𝑁𝐷 detectors with known and fixed positions. The response of this detector network is considered to be a vector-valued random variable , where we assume that each 𝑖 , corresponding to the counts measured by the 𝑖th detector (𝑖 = 1, 2, … , 𝑁𝐷 ), is an independent Poisson-distributed random variable. We are provided with a vector of count data that is measured in the field by the detector network, denoted 𝐃, which represents a realization of . We then seek to compute the posterior density 𝑃 (𝐫, 𝐼 ∣ 𝐃) (henceforth, the posterior), which is the probability that a given (𝐫, 𝐼) is the true source location and intensity, conditioned on the observations 𝐃. The posterior can be expressed directly via Bayes’ formula as a normalized product of a prior density 𝑃0 (𝐫, 𝐼) and a likelihood function 𝐃 (𝐫, 𝐼 ; 𝛴) [11]: 𝐃 (𝐫, 𝐼 ; 𝛴) ⋅ 𝑃0 (𝐫, 𝐼) 𝐼



recorded, while X and (𝐼min , 𝐼max ) represent the limits for the position and intensity of the source. In the experiments described in Section 3 we employ a uniform prior distribution, corresponding to the assumption that we initially only know a wide area where the source is located and nothing more. Such a prior is usually referred to as uninformative since it assumes that all points in the search space are equally probable. 𝐃 (𝐫, 𝐼 ; 𝛴) is called the likelihood function and employs a statistical model of the detector response to approximate the probability that the measurements 𝐃 would be observed if the source is located at 𝐫 with intensity 𝐼. In this model we assume that the detector responses are governed by Poisson counting statistics and depend on the source location and intensity, as well as 𝛴(𝐫 ′ ), the total macroscopic crosssection of materials in the scene at location 𝐫 ′ . Direct evaluation of Eq. (1) is typically impractical due to the difficulty of evaluating the integral normalization term in the denominator. Instead, we can employ a Metropolis-type sampler to draw samples from the posterior via rejection, a process typically referred to as Markov Chain Monte Carlo (MCMC) [12]. This sampler avoids computing the integral in Eq. (1) by evaluating ratios of 𝑃 (𝐫, 𝐼 ∣ 𝐃) at different points in the domain, causing the integral terms to cancel [11,12]. For our use, we deviate slightly from the traditional Metropolis algorithm and instead employ a variation modified to include fixed parameter distributions. This is due to the impact of uncertainty in 𝛴(𝐫 ′ ), which we consider to be a fixed epistemic uncertainty arising from our inability to precisely know the cross-section of materials in the scene. A traditional Metropolis sampler which included sampling of the material crosssections would also ‘‘update’’ the probability distribution for 𝛴(𝐫 ′ ) to reflect the knowledge obtained from 𝐃, whereas the modified sampler holds the uncertainty in 𝛴(𝐫 ′ ) constant and propagates it onto the posterior distribution for 𝐫 and 𝐼. A full description of the localization algorithm is provided in Ref. [13].

2. Methodology

𝑃 (𝐫, 𝐼 ∣ 𝐃) =

)

2

(1)

1 There has been some success in designing parallel Metropolis samplers beyond simply running multiple independent chains, for example Ref. [14], however many of these modifications either achieve only modest parallel utilization or rely on assumptions that are not valid in our application.

∫𝐼 max ∫X 𝐃 (𝐫, 𝐼 ; 𝛴) ⋅ 𝑃0 (𝐫, 𝐼) 𝑑𝐫 𝑑𝐼 min

𝑃0 (𝐫, 𝐼) is a probability distribution which represents our knowledge about the source location and intensity before any measurements are 2

Please cite this article in press as: J. Hite, et al., Localization of a radioactive source in an urban environment using Bayesian Metropolis methods, Nuclear Inst. and Methods in Physics Research, A (2018), https://doi.org/10.1016/j.nima.2018.09.032.

J. Hite et al.

Nuclear Inst. and Methods in Physics Research, A (

Table 1 Model parameters.



2.3. Bayesian uncertainty and material macroscopic cross sections

Parameter

Meaning

𝑑𝑖 𝛥𝑡𝑖 𝜖𝑖int 𝐴𝑖 𝐫𝑖 𝛴(𝐫 ′ )

Counts recorded by detector due to the source (counts) Measurement time (s) Detector intrinsic efficiency (%) Detector face area (m2 ) Detector position vector (m) Total cross-section at 𝐫 ′ (1∕m)

We have described uncertainties in the source position and intensity using a strictly Bayesian interpretation of uncertainty. In this context, the uncertainty is a measure of confidence in the estimates for 𝐫 and 𝐼. We begin with a very high degree of uncertainty in 𝐫 and 𝐼 and gradually reduce the uncertainty as information is accumulated through measurements made by the detector network, with the decrease in uncertainty being weighted by the corresponding uncertainty of each measurement (Eq. (1)). Our method deviates from the traditional formulation in the treatment of the macroscopic cross sections of other objects in the scene, which determine the mean-free path of gammas as they travel from source to detector. Accurate and precise knowledge of the composition and density of buildings and objects in the scene is unavailable, so instead we employ rough estimates based on typical internal dimensions and construction materials [19]. Because these estimates are imprecise, we consider them to be highly uncertain (that is, we have a low degree of belief in their accuracy) and we wish to propagate this uncertainty onto the estimated posterior for the source location and intensity. This differs from the uncertainty in the source location in the sense that the macroscopic cross section uncertainties are fixed measurements from the detector network do not improve our confidence in the estimated macroscopic cross sections. At each step of the sampling process, the macroscopic cross sections are perturbed independently from the current state of the sampler by randomly sampling from a probability distribution that is specified by the user. This process therefore combines the MCMC sampler with a traditional forward Monte Carlo error propagation and produces a posterior distribution that also incorporates the uncertainty in the cross section values. In this context, it is informative to think of uncertainty as a measure of the ‘‘resolution’’ of our prediction for the source location. The input uncertainties we model (macroscopic cross sections, counting uncertainties for a fixed measurement time) result in a fundamental limit to how precisely we can identify the location of the source. With longer measurements and more precise estimates of the macroscopic cross sections the precision of our prediction is increased and we are able to identify a smaller area in which we predict the source is located, ideally up to the limiting case of infinitely long measurements and perfectly known macroscopic cross sections where the posterior becomes a delta function. Separate from that we have the accuracy of said prediction relative to the true source location, which is affected by how well our model of the detector response predicts the actual source location, an effect that is neither uncertain nor random.

The response of the detector to the source is assumed to be a Poissondistributed random variable with mean 𝑑𝑖 (𝐫, 𝐼 ; 𝛴). We also allow for a random contribution from background that is independent of the source, represented by the random variable 𝑖 . We typically assume 𝑖 is also Poisson distributed with a mean 𝑏𝑖 , which is allowed to vary with detector location and time,2 and hence our statistical model for the total detector response 𝑖 is [ ] 𝑖 = Po 𝑑𝑖 (𝐫, 𝐼 ; 𝛴) + 𝑖 [ ] [ ] = Po 𝑑𝑖 (𝐫, 𝐼 ; 𝛴) + Po 𝑏𝑖 with every 𝑖 assumed to be mutually independent. 2.2. Effect of detector orientation Eq. (2) relies on an approximation for the solid angle of the 𝑖th detector, 𝛺𝑖 (𝐫), as viewed from the location of the source: 𝛺𝑖 (𝐫) ≈

)

𝐴𝑖 ‖𝐫 − 𝐫𝑖 ‖22

This assumes that the detector face area exposed to the source, 𝐴𝑖 , is approximately constant with respect to the relative positions of the source and detector. Such an assumption is true for a spherical detector and approximately true for a cubic detector; however, the experiments described in Section 3 used NaI detectors with crystal dimensions of 2′′ × 4′′ × 16′′ . The solid angle subtended by each of these detectors therefore depends strongly on the orientation of the detector in space about its own center and relative to the source location 𝐫. We will first demonstrate the simplest possible approach, where we use a constant face area computed by averaging the face areas of each detector; equivalently we assume that each detector has a face profile whose area is the average face area of that detector. In Section 4 we will see that this produces significant errors in the estimated source location, caused by variations in the measured count rates due to varying geo geometric efficiencies, 𝜖𝑖 ∶= 𝛺𝑖 ∕4𝜋, an effect that cannot be accounted for by this simplistic model. Therefore we extend our model to also account for variations in the detector solid angle using the method of Van Oosterom and Strackee [18], which provides an exact expression for the solid angle subtended by an arbitrarily oriented right triangle as viewed from a specified position in free space. Eq. (2) then becomes ( ) 𝛺 (𝐫) 𝑑𝑖 (𝐫, 𝐼 ; 𝛴) = 𝐼𝛥𝑡𝑖 ⋅ 𝜖𝑖int ⋅ 𝑖 ⋅ exp − 𝛴(𝐫 ′ ) 𝑑𝑠 (3) ∫𝐫→𝐫𝑖 4𝜋

2.4. Relation to maximum likelihood estimation As mentioned in Section 1, a common alternative approach to source localization is to compute the Maximum Likelihood Estimate (MLE), which is defined as the 𝐫, 𝐼 that maximize the likelihood function 𝐃 : (𝐫MLE , 𝐼MLE ) ∶= arg max 𝐃 (𝐫, 𝐼 ; 𝛴)

(4)

𝐫,𝐼

In frequentist inference the MLE is often used to provide an estimate of the model parameters without invoking prior information. From a Bayesian perspective, the MLE is related to the maximum a posteriori (MAP) estimate, defined as the 𝐫, 𝐼 that maximize the posterior given by Eq. (1):

where 𝛺𝑖 (𝐫) is computed by summing the solid angles of each face of the detector exposed to the source, calculated using the expression given in eq. (8) of Ref. [18].

(𝐫MAP , 𝐼MAP ) ∶= arg max 𝑃 (𝐫, 𝐼 ∣ 𝐃) 𝐫,𝐼

2

Note that we are treating the background as known, at least at the detector locations. In our experiments the background distribution was recorded independently in separate measurements and exhibited noticeable spatial variation, but only weak temporal dependence. Ref. [16] demonstrated a method for calibrating background online concurrently with measurements of the source, while we consider the technique described in [17] to be promising for constructing predictive models for background.

= arg max 𝐃 (𝐫, 𝐼 ; 𝛴) ⋅ 𝑃0 (𝐫, 𝐼)

(5)

𝐫,𝐼

As such, the MLE can be interpreted as a special case of the MAP estimate with a constant (i.e., uniform) prior distribution 𝑃0 , since multiplication by a constant does not change the location of the maximum in the parameter space. 3

Please cite this article in press as: J. Hite, et al., Localization of a radioactive source in an urban environment using Bayesian Metropolis methods, Nuclear Inst. and Methods in Physics Research, A (2018), https://doi.org/10.1016/j.nima.2018.09.032.

J. Hite et al.

Nuclear Inst. and Methods in Physics Research, A (

)



Fig. 2. Overhead view of experiment site, the two source locations are shown in orange, while the overall experiment area is marked in red (satellite imagery courtesy of Google Maps) . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

2.6. Algorithm performance

Fig. 1. Sum-of-squares error as a function of source location for experiment 1 described in Section 3. The peaks are clipped to show detail of the error surface.

Finally, let us conclude discussion of the MCMC source localization algorithm by briefly addressing performance. As alluded to at the beginning of this section, MCMC-based algorithms are inherently computationally intensive, particularly when compared to techniques based on deterministic numerical optimization. Our method is no exception and requires a typical run-time of approximately 45 min for each case shown in this report when run on a typical desktop computer, versus the 2 to 3 min we achieved using hybrid optimization algorithms on a similar model geometry in Ref. [20]. This is qualified by the consideration that the computational cost of the MCMC sampler is dominated by the evaluation of the detector response model in Eq. (2). Our current implementation of the detector response model is implemented in Python and designed for flexibility; it is not optimized for performance and we expect that significant improvements are possible.

The output of our source localization algorithm is the full posterior distribution, but for the sake of discussion we typically also identify the mode of the posterior as ‘‘the’’ predicted source location since it is, by definition, the source location we are most confident is the true one considering the count rates we observed. It is clear from Eq. (5) that the mode of the posterior is also the MAP estimate. As noted previously, for our analysis we assume uniformly distributed prior distributions for all parameters, hence the mode of the posterior distribution is also the MLE. The distinction is that our method provides additional information since the full posterior distribution is available, whereas the MLE is only a point estimate. The MLE is also incompatible with other forms of prior distribution, which are of interest in situations where additional information about the source characteristics is available.

3. Description of experiment 2.5. Difficulties with direct numerical optimization In May 2017, we performed a measurement campaign at the Energy Systems Test Complex at ORNL, also known as the site of the former Experimental Gas-Cooled Reactor. The measurement campaign took place outdoors in a cluttered environment chosen to mimic a real-world search in an urban setting. We placed a Cs-137 source with a nominal activity of approximately 37 mCi in two different locations, shown in orange in Fig. 2. These source placements divided the measurement campaign into two separate experiments, which we refer to as experiment 1 and experiment 2. For each experiment, we used a set of 6 networked mobile detectors (described in Section 3.1) to record count rates at 12 locations, both with and without the source present. Measurements with the source and of background were each taken for approximately 30 min realtime, and count rates were low enough that dead-time effects were negligible (typically < 0.1%). The placement of these detectors is shown in Fig. 3, with all detector positions recorded using precision differential GPS. Source-to-detector distances ranged from approximately 20 m to 200 m. Additionally, all detectors were nominally oriented to be facing due north using hand compasses, since the source location is treated as unknown. Manual orientations were error-prone, with substantial deviations from north in many of the detectors; however, the actual orientations were recorded using the detectors’ onboard compasses, mounted such that they record orientation relative to the longest axis of the detector. The actual placements of the detectors in the field and their orientations as reported by the onboard compasses are summarized in Table 2. Each location is given a name denoting the experiment number, measurement set and detector number — for example, measurement 1A-2

The consideration of attenuators significantly complicates the application of many of the alternative techniques for source localization described in Section 1, which generally do not account for variable attenuation in a heterogeneous environment. In particular, the most common approach to compute the MLE defined in Section 2.4 is to perform a direct numerical maximization of the likelihood function,3 but this becomes difficult when attenuators are included due to the non-smooth error surface that results. As an illustration, Fig. 1 plots the error surface in terms of the sum-of-squares error (𝑆𝑆𝐸) for the measurements collected in the first of the experiments described in Section 3. The complexity of this error surface, which includes discontinuities along every object boundary as well as singularities at the detector locations, makes it difficult to provide a direct comparison to methods based on numerical optimization. For instance, in principle the MLE produced by direct numerical optimization should agree with the mode of the posterior density produced by our algorithm (see Section 2.4); however, when we implemented this the search failed to converge in all cases using a variety of optimization algorithms including Nelder–Mead, Broyden–Fletcher–Goldfarb–Shanno (BFGS), and Simulated Annealing. Interested readers should see Ref. [20] for more information regarding the numerical challenges involved, as well as the application of more advanced hybrid optimization techniques to the same model as in Section 2.1. 3 It is actually more common to minimize the log of the likelihood for numerical reasons, but the problem is equivalent because the log function is monotonic.

4 Please cite this article in press as: J. Hite, et al., Localization of a radioactive source in an urban environment using Bayesian Metropolis methods, Nuclear Inst. and Methods in Physics Research, A (2018), https://doi.org/10.1016/j.nima.2018.09.032.

J. Hite et al.

Nuclear Inst. and Methods in Physics Research, A (

)



Fig. 3. Detector and source placements for each experiment (satellite imagery courtesy of Google Maps).

Fig. 4. Illustration of the detectors used for the experiment.

Each detector was housed in a simple wheeled cart (Fig. 4b) so that it could easily be moved into place between measurements. Every cart was equipped with a GPS antenna (the white disc in the bottom right of Fig. 4a) and sensors to record environmental data including temperature, pressure, and humidity, as well as the orientation of the cart. The carts were covered with a Mylar sheet to reflect sunlight and minimize heating of the electronics. Calibrations for every detector were taken at the beginning and end of each experiment using Eu-252, Co-60, Cs-1374 and Ba-133 laboratory test sources, as well as a camera lens containing a significant amount of thorium. These were used to generate energy calibrations for each detector, as well as to provide reference points for the automatic gain stabilization algorithm used by the detectors. Examination of these calibrations after the conclusion of the measurement campaign indicated that they experienced negligible calibration drift over the course of the experiments. When modeling the detectors we used an intrinsic efficiency of 29%, determined from the calibration spectra taken at the beginning and end of each experiment. Since the detectors were all placed facing north, we modeled all detectors as having an face area of 224 cm2 , corresponding to the average face area of the 2′′ × 4′′ × 16′′ crystal in the detector (see Section 2.2). Dwell times were taken to be the live-time recorded by each detector and were all approximately equal to the nominal real-time for each measurement.

refers to the position of detector number 2 (of 6) during the first set of measurements (set A) of experiment 1. 𝑥 and 𝑦 coordinates are reported in meters, relative to 35◦ 56′ 12.1′′ 𝑁 84◦ 16′ 36.8′′ 𝑊 (WGS 84 PseudoMercator/EPSG:3857 coordinate reference), with the positive 𝑥 and 𝑦 directions being east and north, respectively. Orientations are reported as the angle between the longest axis of the detector and magnetic north, with positive angles in the counter-clockwise direction. Table 2 also provides a qualitative description of how obscured the source was from the position of each detector. These are organized into 4 categories, with ‘‘direct’’ indicating detectors that had an unobstructed view of the source, while ‘‘highly occluded’’ indicates detectors with several mean free paths of material between themselves and the source. Detectors designated as ‘‘occluded’’ had approximately 1–2 mean free paths of material between themselves and the source, while ‘‘lightly occluded’’ refers to detectors with lines of sight interrupted by less than 1 mean free path of material. 3.1. Detectors The detectors used for the experiment were provided by ORNL and were comprised of a 2′′ × 4′′ × 16′′ NaI scintillator coupled to an integrated photomultiplier tube and Canberra Osprey multi channel analyzer (MCA) (Fig. 4a). The output of the MCA was fed to an on board computer, which reported the measured spectra back to a central base station over a WiFi link. This wireless link also provided a centralized command and control system, allowing the detectors to be synchronized and monitored from a single location.

4 The Cs-137 calibration source was not the same as the target source used in the experiment.

5 Please cite this article in press as: J. Hite, et al., Localization of a radioactive source in an urban environment using Bayesian Metropolis methods, Nuclear Inst. and Methods in Physics Research, A (2018), https://doi.org/10.1016/j.nima.2018.09.032.

J. Hite et al.

Nuclear Inst. and Methods in Physics Research, A (

)



Table 2 Summary of actual source and detector placement during the experiment, including a qualitative description of the visibility of the source for each detector. Name

𝑥 (m)

𝑦 (m)

Dist. to source (m)

Orientation (deg. N)

Source visibility

Source 1A-1 1A-2 1A-3 1A-4 1A-5 1A-6 1B-1 1B-2 1B-3 1B-4 1B-5 1B-6

114.9 7.84 26.61 62.47 60.08 87.34 130.68 84.50 101.08 149.72 133.66 164.86 163.09

88.9 70.71 36.56 88.61 88.61 63.13 65.92 74.99 112.38 113.99 89.00 60.81 17.05

108.61 102.67 52.46 69.20 37.79 27.87 33.49 27.22 42.82 20.32 57.27 86.49

61 −15 2 18 −2 0 33 19 15 84 21 0

Highly occluded Highly occluded Occluded Occluded Direct Lightly occluded Highly occluded Lightly occluded Lightly occluded Direct Direct Highly occluded

Source 2A-1 2A-2 2A-3 2A-4 2A-5 2A-6 2B-1 2B-2 2B-3 2B-4 2B-5 2B-6

164.1 190.23 206.51 193.30 142.84 153.99 111.22 −9.46 53.72 81.98 142.28 215.31 160.91

100.8 105.89 79.54 65.78 92.36 119.94 119.34 57.86 95.20 63.52 125.97 102.44 66.12

26.65 47.51 45.68 23.56 21.54 55.94 178.83 110.45 90.16 33.19 51.26 34.59

−68 −65 −27 −16 −137 −84 −57 −103 −80 −29 −97 −63

Highly occluded Highly occluded Occluded Lightly occluded Lightly occluded Lightly occluded Highly occluded Lightly occluded Lightly occluded Direct Highly occluded Highly occluded

Fig. 5. Model geometry of the site, colored by the calculated mean free path. The DGPS reference points are also indicated by the black marks. Note that there are some small disagreements with the objects visible in the satellite imagery, this is due to the satellite photos being taken on a different day than that of the experiments.

Fig. 6. Elevation difference in experiment 2.

Material macroscopic cross-sections were assigned on a ‘‘best estimate’’ basis. Based on a photographic survey and measurements derived from the DGPS data we computed dimensions for each object. Next, we consulted a standard handbook to determine material compositions [21] and then used the NIST XCOM database to compute the nominal microscopic cross-sections5 [22]. Finally, the microscopic cross-sections were homogenized over the volume of the object to compute macroscopic cross-sections, accounting for the estimated internal dimensions and voids. These macroscopic cross-sections are visualized in terms of the corresponding mean free paths in Fig. 5, with the color corresponding

3.2. Site characterization In order to use the technique described in Section 2, it was necessary to construct a computer model for the geometry of the scene. This was done by combining differential GPS (DGPS) measurements and georeferencing them with satellite photography of the site. The DGPS system we employed is notionally capable of sub-centimeter accuracy under ideal circumstances; however, for most measurements reported accuracies of 3 to 10 cm were typical. We also took manual measurements of several objects in the field using a tape measure to validate the dimensions computed from the DGPS points. The geometry and DGPS reference points are shown in Fig. 5. When constructing the model, we decided to omit several objects in the scene. This was done considering that the detector response model of Section 2.1 only computes the uncollided flux arriving at the detectors, hence objects that do not lie between a source and a detector have no effect on the model’s estimate of the count rate.

5 Microscopic cross-sections were taken to be the total interaction crosssection at 662 keV as we are only concerned with computing the uncollided flux arriving at the detector. This violates our assumption that we know nothing about the source, however, the selection of this value is not critical since we will be treating them as highly uncertain in our later analysis. Furthermore, we analyzed the total counts acquired by each detector, not the counts under the 662 keV photopeak.

6 Please cite this article in press as: J. Hite, et al., Localization of a radioactive source in an urban environment using Bayesian Metropolis methods, Nuclear Inst. and Methods in Physics Research, A (2018), https://doi.org/10.1016/j.nima.2018.09.032.

J. Hite et al.

Nuclear Inst. and Methods in Physics Research, A (

)



Fig. 7. Foreground and background count rates in both experiments.

to the mean free path in that object, categorized into four groups representing different qualitative levels of attenuation (e.g. black and red are strong attenuators, while yellow and blue are less so). We emphasize that accurate selection of these values is not critical since they only serve as a reference point for assigning the cross section uncertainties. The localization algorithm randomly varies the macroscopic cross section values over a wide range during the sampling process and does not rely on the reference values directly. The model geometry is strictly two-dimensional, hence it does not account for changes in elevation across the scene. The site where the experiment was performed is approximately flat, with the exception of a downward drop-off in the grassy area located in the northeast section of the site. During the second experiment, one of the detectors

was positioned on this hill, resulting in an elevation change of 1 to 2 m relative to the source position (this placement is shown in Fig. 6). This was a deliberate choice, intended to introduce some of the variability expected in a real-world deployment. 4. Results and analysis Table 3 summarizes the measurements recorded during the experiments, while Fig. 7 plots the net foreground and background count rates recorded during both experiments. Total counts refer to the total counts recorded during measurements with the source present, while background counts were accumulated from measurements with the source absent, and net counts are calculated from measurements with 7

Please cite this article in press as: J. Hite, et al., Localization of a radioactive source in an urban environment using Bayesian Metropolis methods, Nuclear Inst. and Methods in Physics Research, A (2018), https://doi.org/10.1016/j.nima.2018.09.032.

J. Hite et al.

Nuclear Inst. and Methods in Physics Research, A (

)



Fig. 8. Marginal posterior density for source location 𝐫 in experiment 1 with 50% relative uncertainty in all cross-sections. Uncertainty: 𝜎𝑥 = 2.01 m and 𝜎𝑦 = 1.88 m.

the source present subtracting the background counts after normalization for live time. Table 3 also shows the ratio of total counts to background, which was approximately 3 to 13 for detectors near the source, while detectors on the perimeter had total to background ratios between 1 and 2. As can be seen in Fig. 7, detectors with a direct lineof-sight to the source recorded the highest number of counts and hence contribute the most to the localization, while the Bayesian formulation guarantees robustness against detectors with lower count rates since their contributions to the localization are weighted by their uncertainty. Figs. 8 and 9 show the marginal posterior densities for the source position from experiments 1 and 2 using the detector response model that does not account for detector orientation (see Section 2.2). These posteriors were generated by drawing 105 samples,6 with convergence being determined by examining both the chain autocorrelation and Gelman–Rubin scores [23]. The measurement data was taken to be the net counts recorded by each detector and did not exploit any spectral features as it would require a priori knowledge of the source to identify the relevant features of interest. Background was treated as a Poisson random field, with mean values at each detector location taken from the measurements made without the source present. All cross-section values were assigned a uniformly distributed uncertainty within ±50% of their respective nominal values. In Figs. 8a and 9a we see that the method is able to identify the source location with reasonable accuracy relative to the scale of the problem. On further inspection in Figs. 8b and 9b we can see that there is a distinct bias in the estimated source position, one which exceeds the posterior uncertainty that arises from the ad hoc estimation of the cross-sections and measurement uncertainties. These biases, which are on the order of 5 to 10 m, are small enough that the estimates are sufficiently accurate to guide a more targeted search, which is our primary objective. However we shall proceed to examine them in more detail in order to account for detector orientation.

vehicle was located in close proximity to several of the detectors as well as the source, hence its movement produced significant changes in the count rates that were not included in the detector response model. Furthermore, the path traveled by the vehicle crossed between two of the detectors and the source, resulting in a large drop in the count rates measured when the vehicle interrupted their view of the source. Fig. 10 shows time series plots of the count rates measured by the three detectors that were affected by the vehicle’s movement. Fig. 10a shows the detector that was situated near the initial position of the vehicle, with a large drop in the count rate visible when the vehicle was passing between the source and detector, followed by an increase in the average count rate after the vehicle had passed due to reduced attenuation. Fig. 10b shows the count rate recorded by the detector located near the final position of the vehicle, with a sharp drop in the count rate as the vehicle moved between the source and the detector, followed by an overall decrease in the count rate due to increased attenuation in the presence of the vehicle. Finally, Fig. 10c shows the detector that was located on the opposite side of the parking lot. Here the change in count rate is smaller than the other two detectors, nevertheless there is a clear trend visible in the recorded count rates. All three cases show statistically significant changes in the observed count rates when compared to the average count rate over the full measurement period and as such we consider them outliers. Preliminary analysis suggested that these count rate fluctuations had an impact on the accuracy of the source localization algorithm applied to experiment 2. Our current implementation of the detector response model does not allow for time-dependent changes in the geometry; instead we elected to discard the data recorded by these detectors. This naturally results in higher uncertainty due to the decreased number of counts recorded by the overall detector network, however the original experimental design included substantial redundancy in the number of detectors in case of just such an abnormality. For the remainder of this article, all results shown for experiment 2 exclude these measurements. Note that we also examined the results when these measurements were excluded, but without incorporating a correction for the detector orientation. This case showed modest improvement compared to those in Fig. 9, with the error improving to approximately 7 m.

4.1. Count rate anomalies during measurements By examining the footage taken from cameras situated throughout the experiment area, we discovered that a large vehicle in the eastern parking lot moved during the measurements for experiment 2. The 6

4.2. Detector orientation As discussed in Section 2.2, the geometric efficiency of the detectors used for the measurements varies significantly with the orientation of

All subsequent posteriors shown were also constructed using 105 samples. 8

Please cite this article in press as: J. Hite, et al., Localization of a radioactive source in an urban environment using Bayesian Metropolis methods, Nuclear Inst. and Methods in Physics Research, A (2018), https://doi.org/10.1016/j.nima.2018.09.032.

J. Hite et al.

Nuclear Inst. and Methods in Physics Research, A (

)



Table 3 Count rates recorded by the detectors. Name

Total live time (s)

Total counts

Background live time (s)

Background counts

Total/Bg

SNRa

1A-1 1A-2 1A-3 1A-4 1A-5 1A-6 1B-1 1B-2 1B-3 1B-4 1B-5 1B-6

1789.13 1789.52 1785.44 1789.72 1770.83 1756.45 1787.34 1764.57 1783.96 1710.37 1781.88 1763.21

2,008,213 1,836,856 2,656,026 1,772,977 5,497,133 7,310,350 2,407,008 6,513,106 3,053,585 15,719,863 3,474,658 2,852,067

1716.42 1790.66 1789.15 1746.49 1790.07 1739.25 1791.92 1789.12 1790.91 1792.64 1790.64 1784.61

1,702,957 1,615,817 1,815,188 1,373,245 1,841,017 1,221,574 1,430,501 1,954,174 1,632,757 1,272,534 1,699,812 2,570,414

1.13 1.14 1.47 1.26 3.02 5.93 1.69 3.38 1.88 12.95 2.05 1.12

175.0 174.8 627.5 308.3 2723.8 5471.1 820.6 3303.2 1119.1 13164.6 1371.1 196.1

2A-1b 2A-2 2A-3 2A-4b 2A-5 2A-6 2B-1 2B-2 2B-3 2B-4 2B-5 2B-6b

1782.81 1786.82 1776.18 1722.84 1738.96 1776.83 1790.52 1789.14 1788.92 1759.72 1789.42 1785.37

3,347,193 2,385,999 4,276,542 13,224,755 11,710,996 3,932,848 1,736,153 1,945,903 1,967,095 7,098,048 2,050,384 2,594,652

1793.62 1789.80 1790.73 1791.80 1788.52 1788.74 1791.23 1791.20 1790.03 1788.51 1794.64 1789.89

1,150,058 1,674,788 1,398,868 1,357,271 2,055,420 1,814,991 1,620,411 1,501,151 1,641,784 2,005,485 983,992 1,702,360

2.93 1.43 3.08 10.1 5.86 2.18 1.07 1.30 1.20 3.60 2.09 1.53

2061.5 552.2 2452.7 10434.1 6870.4 1586.3 91.4 364.6 254.8 3648.3 1079.5 688.0

a

The definition for signal-to-noise ratio in the context of radiation detector measurements is somewhat ambiguous. In Table 3, we define the SNR as √ 𝐶Net ∕ 𝐶Bg ∗ 𝓁Tot ∕𝓁𝐵𝑔 , where 𝐶 denotes counts and 𝓁 the live time. That is, the ratio of the net counts to the standard deviation of the background noise, adjusted for differences in live time between the foreground and background measurements. b Marked entries denote measurements with anomalous count rates (see Section 4.1). Total-to-background ratios and SNRs are adjusted for live time.

Fig. 9. Marginal posterior density for source location 𝐫 in experiment 2 with 50% relative uncertainty in all cross-sections. Uncertainty: 𝜎𝑥 = 5.34 m and 𝜎𝑦 = 3.17 m.

the detector relative to the source. Figs. 11 and 12 show the posterior densities obtained from the MCMC algorithm using the same parameters as the plots in Figs. 8 and 9, but with the correction for detector orientation described in Section 2.2 applied. Observe that the error in the source location is reduced to 2 m to 3 m when the effect of detector orientation is included.

on the desired dwell time. We used this technique to run the MCMC sampler separately for a series of simulated experiments with increasing counting time to examine its effect on the bias in the estimated source location. Fig. 13 plots the Cartesian distance between the mode of the posterior distribution (i.e., the most probable source location) and the true source location versus the count time. There is a clear trend of decreasing bias as the count time is increased, with a count time of 5 min resulting in a decrease in the bias to approximately 4 m. On the other hand count times beyond 5 min only show modest improvements in accuracy of approximately 1 to 2 m, as well as little change in uncertainty. These observations together suggest that a significant portion of the

4.3. Effect of counting time on posterior error All measurements were taken for 30 min live-time, but shorter measurement times can be simulated by assuming the counts recorded follow a Poisson distribution and resampling measurement data based 9

Please cite this article in press as: J. Hite, et al., Localization of a radioactive source in an urban environment using Bayesian Metropolis methods, Nuclear Inst. and Methods in Physics Research, A (2018), https://doi.org/10.1016/j.nima.2018.09.032.

J. Hite et al.

Nuclear Inst. and Methods in Physics Research, A (

)



Fig. 10. Time series of fluctuating count rates recorded by three of the detectors in experiment 2 (blue) versus the average count rate over the entire measurement (green).

Fig. 11. Marginal posterior density for source location 𝐫 in experiment 1 with 50% relative uncertainty in all cross-sections and including detector orientation. Uncertainty: 𝜎𝑥 = 2.78 m and 𝜎𝑦 = 2.45 m.

10 Please cite this article in press as: J. Hite, et al., Localization of a radioactive source in an urban environment using Bayesian Metropolis methods, Nuclear Inst. and Methods in Physics Research, A (2018), https://doi.org/10.1016/j.nima.2018.09.032.

J. Hite et al.

Nuclear Inst. and Methods in Physics Research, A (

)



Fig. 12. Marginal posterior density for source location 𝐫 in experiment 2 with 50% relative uncertainty in all cross-sections and including detector orientation.

causing a deviation in the location predicted by the posterior from the true source location relative to the position of objects in the scene. For an overview of the georeferencing process see Ref. [24]. The second potential source of error is in the reconstruction of the geometry and is also related to the georeferencing of the scene. The georeferencing error is an epistemic effect from interpolation of the reference points and in particular affects the distances computed between locations, however, there is also an aleatory uncertainty arising from imperfect knowledge of the dimension of objects in the scene. The exterior dimensions of objects were determined via a combination of the DGPS points recorded and manual measurements and have an associated uncertainty. Imperfect measurements may have yielded object dimensions that were too large or too small, resulting in errors in the estimated source position. The exact location of vehicles in the Eastern parking lot is especially relevant, as it was difficult to exactly determine their positions and we had to rely on best estimates. Lastly, a portion of the bias almost certainly resulted from physical effects that are not accounted for in the detector response model. The model only computes the uncollided flux arriving at the detector, whereas we used net counts recorded as the measurement data. The net count rate obviously includes counts resulting from gammas that have undergone scattering and the resulting discrepancy would cause bias if its effect violated the assumption of independently distributed error in the statistical model.

Fig. 13. Error in estimated source location versus count time for experiment 1.

bias observed is dependent on the counting time, but also that other independent sources of bias are present.

4.5. Trilateration

4.4. Other sources of error

For comparison, Tables 4 and 5 shows the results of applying trilateration to the data collected during each experiment using the least-squares formulation described in Ref. [5]. This is applied without accounting for the presence of objects in the scene, equivalent to the model described in Sections 2.1 and 2.2 without the exponential term included. The resulting least-squares problems were then solved using an implementation of the Levenberg–Marquardt algorithm, with error estimates computed via the standard method of linearizing the error surface. Table 4 shows the predicted source location without considering detector orientation, while the results in Table 5 are calculated using geometric efficiencies that have been adjusted for orientation. Run times for these calculations are negligible (less than a second), but the resulting predictions of the source location are poorer than the results using MCMC.

We further identify several possible candidates as the cause(s) of the remaining bias. The first is in the georeferencing of the satellite photography with the DGPS control points taken. This process proceeds by referencing different locations in the overhead imagery with known coordinates, which were recorded using the DGPS system described in Section 3.2. An interpolating function is generated using these reference locations to associate geographic coordinates with every point in the image, followed by a transformation to align the image with the new coordinate system. Since this transformation is generally nonlinear (and non-affine) it can introduce non-uniform errors, particularly when computing real-world distances between different locations in the image. We hypothesize that small errors in the relative positions of objects in the scene may have an effect on the estimated source location, 11

Please cite this article in press as: J. Hite, et al., Localization of a radioactive source in an urban environment using Bayesian Metropolis methods, Nuclear Inst. and Methods in Physics Research, A (2018), https://doi.org/10.1016/j.nima.2018.09.032.

J. Hite et al.

Nuclear Inst. and Methods in Physics Research, A (

Table 4 Predicted source locations using simple trilateration with error estimates without accounting for detector orientation. Experiment

𝑥 (m)

𝑦 (m)

𝜎𝑥 (m)

𝜎𝑦 (m)

Residual norm (m)

1 2

148.9 176.0

83.5 76.3

8.3 9.1

1.8 21.0

34.4 13.3

𝑥 (m)

𝑦 (m)

𝜎𝑥 (m)

𝜎𝑦 (m)

Residual norm (m)

1 2

110.3 145.7

76.3 119.4

9.1 4.6

21 6.0

13.3 26.2



Acknowledgments This work was funded by the NNSA Office of Defense Nuclear Nonproliferation R&D through the Consortium for Nonproliferation Enabling Capabilities (CNEC) and the experiments described in this article were performed in cooperation with Oak Ridge National Laboratory.

Table 5 Predicted source locations using simple trilateration with error estimates including detector orientation. All distances are in meters. Experiment

)

References [1] J. Hite, J. Mattingly, K. Schmidt, R. Ştefănescu, R. Smith, Bayesian Metropolis techniques for source localization (Poster), in: IEEE Symposium on Radiation Measurements and Applications, Berkeley, CA, 2016. [2] J. Hite, J. Mattingly, K. Schmidt, R. Ştefănescu, R. Smith, Bayesian metropolis methods applied to sensor networks for radiation source localization, in: IEEE International Conference on Multisensor Fusion and Integration for Intelligent Systems (MFI), IEEE, 2016, pp. 389–393. [3] N. Rao, M. Shankar, J. Chin, D. Yau, S. Srivathsan, S. Iyengar, Y. Yang, J. Hou, Identification of low-level point radiation sources using a sensor network, Int. Conf. Inf. Proc. Sensor Netw. 7 (3) (2008) 1–35. [4] M. Morelande, B. Ristic, Radiological source detection and localisation using Bayesian techniques, IEEE Trans. Signal Process. 57 (11) (2009) 4220–4231. [5] Y. Liu, Z. Yang, X. Wang, L. Jian, Location, localization, and localizability, J. Comput. Sci. Tech. 25 (2) (2010) 274–297. [6] R. Stolsmark, E. Tossebro, Uncertainty in trilateration, in: Proceedings of the 1st International Conference on Sensor Networks, 2012, pp. 237–241. [7] M. Morelande, B. Ristic, A. Gunatilaka, Detection and parameter estimation of multiple radioactive sources, in: FUSION 2007 - 10th International Conference on Information Fusion, 2007. [8] Z. Liu, Mobile radiation sensor networks for source detection in a fluctuating background using geo-tagged count rate data (M.s., ), University of Illinois at Urbana-Champaign, 2016. [9] J. Howse, L. Ticknor, K. Muske, Least squares estimation techniques for position tracking of radioactive sources, Chem. Eng. 37 (2001) 1–15. [10] N. Rao, S. Sen, N. Prins, D. Cooper, R. Ledoux, J. Costales, K. Kamieniecki, S. Korbly, J. Thompson, J. Batcheler, R. Brooks, C. Wu, Network algorithms for detection of radiation sources, Nucl. Instrum. Methods Phys. Res., Sect. A 784 (2015) 326–331. [11] R. Smith, Uncertainty Quantification: Theory, Implementation, and Applications, Society for Industrial and Applied Mathematics, Philadelphia, 2013. [12] W. Hastings, Monte Carlo sampling methods using Markov chains and their applications, 1970. [13] J. Hite, J. Mattingly, Bayesian Metropolis methods for source localization in an urban environment, Radiat. Phys. Chem. Press and Corrected Proofs (2018) URL https://doi-org.prox.lib.ncsu.edu/10.1016/j.radphyschem.2018.06.024. [14] P. Jacob, C. Robert, M. Smith, Using parallel computation to improve independent Metropolis–Hastings based estimation, J. Comput. Graph. Statist. 20 (3) (2011) 616–635. [15] J. Favorite, K. Bledsoe, D. Ketcheson, Surface and volume integrals of uncollided adjoint fluxes and forward-adjoint flux products, Nucl. Sci. Eng. 163 (2009) 73–84. [16] K. Schmidt, Uncertainty Quantification for Mixed-effects Models with Applications in Nuclear Engineering (Ph.d dissertation), North Carolina State University, 2016. [17] Z. Liu, C. Sullivan, Mobile radiation sensor networks for source detection in a fluctuating background using geo-tagged count rate data, Trans. Nucl. Sci. (2016). [18] A. Van Oosterom, J. Strackee, The solida angle of a plane triangle, IEEE Trans. Biomed. Eng. BME-30 (2) (1983) 125–126. [19] R. Chudley, Building Construction Handbook, Routledge, London New York, 2014. [20] R. Ştefănescu, K. Schmidt, J. Hite, R. Smith, J. Mattingly, Hybrid optimization and Bayesian inference techniques for a non-smooth radiation detection problem, Internat. J. Numer. Methods Engrg. (2017) 1–36. [21] R. McConn, C. Gesh, R. Pagh, R. Rucker, R. Williams III, Compendium of material composition data for radiation transport modeling, Tech. rep.,, Pacific Northwest National Laboratory (PNNL), Richland, WA (US), 2011. [22] M. Berger, J. Hubbell, S. Seltzer, J. Chang, J. Coursey, R. Sukumar, D. Zucker, K. Olsen, XCOM: photon cross sections database (NIST Standard Reference Database 8), Online database, National Institute of Standards and Technology, 2010. [23] A. Gelman, D. Rubin, lnference from iterative simulation using multiple sequences, Statist. Sci. 7 (4) (1992) 457–472. [24] A. Hackeloeer, K. Klasing, J. Krisp, L. Meng, Georeferencing: a review of methods and applications, Ann. GIS 20 (1) (2014) 61–69.

5. Conclusions In this article, we presented a method for estimating the location of an unknown source of radiation in an urban environment using measurements from a network of distributed radiation detectors. By employing a Bayesian framework based on a specialized Markov-chain Monte Carlo sampler, we calibrate a simplified statistical model for the response of the detector network based on the count rates measured in the field. This results in a posterior probability density for the estimated source location which incorporates all measurements, representing the best estimate of the source position while also describing the uncertainty associated with counting statistics and the composition of objects in the scene. To test this method we performed a pair of experiments mimicking a real-world search scenario at ORNL in the area surrounding the Energy Systems Test Complex. This area features several buildings and large objects as well as open spaces, providing a good surrogate for a typical urban scene. We placed a 37 mCi source at two different positions and measured the detector response at 12 locations per source placement using 6 mobile detectors, each with a 2′′ × 4′′ × 16′′ NaI scintillator. We also performed a survey of the site, which we later used to reconstruct the geometry of the scene for use with the detector response model. The results of the MCMC sampler for both experiments produce estimated source locations that are sufficiently close to guide a targeted follow-up search, however, they consistently show a bias of 5 to 10 m for a count time of 30 min in both experiments. By extending our model for the detector response to account for the orientation of the detectors as well as excluding some anomalous count rates recorded in the second experiment we are able to reduce the bias to be on the order of 2 m. The demonstration of the importance of accounting for detector orientation serves to highlight the utility of the Bayesian Metropolis localization methodology. Simple intuition would likely lead many readers to conclude that the localization errors in Figs. 8 and 9 arise from the large uncertainties in the material cross sections. Access to the full posterior densities allowed us to recognize that the localization errors exceeded what is expected from cross section uncertainties and subsequently to recognize the importance of accounting for detector orientation, as well as to identify count rate anomalies caused by a vehicle moving through the scene in experiment 2.

12 Please cite this article in press as: J. Hite, et al., Localization of a radioactive source in an urban environment using Bayesian Metropolis methods, Nuclear Inst. and Methods in Physics Research, A (2018), https://doi.org/10.1016/j.nima.2018.09.032.