Metrological approach of γ-emitting radionuclides identification at low statistics: application of sparse spectral unmixing to scintillation detectors

The problem addressed in the following is a full-spectrum analysis of low-statistics γ-spectra based on spectral unmixing in order to estimate counting attached to each γ-emitting radionuclide contributing to the measurement. In γ-spectrometry, an observed spectrum s = [s₁, ..., s_M ] is composed of M channels corresponding to intervals of deposited energy by interacting photons in a detector. Hence, counting in each channel s_m for 1 ⩽ m ⩽ M represents the number of events at a given energy interval. The probability to obtain the counting s_m in channel m is given by the following Poisson distribution:

$$\begin{equation}\mathcal{P}\left({s}_{m}\vert {\lambda }_{m}\right)=\frac{{\lambda }_{m}^{{s}_{m}}}{{s}_{m}!}\enspace {\mathrm{e}}^{-{\lambda }_{m}}\end{equation} \tag{ 1 }$$

where λ_m is the average value of s_m . Let $\mathbf{\Phi }\in {{\mathbb{R}}_{+}}^{M{\times}N}$ a matrix, called hereafter dictionary, containing on its columns the normalized spectral signatures corresponding to the detector response of N γ-emitting radionuclides. It is noteworthy that in the present work, natural Bkg is considered as a γ-emitting radionuclide to be identified. The Poisson parameters λ_m are linear combinations of the individual spectral signatures in channel m related to the γ-emitting radionuclides (with index n)

$$\begin{equation}{\lambda }_{m}=\sum _{n=1}^{N}{\phi }_{m,n}{a}_{n}\quad \forall \enspace m=1,\dots ,M\end{equation} \tag{ 2 }$$

where the expected counting contained in the vector a = [a₁, ..., a_N ] are related to each spectral signature in the dictionary. When the M channels are considered for full-spectrum analysis, the statistical independence of the stochastic measurements yields the following expression for the likelihood:

$$\begin{equation}\mathcal{P}\left(\boldsymbol{s}\vert \boldsymbol{a}\right)=\prod _{m=1}^{M}\frac{{\left({\sum }_{n=1}^{N}{\phi }_{m,n}{a}_{n}\right)}^{{s}_{m}}}{{s}_{m}!}\enspace {\mathrm{e}}^{-\sum _{n=1}^{N}{\phi }_{m,n}{a}_{n}}.\end{equation} \tag{ 3 }$$

The above likelihood represents the conditional probability to measure the spectrum s given the counting vector a . Computing the maximum likelihood estimator $\hat{\boldsymbol{a}}$ of the vector a is obtained by minimizing the negative log-likelihood (cost function) expressed as follows:

$$\begin{equation}\mathcal{L}\left(\boldsymbol{a}\right)=\sum _{m=1}^{M}\left(\sum _{n=1}^{N}{\phi }_{m,n}{a}_{n}\right)-{s}_{m}\enspace \mathrm{log}\left(\sum _{n=1}^{N}{\phi }_{m,n}{a}_{n}\right).\end{equation} \tag{ 4 }$$

The estimated counting are obtained by minimizing this cost function, which is usually written as:

$$\begin{equation}\hat{\boldsymbol{a}}=\text{argmin}{}_{\hspace{1pt}\boldsymbol{a}}\enspace \mathcal{L}\left(\boldsymbol{a}\right).\end{equation} \tag{ 5 }$$

Assuming that $\forall \enspace m=1,\dots ,M;\quad \sum _{n=1}^{N}{\phi }_{m,n}{a}_{n}\ne 0$, the right expression of the cost function $\mathcal{L}\left(\boldsymbol{a}\right)$ is differentiable and its gradient with respect to counting a_r (for each γ-emitting radionuclide r to be identified) is defined as follows:

$$\begin{equation}{\nabla }_{{a}_{r}}\mathcal{L}\left(\boldsymbol{a}\right)=\sum _{m=1}^{M}{\phi }_{m,r}-\frac{{s}_{m}}{{\sum }_{n=1}^{N}{\phi }_{m,n}{a}_{n}}{\phi }_{m,r}.\end{equation} \tag{ 6 }$$

A classical procedure to minimize a cost function is given by the gradient descent method, which is a first-order iterative optimization technique. In short, this method consists in decreasing the cost function according to a direction given by its gradient. The update of the unknown parameters to be estimated according to a gradient-descent scheme is:

$$\begin{equation}\forall \enspace r=1,\dots ,N;\quad {\hat{a}}_{r}^{\left(\ell +1\right)}={\hat{a}}_{r}^{\left(\ell \right)}-{\eta }_{r}^{\left(\ell \right)}\cdot {\nabla }_{{a}_{r}}\mathcal{L}\left({\hat{\boldsymbol{a}}}^{\left(\ell \right)}\right)\end{equation} \tag{ 7 }$$

where ${\nabla }_{{a}_{r}}\mathcal{L}\left({\hat{\boldsymbol{a}}}^{\left(\ell \right)}\right)$ and ${\eta }_{r}^{\left(\ell \right)}$ are respectively the gradient of the cost function and the gradient step at the iteration ℓ for the rth radionuclide. In order to obtain an algorithm that can be implemented in embedded digital systems, the solution proposed in [4] defines the gradient step at each iteration as follows:

$$\begin{equation}\forall \enspace r=1,\dots ,N;\quad {\eta }_{r}^{\left(\ell \right)}=\frac{{\hat{a}}_{r}^{\left(\ell \right)}}{{\sum }_{m=1}^{M}{\phi }_{m,r}}.\end{equation} \tag{ 8 }$$

