Bayesian inference of 1D activity profiles from segmented gamma scanning of a heterogeneous radioactive waste drum
BB AYESIAN INFERENCE OF ACTIVITY PROFILES FROMSEGMENTED GAMMA SCANNING OF A HETEROGENEOUSRADIOACTIVE WASTE DRUM
Eric Laloy ∗ Bart Rogiers An Bielen Sven Boden
Institute for Environment, Health and Safety, Belgian Nuclear research Centre (SCK CEN)January 7, 2021 A BSTRACT
We present a Bayesian approach to probabilistically infer vertical activity profiles within a radioactivewaste drum from segmented gamma scanning (SGS) measurements. Our approach resorts to Markovchain Monte Carlo (MCMC) sampling using the state-of-the-art Hamiltonian Monte Carlo (HMC)technique and accounts for two important sources of uncertainty: the measurement uncertaintyand the uncertainty in the source distribution within the drum. In addition, our efficiency modelsimulates the contributions of all considered segments to each count measurement. Our approachis first demonstrated with a synthetic example, after which it is used to resolve the vertical activitydistribution of 5 nuclides in a real waste package.
Segmented gamma scanning (SGS) is a longstanding but still actively researched technique to assess the activity ofradioactive waste drums in a non-destructive way [e.g., 2, 15, 20, 8, and references therein]. Classical SGS interpretationmethods either sum up all gamma spectra to derive representative activities for the whole drum or treat the measuredcounts associated with each segment individually. Efficiency calibration can be performed in a classical way usingreference point sources that are placed in a waste drum matrix in a special way to simulate source homogeneity.Nowadays, Monte Carlo (MC) modeling techniques such as implemented in the commercially available ISOCS software[ISOCS: In Situ Object Counting System, 16] are commonly used to determine detector’s efficiency. Though offeringsignificant geometrical flexibility, the ISOCS variants cannot always fully account for the spatial distribution of activitieswithin the drum. Furthermore, the uncertainty in the estimated activities by radionuclide quantification techniques isoften calculated in a relatively simple way, whether based on first-order Taylor expansion or MC error propagation,ignoring for instance drawbacks related to non-physical values, non-Gaussian distributions or non-linear expressions,or the omission of, e.g., systematic uncertainties like that of the efficiency [e.g., 14]. Bayesian inference, also calledBayesian data inversion, is well suited to assess the uncertainty of quantities that result from the application of amathematical model [see the book by 10, for a comprehensive description of Bayesian data analysis]. In our opinion,Bayesian inference is particularly attractive for radiological characterization for the following reasons.• We are not simplifying the mathematical model, as is the case for first-order Taylor expansion, and can workwith any type of distribution (this is the same for MC error propagation)• We can account for prior information and expert opinion. Also MC error propagation allows this. Thedifference is however, that if there is information in the data, the Bayesian inference will update those priorbeliefs into posterior beliefs, while in the MC error propagation, there is no way to do that. Also, for theBayesian inference, we can express our prior beliefs on the quantities of interest (e.g. activities), while for theMC error propagation there is no way to do that. ∗ Corresponding Author, [email protected] a r X i v : . [ phy s i c s . d a t a - a n ] J a n We do in principle not need to give special attention to performance characteristics like decision threshold,detection limit and minimum detectable activities (MDA). The Bayesian inference will, based on our priorbeliefs, always provide the possible range of outcomes, even for cases where a classic workflow would haveresulted in a reported MDA.• With this approach it is straightforward to deal with situations for which multiple observations inform aboutthe same quantities and/or processes. For instance, when one or more inferred radionuclide activities areassociated with multiple energy peaks in gamma spectrometry. Here for a given nuclide, Bayesian analysiswill provide a unique posterior activity distribution that is consistent with all the energy peaks. In contrast,this is something neither conventional first-order Taylor expansion nor MC uncertainty propagation can copewith. Indeed, in this case the uncertainty propagation would deliver separate activity distributions for eachpeak, with no straightforward way to combine them. Other situations of multiple observations informing aboutthe same quantity arise when different radiological measurement techniques sensing the same radionuclideactivities and/or relying on the same efficiency model are combined.The pioneer conference paper by
Clément et al. [6] introduced the use of Bayesian analysis in the radioactive wastecharacterization community, using a standard random walk Metropolis (RWM) Markov chain Monte Carlo (MCMC)algorithm [see 10, for details about MCMC]. In this work, we propose a fully Bayesian approach to infer the verticaldistribution of the activities of the measured nuclides and to quantify the corresponding uncertainty. Our approach relieson state-of-the-art MCMC sampling using the Hamiltonian Monte Carlo (HMC) technique [19, 4] and accounts for thenet and background count uncertainties (i.e., Poisson statistics) together with the uncertainty in the detector’s efficienciescaused by the uncertainty in the source distribution. Furthermore, our efficiency model accounts for the contributions ofall segments to each count measurement given the assumption of a constant matrix density and corresponding gammaray attenuation. To the best of our knowledge, this is the first application of Bayesian inference to estimation anduncertainty quantification of spatially-distributed activity from gamma spectrometry.The remainder of this paper is organized as follows. Section 2 describes the considered waste drum, the used SGSapparatus and the measured data, together with the considered model and the MCMC-based Bayesian inversion. Section3 then details the results of our high-dimensional MCMC sampling for spatially-distributed radionuclides’ quantificationfor both a synthetic and a real case, before section 4 provides some discussion and outlines some future developments.Finally, section 5 summarizes our main findings and provides a short conclusion.
The non-destructive technique used in this study to identify and quantify the gamma-emitters present in a radioactivewaste drum is segmented gamma scanning (SGS), which scans the drum slice by slice (segment), while the wastedrum is rotating. We have used the so-called 3AX SGS [5]. The 3AX device can scan in the horizontal plane whilethe drum is rotating. Here, we only used the rotation capability and varied the height of the 3AX detector alongthe vertical axis of the drum. A 3AX SGS is made of three parts: the mechanics, the detector system and the dataacquisition and processing unit. The mechanical system consists of an evaluation unit for the turntable and two trolleyspositioned respectively on the left and right hand side of the turntable opposite to one another. The right trolley holdsthe detector-collimator system while the left one can be equipped with a transmission source (which has not beenused here). Both trolleys can move horizontally towards and from the drum. The SGS is equipped with a high puritygermanium detector (HPGe) connected to a cooling unit using liquid nitrogen. This detector is fixed in a lead collimator,enabling it to “see" only a fraction of the waste package. The slice of a drum that can be seen by the detector is definedby the aperture of the collimator and the distance between detector and waste package. The output of the detectoris amplified by a preamplifier and a spectroscopy linear amplifier. This pulse is converted into a digital number bythe analogue-to-digital converter (ADC) for further computer-based processing, and is in our case combined with themultichannel analyzer (MCA). The used detector was calibrated by applying a custom correction function to anothersimilar detector that is ISOCS/LabSOCS-calibrated [17]. Data acquisition is performed and controlled by the mainprogram which communicates with the Genie-2000 [12] software for conventional activity estimation.Note that for nuclear fuel measurements, better activity estimates could be obtained by using a dedicated isotopiccomposition software [e.g., FRAM, MGA 23, 18] instead of Genie-2000. However, here estimation accuracy is notas critical as usual because these values only serve to define a weakly informative prior distribution for the MCMCinference (see section 2.5). 2 .2 Waste package and setup characteristics
We used the 3AX scanning results of a real non-conditioned radioactive waste package. The latter, a galvanized200-liter drum (height: 890 mm; diameter: 580 mm; wall thickness: 1.25 mm), contains waste originating from thedecommissioning of a MOX (Mixed Oxide Fuel) glove box. The matrix of the waste mainly consists of metals such ascarbon steel ( ∼
70 wt%) and lead ( ∼ ∼
22 wt%) is present. The drum was filled up completely with a net weight of54 kg. This roughly corresponds to an average matrix density of 0.27 g/cm . The drum dimensions, filling degree,matrix composition and matrix density distribution are considered fixed for this study and the associated uncertaintiesare therefore not taken into account. Instead, we focus on one of the main sources of uncertainties for this kind ofmeasurements: the source distribution uncertainty (see section 2.3 for details).Considering 20 individual measurements starting from the bottom of the drum and raising the detector 46.8 mm foreach next measurement, the drum can be discretized into 20 horizontal segments (Figure 1): two segments of 23.4 mmheight for bottom and top measurement (S1 & S20) and 18 segments of 46.8 mm height (S2 − S19). Measurementsetup details are as follows:• The detector is surrounded by a rectangular lead slit collimator with a 1 mm copper coating (5 ◦ opening; 18mm x 110 mm; 100 mm depth);• The 200-liter drum is placed on a turning table and is continuously rotated during the measurements.• The measured 20 segments result in 20 individual spectra (M1 up to M20) and measurement time is 300 s. Forthe first measurement, the bottom of the drum is placed in front of the middle of the detector (M1). For eachfollowing measurement the detector is raised by 46.8 mm. In this way, the top of the drum is placed in front ofthe middle of the detector for the 20 th measurement (M20, see Figure 1).Figure 1: Vertical cross sections of the 3AX measuring setup for a 200-L drum. (a) 20 positions of the detector (M1 upto M20, blue) indicating the solid angle umbra (blue) for each measurement and the penumbra for measurement M11.The contents of the 200-L drum is discretized into 20 cylindrical volumes (segments) (red): S2 up to S19 in the middleof the drum with a height of 46.8 mm and S1 (bottom) & S20 (top) with a height of 23.4 mm. (b) Homogeneous (grayrectangles) and point (red dots) source distribution within each segment.Table 1 lists the identified nuclides in the drum and their associated energy peaks.3able 1: Detected nuclides and associated photon energies.Nuclide Energy [keV]Am-241 125.3662.4722.01Pu-238 152.72766.39Pu-239 129.30345.01375.05413.71451.48Pu-240 160.31Pu-241 148.57 To account for the effect of the uncertainty associated with the source distribution within the drum on the detectorefficiencies, we modeled the whole measurement system using the complex cylinder model of the Geometry ComposerV4.3 library of the ISOCS/LabSOCS software by
ISOCS/LabSOCS - Mirion Technologies (Canberra), Inc. [17]. Takingthe drum rotation during the measurement into account, we modeled two extreme source distributions (Figure 1):• Maximum efficiency: a homogeneous source distribution (hypothesis h ) within each segment (Figure 1a).• Minimum efficiency: a point source (hypothesis p ) placed on the drum axes on the border of a segment (on thebottom of the segment for S1 up to S10 and on the top of the segment for S11 up to S20, see Figure 1b)Since there is symmetry between measurements M1-M10 and measurements M11-M20 (Figure 2), the efficiencycalculations for the lower part of the drum could be duplicated to the upper part of the drum. We designed an individualmodel for every individual measurement at the lower part of the drum (M1 up to M10), placing the source stepwise insegments S1 up to S20 for the homogeneous (hypothesis h ) and point source (hypothesis p ) scenarios, respectively. Thismeans that we designed a total of 400 models (10 detector locations ×
20 source locations × h and p assumptions, the detector’s efficiency was thus calculated for each configuration ofsource location (S, the source is located in a given segment while the other 19 segments do not contain any source),detector location (M), and photon energy (among the 12 energies listed in Table 1). These efficiencies were thenencapsulated into 2 n source × n peaks × n seg
3D arrays, E h and E p , for the h and p assumptions, respectively. Here n source = n seg = 20 and n peaks = 12 . By using a λ ∈ [0 , coefficient, one can balance the efficiency (cid:15) ( e ) associated with a given detection at energy e between the h and p assumptions for a given segment : (cid:15) ( e ) = λ(cid:15) h ( e ) + (1 − λ ) (cid:15) p ( e ) . This gives rise to our followingcount simulation model for a given count rate, c rate c ratei,k = n seg (cid:88) l =1 a j,k =1 ··· n k ,l P γ j,k (cid:104) λ i (cid:15) hi,l ( e k ) + (1 − λ i ) (cid:15) pi,l ( e k ) (cid:105) , (1)where the subscript i indexes the detector location, i = 1 , · · · , , the subscript j indexes the considered nuclide, j = 1 , · · · , , the subscript k indexes the considered energy peak, k = 1 , · · · , , a j,k =1 ··· n k ,l denotes the activity ofnuclide j at segment location l that is associated with the k = 1 · · · n k enery peaks, and P γ j,k is the emission probabilityof nuclide j at energy k .The simulated gross count at detector location i for the energy k is then calculated as c grossi,k = c ratei,k t m + b i,k , (2)where t m is the measurement time and b i,k is the background continuum count at detector location i for the energy k .4 .5 Bayesian inference A common representation of the forward problem is d = F ( θ ) + e , (3)where d = [ d , . . . , d N d ] ∈ R N d , N d ≥ is the measurement data, F ( θ ) is a deterministic forward model withparameters θ and the noise term e lumps all sources of errors.In the Bayesian paradigm, parameters in θ are viewed as random variables with a posterior probability density function(pdf), p ( θ | d ) , given by p ( θ | d ) = p ( θ ) p ( d | θ ) p ( d ) ∝ p ( θ ) L ( θ | d ) , (4)where L ( θ | d ) ≡ p ( d | θ ) signifies the likelihood function of θ . The normalization factor p ( d ) = (cid:82) p ( θ ) p ( d | θ ) d θ isnot required for parameter inference when the parameter dimensionality is fixed. In the remainder of this paper, we willthus focus on the unnormalized density p ( θ | d ) ∝ p ( θ ) L ( θ | d ) .If we assume d to follow a Poisson process, which is the norm for count data, L ( θ | d ) can be written as L ( θ | d ) = N d (cid:89) i =1 exp (cid:16) − ˜ d i (cid:17) ˜ d id d ! , (5)where ˜ d = (cid:104) ˜ d , . . . , ˜ d N d (cid:105) = F ( θ ) contains the simulated responses.Here the vector of inferred variables, θ , consists of the 100 activities, a (5 nuclides ×
20 segments), 237 backgroundcontinuum counts b (12 energy peaks ×
20 segments less the 3 measurements for which neither a background continuumcount nor a net count could be measured) and 20 λ coefficients. As of the d vector, it comprises N d = 237 gross countsobtained from the sum of the 237 net counts and 237 background continuum counts. In addition, our derived posteriordistribution is not only conditioned to d but also to the chosen E h and E p arrays. This results in assessment of the p (cid:16) a , b , λ | d , E h , E p (cid:17) distribution p (cid:16) a , b , λ | d , E h , E p (cid:17) ∝ L (cid:16) a , b , λ | d , E h , E p (cid:17) p ( a , b , λ ) . (6)The marginal posterior distribution of a given quantity, say a , is obtained by integrating the posterior distribution overall other inferred variables p (cid:16) a | d , E h , E p (cid:17) = (cid:90) (cid:90) (cid:90) p (cid:16) a ∼ , b , λ | d , E h , E p (cid:17) d a ∼ d b d λ , (7)where the a ∼ vector contain all elements of a but a .We assume the prior distributions for a , b and λ to be independent, which gives p ( a , b , λ ) = p ( a ) p ( b ) p ( λ ) . (8)For the activities, a , we have some prior knowledge in the form of estimates, ˆ a , derived by application of the Genie-2000 [12] procedure to each segment separately. Assuming that these estimates have an accuracy of ± log (ˆ a ) ± . Further assuming that log ( a ) is normallydistributed and that log (ˆ a ) ± provides a 99% uncertainty interval induces a standard deviation, σ a , for log ( a ) of σ a = 13 log (10) where the log (10) term accounts for the change of basis between log ( · ) and log ( · ) . This formsthe rationale for using a lognormal prior for a : p (log ( a )) = N (cid:0) log (ˆ a ) , σ a I (cid:1) . As demonstrated in section 3 this is aweakly informative prior from which the posterior distribution will easily depart if needed to appropriately fit the countdata. Note that the values in ˆ a are likely to be overestimations as interpreting each segment separately implies that allthe measured counts at a given detector location are assigned to the activity of the corresponding segment, which isincorrect.Cases for which the detected net count(s) for a given (set of) measurement(s) { i, k = 1 , · · · , n k } , are too smallcompared to the derived background continuum count(s) for the Genie-2000 procedure to estimate an activity value,are easily included in the analysis as for non-influential activity value(s), the Bayesian approach will simply returnthe prior distribution. Overall, these “non-detects" concerned 43 out of the 100 inferred activities, including the 20Pu-240 activities. In these cases the Genie-2000 procedure returns a calculated minimum detectable amount (MDA)5hich means that in absence of any other information, the “true" value could be anything between zero and somethingsomewhat close to the calculated MDA. For those activities, a mda , we decided to also use a lognormal prior but withmean parameter set to MDA/2: p (log ( a mda )) = N (cid:0) log (0 . MDA ) , σ a I (cid:1) . Indeed, some practitioners use half theMDAs as working values. This is arguably subjective and another choice could be to use an uniform prior between zeroand some upper bound related to the MDA, say: p ( a mda ) = U ( , MDA ) . Note however that no matter the specifiedupper bound, the inference should point out automatically that activities above the MDA are inconsistent with theobserved counts. This is further discussed in section 4.Regarding the 237-dimensional background continuum count vector, b , the obvious choice is to set p ( b ) as the productof independent Poisson prior distributions with shape parameter equal to the (rounded) measured background continuumcounts: P ois ( b ) . However, for technical reasons linked to the used software (see next section) we cannot use aPoisson prior (in short, Poisson priors cannot be assigned to real variables). Instead we use for p ( b ) an uncorrelatedand independent normal prior distribution with mean and variance vectors both equal to the measured count values: p ( b ) = N ( b , C b ) with C b a diagonal matrix with the b , · · · , b values as diagonal elements. A N ( x, x ) distributionprovides an increasingly accurate approximation to P ois ( x ) as x increases [see, e.g., 3]. The approximation iscommonly deemed excellent for x > and reasonably accurate for x > [24]. Lastly, the prior distribution forthe 20-dimensional λ vector is taken as a bounded uniform distribution: U ( , ) .As no analytical solution of the 357-dimensional p (cid:16) a , b , λ | d , E h , E p (cid:17) distribution is available, we sample from p (cid:16) a , b , λ | d , E h , E p (cid:17) by MCMC simulation [see 10] using a state-of-the-art implementation of the HMC sampler [see19, 4, for an extensive description]. Convergence of the MCMC to the posterior target is monitored by means of thepotential scale reduction factor, ˆ R [9, 10], using four independent Markov chains evolved in parallel (see next section).The ˆ R statistic compares for each parameter of interest the average within-chain variance to the variance of all theMarkov chains mixed together. The closer the values of these two variances, the closer to one the value of the ˆ R diagnostic. Values of ˆ R jointly smaller than 1.02 for all sampled variables are deemed to indicate convergence to alimiting distribution. We used the open-source greta package [13] to perform the HMC-based MCMC sampling.The greta package is an R[22] interface to some of the MCMC sampling algorithms implemented in the Tensorflow-probability package [TFP,7] which itself relies on the Tensorflow (TF) machine learning platform [1]. The standard language to interact withTFP is Python [21] and the TFP-Python syntax is somewhat complicated. The main advantage of greta is to make itstraightfoward to build probabilistic models by constructing a directed acyclic graph and to perform MCMC samplingwith TFP, using a seamless R-like syntax. The most useful MCMC sampler available through greta and used herein isHMC [19]. The HMC sampler has proven to be quite efficient when the gradient of the likelihood or objective functioncan be calculated by applying the chain rule, using an automatic differentiation (AD) technique such as implemented inTFP. Not all numerical models are suited to AD (for instance, numerical solvers of partial differential equations are not)but numerical models used for routine interpretation of radiological characterization measurements typically are. TheTFP MCMC algorithms can be run in parallel on both CPUs and GPUs which makes greta to be quite computationallyefficient. In this study, we ran four separate HMC trials in parallel over 4 CPUs.
The four parallel HMC runs were performed for a total of 30,000 warmup iterations each, after which 10,000 posteriorsamples were collected per run. The ˆ R convergence criterion computed using these 4 Markov chains was below 1.02 forall of the 357 sampled variables, indicating official convergence to the posterior target by the four runs. The collected40,000 samples were thus used to approximate p (cid:16) a , b , λ | d , E h , E p (cid:17) . On the used 4-core notebook (equipped with Inteli7-6820HQ CPU @ 2.70GHz), achieving these 40,000 iterations per/run incurred a computational time of about 2 hoursin total. Before proceeding with the inversion of the real data, we study in this section whether our proposed approach is ableto correctly recover the activities for a synthetic problem for which the true values are known. We used the samedrum, efficiency data and count simulation model as for the real case, but considered only one nuclide and energy peak:Pu-239 at 413.71 keV. The “true" λ was set to 1 in every segment except for segment 16 where it was set to 0. The6u-239 activity was set to 1 × Bq in every segment except for segments 5 and 16 were it was set to 1 × Bq. Theresulting i = 1 , · · · , simulated net counts, c nettrue,i , were used as mean parameters of 20 Poisson distributions fromwhich the 20 observed nets counts were drawn: c netobs,i ∝ P ois (cid:0) c nettrue,i (cid:1) . The 20 observed background continuum counts, b obs,i , were also sampled from a Poisson distribution with a mean parameter inspired from the real measurements: c backobs,i ∝ P ois (100) . Note that due to the added errors in the observed net counts, no perfect inverse solution existswhich is required for Bayesian inference. The same normal p ( b ) prior distribution and uniform p ( λ ) prior distributionas for the real case were used. We defined p ( a ) = N (cid:0) log (ˆ a synth ) , σ a I (cid:1) with the 20 elements of the ˆ a synth vectorbeing × Bq. This way, the prior mean of the Pu-239 log-activity is (i) one order of magnitude too large for the 18segments where the “true" Pu-239 activity of the waste ( × Bq) is small compared to the background activity, and(ii) two order of magnitude too large for the 2 segments where the “true" Pu-239 activity is large ( × Bq). Inaddition, setting the standard deviation of the log-prior activity to σ a (that is, ) induces a wide and thus quiteuncertain prior.Figure 2 presents the resulting posterior parameter distribution. Note that segment numbering goes from bottom(segment 1) to top (segment 20). With respect to the activities (Figure 2a), it is seen that the MCMC sampling recoversthe 2 large true activities in segments 5 and 16 very accurately, with a high degree of certainty. This even though themode of the lognormal prior is way off. The posterior log-activity distributions in the other segments (all but 5 and 16)often peak near the true values and are somewhat asymetric with a larger probability mass at the right side. This isbecause here the MCMC sampling finds a balance between simulating sufficiently small gross counts and honoring theprior. This nicely illustrates how Bayesian analysis handles MDAs (non-detects). If the range of values spanned by theprior distribution is sufficiently low, then it simply returns the prior distribution. Otherwise, the posterior distributionforms an intermediate distribution between the prior and the likelihood (equation (4)). The posterior backgroundcontinuum count distributions globally peak around the true values, with slight deviations from their associated priors(Figure 2b). As expected, the posterior λ distributions are the same as their prior counterparts for the low activities(Figure 2c). For the 2 segments where the true activity is large, the inference put most of the probability mass towardsthe true boundary values of λ , λ = 1 for segment 5 and λ = 0 for segment 16. Overall these results demonstrate thatour approach can correctly infer the joint posterior distribution of a , b and λ .7igure 2: Prior and posterior (a) activity, (b) background continuum count and (c) λ distributions for the synthetic testcase. The black dots denote the true values used to generate the gross count data. Note that the x -axis of subplot (a) hasa base 10 logarithmic scale. For visual convenience, the background continuum counts have been standardized using thenormal prior parameters. For the λ variable, we show a posterior histogram instead of a kernel density estimate. This isbecause when applied to a bounded data sample, kernel density smoothing tends to create artifacts near the bounds. Figure 3 depict the marginal posterior distributions of the inferred Am-241, Pu-238, Pu-239, Pu-240 and Pu-241activities: p (cid:16) a | d , E h , E p (cid:17) . It is observed that the “multi-energy" nuclides, that is, Am-241, Pu-238 and Pu-239 whichare associated with three (Am-241), two (Pu-238) and five (Pu-239) energy peaks, respectively (Table 1), are globallywell resolved. This is indicated by a posterior distribution that is narrower than the prior distribution. Strikingly, for8ll 100 activities but that of Pu-239 in segment 16, the posterior mode is lower than the prior mean of the log-activity( log (ˆ a ) or log (0 . MDA ) , see section 2.5). As written above, this is because processing each segment individually,as done to define p ( a ) , is likely to overestimate the actual activities as it overlooks the fact that the detected countsat a given detector location are not only caused by the activity in the considered segment but also by the activities inthe neighboring segments. In contrast, for a given nuclide our net count rate simulation model (see equation (1) andassociated text) accounts for the contributions from all activities to each measurement. For Pu-239 at segment 16 only(Figure 3), the posterior mode of the log-activity is above the selected prior mean: about 3.9 × Bq against 1.13 × Bq. Note that since this activity is among those for which the Genie-2000 procedure returned a MDA value, herethe posterior mode is about two times the MDA (prior mean of the log-activity was taken as log (
M DA/ ).The “single-energy" Pu-240 and Pu-241 nuclides are relatively well resolved though, perhaps not surprisingly, someare more uncertain than the multi-energy nuclides. This is especially the case for Pu-240 at segment 1 and Pu-241 atsegments 13 −
16 and 20. Overall, at these locations the 95% uncertainty intervals for Pu-240 and Pu-241 cover aboutone order of magnitude.Table 2 lists the corresponding 95% posterior uncertainty (or credible) intervals of the total activity of each nuclide overthe whole drum. For this package, the Pu-241 nuclide has by far the largest activity, one to two orders of magnitudehigher than for the other nuclides. As written earlier, the measured net count for the 20 Pu-240 activities is zero and, inthis case, the prior distribution was set to N (cid:0) log (0 . MDA ) , σ a I (cid:1) . Since the range spanned by the posterior Pu-240log-activities is systematically lower than the range spanned by the prior Pu-240 log-activities, our results show that forthe considered case study setting the Pu-240 non-detects to half the MDA is overly conservative.Table 2: 95% posterior uncertainty intervals [10 Bq] of the 5 nuclides’ activities across the whole drum.Nuclide 95% uncertainty [10 Bq]Am-241 . − . Pu-238 . − . Pu-239 . − . Pu-240 . − . Pu-241 . − . x -axis.Figure 4 provides more insights into the posterior activity distribution by displaying the (Pearson) linear correlationcoefficients between (a) the inferred Am-241 activities in each segment (Figure 4a) and (b) the 5 nuclides across allsegments (Figure 4b). The Am-241 correlations between segments (Figure 4a) are representative of the between-segment correlations of the other nuclides. The following pattern is observed. For a given nuclide, there is a negativecorrelation between direct neighbor (or lag-1) activites. This makes sense as direct neighbors contribute more thansegments further away to the simulated count at a considered location. Hence, the MCMC inference can counterbalancethe effect of increasing the activity in a segment by decreasing the activity in a direct neighbor location. Lag-2 activities(that is, activities separated by one segment) then correspondingly show a weak positive correlation while correlationsfade away from lag-3 on. Moreover, looking at the correlations between the 5 nuclides across all segments mergedtogether (Figure 4b) reveals only one important correlation: the correlation coefficient between Pu-238 and Pu-241 is0.85.To complete this analysis of the posterior activity distribution, Figure 5a displays the total posterior activity through theconsidered drum, from bottom (segment 1) to top (segment 20). Globally, the total posterior activity decreases frombottom (segment 1) to top (segment 20), from about . × Bq to less than × Bq.10igure 4: Posterior (Pearson) correlations corresponding to every pair of inferred activities.11igure 5: Posterior distributions of (a) the total activity and (b) the λ parameter within each segment. The chosenuniform prior For λ is also shown in subplot (b). No prior distribution is shown for the total activity as no closed formexpression is available for the sum of lognormal distributions. Segment numbering goes from bottom (1) to top (20).Note that the x -axis of subplot (a) has a base 10 logarithmic scale. Note that for the λ variable, we show a posteriorhistogram instead of a kernel density estimate. This is because when applied to a bounded data sample, kernel densitysmoothing tends to create artifacts near the bounds. 12he marginal posterior background continuum activities, p (cid:16) b | d , E h , E p (cid:17) , are depicted in Figure 6. The posteriordistributions are generally close to their associated prior distribution but some depart from it. Consistently with ourBayesian framework, these deviations are needed to maximize the posterior density. We would like to stress that the prior, p ( b ) , is based on a measured background continuum count that is considered to be a realization of a Poisson distributionof which the shape parameter is unknown. By setting p ( b ) = P ois ( b ) (or p ( b ) = N ( b , C b ) ≈ P ois ( b ) as doneherein) one assumes that the measured background continuum count is equal to the mean of its underlying distribution,which is obviously not necessarily the case. Some deviations from p ( b ) should thus come at no surprise. If deemednecessary, an alternative solution to enforce a tight closeness between N ( b , C b ) ≈ P ois ( b ) and p (cid:16) b | d , E h , E p (cid:17) isdiscussed in section 4.Figure 6: Prior and posterior distributions of the inferred background continuum counts. Each facet’s label mentions theconsidered combination of nuclide and energy peak (keV). For visual convenience, the counts have been standardizedusing the normal prior parameters and only 7 out of the considered 12 peaks are shown (the following energy peaks areleft out of the plot: Am-241 & 662.24 keV, Pu-238 & 152.72 keV, Pu-239 & 129.3 keV, and Pu-239 & 375.05 keV).Segment numbering goes from bottom (1) to top (20).Figure 5b presents the marginal posterior densities of the elements of the λ vector, which is our parameter that accountsfor the uncertainty in the efficiencies. For segments 2 to 19 the posterior probability mass is mostly concentrated on theupper bound of 1, meaning that the MCMC is strongly favoring the hypothesis of a homogeneous radioactive source13istribution in these locations. In contrast, for the bottom (1) and top (20) segments every intermediate state between afully homogeneous source distribution and a fully point source distribution has non-zero probability. it is importantto note that our efficiency model (see section 2.3) might be biased, especially for the top segment. Hence, a smallerfilling height than the drum top might cause the large uncertainty in the source distribution for top segment 20 and thelong-tail uncertainty distribution in the underlying segment 19.The fit of the posterior gross count simulations, F (cid:16) a , b , λ , E h , E p (cid:17) to the measured gross counts, d is depicted inFigure 7 using a base 10 logarithmic scale to make the discrepancies between (very) small measured and simulatedcounts visible. The vertical bars denote the 95% posterior uncertainty intervals associated with the simulated counts.The observed counts are mostly well fitted and the 95% uncertainty intervals are relatively tight and most of the timeinclude the 1:1 line. Furthermore, the most important discrepancies, in relative terms, are for some rather small grosscounts, below 10 −
15. This is consistent with both our expectations and the Poisson count statistics.Figure 7: Measured gross counts against a posteriori simulated gross counts. The green points represent the maximumlikelihood solution (solution with the largest Poisson likelihood) among the collected 40,000 posterior samples, and thevertical bars denote the 95% posterior uncertainty intervals. A base 10 logarithmic scale is used to make the deviationsbetween the smallest observed and simulated gross counts visible.
We have shown that vertical activity profiles can be probabilistically inferred from SGS data accounting for the effectof the source distribution uncertainty on the detector efficiency. The following points and/or potential improvementsnevertheless deserve further attention.As for any prior distribution, the choice of the priors for the activities flagged as MDAs by the standard analysisprocedure [12] applied to each segment separately, a mda , is a relatively subjective expert-based decision. Note, however,14hat if the prior is sufficiently weakly informative or wide and the information content of the count measurement datais large enough (sufficiently large net counts compared to the corresponding background counts and not too manymissing values), then the actual shape of the prior should not substantially influence the posterior results. An additionalMCMC run with using p ( a mda ) = U ( , MDA ) instead of p (log ( a mda )) = N (cid:0) log (0 . MDA ) , σ a I (cid:1) shows that,unsurprisingly, p ( a mda ) = U ( , MDA ) leads to posterior distributions for the elements of a mda that are close to theirassociated uniform prior (Figure 8). For the non-MDAs, the two posterior distributions are very similar (Figure 8).In some relatively few cases, the marginal posterior background continuum activities, p (cid:16) b | d , E h , E p (cid:17) , differ largelyfrom the used prior distribution for b , N ( b , C b ) ≈ P ois ( b ) . This is probably due to some small inconsistenciesbetween our model and the real drum. Notwithstanding, if having p (cid:16) b | d , E h , E p (cid:17) closely approximating N ( b , C b ) ≈ P ois ( b ) is judged necessary, then setting p ( b ) = N ( b , s b C b ) with s b < can do the job. Indeed, we found forthe considered case study that using s b = 0 . allows for p (cid:16) b | d , E h , E p (cid:17) to closely reproduce N ( b , C b ) while onlyslightly changing p (cid:16) a | d , E h , E p (cid:17) . However, this comes at the expense of a less good fit of the simulated gross counts, F (cid:16) a , b , λ , E h , E p (cid:17) to the measured gross counts, d (not shown).This study only considers the count measurement uncertainty and the uncertainty of the source distribution (albeit ina simplified way), the latter being one important source of uncertainty in the detector’s efficiency. Nonetheless, thematrix density distribution of the package is assumed to be homogeneous. Since the matrix-related uncertainties canbe relatively large and thereby impacting the detector’s efficiency too, more work is needed to account for it withinthe inference. This could be approached by using a Monte Carlo particle transport code to simulate the impact ofvarious matrix and source configurations on the efficiency of the considered detector. The so-derived dataset couldthen be used to construct a computationally cheap nonlinear regression model (also called proxy or metamodel) thatwould predict the simulated efficiencies from the sampled matrix and source parameters, for each considered nuclideand segment location. Both matrix-related and source distribution uncertainty could thus be accounted for within theBayesian inference by sampling the surrogate model parameters along with the other variables. This warrants furtherinvestigations.More generally, an extended Bayesian approach based on the one presented herein could also be used to infer spatialactivity distributions in 1D, 2D and 3D within a drum waste package from combined measurements from differenttechniques. This will be investigated in future work, in the framework of the MICADO EU project (see section 6).15igure 8: Posterior activity distributions obtained when using (1) p (log ( a mda )) = N (cid:0) log (0 . MDA ) , σ a I (cid:1) (Posterior1) and (2) p ( a mda ) = U ( , MDA ) . For visual convenience, the lower limit of the x -axis is thresolded to × Bq.Also, note the base 10 logarithmic scale of the x -axis. We propose a Bayesian approach to probabilistically infer vertical activity distributions within a radioactive wastedrum from segmented gamma scanning (SGS) measurements. Our approach relies on state-of-the-art Markov chainMonte Carlo (MCMC) sampling using the Hamiltonian Monte Carlo (HMC) technique and accounts for the net andbackground count measurement uncertainty together with the uncertainty in the detector efficiencies caused by theuncertainty in the source distribution within the drum. Furthermore, our efficiency model accounts for the contributionsof all considered segments to each count measurement. For the considered case study, our approach is used to resolvethe vertical activity distribution of 5 nuclides over 20 locations and produce sound uncertainty estimates. Main topics forfuture research include how to incorporate the impact of the uncertainty in matrix density distribution on the modeledefficiencies, to extend our setup from 1D to 2D and/or 3D, and to combine measurements from different techniques forbetter radionuclides’ quantification. 16
Acknowledgments
This work received funding by the EU project MICADO Consortium. 2018: “Horizon 2020, Nfrp 2018-10, Proposal847641: Measurement and Instrumentation for Cleaning and Decommissioning Operations".
References [1] Abadi, M., Agarwal, A., Barham, P., et al. TensorFlow: Large-Scale Machine Learning on HeterogeneousDistributed Systems. ArXiv e-prints, March 2016,
StatisticalScience,7 , 457–472.[10] Gelman, A., Carlin, J.B., Stern, H.S., Dunson, D.B., Vehtari A., and Rubin, D.B., 2014. Bayesian Data Analysis,Third Edition, Chapman & Hall/CRC.[11] Gelman, A., 2015. https://statmodeling.stat.columbia.edu/2015/08/25/can-you-change-your-bayesian-prior/ .[12] Genie 2000 Spectroscopy Software. Customisation Tools Manual 9230847F 3/01. Mirion Technologies (Canberra),Inc. .[13] Golding, N., 2019. greta: simple and scalable statistical modelling in R. J. Open Source Softw., 4, p. 1601.[14] Kirkpatrick, J. M., Venkataraman, R., Young, B. M. 2013. Minimum detectable activity, sys-tematic uncertainties, and the ISO 11929 standard. J. Radioanal. Nucl. Chem., 296, 1005–1010.https://link.springer.com/article/10.1007/s10967-012-2083-5.[15] Krings, T., Mauerhofer, E. 2011. Reconstruction of the activity of point sources for the accurate char-acterization of nuclear waste drums by segmented gamma scanning. Appl. Radiat. Isot., 69(6), 880–889.https://doi.org/10.1016/j.apradiso.2011.02.009.[16] The ISOCS (In Situ Object Counting System) Calibration Software . .[17] ISOCS/LabSOCS Detector Characterization Report. CSO .[19] Neal R (2011). MCMC Using Hamiltonian Dynamics. In S Brooks, A Gelman, GL Jones,XL Meng (eds.),Handbook of Markov Chain Monte Carlo, pp. 116–162. Chapman & Hall/CRC.[20] Patra, S., Agarwal, C. 2019. Segmented gamma-ray assay of large volume radioactive waste drums containingplutonium lumps. Appl. Radiat. Isot., 153, 108827. https://doi.org/10.1016/j.apradiso.2019.108827.[21] Python Software Foundation. 2020. .[22] R Core Team. 2020. R: A language and environment for statistical computing. R Foundation for StatisticalComputing, Vienna, Austria. .1723] Sampson, T.E., Kelly, T.E., Vo, D.T. 2003. Application guide to gamma-ray isotopic analysis using the FRAMsoftware, LA-14018, Los Alamos.[24] SOCR/UCLA. http://wiki.stat.ucla.edu/socr/index.php/AP_Statistics_Curriculum_2007_Limits_Norm2Poissonhttp://wiki.stat.ucla.edu/socr/index.php/AP_Statistics_Curriculum_2007_Limits_Norm2Poisson