Research Article | Open Access
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. https://doi.org/10.5402/2011/189735
Comparison of Two Approaches for Detection and Estimation of Radioactive Sources
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.
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 ). 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.  applied a nonlinear, recursive least-squares algorithm for tracking the positions of a moving radioactive point source. Stephens Jr. and Peurrung  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.  and Brennan et al.  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.  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.  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: 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 , the th point source of gamma radiation (where ) is parameterised by (i)its location in a Cartesian coordinate system with , where is some large region that is assumed to contain all the sources, (ii)its intensity rate (source strength), .
The parameter vector of source is thus given by where denotes the matrix transpose. The source parameter vectors are collected into a stacked vector , 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: 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: (). In (2.4), is the binomial rate, and, and are, respectively, the minimum and maximum number of sources.
Two alternative functional forms are used for the prior distribution of (). The first of the two alternatives assigns the prior for the source strength to be a gamma distribution: 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: where is the probability that the source is active (), is the Dirac delta function, and is an a priori upper bound on the expected source intensity rate. In (2.6), denotes the indicator (characteristic) function for set , with if and if .
Finally, the prior distribution for the source location is assigned to be uniformly distributed over the region that is assumed a priori to contain the source: where 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 ), implying that the probability that a gamma radiation detector registers counts ( 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 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 () be a measured count from the GM counter, taken at (receptor) location . The following assumptions are made: (1) the GM counter has a uniform directional response; (2) the radiation sources are point sources; (3) there is negligible attenuation of gamma radiation due to air; (4) the radiation measurements are independently distributed, and, (5) the exposure time (or the sampling interval) is constant for all measurements.
2.2.1. Measurement Model: Known Background Rate
The likelihood function is simply the joint density of the measurement vector 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 ; Tsoulfanidis ) as follows: where is the mean radiation count for the th sensor location: with 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 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 (). The mean signal count rate can be described by the radiation data and by the quantitative model (cf. (2.12)), so 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), is the uncertainty in the estimate of the mean signal count rate, and 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 where is the composite error. If the expectation value of the composite error is zero and the variance is assumed to be given by , then application of the principle of maximum entropy (see, Jaynes ) provides the following explicit form for the likelihood function (assuming that at least one radiological source is present, so ): 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 where is the variance of the uncertainty in the estimate of the mean signal count rate (obtained from the data ) and is the variance of the model error . In this paper, the model error standard deviation (or square root of the variance) is specified as .
To complete the specification for the likelihood function in (2.15), we need to provide the estimate along with the variance in this estimate. To this purpose, we assume that we have available a measurement of background radiation (obtained in the absence of any point radiation sources) given to us as counts obtained over an exposure time . Now, suppose that a radiation sensor at a fixed receptor location obtains a measurement of counts in an exposure time . We will now show how to use this information to obtain an estimate for the mean signal count rate , as well as to characterize the uncertainty in this estimate in terms of the standard deviation .
To accomplish this, we utilize a Bayesian solution for the extraction of weak signals in strong background interference described by Loredo . 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 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 ) which has the following form (conditioned on a known mean background rate ): 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 ): since does not have any bearing on the determination of . The prior is assigned the nonsinformative Jeffreys prior so and the likelihood function , so
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: 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 The marginal probability distribution for the mean signal count rate can be derived by marginalizing with respect to the (unknown) mean background count rate: which on using (2.22) and performing the integration finally gives (on normalization of the probability density function) where 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: The second moment of about zero can be evaluated as from which the uncertainty in the estimate for given in (2.26) can be determined as 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 ). 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 , 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 ): which on insertion of the normalized form of (2.20) for and evaluation of the integral gives 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 and the likelihood function of (2.30) valid for can be combined to give the following general likelihood function for the measurement model assumed in this section: 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): The minimum mean-squared error (MMSE) estimate of is then the posterior expectation,
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: where is a “particle” set consisting of samples drawn from the importance density and is the set of importance weights. The sum of these weights is equal to one (namely, ).
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 , denote the posterior PDF at stage . A series of can be constructed by setting where the intermediate likelihood at stage is with , such that , for and . Assume that a random sample from 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 , the sample is drawn directly from the prior with .
The performance of the procedure depends somewhat on the number of steps and the expansion factors . The expansion factors were selected in line 8 of Algorithm 1 using an adaptive scheme proposed in Musso et al. . 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 to give samples that have uniform weights. This resampling is undertaken using the systematic resampling algorithm (see, Kitagawa ), which is summarized in the form of a pseudocode given in Table 3.2 of Ristic et al. . 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. ).
The case 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 (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 (, Chapter 7). Parameter is a tuning parameter, which needs to be determined by field trials and calibration.
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  page 225]); Rissanen ): where is the maximum likelihood estimate of assuming sources are present and is the dimension of under the hypothesis . Here, denotes a set of candidate source numbers up to some maximum . While MDL is derived based on the maximum likelihood estimate (see (3.6)), we will use the MDL in conjunction with the PC by simply replacing with the PC estimate in (3.6), so Finally, the estimate of the number of sources is then determined as follows:
To estimate the number of sources using the MDL criterion, we first run the importance sampling using progressive correction for all values of . 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 (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 when .
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 , ). 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 . 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 . To this purpose, let () 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 radiological point source and, in so doing, changes the parameter space dimension by .
In the first category of moves (which are dimension-conserving), the parameter vector ( fixed) is partitioned as follows: where and . The parameters collected together in are linearly related to the radiation “measurements” (cf. (2.12)). For these parameters, we apply a Gibbs sampler for the update (propagation) move involving drawing samples directly from the univariate full conditional distribution for (which happens to be a truncated Bernoulli-Gaussian distribution in this case). The remaining parameters in are related nonlinearly to the radiological measurements (cf. (2.12)). Owing to the fact that the univariate conditional posterior distribution for these parameters cannot be determined analytically, we apply a Metropolis-Hastings (M-H) sampler to update the parameters in . Towards this purpose, we update () by generating a new value from a proposal transition kernel that is chosen here to be a Gaussian mixture, with each component of the mixture having a mean and different variances . 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 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 , the probabilities for the birth () and death () moves are specified as follows: ensuring that the probability of a jump move (either birth or death) lies in at each iteration. We note that for , and for , .
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 and which define the prior distribution , as well as a value for (maximum number of iteration steps to take). (2)Set the iteration counter and choose an initial state for the Markov chain by sampling from . (3)Starting from , conduct the following sequence of moves to update the state vector to : where and denote some intermediate transition states between iterations and . (4)Change the counter from to , and return to step 3 until a maximum number of steps () has been taken.
In Step 3, denotes a jump move involving either the birth (addition) of a radiological point source at a random locations or the death (removal) of an existing radiological point source with probabilities and , respectively; denotes the update of the source intensity rates using a Gibbs sampler; and, denotes the update of the source locations using an M-H sampler.
3.2.1. Simulated Annealing
To improve convergence of the Markov chain, we have implemented a simulated annealing procedure. This procedure is very similar to the progressive correction that is used in conjunction with the importance sampling (as described in Section 3.1). In a nutshell, simulated annealing mimics the quasistationary evolution of a system through a sequence of equilibrium samplings, involving an infinitely slow switching between the initial and final states. More concretely, let us consider an ensemble of different source distributions (system) that has been randomly drawn from a modified posterior having the following form: , where is a “coolness” parameter. These samples will be labelled (). When (initial state), the likelihood function is switched off and the modified posterior distribution reduces exactly to the prior distribution. On the other hand, when (final state) the modified posterior distribution is exactly the posterior that we wish to sample from. In between these two extremes, with , 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 (), so implies that . The modified posterior () corresponds to “heating” the posterior distribution to a temperature . 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 (infinite temperature) by randomly drawing samples of source distributions () from the prior distribution . Given an ensemble of samples that has achieved equilibrium (at temperature ) with respect to the modified posterior , an ensemble of samples that is consistent with (at the reduced temperature , ) can be obtained by using the weighted resampling method (see, Gamerman and Lopes ) applied to (). To this purpose, each sample is assigned a normalized importance weight as follows: for . 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 with weights . This resampling involves mapping the random discrete measure defined by into the new random discrete measure defined by having uniform weights (namely, 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 , but in this paper we use the systematic resampling scheme described by Kitagawa  (namely, the same resampling scheme that is used in the PC algorithm described in Section 3.1). This scheme requires time to execute and minimizes the Monte Carlo variation in the resample. For completeness, we briefly describe this resampling scheme as applied to the current problem. A vector of copies of the samples () is obtained by computing a vector of the cumulative sums of , generating a random draw from the uniform distribution , and determining ( from where denotes the “integer part of”. After resampling, the weights for each member of the resample are reset to (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 (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 from to . This sequence defines the annealing schedule for . 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 where the stochastic sampling scheme is exploring the prior) the step size is chosen so that where are the weights computed in accordance to (3.12) and is a small constant (). 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 () with for and , that define the annealing schedule need not sum to (namely, in general).
When , the annealing phase is complete (corresponding, as such, to the burn-in phase of the algorithm) and probabilistic exploration of the hypothesis space proceeds using the RJMCMC algorithm described at the end of Section 3.2. In the application of this algorithm following simulated annealing, Step 2 is modified as follows: instead of sampling from the prior distribution , these initial states are simply chosen to coincide with the samples from the ensemble of source distributions 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 without saturating , and it has a fairly flat response () from 0.1 to above 1 MeV (see ; page 7). Measured dose rate data were recorded in Sv . We converted these dose rate data into raw count measurements by multiplication by a conversion factor (which was 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 (Cs) and one cobalt (Co) 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 , 392.28, and 98.07 Sv , for Sources 1, 2, and 3, respectively. We will use these values as the “true” source strengths in our source estimations.
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 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.
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 counts .
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: 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 and a scale parameter Sv . 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: , where denotes Cartesian product.
The PC algorithm is initialized by drawing a random sample in -dimensional parameter space from . The number of samples in the PC algorithm is set to .
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 :
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: , , , , Sv m2 , and 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 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 . After termination of the annealing at , a further iterations of the RJMCMC procedure were applied to each of the 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 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: where the (marginal) posterior probability for the number of sources is estimated from the samples as follows: Here, is the number of samples drawn from a sampled realization of the Markov chain. Similarly, using the Markov chain samples drawn from the posterior distribution, the posterior means of the various source parameters (various components of ; see (2.2)) can be estimated by 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 ).3 In this approach, can be evaluated as follows: 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: where is determined according to (3.14); and, the annealing schedule so defined provides the set of quadrature points .
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 ): 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 ).
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 and 48 s. In the application of importance sampling using PC, the mean background count rate was assumed to be exactly known ( counts was used in the algorithm). For the application of RJMCMC used in conjunction with simulated annealing, the first 24 s ( s) of a measurement of the background counts at receptor location 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 , 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 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 and 48 s), the minimum value of was found to correspond to , 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 (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 s (upper portion of Table 3) and for 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 and 48 s.
(a) Source 1
(b) Source 2
(c) Source 3
(a) Source 1
(b) Source 2
(c) Source 3
For the application with 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 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 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 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 s). The comparable result for 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 to 2 (or from to 1 and from to 0) do not occur. Furthermore, birth moves involving transitions from to 3 do not occur, either. On the other hand, dimension-changing moves (births or deaths) involving transitions from to 4 (and vice versa) are seen to occur. The residence time for occupation of a state with is seen to be significantly longer than that for a state with . 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 (), which corresponds to the correct number of sources for this test. More specifically, 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 sources (and not sources owing to the fact that 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 is now favored with a probability greater than about 0.995.
Given the fact that our best estimate for the number of sources is , 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 and 48 s.
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 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 and 48 s. Furthermore, the intensity rates of the sources are estimated quite well for both 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 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 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 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 s. Again, the comparable results for 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.
(a) Source 1
(b) Source 2
(c) Source 3
Finally, the information gain (see (5.8)) obtained from the radiation measurements for and 48 s (for test case 3) was found to be and 60.1 natural units (nits), respectively.4 This implies that the information contained in the radiation measurements with exposure times 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 and , 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 s provided about 7 nits more information gain than that using radiation data measured with an exposure time of s.
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 did not adequately characterize the uncertainty in the source parameter estimates for the radiation measurements obtained with an exposure time 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 ). In fact, as discussed in Jaynes , 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. ).
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.
- 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.
- 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).
- Actually, the technique of thermodynamic integration dates back to Kirkwood  who first used it within the context of the statistical mechanics of fluid mixtures to compute free energy differences.
- 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)).
- W. K. H. Panofsky, “Nuclear proliferation risks, new and old,” Issues in Science and Technology, vol. 19, no. 4, pp. 73–75, 2003.
- 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.
- 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.
- 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.
- 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.
- 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.
- 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.
- 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.
- 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.
- N. Tsoulfanidis, Measurement and Detection of Radiation, Taylor & Francis, 1995.
- A. Martin and S. A. Harbison, An Introduction to Radiation Protection, Chapman & Hall, 3rd edition, 1986.
- E. T. Jaynes, Probability Theory: The Logic of Science, Cambridge University Press, Cambridge, UK, 2003.
- 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.
- 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.
- 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.
- B. Ristic, S. Arulampalam, and N. Gordon, Beyond the Kalman Filter: Particle Filters for Tracking Applications, Artech House, 2004.
- R. Mahler, Statistical Multisource Multitarget Information Fusion, Artech House, 2007.
- S. M. Kay, Statistical Signal Processing: Detection Theory, Prentice-Hall, 1998.
- J. Rissanen, “Modelling by shortest data description,” Automatica, vol. 14, pp. 465–471, 1978.
- P. J. Green, “Reversible jump Markov chain Monte Carlo computation and Bayesian model determination,” Biometrika, vol. 82, no. 4, pp. 711–732, 1995.
- 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.
- 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.
- D. Gamerman and H. F. Lopes, Markov Chain Monte Carlo: Stochastic Simulation for Bayesian Inference, Chapman & Hall, 2nd edition, 2006.
- AN/PDR-77 User’s Guide, Canberra Industries, Meriden, Conn, USA, 2000.
- 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.
- T. M. Cover and J. A. Thomas, Elements of Information Theory, John Wiley & Sons, New York, NY, USA, 1991.
- 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.
- J. G. Kirkwood, “Statistical mechanics of fluid mixtures,” The Journal of Chemical Physics, vol. 3, no. 5, pp. 300–313, 1935.
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.