By including the above gradient step in expression (7) and considering that the spectral signatures are normalized (i.e. $\sum _{m=1}^{M}{\phi }_{m,r}=1$), the following multiplicative update is obtained for the estimation of the counting a :

$$\begin{equation}\forall \enspace r=1,\dots ,N;\quad {\hat{a}}_{r}^{\left(\ell +1\right)}={\hat{a}}_{r}^{\left(\ell \right)}\sum _{m=1}^{M}\frac{{s}_{m}}{{\sum }_{n=1}^{N}{\phi }_{m,n}{\hat{a}}_{n}^{\left(\ell \right)}}{\phi }_{m,r}.\end{equation} \tag{ 9 }$$

The first-order optimality condition is obtained when $\nabla \mathcal{L}\left(\hat{\boldsymbol{a}}\right)=0$, which is equivalent to:

$$\begin{equation}\forall \enspace r=1,\dots ,N;\quad \sum _{m=1}^{M}\frac{{s}_{m}}{{\sum }_{n=1}^{N}{\phi }_{m,n}{\hat{a}}_{n}^{\left(\ell \right)}}{\phi }_{m,r}=1.\end{equation} \tag{ 10 }$$

It can be pointed out that the multiplicative update rule leads to estimated counting in expression (9) that cannot be negative when the mixing weights ${\hat{\boldsymbol{a}}}^{\left(\ell \right)}$ are non-negative. Thus, by initializing counting with non-negative values, the multiplicative update rule implicitly includes non-negativity constraint avoiding the assessment of non-physical negative values. The multiplicative update algorithm, called hereafter non-negative Poisson unmixing (NNPU), allows a full-spectrum analysis without the use of region-of-interest related to each γ-emitting radionuclides to be identified. The pseudo-code of the NNPU algorithm is given in [5].

The problem of decision-making was first investigated in a previous work [5] in the case of anomaly detection of γ-emitting radionuclides in natural Bkg radiation when using a 3'' × 3'' NaI(Tl) detector. At low statistics, the problem of overfitting has a direct impact on false positive alarms and thus on the robustness of decision-making. This issue was addressed by computing DT when only natural Bkg is present by using MC calculations [5] for several γ-emitting radionuclides generally used for RPM testing (⁵⁷Co, ⁶⁰Co, ¹³³Ba, ¹³⁷Cs, ²⁴¹Am) [7] according to a false positive rate equal to 0.1% (probability of type I error) [6]. This first study figures out that this solution was not optimal to account for the significant variation of DT values when γ-emitting radionuclides are mixed with natural Bkg. In order to overcome this overfitting problem, which is also sensitive to overlapping between spectral signatures when using scintillation detectors, decision-making is implemented in the following using a sparsity constraint combined to the NNPU algorithm in order to obtain an optimized solution of identified radionuclides.

Sparsity was previously studied for spectral unmixing of γ-spectra obtained with HPGe detectors [8] using the Chambolle–Pock algorithm [9] for the unmixing technique. The aim of the sparsity constraint is to determine a minimum subset of identified radionuclides from the initial full dictionary Φ that corresponds to the best fit of a measured γ-spectrum. The optimization problem for spectral unmixing expressed in (5) combined with the sparsity constraint can be rewritten as follows:

$$\begin{equation}\hat{\boldsymbol{a}}=\text{argmin}{}_{\hspace{1pt}\boldsymbol{a}}\enspace \mathcal{L}\left(\boldsymbol{a}\right)\enspace \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\enspace \mathrm{t}\mathrm{o}\enspace {\Vert}\boldsymbol{a}\vert {\vert }_{0}{\leqslant}K\end{equation} \tag{ 11 }$$

where || a ||₀ is the number of non-zero elements in the vector a and K is the number of identified radionuclides that are actually present in the measured spectrum. In order to solve this sparsity constraint, the $\mathcal{P}$-OMP algorithm (Poisson-based orthogonal matching pursuit) described in [8] was implemented to determine the most likely spectral signatures that maximize the likelihood function (3).

The principle of the $\mathcal{P}$-OMP algorithm consists in a sequential selection of spectral signatures to be included in an optimized downsized dictionary. To that end, the $\mathcal{P}$-OMP algorithm first selects a new spectral signature (from the initial full dictionary Φ) that maximizes the likelihood computed with the NNPU algorithm. An hypothesis test is subsequently performed using a likelihood ratio test (LRT) usually applied for model selection [10]. For that purpose, the statistical deviance between the models obtained with two dictionaries Φ₁ and Φ₀ (Φ₁ corresponds to Φ₀ plus the selected spectrum signature) is defined by the difference of their respective negative log-likelihoods $D=-2\left({L}_{1}-{L}_{0}\right)$ computed with the NNPU algorithm. The selection process of a new spectral signature is stopped (the null hypothesis corresponding to Φ₀ is accepted) when the statistical deviance D is lower than a DT value depending on the expected false positive ratio (probability of type I error). As a result, Φ₀ corresponds to the optimized solution of identified radionuclides with their respective estimated counting given by the NNPU algorithm. On the contrary, the selection process goes on when the null hypothesis corresponding to Φ₀ is rejected according to an expected significance level (the statistical deviance D is greater than the DT value). The $\mathcal{P}$-OMP algorithm based on the sparsity constraint is described in table 1 in the case of anomaly detection of γ-emitting radionuclides in natural Bkg radiation.

Table 1. Pseudocode of the $\mathcal{P}$-OMP algorithm in the case of anomaly detection of γ-emitting radionuclides in natural Bkg radiation.

Input

