A novel approach to the localization and the estimate of radioactivity in contaminated waste packages via imaging techniques
Michele Pastena, Bastian Weinhorst, Günter Kanisch, Edgar Perez Lezama, Johannes Ratdke, Tim Thomas
AA novel approach to the localization and the estimate of radioactivity incontaminated waste packages via imaging techniques
Michele Pastena a,b, ∗ , Bastian Weinhorst a,b , Günter Kanisch c , Edgar Pérez Lezama a,b , Johannes Radtke a,b , TimThomas a,b, a Safetec Entsorgungs- und Sicherheitstechnik GmbH b Kurpfalzring 98a, D-69123 Heidelberg, Germany c Wittland 44g, D-22589 Hamburg, Germany
Abstract
Dismantling nuclear power plants entails the production of a large amount of contaminated (or potentiallycontaminated) waste that must be disposed according to national and international regulations. A large partof the end products needs to be stored in special repositories, but a significant part of it is slightly contaminated ornot contaminated at all, making it possible to free release it. One possible approach to free release measurementsuses Large Clearance Monitors, chambers surrounded by plastic scintillation detectors that can measure up to1000kg of waste. Due to the composite nature of the detection system in a Large Clearance Monitor, it is easy toimagine that one can apply 3D imaging algorithms to localize radioactive sources inside a waste package. In thiswork we will show how a special algorithm that maximizes the conditional informational entropy allows decisionsabout the clearance of portions of the sample.
Keywords:
Clearance measurements, Conditional Entropy, nuclear waste imaging
Introduction
One major path for clearance of large amount ofwaste arising from the decommission of nuclear powerplants is based on measurements with Large ClearanceMonitors (LCM), devices that utilize multi-detectorplastic scintillation arrays (LCMs can consist of up to 24plastic scintillation detectors) to offer both, extremelylow minimum detectable activities (MDA) and shortmeasurement times for waste packages weighing up toone metric ton. LCMs represent in general the mostefficient way to measure big and heavy samples. Suchclearance measurements are performed also by meansgermanium detectors (or arrays of germanium detectors)with beam collimators both for waste packages (rotatingthe sample and measuring its emitted radiation in severaldifferent directions) and for whole rooms or buildings(measuring the radiation emitted by the surface of aroom by means of a collimated static detector). Ithas to be pointed out though that germanium detectorsare slow and inefficient as each single measurementrequires 25-30 minutes and furthermore such devicesrequire delicate cooling systems with liquid nitrogen. Inorder to handle the large amount of waste produced ∗ Corresponding author Tel. +49 (0) 6221 / 65175138
Email address: [email protected] (MichelePastena) during the decommissioning of a nuclear power plant,the measurement for the release procedure must besufficiently short. For this reason, Large ClearanceMonitors (LCM) are used for this task.However, a major problem with these devices isto obtain reliable information of so-called hotspots,i.e. contaminated or activated material in the waste-package. While the π -setup of the detectors allows acrude localization of one hotspot (Mirion TechnologiesRTM644Inc User Manual), no analysis has beenconducted so far to test the uncertainty of the estimatedactivity. Moreover, the localization of hotspots areroutinely implemented in a sufficiently high radioactiveenvironment [1, 2] by employing high-purity germaniumdetectors, but the low radioactivity needed for clearancemeasurements, presents a different problem. Rather thanposing the question where the activity is likely to belocated for signals which can be clearly distinguishedfrom the background, for clearance measurements allpossible geometrical and hardware uncertainties have tobe taken into account for the minimum detection limitand thus the possible activity distributions have to bederived according to ISO 11929. In order to solve theproblem of distinguishing the background radiation fromthe one produced by a sample with a contaminatedhotspot, we propose a new method that resembles theimaging techniques used in medical physics.In this paper we want to present a novel approach a r X i v : . [ phy s i c s . d a t a - a n ] J a n o the localization and the estimate of radioactivityin contaminated waste packages. The new procedureuses an imaging algorithm to reconstruct accurateimages of the radioactivity distribution in a given wastepackage, allowing fast and reliable decisions for clearanceof potentially contaminated or activated material.The proposed method relies on the maximization ofinformational entropy. This algorithm is used in themedical field for the imaging of tissues in the humanbody by means of PET/SPECT scans.In the section 1 we will summarize the mostcommon techniques used to determine the activityof low activated material that can be potentiallyfree released (clearance measurements) mentioning themodels and techniques prescribed by the internationalstandards. In the sections from 2 to 4 we will presentthe Conditional Entropy Maximization, an imagingtechnique widely used in medical physics and we willshow the changes needed to adapt it to the purposeof clearance measurements. In section 5 we willpresent the mathematical machinery needed to make theaforementioned algorithm compliant with the regulatorynorms regarding the determination of the characteristiclimits (decision threshold, detection limit and limitsof the coverage interval) for measurements of ionizingradiation provided by the international standard. In thefinal section 6 we will apply the Conditional EntropyMaximization to some simulated samples. This researchdid not receive any specific grant from funding agenciesin the public, commercial, or not-for-profit sectors.
1. State of art of the clearance measurements
Usually a LCM consists of a chamber encompassed by agamma ray detection system distributed over the wholesolid angle around the sample, consisting of 24 polyvinyltoluene (PVT)-based large volume plastic scintillationdetectors. In order to reduce the background radiation,massive shielding (i.e. 130 mm steel) surrounds thechamber. The activity of the waste is then estimatedusing a combined measurement of the 24 count rates.This way the LCM can measure up to 1000 kg of waste inone minute. The LCM is also equipped with a conveyorsystem for the waste packages.One major drawback of the plastic scintillation detectorsused in the LCMs is the poor energy resolution thatmakes the identification of the radionuclides in the wastealmost impossible. This means that gamma spectracannot be identified, but only the total count in agiven time (usually 60 s) is measured for a really widerange of energies. The estimate of the activity dependsthen on the knowledge of the radionuclide inventory,also known as nuclide vector (NV). Furthermore, thelacking identification of gamma lines does not allow todiscriminate between natural radio nuclides and those
Figure 1: Reproduction of the chamber of a Large ClearanceMonitor produced during the operation of the nuclear powerplant.To obtain the activity of a waste package, the totalcount rate and a calibration factor w are needed. Thefactor w = 1 /ε is the inverse of the detection efficiencyof the whole system (waste and measurement device)i.e. the ability of detecting the photons emitted by theradioactive sources, and depends on the composition ofthe nuclide inventory. However, the absolute efficiencydepends not only on the detection properties of thedetector (also known as detector characterization), butis also determined by the waste package, material, andactivity distribution (referred to as geometry). The aforementioned calibration factor w is one ofthe quantities needed to estimate the activity of awaste package. A physical model of the connectionbetween the activity and the other physical observables(count of gamma photons, shielding) is suggested bythe international standard ISO 11929 for each typeof measurement and environment (high activity, lowactivity, clearance measurement, etc.). The reason ofthe importance of this international standard lays onthe fact that the law in the field of nuclear safetyprescribes the adoption of state-of-art models, techniquesand technologies that make the handling of radioactivesources as safe as possible. The duty of the ISO 11929standard is then to summarize all the most recent andefficient models and prescriptions to be adopted. Withreference to the ISO 11929:2019 standard, from now onwe will take care exclusively of clearance measurements2i.e. measurements of samples with such a low activitythat they are potentially free releasable). Once we knowthe calibration factor w of the measurement device thatwe are using (both plastic and germanium detectors), thesimplified model to evaluate the activity a is a = ( r g − r · x − x + x ) w (1)with• r g gross count rate measured with the sample in thechamber,• r background count rate measured with emptychamber,• x correction of the background count rate for itsvariability due to work activities near the device forclearance measurements,• x correction for natural radionuclides, i.e. K andthe Ra and T h decay series, in the material tobe measured,• x correction for shielding of the background by thematerial to be measured.It has to be remarked that in the present work we willbe analyzing only Co sources that emit two photonsper decay with a probability close to 1. In general itcan happen that the nuclide inventory (i.e. the typeand relative percentage of radionuclides) is much morecomplex, requiring more attention in the formulationand solution of the problem. In order to lower theexternal dose of ionizing radiation, international normsprescribe to fully characterize the activity distributionof a sample (e.g. a waste package) to be analyzed.The core idea is that it is not sufficient that themeasurand is below the prescribed threshold, but thatwe also need a high statistical confidence (typically 95%).This can be achieved either by comparing case by casethe whole statistical distribution with the prescribedthresholds or by calculating the characteristic values (inparticular the upper limit of the coverage interval and thedetection limit). The second option is clearly the moststraightforward, and to this end we need the uncertaintyand all the statistical characteristic values according to[3]. The uncertainty for the activity (1) is given by u ( a ) = w (cid:20) r g t g + x r t + r x u ( x )+ u ( x ) + u ( x ) (cid:3) + a u rel ( w ) with u rel ( w ) = u ( w ) /w relative uncertainty of w .This formula is the starting point to calculate all thecharacteristic limits of the activity distribution, as it is explained in Appendix A. Indeed, the uncertainty for anassumed true value ˜ a is then ˜ u (˜ a ) = w (cid:20) t g (cid:18) ˜ aw + n t x + x + x (cid:19) + x r t + r · u ( x ) + u ( x ) + u ( x ) (cid:3) + ˜ a u rel ( w ) and thus we can calculate the following characteristicvalues for the simplified evaluation model:• decision threshold a ∗ = k − a ˜ u (0) = k − α w (cid:26) t g (cid:18) n t x + x + x (cid:19) + x r t + r u ( x ) + u ( x ) + u ( x ) (cid:9) / where k − α is the (1 − α ) -quantile of a normaldistribution (typically α = 0 . ),• detection limit: it is obtained by solving therecursive equation a = a ∗ + k − β ˜ u (cid:0) a (cid:1) (2)where k − β is the (1 − β ) -quantile of a normaldistribution (typically β = 0 . ),• upper and lower limit of the coverage interval a (cid:47) = a − k p u ( a ) ; a (cid:46) = a + k q u ( a ) where p = ω (1 − γ/ , q = 1 − ωγ/ , with ω = Φ ( y/u ( y )) probability corresponding to y/u ( y ) (assuming that the activity is normally distributed),and γ = 0 . .These definitions hold in the case of normally distributedactivities. In the most general case, we do not needto specify a probability distribution as it is discussedin Appendix A. As it has been pointed out, thesevalues allow us to make decisions about the riskiness ofcontaminated and nuclear waste in the case of clearancemeasurements.
2. Conditional Entropy Maximization forcontaminated waste
In the field of medical physics, reconstructing images ofthe inner parts of a human body for clinical purposeshas always been a challenge. The Single PhotonEmission Tomography (SPECT) is mostly used to detectfunctional or metabolic images of a tissue [4] and hasanalogies with the characterization of radioactive waste3ackages. The core idea of the SPECT in medical physicsis to inject a radioactive substance (typically emittingX rays) in living tissues and measure the count rate ofphotons out of the body along several different directions.This is achieved by encompassing the body in a 360° (in aplane) or a π (in the space) detection system: this allowsthe reconstruction of the tissue itself. A similar imagingtechnique is the Positron Emission Tomography (PET)in which, instead, a positron-electron pair annihilatesand emits two photons in opposite directions [5]; in thiscase the image is reconstructed via measurements of thecoincidence count rate. A plethora of imaging algorithmshas been proposed [6] in the past decades; here we willfocus on the Conditional Entropy Maximization (CEM),in which the image of a radioactive or contaminatedsource is reconstructed by means of the maximizationof the informational entropy.Entropy maximization is a general technique to solveproblems whose analytic solution is not possible orcomputationally complicated. This approach requiresa straightforward and unambiguous definition of theentropy [7, 8]. Entropy maximization overcomes issuesof other methods and furthermore allows to bring priorinformation into the calculation.The general idea behind the entropy maximization is todefine a space Y of the measured observable y ∈ Y (asingle quantity or a vector of quantities) and a parameterspace Λ of the observable λ ∈ Λ to be estimated. Wecan then define a probability density function P ( y ; λ ) over the parameter space that is in general difficultto maximize (i.e. we cannot find in general the mostlikely λ that realizes the observed y ). Therefore, weintroduce a fictitious space X and a mapping x → y ,and thus we move our problem of finding the parametervalue (or values) to the problem of maximizing anotherfunctional (in this case the conditional entropy) definedon a different space. The task is then to assign an entropyfunction that is compatible with the prior informationabout the system and the measurement process, and thenmaximize it with a constraint given by the outcome ofthe measurement.We will now show how this general approach can beapplied to the reconstruction of an image by measuringthe count of emitted photons (or coincidence counts).The measurement is performed by enclosing the samplein a chamber surrounded by detectors (in a plane orin the space) that can measure the number of emittedphotons in different directions (see Fig. 2). Themeasured count rates y j ( j = 1 , . . . , M with M totalnumber of detectors) can be modeled as a Poissonvariable and the sample as a matrix of N pixels (or cells),each with activity λ i . The goal is thus to properly modelthe emission and measurement process in order to definethe probability (or likelihood) function P ( y ; λ ) . Firstof all, we can define a matrix p ij of the ratios between the intensity detected by detector j and the decay rateat cell (i.e. pixel) i . In the following we will refer tothese quantities as efficiencies. The mean value of themeasured count is then given by µ j = t m N (cid:88) i =1 λ i p ij (3)with t m time of the measurement. If we assumethe measured counts to be Poisson distributed,then the conditional probability that a realization y = ( y , . . . , y M ) of the measurand is obtained given aset of λ i is the joint distribution of M Poisson variables P ( y | λ ) = M (cid:89) j =1 e − µ j µ y j j y j ! that takes the following form by inserting the (3) P ( y | λ ) = M (cid:89) j =1 exp( − t m (cid:80) Ni =1 λ i p ij )( t m (cid:80) Ni =1 λ i p ij ) y j y j ! . The conditional entropy of a measurement y ∈ Y givena distribution λ ∈ Λ of the parameter to be estimated isthen given by [9, 10, 11] H ( y | λ ) = − (cid:88) λ P ( λ ) (cid:34)(cid:88) y P ( y | λ ) log P ( y | λ ) (cid:35) (4)with P ( λ ) prior distribution function. The conditionalentropy (4) measures the information about theknowledge of our system, and the realization ˆ λ = (cid:16) ˆ λ , . . . , ˆ λ N (cid:17) that maximizes (4) is the one mostlikely to be realized. Following [9], the maximizationof H ( y | λ ) with respect to the parameter (or the set ofparameters) λ yields P (cid:16) y | ˆ λ (cid:17) (cid:20) λ i ∂P∂λ i log P ( y | λ )+ λ i t m P ( λ i ) (1 + log P ( y | λ )) M (cid:88) j =1 y j p ij (cid:80) Ni =1 λ i p ij − λ i t m (1 + log P ( y | λ )) M (cid:88) j =1 p ij λ i =ˆ λ i = 0 (5)where ˆ λ i are the values of the parameters that maximizethe conditional entropy. Equation (5) is verified when4 igure 2: Model of the measurement process and of the sample.In the top picture it is possible to see a section of the sample andof the LCM showing how the emission and detection is modeled inorder to properly define the conditional entropy and the detectionefficiency. In the bottom picture an illustrative 3D model of theLCM and of the measurement process. the second factor vanishes ( P ( y | λ ) (cid:54) = 0 ), leading to ˆ λ i = ˆ λ i P (cid:16) ˆ λ i (cid:17) (1 + log P ( y | λ i )) t m (cid:80) Mj =1 p ij × (cid:34) ∂P∂λ i (cid:12)(cid:12)(cid:12)(cid:12) λ i =ˆ λ i log P (cid:16) y | ˆ λ i (cid:17) + P (cid:16) ˆ λ i (cid:17) (cid:16) P (cid:16) y | ˆ λ i (cid:17)(cid:17) t m M (cid:88) j =1 y j p ij (cid:80) Ni =1 ˆ λ i p ij . (6)This equation has no closed solution, thus we need tosolve it numerically. In particular, given the special formof (6), an iterative solution can be achieved by settingthe right-hand-side as the estimate for the ( k + 1) -thiterative step of the left-hand-side: ˆ λ k +1 i = ˆ λ ki P (cid:16) ˆ λ ki (cid:17) (cid:0) P (cid:0) y | λ ki (cid:1)(cid:1) t m (cid:80) Mj =1 p ij (7) × (cid:34) ∂P∂λ i (cid:12)(cid:12)(cid:12)(cid:12) λ i =ˆ λ ki log P (cid:16) y | ˆ λ ki (cid:17) + P (cid:16) ˆ λ ki (cid:17) (cid:16) P (cid:16) y | ˆ λ ki (cid:17)(cid:17) t m M (cid:88) j =1 y j p ij (cid:80) Ni =1 ˆ λ ki p ij . Considering that P (cid:16) y | ˆ λ (cid:17) is a probability (and thusnon-negative) and that in a PET or SPECT systemit is much smaller than 1, we can then use theapproximation P (cid:0) y | λ ki (cid:1) ≈ log P (cid:0) y | λ ki (cid:1) . Thisleads to the formula for the estimate of the activity ofeach pixel of a waste package: ˆ λ k +1 i = ˆ λ ki P (cid:16) ˆ λ ki (cid:17) t m (cid:80) Mj =1 p ij (8) × ∂P∂λ i (cid:12)(cid:12)(cid:12)(cid:12) λ i =ˆ λ ki + P (cid:16) ˆ λ ki (cid:17) t m M (cid:88) j =1 y j p ij (cid:80) Ni =1 ˆ λ ki p ij . This means that we can estimate the activity and theposition of radioactive sources if we choose λ i to be theactivity of pixel i .As it is clear from (8), the role of the prior iscrucial. What one needs is a function that encompassesall the prior information about the system and themeasurement. In our case, we have only considered thecase of a uniform prior P ( λ ) = const. = 1 , meaningthat we assume to have no prior knowledge: this choicesimplifies a lot the formula (8) and reduces our algorithmto the Maximization of the Likelihood function [7, 8]: ˆ λ k +1 i = ˆ λ i (cid:80) Mj =1 p ij M (cid:88) j =1 y j p ij (cid:80) Ni =1 ˆ λ ki p ij .
3. Efficiencies
As it has been shown in section 2, a crucial role is playedby the efficiencies p ij , i.e. the probability that a decayhappening at cell (pixel) i emits a photon that is detectedby detector j . To determine the efficiencies, the followingproperties of the system are considered:• the sample (composition, geometry, density, nuclearcross section),• the detectors (type, position, geometry, tally),• the relative position of the pixel and the detectorand the angle of view,• the world around the sample and the detectors (thesample carrier, the air between the sample andthe detector, the external shielding material, andwhatever is not directly either the sample or thedetector),• the source (the isotopic composition, the emissionprobability and energy of the photons),Given the large number of parameters, this calculationis highly nontrivial and is practically impossibleanalytically: it is then necessary to perform Monte Carlosimulations. In the following, we will first describe howMonte Carlo N-particle simulations (MCNP) are usedto calculate the energy transported from source insidethe sample to one of the detectors. Then we will showhow, starting from the properties of the detectors (inour case plastic scintillation detectors), we can transformthat into an estimate of the photon count.The procedure to calculate the efficiencies p ij is then thefollowing:1. the sample and the detection are modeled in theframework of MCNP;2. a standard radioactive source is placed in each cellof the sample;3. a Monte Carlo simulation of the propagation ofgamma photons through the system is performed;4. the energy deposited on the detectors is convertedinto a count of detected photons;5. the efficiency is calculated as the ratio of detectedintensity and decay rate. MCNP® is a general-purpose, continuous-energy,generalized-geometry, time-dependent, Monte Carloradiation-transport code designed to track many particletypes over broad ranges of energies [14].Specific areas of application include, e.g. radiationprotection and dosimetry, radiation shielding,radiography, medical physics, nuclear criticality safety,detector design and analysis, nuclear oil well logging,accelerator target design, fission and fusion reactordesign, decontamination and decommissioning. Thecode treats an arbitrary three-dimensional configurationof materials in geometric cells bounded by first- andsecond-degree surfaces.The transport is based on point-wise cross-sectiondata. For photons, the code accounts for incoherentand coherent scattering, the possibility of fluorescentemission after photoelectric absorption, absorption inpair production with local emission of annihilationradiation, and bremsstrahlung.The code provides easy-to-use and very versatile generalsources and flexible tallies. In the specific case ofthe calculation of efficiencies for detectors in a LargeClearance Monitor the sources are modeled by a isotropicpoint source located in the center of each pixel. Theenergy distribution of the source photons reflects thedecay gamma energies for the respective radio nuclide.The code then transports the particles through thegeometry and a pulse-height tally which records theenergy deposited in the detector volume by each primaryparticle and its secondary particles. The spectrum of thedeposited energies is calculated on a very fine resolution(~3000 bins) between 0 eV and 3 keV.
SimPS (Simulated Plastic Scintillator) is a softwarefor plastic scintillation detectors to calculate detectionprobabilities of photons with specific energies whichinteract with the detector. Moreover, it takes intoaccount the emission probability, i.e. the number ofphotons emitted per decay of the radioactive source.In order to do so, a detector characterization andcalibration needs to be performed by estimating theenergy dependence of the detection probability. Thisprobability is mainly influenced by the discriminatorvoltage of the detector and can be modeled by means of acumulative distribution function of a Pareto-distribution(9) whose parameters x m and α are estimated by fittingcalibration measurements F ( x ) = (cid:40) − x αm x , x ≥ x m , x < x m . (9)6sing this detector characterization and the energydeposition spectrum calculated by MCNP for eachdetector-pixel pair, it is possible to estimate theefficiencies p ij used in the CEM algorithm.
4. Application of the CEM to clearancemeasurements
As it has been already pointed out, the analysis ofa potentially contaminated sample is subject to thecompliance with the international guidelines set by ISOand adopted by the national governments [15]. As ithas been reported in section 1.1, the idea behind thestandard is that not only the best estimate of the activityis needed to make a decision about a sample, but also itsdistribution has to be fully characterized by calculatingall the statistical characteristic values, i.e. the coverageinterval, the detection limit and the decision threshold.In order to achieve this for each measured sample, weneed an evaluation model to estimate the activity and itsuncertainty once we have performed the measurements asit has been already explained.It is then crucial to show if and how the method describedin the previous section is able to estimate activity valuesand their uncertainty. The norm prescribes to use analgebraic model for the answer to this question, and thusit might seem that the CEM method is not compliantwith the ISO 11929:2019. In fact, the CEM method is notonly an imaging technique, but also a way to estimate theactivity (and the respective uncertainty) of any portionof a sample, making it suitable to our purpose.Looking back at the general formulation of the CEMmethod, we can fix the general parameter λ to be theactivity of a cell of the sample and this allows theCEM to produce an estimate of the activity of eachcell: this is a good result, but still incomplete as theCEM itself provides no direct method to characterize thedistribution of the activity. We will address this problemin the following section 5; for now it is enough to knowthat the CEM algorithm can be opportunely integratedto fully characterize the activity distribution of each cellof a measured sample. This puts us in the condition to becompliant with the ISO 11929:2019, since we have both amodel of evaluation of the activity (resp. activities of thesingle cells or portions of the sample) and an estimateof the uncertainty (resp. uncertainty for each cell orportions of the sample). In the following section 4.2 wewill show our proposal to calculate all the characteristicvalues of the activity distribution for single cells and forwhatever portion of the sample, showing also that thisimproves a lot the limits of the coverage interval of theactivity of the whole sample. The general CEM method presented in section 2has a direct and easy application for PET imagingsystems. Despite the application to SPECT systemsis straightforward as well, there are some technicaldifficulties that make it tricky. In fact, in a PET system,given the nature of the emission process, the devicemeasures coincidence count rates and thus can directlymeasure the net count rate of pair annihilations. Onthe other hand, in a SPECT system (as it is in our case),the detectors measure also the background radiation (dueto natural radioactive sources) and thus one can eithersubtract the background from the gross count rate, orslightly change the formula (8) in order to use the grosscount rate instead. Using the gross count rate is of helpbecause in our case the count due to the sample canbe similar or even lower than the count of backgroundphotons. Moreover, using the gross count rate will makealso the calculation of the uncertainties easier, as it willbe shown in the following. The most straightforward wayof circumventing the problem is to substitute the meancount (3) with the corresponding expression of the grossnumber of emitted photons according to the evaluationmodel (1) of the ISO 11929:2019 µ j = t m N (cid:88) i =1 λ i p ij −→ t m (cid:32) x r ,j − x + N (cid:88) i =1 λ i p ij (cid:33) . (10)This modifies the formula (6) leading to ˆ λ i = ˆ λ i P (cid:16) ˆ λ i (cid:17) (1 + log P ( y | λ i )) t m (cid:80) Mj =1 p ij (cid:34) ∂P∂λ i (cid:12)(cid:12)(cid:12)(cid:12) λ i =ˆ λ i log P (cid:16) y | ˆ λ i (cid:17) ++ P (cid:16) ˆ λ i (cid:17) (cid:16) P (cid:16) y | ˆ λ i (cid:17)(cid:17) t m × M (cid:88) j =1 y j p ij r ,j x − x + (cid:80) Ni =1 ˆ λ i p ij where y = ( y , . . . , y N ) are now the gross count ratesof the N detectors. It has to be pointed out thatthe substitution (10) is also needed to calculate thecharacteristic statistical limits of the activity distributionas it will be shown in the next section. The recursiveformula (8) for the estimate of the activity becomes then ˆ λ k +1 i = ˆ λ ki P (cid:16) ˆ λ ki (cid:17) t m (cid:80) Mj =1 p ij (cid:34) ∂P∂λ i (cid:12)(cid:12)(cid:12)(cid:12) λ i =ˆ λ ki + P (cid:16) ˆ λ ki (cid:17) t m M (cid:88) j =1 y j p ij r ,j x − x + (cid:80) Ni =1 ˆ λ ki p ij . (11)7s a direct consequence, also the maximum likelihoodformula (i.e. with constant prior) is modified: ˆ λ k +1 i = ˆ λ ki (cid:80) Mj =1 p ij M (cid:88) j =1 y j p ij r ,j x − x + (cid:80) Ni =1 ˆ λ ki p ij .
5. Calculation of the characteristic limits
The first step to calculate the characteristic statisticallimits of the activity distribution is to estimate theuncertainty of the activities calculated by means of theCEM algorithm. As we have already pointed out, theCEM algorithm itself does not provide a direct wayto estimate the uncertainty, but this problem can beeasily circumvented by propagating the uncertainty ofthe measurement and of any other input parameter(according to their statistical distribution). Themeasured values y i are count rates of photons, and thusthey are Poisson distributed. The idea is then to generate P Poisson random sets of data y η = ( y η , . . . , y ηM ) , η = 1 , . . . , P (or repeating the measurement P times)and perform the CEM calculation for each of thisrandomly generated (or measured) set of data. Thiswill provide P different activity distributions λ η =( λ η , . . . , λ ηn ) and then we can use the obtained valuesto calculate the activity and the uncertainty for eachportion of the sample (in the most simple case for eachsingle cell) as the arithmetic mean and the standarddeviation of the values produced by the CEM algorithm: λ i = 1 P P (cid:88) η =1 λ ηi (12) u i ( λ i ) = 1 P − P (cid:88) η =1 ( λ ηi − λ i ) . (13)This procedure can be repeated if we have any otherknown source of uncertainty for our count rates. Inparticular, looking at (11), it is necessary to propagatealso the uncertainty of the background r , of theefficiencies p ij (that are normally distributed), of thefactor x and of the shielding x .Each cell will then have a distribution of estimatedactivities (cid:0) λ i , . . . , λ Pi (cid:1) that we can sort such that λ ηi ≤ λ η +1 i . Then we can assign to each value λ ηi a cumulative probability η/P and build a discretecumulative distribution F ( λ i | y η ) .We now have all the ingredients to fully characterize thestatistical activity distributions of each portion of thesample. This can be done by following the guidelines ofthe ISO 11929:2019, i.e. by estimating the uncertaintyfor specific values of the assumed true value. We needthen to calculate the outcome of a measurement ˜ y j when the activities λ i assume some special values ˜ λ i . This canbe done following (10): ˜ y j = x r ,j − x + N (cid:88) i =1 ˜ λ i p ij . Then the ˜ y j have to be varied around the mean andgiven as input to the CEM formula (11) to estimate thecharacteristic limits pixel by pixel ˆ λ k +1 i = ˆ λ ki P (cid:16) ˆ λ ki (cid:17) t m (cid:80) Mj =1 p ij (cid:34) ∂P∂λ i (cid:12)(cid:12)(cid:12)(cid:12) λ i =ˆ λ ki (14) + P (cid:16) ˆ λ ki (cid:17) t m M (cid:88) j =1 (cid:16) x r ,j − x + (cid:80) Ni =1 ˜ λ i p ij (cid:17) p ij r x − x + (cid:80) Ni =1 ˆ λ ki p ij . The mean and the standard deviation are then calculatedusing (12) and (13).Formula (14) can be seen as the bridge between the CEMalgorithm (i.e. the evaluation model for all the cells)and the statistical characteristic limits (for single cellsor blocks). In the following we will indeed address theproblem not only for single cells, but also for blocks, thatcan be reduced to single-cell-blocks. To this end we cansplit the sample itself into two parts B (the part ofinterest, of which we want to calculate a characteristiclimit) and B (the rest), such that we can split the sum (cid:80) Ni =1 λ i p ij in (14) into two parts as follows N (cid:88) i =1 λ i p ij = (cid:88) i ∈ B λ i p ij + (cid:88) i ∈ B λ i p ij where a block can consist even of one single cell or allthe cells. In order to keep the pattern of the estimatedactivity distribution fixed in B and B respectivelywhile varying the assumed true value we define thespecific efficiency for the two blocks p B ,j = (cid:80) i ∈ B λ i p ij (cid:80) i ∈ B λ i ; p B ,j = (cid:80) i ∈ B λ i p ij (cid:80) i ∈ B λ i . With these constrains the iterative CEM formula for theactivity in B can be written as follows ˆ λ k +1 B = ˆ λ kB P (cid:16) ˆ λ kB (cid:17) t m (cid:80) Mj =1 p B ,j × (cid:34) ∂P∂λ B (cid:12)(cid:12)(cid:12)(cid:12) λ B =ˆ λ kB + P (cid:16) ˆ λ kB (cid:17) t m × M (cid:88) j =1 (cid:16) x r ,j − x + ˜ λ B p B ,j + ˜ λ B p B ,j (cid:17) p B ,j r ,j x − x + ˜ λ B p B ,j + ˆ λ kB p B ,j . (15)8ssuming a true value for ˜ λ B , we can then samplethe input parameters and estimate the characteristicvalues. It has to be remarked that the numerator andthe denominator of (15) have to be sampled separatelyas they represent measurements of the backgroundthat happen respectively with (numerator) and without(denominator) the sample (i.e. the block of interest) inthe chamber. Following the prescriptions of the ISO 11929:2019, if wewant to calculate the decision threshold of the activityof a sample, we need to estimate the uncertainty ofa vanishing true value of the activity. The decisionthreshold can be defined as the minimum value thatcan be distinguished from the background with a givenconfidence and in this case the definition of backgrounditself is crucial. In order to calculate the decisionthreshold for a block of pixels (or eventually for a singlepixel), we need to consider the radiation emitted by therest of the sample as part of the background.In order to calculate the decision threshold for B , weneed to set its total activity ˜ λ B = 0 while the activitiesof the cells in B are left unchanged (i.e. the valueestimated by the CEM) and are taken into account aspart of the background radiation. This can be obtainedonly by setting all the ˜ λ i in B to zero (since activities arenon negative), leading to the following iterative formulafor the decision threshold λ ∗ B of block B that followsfrom (15) ˆ λ k +1 B = ˆ λ kB P (cid:16) ˆ λ kB (cid:17) t m (cid:80) Mj =1 p B ,j (cid:34) ∂P∂λ B (cid:12)(cid:12)(cid:12)(cid:12) λ B =ˆ λ kB + P (cid:16) ˆ λ kB (cid:17) t m M (cid:88) j =1 (cid:16) x r ,j − x + ˜ λ B p B ,j (cid:17) p B ,j r ,j x − x + ˜ λ B p B ,j + ˆ λ kB p B ,j . By sampling the input parameters of this formula aroundthe fixed true values ( r , x , x , ˜ λ B , p B ,j , p B ,j ), we getmany values of λ ηB that are distributed according to F (cid:16) ˜ λ B = 0 | y η (cid:17) as described in the previous paragraph.The (1 − α ) -quantile (typically with α = 0 . ) of thisdistribution is the decision thresholt λ ∗ B . The detection limit λ B of B is the smallest true valueof the activity in the block, for which the probability ofa false negative (as it is explained in section AppendixA) does not exceed the specified probability β (typically0.05).By applying the iterative formula (15), the detectionlimit λ B is the value of ˜ λ B for which the β -quantile of the distribution function F (cid:16) ˜ λ B | y η (cid:17) is equal to thedecision threshold. This value can to be calculatediteratively by applying root-finding methods. The limits of the coverage interval ( λ (cid:47)B , λ (cid:46)B ) are definedin such a way that the coverage interval contains the truevalue of the measurand within the specified probability − γ (typically with γ = 0 . ). In the following weuse the definition of the probabilistic symmetric coverageinterval which defines the limits of the coverage intervalas the γ/ -quantile ( λ (cid:47)B ) and (1 − γ/ -quantile ( λ (cid:46)B )of the distribution function F ( λ B | y η ) .
6. Examples
In the following we will show the accuracy of theproposed method by simulating measurements withMCNP and reconstructing the sources by means ofthe CEM method. The two examples that follow willsummarize all the things that have been discussed sofar and will show that the method presented is safe andreliable enough to handle waste packages from nuclearfacilities.The first simulated sample is an iron block of size 120x 80 x 23.5 cm (x, y and z dimension) and density1.681 g/cm carried by a steel sample-carrier-box andan aluminum pallet. A point source of Co of activity5 kBq is placed right in the middle of the sample asit is shown in Fig. 3 while the CEM-reconstructeddistribution is in Fig. 4. For the reconstruction, wehave sampled the input parameters 10000 times aroundthe measured (or in this case simulated) values accordingto their statistical distribution. For each sampled setof data, the recursive CEM formula (11) has been usedwith k = 10 to guarantee that the calculations havesatisfactorily converged. The activity for each pixel isthen calculated as the arithmetic mean of all the 10000runs of the algorithm. The characteristic values havebeen calculated according to the prescriptions layed outin the previous section. For the whole sample, the totalactivity is λ tot = 4890 Bq with a standard deviation of σ tot = 116 Bq . The lower and upper limit of the 95%symmetric coverage interval interval are λ (cid:47)tot = 4695 Bq and λ (cid:46)tot = 5077 Bq . The decision threshold is λ (cid:63)tot =60 Bq and the detection limit λ tot = 124 Bq . The integralefficiency of the whole detection system is 21.2% for thatactivity distribution. The characteristic values of thissample are summarized in Table 1.The second simulated sample consists of the same ironblock, with two Co point sources of 5 kBq each, locatedat the top "right” and at the bottom "left” (from theview of the experimenter) corners (see Fig. 5). In Fig.9 igure 3: Simulated sample n.1. Uniform iron block of density1.681 g/cm with one simulated Co source in the middle ofactivity 5 kBq . The three plots represent the three layers of thesample from bottom to top. Figure 4: Simulated sample n.1. This figure shows the CEM-reconstructed activity distribution of the sample shown in Fig. 3 Table 1: Table with all the characteristic values of the simulatedsample n. 1 (activity distribution of Fig. 3)
6, the reconstructed activity distribution is shown. Asbefore the same number of data sets (10000) and CEMiteration steps ( k = 10 ) were used for the calculationand the values shown are the arithmetic means for eachpixel. For the whole sample, the total estimated activityis λ tot = 9954 Bq with a standard deviation of σ tot =163 Bq . The lower and upper limit of the 95% symmetriccoverage interval interval are λ (cid:47)tot = 9706 Bq and λ (cid:46)tot =10241 Bq . The decision threshold is λ (cid:63)tot = 42 Bq and thedetection limit λ tot = 90 Bq . The integral efficiency overall detectors was calculated to be 31.6% for that activitydistribution. The characteristic values of this sample aresummarized in Table 2.The sample was split into two blocks (as in Fig. 7where they are marked in red and green) which wereanalyzed separately. The activity in Block 1 (red) islocated in the bottom of the sample whereas the activityin Block 2 (green) is located in the top. Due to thefact that the lower source is better shielded (conveyorbelt, aluminum pallet, sample carrier box), the blockcontaining this source will clearly have a lower detectionefficiency: indeed it has a total efficiency of 27.1% whilethe other block has an efficiency of 36.4%. Due to thedifference in efficiency for block 1 and 2, the decisionthresholds were estimated to be respectively λ (cid:63) = 71 Bq and λ (cid:63) = 64 Bq . These two values are larger thanthe decision threshold for the whole sample due to thefact that, when calculating the characteristic limits for ablock, the radiation generated by the other one needs tobe considered as part of the background. The detectionlimits are λ = 163 Bq and λ = 150 Bq . The meanactivity in block 1 was calculated to be λ = 5124 Bq with a standard deviation of σ = 284 Bq , whereas the forblock 2 a mean of λ = 4830 Bq and standard deviationof σ = 241 Bq was estimated.The coverage intervals andall other values are summarized in Table 2.A convergence study of the CEM has also beenperformed. If we consider the ( k + 1) -th step of theiterative CEM calculation of formula (11), we can definethe deviation of cell i from the previous iterative step Λ ki = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) λ k +1 i − λ ki λ ki (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . Figure 5: Simulated sample n.2. Uniform iron block of density1.681 g/cm with two simulated Co sources of activity 5 kBq each. The three plots represent the three layers of the sample frombottom to top. igure 6: Simulated sample n.2. This figure shows the CEM-reconstructed activity distribution of the sample shown in Fig. 5 Figure 7: Splitting of sample n.2. Each of the three layers are splitequally. The two blocks (red and green) are such that they haveapproximately the same number of cells Total Block 1 Block 2Activity (Bq) 9954 5124 4830Uncertainty (Bq) 163 284 241Lower limit of coverage (Bq) 9706 4743 4383Upper limit of coverage (Bq) 10241 5668 5159Decision threshold (Bq) 42 71 64Detection limit (Bq) 90 163 150Efficiency (%) 31.6 27.1 36.4
Table 2: Table with all the characteristic values of the simulatedsample n. 2 (activity distribution of Fig. 5) igure 8: Plot in logarithmic scale of the average deviation as afunction of the iterative step k . The first picture on top shows theplot for a single point source in the middle of the sample (as inFig. 3); the second pictures shows the plot for two point sourcesplaced at opposite corners of the sample as in Fig. 5 If we average this quantity over all the cells, we have anindicator Σ k = 1 N N (cid:88) i =1 Λ ki of how much the algorithm has converged. In Fig. 8 wehave plotted Σ k for the two cases of a point source in themiddle (as in Fig. 3) and of two point sources in oppositecorners (as in Fig. 5). Conclusion
In this paper we have presented an innovative approachto clearance measurements of waste packages fromnuclear facilities. The new procedure, borrowed frommedical physics, does not adopt algebraic evaluationmodels for the estimate of the activity and goes onestep further: the maximization of the conditional entropyfunctional allows the estimate of the activity distribution, cell by cell (as it happens for the image reconstruction).The presented CEM method can be also adapted to theinternational ISO 11929 standards by properly definingthe characteristic limits for the measured activities.Moreover, it allows to lower the characteristic limitsfor the free release of the contaminated waste, e.g. bysplitting the sample in two or more parts: if a non-releasable sample containing radioactive hotspots can besplit into portions whose hotspots do not exceed the legalthresholds separately, then we have achieved the goal ofreducing the amount of waste that needs to be stored inspecial repositories with consequent saving of money andresources.In this first study, we have then shown that theCEM method reconstructs faithfully the radioactivedistribution of a simulated metal block with respectivelyone and two point sources. The calculation of thecharacteristic limits of the whole block and of partsof the block shows that we can significantly increasethe amount of releasable waste e.g. by splitting thesample. It was shown that the calculated activitymatches very well the original activity, easily withinthe coverage interval. Furthermore, the localizationand magnitude of the activity was found to be ingood agreement with the activity placed in the sampleboth for a single hotspot and for two hotspots. Thedescribed procedure to calculate characteristic limitsprovided reasonable estimates of the decision threshold,of the detection limit as well as of the coverage intervalsfor both the whole sample and for parts of it. It wasshown that the characteristic limits strongly depend onthe activity distribution, making clear the advantage ofthis approach. Also we have found that the recursivealgorithm requires a finite number of iterations toconverge within a given accuracy.Further studies will include a quantitative analysis of theagreement of the reconstructed image with the originalsample, as well as a general rule to fix the number ofiterations needed to converge significantly. Moreover wewill report the results and conduct further analysis ofreal samples. Another future extension of the methodwill include techniques to improve the spatial resolution.
Acknowledgment
We are thankful to prof. Rolf Michel for supporting ourwork and encouraging the writing of the present paper.
References [1] G. Amoyal, V. Schoepff, F. Carrel, M. Michel,N. Blanc de Lanaute, J. C. Angélique, Developmentof a hybrid gamma camera based on timepix3 fornuclear industry applications, Nucl. Inst. and Meth.13n Phys. Res. Sec. A In Press, Journal Pre-proof(2020).[2] M. Duerr, K. Krycki, B. Hansmann, T. Hansmann,A. Havenith, M. Fritzsche, D. Pasler, T. Hartmann,Advanced sectorial gamma scanning for theradiological characterization of radioactive wastepackages, Atw. Internat. Zeitsch. fuer Kernen. 64 (3)(2019) 160.[3] Evaluation of measurement data - guide to theexpression of uncertainty in measurement. s.l.,BIPM, IEC, IFCC, ILAC, ISO, IUPAC, IUPAP andOIML JCGM 100 (2008).[4] Y.-M. Zhu, S. M. Cochoff, Likelihood maximizationapproach to image registration, IEEE Trans. Im.Proc. 11 (12) (2002) p. 1417. doi:10.1109/TIP.2002.806240 .[5] A. C. Kak, M. Slaney, Principles of ComputerizedTomographic Imaging, 2001.[6] J. B. A. Maintz, M. A. Viergever, A surveyof medical image registration, Med. Image Anal.2 (1) (1998) p. 1. doi:10.1016/S1361-8415(01)80026-8 .[7] B. R. Frieden, Restoring with maximum likelihoodand maximum entropy, J. Opt. Soc. Am. 62 (1972)p. 511. doi:10.1364/JOSA.62.000511 .[8] B. R. Frieden, D. C. Wells, Restoring with maximumentropy: Poisson sources and backgrounds, J.Opt.Soc. Am. 68 (1) (1978) p. 93. doi:10.1364/JOSA.68.000093 .[9] P. P. Mondal, K. Rajan, Image reconstruction byconditional entropy maximisation for pet system,IEEE Proc. Vis. Im. Sign. Proc. 151 (5) (2004) p.345. doi:10.1049/ip-vis:20040717 .[10] C. L. Byrne, Likelihood maximization for list-modeemission tomographic image reconstruction, IEEETrans. Med. Imaging 20 (10) (2001) 1084. doi:10.1109/42.959305 .[11] C. L. Byrne, Iterative image reconstructionalgorithms based on crossentropy minimization,IEEE Trans. Image Process 2 (1) (1993) 96. doi:10.1109/83.210869 .[12] H. Zhu, H. Shu, J. Zhou, X. Dai, L. Luo,Conditional entropy maximization for pet imagereconstruction using adaptive mesh model, Comp.Med. Im. and Gr. 31 (2007) p. 166. doi:10.1016/j.compmedimag.2007.01.001 .[13] A. Lopez, R. Molina, A. K. Katsaggelos, Spectimage reconstruction using compound models, 2001IEEE Int. Conf. Ac. Proc. (2001) p. 1909 doi:10.1109/ICASSP.2001.941318 . [14] C. J. Werner, J. S. Bull, C. J. Solomon, F. B. Brown,G. W. McKinney, M. E. Rising, D. A. Dixon, R. L.Martz, H. G. Hughes, L. J. Cox, A. J. Zukaitis, J. C.Armstrong, R. A. Forster, L. Casswell, Mcnp version6.2 release notes (2 2018). doi:10.2172/1419730 .[15] Iso 11929:2019: Determination of characteristiclimits (decision threshold, detection limit, and limitsof the confidence interval) for measurements ofionizing radiation - fundamentals and applications,International Organization for Standardization,Geneva (2019).[16] K. Weise, G. Kanisch, R. Michel, M. Schläger,D. Schrammel, M. Täschner, Characteristic valuesin measurement of ionizing radiation - material for acritical discussion on fundamentals and alternatives,Fachver. Str. AK Sigma 167 (2013).
Appendix A. Definitions of the characteristicvalues
The formulae in section 1.1 refer to the calculation ofcharacteristic values in case of a specific evaluation modeland under the assumption of a symmetric probabilitydistribution of the measurands. There are some casesthough in which these assumptions don’t hold or it iseven analytically impossible to determine the probabilitydistribution function and thus only a Monte Carlosampling of the measurand is possible [16]. If we definethe probability distribution f Y ( y | a ) that the measurand Y takes the value y given the set of conditions a and a similar distribution f Y ( y | a , y ≥ of a positivemeasurand, we can build the cumulative distributionfunction F Y ( y | a ) = P ( Y < y | a ) = (cid:90) y −∞ f Y ( η | a ) dη and also I = (cid:90) ∞ f Y ( η | a ; y ≥ dη. With these definitions, we can generalize the expressionsfor all the characteristic limits of the activity distributiongiven in section 1.1:• the mean or expectation value ˆ y is the best estimatefor the measurand ˆ y = E ( Y | a , y ≥
0) = I /I ; I = (cid:90) ∞ f Y ( y | a ) dy • the standard uncertainty u (ˆ y ) of the measurandassociated with the best estimate ˆ yu (ˆ y ) = (cid:114) E (cid:16) ( Y − ˆ y ) | a , y ≥ (cid:17) = (cid:114) I I − ˆ y I = (cid:90) ∞ y f Y ( y | a ) dy