Spectral signatures dictionary Φ in which the first one represents the natural Bkg

Observed spectrum: $\boldsymbol{s}=\left[{s}_{1},\dots ,{s}_{m}\right]$

Decision threshold DT depending on the expected false positive rate

Initialization

Identified radionuclide index I0 = [1] (Bkg)

Radionuclide indices I = [2, ..., N] for sparsity selection

Estimate $\hat{\boldsymbol{a}}$ with Φ[I0] using NNPU (Φ0 = Φ[I0])

Compute L0 using the cost function (4)

Set Flag = 1

While Flag = 1 and I ≠ ∅ do

For i ∈ I do

Build the model to test: Itest ← I0 ∪ i

Estimate ${{\hat{\boldsymbol{a}}}^{\left(\boldsymbol{i}\right)}}_{\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}}$ with Φ[Itest] using NNPU

Compute ${{L}^{\left(i\right)}}_{\text{test}}$ using the cost function (4)

End for

Find j such as ${L}_{1}={{L}^{\left(j\right)}}_{\text{test}}$ be the lowest ${{L}^{\left(i\right)}}_{\text{test}}$ ∀ i

If $-2\left({L}_{1}-{L}_{0}\right){ >}\mathrm{D}\mathrm{T}$ do

Update present radionuclide indices: I0 ← I0  ∪ j

Update radionuclide indices to be tested: I ← I\j

Set L0 = L1 and $\hat{\boldsymbol{a}}={{\hat{\boldsymbol{a}}}^{\left(\boldsymbol{j}\right)}}_{\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}}$ (Φ0 = Φ[I0])

Else

Set Flag = 0 (Φ0 not rejected)

End if

End while

Output

Downsized dictionary of identified radionuclides: Φ0

Estimated counting: $\hat{\boldsymbol{a}}\left[{\mathbf{\Phi }}_{0}\right]$

For nested models and under null hypothesis (i.e. Φ₀ is the true model), the LRT statistics asymptotically follows a Chi-square distribution [11]. Thus, the DT value in the $\mathcal{P}$-OMP algorithm can be determined from a Chi-square distribution with one degree of freedom (${\chi }^{2}\left(1\right)$) and a significance level (probability of rejecting the null hypothesis) corresponding to the expected false positive rate. This assumption is investigated in the following sections in the case of the application of the $\mathcal{P}$-OMP algorithm for decision-making when using scintillation detectors.

In order to provide a full metrological tool, the standard deviations related to the estimated counting a can be assessed by means of the Fisher information matrix $\boldsymbol{I}\left(\boldsymbol{a}\right)$ [12]. The elements of this matrix are defined as

$$\begin{align}\hfill I{\left(\boldsymbol{a}\right)}_{i,j}& =E\left[\frac{{\partial }^{2}\mathcal{L}\left(\boldsymbol{a}\right)}{\partial {a}_{i}\partial {a}_{j}}\right]\hfill \\ \hfill & =\sum _{m=1}^{M}E\left[\frac{{s}_{m}}{{\left({\sum }_{n=1}^{N}{\phi }_{m,n}{a}_{n}\right)}^{2}}{\phi }_{m,i}{\phi }_{m,j}\right]\hfill \end{align} \tag{ 12 }$$

where E stands for the expectation. Considering that ϕ_m,i , ϕ_m,j and $\sum _{n=1}^{N}{\phi }_{m,n}{a}_{n}$ are deterministic, it yields:

$$\begin{equation}E\left[\frac{{\partial }^{2}\mathcal{L}\left(\boldsymbol{a}\right)}{\partial {a}_{i}\partial {a}_{j}}\right]=\sum _{m=1}^{M}E\left[{s}_{m}\right]\frac{{\phi }_{m,i}{\phi }_{m,j}}{{\left({\sum }_{n=1}^{N}{\phi }_{m,n}{a}_{n}\right)}^{2}}.\end{equation} \tag{ 13 }$$

Finally, since $E\left[{s}_{m}\right]=\sum _{n=1}^{N}{\phi }_{m,n}{a}_{n}$, the elements of the Fisher information matrix have the following expression:

$$\begin{equation}I{\left(\boldsymbol{a}\right)}_{i,j}=\sum _{m=1}^{M}\frac{{\phi }_{m,i}{\phi }_{m,j}}{{\sum }_{n=1}^{N}{\phi }_{m,n}{a}_{n}}.\end{equation} \tag{ 14 }$$

According to the theory of MLE, the asymptotic distribution of the estimated vector $\hat{\boldsymbol{a}}$ is Gaussian and its covariance matrix $\mathrm{var}\left(\hat{\boldsymbol{a}}\right)$ is given by the inverse Fisher information matrix as follows:

$$\begin{equation}\mathrm{var}\left(\hat{\boldsymbol{a}}\right)\approx {\left[\boldsymbol{I}\left(\hat{\boldsymbol{a}}\right)\right]}^{-1}.\end{equation} \tag{ 15 }$$

The standard deviations (or uncertainties) of the estimated counting $\hat{\boldsymbol{a}}$ can thus be assessed by calculating the square root of the diagonal elements of the inverse Fisher information matrix. The Gaussian approximation validity of the counting estimator will be investigated in the following.

In this section, the DT value for the statistical deviance test was investigated according to the dictionary size when the measured spectrum contains only natural Bkg. For that purpose, the statistical deviance distribution obtained in the $\mathcal{P}$-OMP algorithm was simulated according to increasing dictionary sizes by successively adding a different γ-emitter spectral signature (see table 2). For each configuration, the DT value for a significance level of 0.1% was determined from MC calculations of 5 × 10⁵ Bkg simulations and for two mean counting defined to 1250 and 10 000.

The results displayed in table 2 are different from the DT value of 10.83 given by the ${\chi }^{2}\left(1\right)$ distribution for a significance level of 0.1%. This DT value is only valid for a dictionary size of 3 spectral signatures (natural Bkg, ²⁴¹Am and ¹³³Ba in table 2). In addition, it has to be noted that the DT value obtained for a dictionary size equal to 2 corresponds to a significance level equal to 0.2% from the ${\chi }^{2}\left(1\right)$ distribution. This result can be interpreted as a consequence of the non-negativity constraint. Indeed, as shown in figure 1 that displays the distribution of counting estimation for ²⁴¹Am with the NNPU algorithm in the case of a dictionary size equal to 2 (natural Bkg and ²⁴¹Am), half of the whole simulation results is close to zero. Without the non-negativity constraint, the distribution of estimated counting of ²⁴¹Am, which is not present in simulated spectra, would be centred on zero. In that case, the marginal values leading to a rejection of the null hypothesis with the LRT, would be equally composed of positive and negative values. Thus, because the non-negativity constraint forces negative counting to be close to zero, the DT value determined from the ${\chi }^{2}\left(1\right)$ distribution yields a false positive rate two times lower than expected. In the case of an initial dictionary size equal to 2, the DT value in the $\mathcal{P}$-OMP algorithm can be defined using a significance level twice the expected false positive rate when using the ${\chi }^{2}\left(1\right)$ distribution.

Zoom In Zoom Out Reset image size

In the case of initial dictionary sizes greater than 2, several γ-emitting radionuclides are available for the sparsity selection in the $\mathcal{P}$-OMP algorithm. As already mentioned, the chosen spectral signature for the statistical deviance test corresponds to the lowest value of the cost function. Based on the result in table 2, an empirical relation was defined between the DT value for an expected false positive rate in the $\mathcal{P}$-OMP algorithm and the significance level given by the ${\chi }^{2}\left(1\right)$ distribution. This relation is an extension of the case of an initial dictionary of two spectral signatures taking into account the number of spectral signatures available for the sparsity selection procedure. In order to match the results in table 2, the DT values in the $\mathcal{P}$-OMP algorithm can be obtained using the ${\chi }^{2}\left(1\right)$ distribution based on a significance level equal to twice the expected false positive rate divided by the number of radionuclides available for the selection procedure.

A comparison is presented in figure 2 between the empirical relation and the simulated DT values in table 2. In the case of a Bkg mean counting equal to 10⁴, both results are consistent. For a Bkg mean counting of 1250, the empirical relation based on the ${\chi }^{2}\left(1\right)$ distribution slightly overestimates the DT values for a number of available radionuclides for sparsity selection greater than 6 in the $\mathcal{P}$-OMP algorithm. This result could be due to a problem of convergence of LRT statistics to a Chi-square distribution at low statistics.

Zoom In Zoom Out Reset image size

The above results confirm that the DT value in the $\mathcal{P}$-OMP algorithm cannot be directly determined from the ${\chi }^{2}\left(1\right)$ distribution. In order to approximately maintain a constant false positive rate, it is mandatory to consider the increase of the DT value with the number of available radionuclides for the sparsity selection in the $\mathcal{P}$-OMP algorithm. This empirical adjustment of the significance level is similar to the Bonferroni correction applied in the case of multiple hypothesis testing [13]. For the sparsity constraint, this correction is performed to each individual LRT to ensure the expected significance level (i.e. 0.1% in the present study) for the full size of the dictionary.

In the following sections, the investigation of the DT calculation based on the Bonferroni correction is extended to the case of mixtures of γ-emitting radionuclides with natural Bkg.

3.2.1. Application of the sparsity constraint in the case of simple mixtures

3.2.2. Application of the sparsity constraint in the case of complex mixtures

The evolution of the false positive rate in the case of more complex mixtures was studied by successively adding a radionuclide to be identified in simulated spectra. For that purpose, γ-emitters used for RPM testing were included according to the following order with a mean counting equal to 1000: ²⁴¹Am, ¹³³Ba, ¹³⁷Cs, ⁶⁰Co and ⁵⁷Co. The evolution of false positive rates (2 × 10⁵ MC calculations) obtained according to each mixture configuration is displayed in figure 3. These results show that the false alarm rate remains robust for each configuration up to the following mixture: natural Bkg, ²⁴¹Am, ¹³³Ba, ¹³⁷Cs and ⁶⁰Co. The last case corresponding to the addition of ⁵⁷Co to the previous mixture (five radionuclides to be identified) yields an important increase of false alarms due to systematic identification of ¹⁵²Eu (see figure 3), which is actually not present in simulated spectra.

Zoom In Zoom Out Reset image size

An analysis of the order of the radionuclide sparsity selection in the $\mathcal{P}$-OMP algorithm revealed that ¹⁵²Eu was never the last γ-emitter to be included in the downsized dictionary. However, it can be observed in figure 4 presenting the distribution of estimated counting for ¹⁵²Eu that half of the values is close to zero. This distribution is similar to the case of ²⁴¹Am in figure 1 when this radionuclide is not present in simulated spectra. Thus, a possible interpretation of the systematic false identification is that in the case of a complex mixture covering a large range of γ-photon energies such as in the last test configuration, ¹⁵²Eu is systematically selected because its spectral signature well explains the fit residuals in the first steps of the sparsity selection. But due to the selection of other present radionuclides performed afterwards, most of estimated counting of ¹⁵²Eu correspond to low values close to zero as shown in figure 4.

Zoom In Zoom Out Reset image size

Based on the above analysis, the $\mathcal{P}$-OMP algorithm was improved to overcome the problem of systematic detection of ¹⁵²Eu in the case of complex mixtures. Once the downsized dictionary is obtained, the LRT is applied to the radionuclide having the lowest estimated counting. In order to keep unchanged the expected false positive rate with a significance level of 0.1%, the applied DT value corresponds to the number of unidentified radionuclides remaining for the sparsity selection. The calculation of the false alarm rate was reconsidered using the modified $\mathcal{P}$-OMP algorithm in the case of the last mixture in figure 3: natural Bkg, ²⁴¹Am, ¹³³Ba, ¹³⁷Cs, ⁶⁰Co and ⁵⁷Co. The new value of the false positive rate equal to 0.011% demonstrates that the simple modification in the $\mathcal{P}$-OMP algorithm solved the problem of systematic identification of ¹⁵²Eu. The false positive rate with the modified $\mathcal{P}$-OMP algorithm remains robust when ¹³⁴Cs, ²⁰⁷Bi and ⁸⁸Y are also successively added in simulated spectra. The specific case of ¹⁵²Eu raises the general problem of radionuclide spectra having strong overlapping with a measured spectrum when using scintillation detectors.

3.2.3. Detection limits

The influence of the sparsity constraint on detection limits (DL) was studied through a comparison with DL values determined from individual DT values for each radionuclide using the NNPU algorithm [5]. In the present work, DL values corresponds to mean counting related to a radionuclide for an expected false negative rate of 0.1% (probability of type II error) and using the DT values given in table 3 corresponding to an expected false positive rate of 0.1%. The DL values were determined by performing 10⁵ MC calculations with simulated spectra composed of natural Bkg and the radionuclide of interest. The results given in table 4 correspond to radionuclides usually applied for RPM testing (⁵⁷Co, ⁶⁰Co, ¹³³Ba, ¹³⁷Cs, ²⁴¹Am) and ¹⁵²Eu. No change was observed on DL values obtained with the $\mathcal{P}$-OMP algorithm except for ¹⁵²Eu in which case the value is slightly lower. These results demonstrate that decision-making based on the same DT value for all radionuclides in the selection process of the $\mathcal{P}$-OMP algorithm does not reduce the detection performances compared with those obtained using individual DT values in the previous work [5].

Table 4. Comparison of DL obtained respectively with the $\mathcal{P}$-OMP and the NNPU algorithms.

Radionuclides	$\mathcal{P}$-OMP	NNPU
241Am	60	60
133Ba	160	160
57Co	110	110
60Co	115	115
137Cs	85	85
152Eu	260	280

As previously shown, the sparsity constraint is an efficient technique to obtain the downsized dictionary of present radionuclides in a spectrum when using spectral unmixing. In the present section, the impact of the sparsity constraint on estimated counting and potential biases are addressed as well as the calculation of standard deviations based on the Fisher matrix. For that purpose, counting estimation obtained respectively with the $\mathcal{P}$-OMP and NNPU algorithms were compared by studying three different mixtures containing natural Bkg (mean counting equal to 1250). The dictionary of spectral signatures remains unchanged. For each mixture, the standard deviations of estimated counting were compared with the values determined with the Fisher matrix computed using the expected mean counting of the radionuclides to be identified. The results were obtained by means of 10⁵ MC calculations.

In order to investigate the case of low overlapping with other spectral signatures in the dictionary, a simple mixture comprising ²⁴¹Am and natural Bkg was first considered. The mean counting of ²⁴¹Am was set to 60 corresponding to the DL value for this radionuclide in the case of a false negative rate of 0.1% (see table 4). The average values of estimated counting are presented in table 5. No significant bias is observed with the expected ²⁴¹Am mean counting for both the $\mathcal{P}$-OMP and NNPU algorithms. On the contrary, a significant deviation was obtained on counting estimation of natural Bkg with the NNPU algorithm. Due to overfitting effect, this bias is approximately equal to the sum of average values of counting estimation of non-present radionuclides. This problem is strongly reduced when the $\mathcal{P}$-OMP algorithm is applied because the sparsity constraint forces counting to be equal to zero for non-identified radionuclides. As a result, no significant bias is observed for counting estimation of natural Bkg with the $\mathcal{P}$-OMP algorithm. The residual average values given for non-present radionuclides in table 5 corresponds to counting related to false positive events.

Table 5. Comparison of average values of estimated counting with the $\mathcal{P}$-OMP and NNPU algorithms in the case of a mixture of natural Bkg and ²⁴¹Am.

Radionuclides	$\mathcal{P}$-OMP	NNPU	Expected mean counting
Bkg	1250.0	1205.6	1250
241Am	59.97	59.90	60
133Ba	0.010	8.76	0
207Bi	0.012	6.78	0
57Co	0.005	6.62	0
60Co	0.011	3.87	0
134Cs	0.008	5.07	0
137Cs	0.004	3.22	0
152Eu	0.019	5.60	0
88Y	0.005	3.69	0

The standard deviations for the estimated counting given by both algorithms are presented in table 6. They are also compared to the values computed with the Fisher matrix using the expected mean counting. The results are consistent for the values related to ²⁴¹Am. On the contrary, the standard deviation obtained in the case of the NNPU algorithm for natural Bkg is significantly greater than the value determined from counting estimated with the $\mathcal{P}$-OMP algorithm. This effect is probably due to the fact that the NNPU algorithm is more sensitive to overfitting yielding larger variability on estimated counting of natural Bkg. For both natural Bkg and ²⁴¹Am, the standard deviations obtained from the Fisher matrix are consistent with the values assessed from the distributions of estimated counting given by the $\mathcal{P}$-OMP algorithm.

Table 6. Comparison of standard deviations of estimated counting with the $\mathcal{P}$-OMP and NNPU algorithms and the values computed with the Fisher matrix in the case of a mixture of natural Bkg and ²⁴¹Am.

Radionuclides	$\mathcal{P}$-OMP	NNPU	Fisher
Bkg	36.31	43.99	36.14
241Am	10.92	10.82	10.78
133Ba	0.91	12.55	—
207Bi	0.99	10.47	—
57Co	0.55	9.10	—
60Co	0.77	6.54	—
134Cs	0.71	8.46	—
137Cs	0.39	5.19	—
152Eu	1.56	12.84	—
88Y	0.51	6.01	—

In the second mixture case, ²⁴¹Am was replaced with ¹⁵²Eu for its large range of γ-photon energies (120 keV to 1410 keV) yielding higher overlapping with other spectral signatures in the dictionary. The mean counting of ¹⁵²Eu was set to 260 corresponding to the DL value for this radionuclide in the case of a false negative rate of 0.1% (see table 4). The average values of estimated counting are displayed in table 7. Contrary to the previous tested mixture, significant biases with expected mean counting was observed for both natural Bkg and the radionuclide to be identified (¹⁵²Eu). As previously shown, the application of the $\mathcal{P}$-OMP algorithm corrects for these deviations and yields standard deviations of estimated counting consistent with those computed with the Fisher matrix (see table 8). A comparison between estimated and simulated spectra is presented in figure 5 corresponding to a mixture of natural Bkg and ¹⁵²Eu with respective mean counting equal to 1250 and 260.

Table 7. As in table 5, with a mixture of natural Bkg and ¹⁵²Eu.

Radionuclides	$\mathcal{P}$-OMP	NNPU	Expected mean counting
Bkg	1249.92	1218.95	1250
241Am	0.004	2.93	0
133Ba	0.013	12.51	0
207Bi	0.005	7.37	0
57Co	0.017	12.76	0
60Co	0.010	6.35	0
134Cs	0.013	7.62	0
137Cs	0.003	3.90	0
152Eu	259.92	233.09	260
88Y	0.008	4.53	0

Table 8. As in table 6, with a mixture of natural Bkg and ¹⁵²Eu.

Radionuclides	$\mathcal{P}$-OMP	NNPU	Fisher
Bkg	51.81	54.24	51.48
241Am	0.36	4.38	—
133Ba	1.12	16.23	—
207Bi	0.60	11.43	—
57Co	1.16	14.43	—
60Co	0.80	9.25	—
134Cs	0.99	11.25	—
137Cs	0.39	6.02	—
152Eu	41.26	44.86	40.74
88Y	0.63	7.01	—

Zoom In Zoom Out Reset image size

The third configuration corresponds to a more complex mixture comprising three radionuclides used for RPM testing and ¹⁵²Eu with the respective expected mean counting: 400 (⁶⁰Co), 500 (¹³³Ba), 250 (²⁴¹Am) and 500 (¹⁵²Eu). In this case the expected mean counting were chosen greater than the DL values mentioned in table 4 to take into account their increase due to the mixture. As previously observed for simple mixtures, the most significant biases obtained with the NNPU algorithm in table 9 correspond to natural Bkg and ¹⁵²Eu having both spectral signatures with a large range of γ-photon energies. The $\mathcal{P}$-OMP algorithm solved the problem of estimation biases of the NNPU algorithm for each radionuclide to be identified. In addition, the standard deviations calculated with the Fisher matrix in table 10 are consistent with those determined from the distribution of estimated counting with the $\mathcal{P}$-OMP algorithm.

Table 9. As in table 5, with a mixture of natural Bkg, ⁶⁰Co, ¹³³Ba, ²⁴¹Am and ¹⁵²Eu.

Radionuclides	$\mathcal{P}$-OMP	NNPU	Expected mean counting
Bkg	1249.74	1230.51	1250
241Am	249.95	250.05	250
133Ba	500.02	504.44	500
207Bi	0.017	10.47	0
57Co	0.028	13.96	0
60Co	400.19	398.01	400
134Cs	0.011	9.51	0
137Cs	0.005	5.31	0
152Eu	500.02	471.77	500
88Y	0.004	5.87	0

Table 10. As in table 6, with a mixture of natural Bkg, ⁶⁰Co, ¹³³Ba, ²⁴¹Am and ¹⁵²Eu.

Radionuclides	$\mathcal{P}$-OMP	NNPU	Fisher
Bkg	68.053	68.71	67.85
241Am	19.02	18.95	19.01
133Ba	40.95	41.81	40.99
207Bi	1.38	15.73	—
57Co	1.92	18.14	—
60Co	35.8	37.82	35.76
134Cs	1.08	14.56	—
137Cs	0.56	8.17	—
152Eu	61.18	69.33	60.81
88Y	0.5	9.11	—

It can be drawn from these results obtained in the case of three different mixtures that the determination of a downsized dictionary of identified radionuclides based on the sparsity constraint, is important to estimate reliable counting. For metrological applications, it also enables the calculations of standard deviations that are consistent with those assessed from the distributions of estimated counting given by the $\mathcal{P}$-OMP algorithm.

In the previous section, it was shown that the standard deviations of estimated counting provided by the $\mathcal{P}$-OMP algorithm are consistent with the values determined with the Fisher information matrix computed with expected mean counting. In this part, the Gaussian approximation mentioned in 2.4 was investigated when the Fisher matrix is directly assessed with estimated counting given by the $\mathcal{P}$-OMP algorithm. For that purpose, the coverage probability related to standard deviations was analysed for three values of the coverage factor k with k = 1–3. The differences between expected mean counting and estimated values were classified according to coverage probabilities in order to be compared with a normal distribution (68.3% for k = 1, 95.5% for k = 2 and 99.7% for k = 3) using the following expression for each identified radionuclide i (with a_i the expected mean counting):

$$\begin{equation}-k\sqrt{{{\left[\boldsymbol{I}\left(\hat{\boldsymbol{a}}\right)\right]}^{-1}}_{ii}}{\leqslant}{\hat{a}}_{i}-{a}_{i}{\leqslant}+k\sqrt{{{\left[\boldsymbol{I}\left(\hat{\boldsymbol{a}}\right)\right]}^{-1}}_{ii}}.\end{equation} \tag{ 16 }$$

The analysis on coverage intervals was carried out by MC calculations with the same mixture cases presented in the previous section 3.3. The results for the two simple mixtures comprising respectively ²⁴¹Am and ¹⁵²Eu are presented in table 11. The coverage intervals obtained in the case of the complex mixture (⁶⁰Co, ¹³³Ba, ¹⁵²Eu, ²⁴¹Am and natural Bkg) are displayed in table 12. Considering the results in tables 11 and 12, it can be concluded that the Gaussian approximation of the standard deviations calculated with the Fisher matrix is verified when using sparse spectral unmixing. Thus, reliable standard deviations with coverage intervals corresponding to those expected from a normal distribution can be obtained when the Fisher matrix is directly calculated with estimated counting. This result is important for metrological application of sparse spectral unmixing.

In this section, the DT calculation based on the Bonferroni correction established for the sparsity constraint in the case of a 3'' × 3'' NaI(Tl) detector was verified for a dictionary of spectral signatures obtained with a 3'' × 3'' EJ240 plastic detector. This study was carried out in the same manner as described in 3.1 for an expected false positive rate equal to 0.1%.

The LRT statistics was first investigated using simulated spectra of natural Bkg according to increasing dictionary sizes by including successively a different γ-emitter in the order given in table 13. Each DT value associated to each dictionary size was determined by performing 5 × 10⁵ MC calculations. The comparison between simulated and calculated DT values is presented in table 13. It can be drawn from these results that the DT calculation based on the Bonferroni correction provides a good approximation of the evolution of the DT value according to the number of radionuclides available for the sparsity selection.

The robustness of the sparsity constraint was also tested in the case of simple mixtures composed of natural Bkg and respectively ¹³³Ba, ⁶⁰Co and ⁵⁷Co according to expected mean counting equal to 1000 and 10 000. The evolution of false positive rates for the three mixtures is reported in table 14. As already observed with the NaI(Tl) detector, no visible trend was obtained on false positive rates when a radionuclide is mixed with natural Bkg according to increasing counting.

The robustness of decision-making for more complex mixtures was performed using the same procedure described in section 3.2.2. For each mixture composed of natural Bkg, simulated spectra were obtained by successively including an additional γ-emitter used for RPM testing according to the following order: ²⁴¹Am, ¹³³Ba, ¹³⁷Cs, ⁶⁰Co and ⁵⁷Co. For each radionuclide, the mean counting was equal to 1200. The evolution of the false positive rates according to increasing dictionary size is displayed in figure 7. It can be observed that the false positive rates remain robust until the addition of ⁵⁷Co in simulated spectra. The last configuration yields an increase of the false alarm rate to about 0.3%. This result is mainly due to false identification of ^99mTc (not present in simulated spectra).

Zoom In Zoom Out Reset image size

In order to figure out the correlations between estimated counting of ⁵⁷Co and ²⁴¹Am and false positive events of ^99mTc, a classification of those results according to three groups are respectively depicted in figures 8 and 9.

• The first group (green set composed of 34 values) in figure 8 represents false positive events with estimated counting of ^99mTc lower than 250. This cluster corresponds to estimated counting of ⁵⁷Co slightly lower than the expected value. In the case of ²⁴¹Am in figure 9, the average value of the plots (equal to about 1230) is slightly greater than the expected mean counting (equal to 1200).
• The second group (blue set composed of 493 values) in figure 8 represents false positive events when ^99mTc is systematically selected instead of ⁵⁷Co (estimated counting equal to zero). In figure 9, this group is correlated to estimated counting of ²⁴¹Am (average value approximately equal to 834) significantly lower than the expected value.
• The last group (red set composed of 99 values) depicts estimated counting of ⁵⁷Co (mainly lower than to 400) significantly lower than the expected value. Similarly to the blue plots, estimated counting of ^99mTc (average value close to 1800) are higher than expected values of present radionuclides. This last group is also correlated to estimated counting of ²⁴¹Am (average value approximately equal to 800) significantly lower than the expected value (see figure 9).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The interpretation of these results is different from the case of ¹⁵²Eu with the NaI(Tl) detector described in section 3.2.2. Indeed, as shown in figures 8 and 9, the false positive events of ^99mTc are mainly represented by the two last groups (blue and red plots) showing significant values of estimated counting of ^99mTc. These two groups are also correlated to underestimated counting of both ⁵⁷Co and ²⁴¹Am or to the non-identification of ⁵⁷Co. Contrary to the case of ¹⁵²Eu, the modified $\mathcal{P}$-OMP algorithm cannot avoid these false positive detections of ^99mTc. The correlation observed with both ⁵⁷Co and ²⁴¹Am estimated counting leads to the possible interpretation that the problem comes from the strong overlapping existing between the spectral signatures of ⁵⁷Co, ²⁴¹Am and ^99mTc. The mixture of five radionuclides and natural Bkg investigated here has pointed out a potential limitation of sparse spectral unmixing when applied to plastic detectors.

Considering the results obtained in the previous section, the impact of spectral signature overlapping on counting estimation was studied with mixtures containing ^99mTc and ⁵⁷Co. The comparison of counting estimation obtained with the $\mathcal{P}$-OMP and NNPU algorithms was implemented by performing 2 × 10⁵ MC calculations.

The first studied case corresponds to a mixture composed of ^99mTc and ⁵⁷Co with natural Bkg with the respective expected mean counting 1300, 1000 and 1500. The comparison results are given in table 15. As already reported in the case of the NaI(Tl) detector, the counting estimation bias observed on natural Bkg when the NNPU algorithm is applied, vanished with the $\mathcal{P}$-OMP algorithm. However, a slight bias still exists on the counting estimation of ^99mTc and ⁵⁷Co that highlights the problem related to overlapping of spectral signatures mentioned in the previous section. In addition, it can be pointed out in table 15 that the average value of ²⁴¹Am is greater than those obtained with the other radionuclides non-present in simulated spectra. This result comes from the 1185 false positive events of ²⁴¹Am (2 × 10⁵ MC calculations) obtained with the $\mathcal{P}$-OMP algorithm and yielding an increase of the total false positive rate to a value equal to about 0.7%. These results can be linked to the standard deviations in table 16 computed with the Fisher matrix using expected mean counting. The comparison with the standard deviations obtained from the distribution of estimated counting, indicates that the values given by the Fisher matrix are underestimated for ^99mTc and ⁵⁷Co.

Table 15. Comparison of average values of estimated counting with the $\mathcal{P}$-OMP and NNPU algorithms in the case of a mixture of natural Bkg with ^99mTc and ⁵⁷Co.

Radionuclides	$\mathcal{P}$-OMP	NNPU	Expected mean counting
Bkg	1499.8	1364.8	1500
241Am	2.2	32.9	0
57Co	996	993.9	1000
60Co	0.03	19.8	0
137Cs	0.03	15.7	0
207Bi	0.05	20.6	0
133Ba	0.05	43.9	0
99mTc	1301.7	1309	1300

Table 16. Comparison of standard deviations of estimated counting with the $\mathcal{P}$-OMP and NNPU algorithms and the values computed with the Fisher matrix in the case of mixture of natural Bkg in the case of a mixture of natural Bkg with ^99mTc and ⁵⁷Co.

Radionuclides	$\mathcal{P}$-OMP	NNPU	Fisher
Bkg	54.5	117.1	53.4
241Am	28.6	54.5	—
57Co	208.9	187	172.1
60Co	1.8	24.1	—
137Cs	1.9	22.1	—
207Bi	3	32 .1	—
133Ba	2.9	45.6	—
99mTc	234.4	234.3	185.5

In table 17, the coverage intervals associated with the standard deviations were assessed using directly estimated counting in the Fisher matrix in the same way as described in section 3.4. It is noteworthy that, despite the biases on standard deviations reported in table 16, the coverage probabilities are still close to those expected from a normal distribution. The difference is more significant for a coverage factor k equal to 3.

Having specifically tested the problem of a mixture of ^99mTc and ⁵⁷Co on counting estimation, a different γ-emitter mixture was addressed with the following composition: natural Bkg with ²⁴¹Am, ¹³³Ba and ^99mTc (without ⁵⁷Co) according to the following mean counting 1500, 600, 700 and 1000. As presented in table 18, contrary to the previous case, the absence of ⁵⁷Co in simulated spectra leads to better counting estimation when the $\mathcal{P}$-OMP algorithm is applied. It is noteworthy that the significant bias on counting estimation of ⁵⁷Co with the NNPU algorithm is also strongly reduced. The false positive rate is equal to 0.088%, which is in good agreement with the expected false positive rate. It is noteworthy that no systematic false event of ⁵⁷Co was observed despite the correlations between ⁵⁷Co, ^99mTc and ²⁴¹Am.

Table 18. As in table 15, with a mixture of natural Bkg with ²⁴¹Am, ¹³³Ba and ^99mTc.

Radionuclides	$\mathcal{P}$-OMP	NNPU	Expected mean counting
Bkg	1499.7	1383.9	1500
241Am	600.4	623.9	600
57Co	0.3	58.3	0
60Co	0.03	17.5	0
137Cs	0.03	14.5	0
207Bi	0.05	19	0
133Ba	700	740	700
99mTc	999.4	943.3	1000

Regarding the standard deviations reported in table 19 related to counting results in table 18, a good agreement was obtained between calculated and simulated values when using the $\mathcal{P}$-OMP algorithm compared to the previous mixture involving both ^99mTc and ⁵⁷Co. As reported in table 20, a better agreement is also obtained on the approximation of coverage intervals with a normal distribution when the Fisher matrix is calculated using directly estimated counting.

Table 19. As in table 16, with a mixture of natural Bkg with ²⁴¹Am, ¹³³Ba, ^99mTc.

Radionuclides	$\mathcal{P}$-OMP	NNPU	Fisher
Bkg	75.9	132.6	75.3
241Am	82.4	89.5	81.9
57Co	12.9	104.2	—
60Co	2	24.1	—
137Cs	2	21.6	—
207Bi	3.2	31.8	—
133Ba	61.6	71.6	61.6
99mTc	97.8	179.4	95.4