Temporal energy window

In this subsection, the temporal energy window is proposed for sample processing of time-series detector response data. Samples were partitioned into multiple segments with consideration of the sample type.

Urban landscapes teem with roads and structures composed of natural and man-made substances. Natural Occurring Radioactive Materials (NORMs) are inherent in these substances, with concentrations differing across materials. Predominantly, NORM comprises isotopes, such as ⁴⁰K, ²³⁸U, and ²³²Th, along with the radioactive daughter products of the latter two, commonly denoted as KUT [26]. As detectors navigate the search zone, particularly when radioactive substances are unmonitored, the makeup of the neighboring structures and their ensuing radioactive signatures shift with each locale [27]. Consequently, the cumulative gamma photon count rate and spectra recorded by detectors might demonstrate notable fluctuations [28]. Adding to the complexity, illicit radioactive substances may be concealed, leading to attenuated detection signals that are challenging to identify. The interplay between gamma photons and diverse substances, mediated by various physical processes, amplifies the dynamism of the observed radiation background signal. Hence, the time-series detector response data acquired in urban settings are profoundly shaped by ambient conditions, often overshadowing the distinctive peaks that mark the presence of radioactive materials in the energy spectra.

To counteract the effects of variable backgrounds, enhance the integrity and dependability of the energy spectrum data, and streamline subsequent data processing, we segmented the time-series detector response data using a temporal energy window. This strategy primarily underscores the features of faint radioactive materials.

In an urban setting, the detection system operates under two potential conditions: with or without the presence of an auxiliary radiation source, which is contextualized against the background radiation. Consequently, the time-series detector response dataset encompasses active and passive samples. An active sample pertains to the detector response data captured in the presence of a radioactive source, while a passive sample relates to data collected in an environment devoid of any radioactive source.

The time-series detector response dataset is defined as $T$. $T=\left\{\left(S_1\mathrm{,}{y}_1\mathrm{,}lo{c}_1\right)\mathrm{,}\left(S_2\mathrm{,}{y}_2\mathrm{,}lo{c}_2\right)\mathrm{,}\dots \mathrm{,}\left(S_M\mathrm{,}{y}_M\mathrm{,}lo{c}_M\right)\right\}$ Si is the matrix of time-series detector response data and defined by Eq.1. Furthermore, M denotes the number of samples in the dataset. i=1,2,…, M. $S_i=\left[t^i\mathrm{,}{e}^i\right]=\left[\begin{array}{cc}{t}_1^i& e_1^i\\ t_2^i& e_2^i\\ ⋮& ⋮\\ t_m^i& e_m^i\end{array}\right]$ (1) tⁱ denotes the time series, $t_j^i\in {\mathbb{Q}}^{+}$, and 𝒆i denotes the energy value recorded by the detector over time, and $e_j^i\in {\mathbb{Q}}^{+}$. m denotes the length of 𝑺i, with j=1,2,...,m.

When 𝑺i is an active sample, $y_i\in \mathcal{Y}$ is the class label of 𝑺i. Furthermore, $lo{c}_i\in {\mathbb{Q}}^{+}$ denotes the time point when the detector is closest to the radioactive source position during the movement, $t_1^i<lo{c}_i<t_m^i$. When 𝑺i is a passive sample, yi=0 and loci=0.

The temporal energy window is proposed for sample processing of time-series detector response data. The quantity and length of a temporal energy window were defined as $w_{\text{q}}\in {\mathbb{Z}}^{+}$ and $w_{\text{l}}\in {\mathbb{Z}}^{+}$. The key of utilizing a temporal energy window for processing time-series detector response data lies in determining the temporal origin of the window, which refers to the initial point from which the temporal energy window conducts the partition task, thereby determining the position of the window within the time series. The temporal origin of an energy window is defined as: tj’ and j’ is calculated by Eq. 2. $\begin{array}{c}{j}^{\prime }=\left\{\lfloor \mathrm{0,}{w}_{\text{q}}-1\rfloor \times \lfloor \frac{m}{2\times w_{\text{q}}}\rfloor \right\}\\ \times \overline{b}+\left\{\underset{j^{\prime }\in [\mathrm{1,}m]}{\text{argmin}}\left|lo{c}_i-t_{j^{\prime }}\right|-\lfloor \frac{1}{2}\times w_l\rfloor \right\}\times b\end{array}$ (2) In Eq. 2, b denotes a Boolean variable. Assuming that the present sample is active, b is true; whereas if the present sample is passive, then b is false. Obviously, tj’ of a passive sample is distributed evenly at several different locations in the time axis, while tj’ of an active sample is fixed due to the demand for obtaining energy fragments as close to the source as possible.

The segmented time-series detector response dataset processed by a temporal energy window is denoted as T_seg and represented by Eq. 3. Samples in T_seg may have been expanded in comparison to T, which is dependent on w_q. $\begin{array}{c}{T}_{\text{seg}}=\{\left(S_{seg}^1\mathrm{,}{y}_1\mathrm{,}lo{c}_1\right)\mathrm{,}\left(S_{seg}^2\mathrm{,}{y}_2\mathrm{,}lo{c}_2\right)\mathrm{,}\cdots \mathrm{,}\\ \left(S_{seg}^M\mathrm{,}{y}_M\mathrm{,}lo{c}_M\right)\}\\ =\{\left(S_{{\text{seg}}_1}^1\mathrm{,}{y}_1\mathrm{,}lo{c}_1\right)\mathrm{,}\left(S_{{\text{seg}}_2}^1\mathrm{,}{y}_1\mathrm{,}lo{c}_1\right)\mathrm{,}\\ \cdots \mathrm{,}\left(S_{{\text{seg}}_{w_{\text{q}}}}^1\mathrm{,}{y}_1\mathrm{,}lo{c}_1\right)\mathrm{,}\cdots \mathrm{,}\left(S_{{\text{seg}}_{w_{\text{q}}}}^M\mathrm{,}{y}_M\mathrm{,}lo{c}_M\right)\}\end{array}$ (3) In Eq.3, $S_{seg}^1$ denotes the group of segmented time-series detector response data of 𝑺i and represented by Eq.4, $seg=\left[{\text{seg}}_1\mathrm{,}{\text{seg}}_2\mathrm{,}\dots \mathrm{,}{\text{seg}}_{w_q}\right]$. $S_{seg}^i=\left\{\left[t_{seg}^i\mathrm{,}{e}_{seg}^i\right]\right\}$ (4) Employing $S_{se{g}_k}^i$ to symbolize each individual segmented sample of $S_{seg}^i$. $S_{se{g}_k}^i$ is denoted by Eq.5. Obviously, the length of each $S_{se{g}_k}^i$ is correlated with the length of the temporal energy window wl. $S_{se{g}_k}^i=\left[t_{se{g}_k}^i\mathrm{,}{e}_{se{g}_k}^i\right]=\left[\begin{array}{ccc}{t}_{j^{\prime }}^{\left(i\mathrm{,}k\right)}& ,& e_{j^{\prime }}^{\left(i\mathrm{,}k\right)}\\ t_{j^{\prime }+1}^{\left(i\mathrm{,}k\right)}& ,& e_{j^{\prime }+1}^{\left(i\mathrm{,}k\right)}\\ ⋮& & ⋮\\ t_{j^{\prime }+w_l}^{\left(i\mathrm{,}k\right)}& ,& e_{j^{\prime }+w_l}^{\left(i\mathrm{,}k\right)}\end{array}\right]$ (5)

Peak-ratio spectrum analysis

In this subsection, we delve into the formation mechanism and measurement principle of the gamma-ray instrument spectrum. These are leveraged as the foundation for extracting spectral features. Key features include aggregated counts, peak-to-flat ratios, and peak-to-peak ratios. This type of an approach aids in the analysis and interpretation of the intrinsic significance of energy spectra.

After processing through the temporal energy window, the segmented time-series detector response data is converted into an energy spectrum format, easing the subsequent feature extraction. The energy spectrum provides a distribution curve mapping the count rate against particle energy, a pivotal tool in detecting and identifying radioactive nuclear materials.

For the context of this study, the relative distance between the detector and radiation source is in constant flux due to the detector’s movement. It is essential to underline that this study primarily focuses on scenarios with static radiation sources. Dynamics, such as the continuous movement of the source or its dissolution in water, have not been contemplated. Owing to the finite number of photon counts within the full energy peak, statistical fluctuations become pronounced. Consequently, the channel with the peak counts might not align with the expected value of a Gaussian distribution [29, 30]. To mitigate the effects of these statistical fluctuations, spectral data are reorganized into multiple bins along the energy axis. Each bin encompasses an energy range, and counts within this range are consolidated to create a novel feature vector.

The transformed spectrum dataset is denoted as T_spe. $T_{\text{spe}}=\left\{\left(x_{seg}^1\mathrm{,}{y}_1\right)\mathrm{,}\left(x_{seg}^2\mathrm{,}{y}_2\right)\mathrm{,}\cdots \mathrm{,}\left(x_{seg}^M\mathrm{,}{y}_M\right)\right\}$

$x_{seg}^i$ represents the transformed spectrum data from $S_{seg}^i$, i=1,2,⋅s,M and $seg=\left\{se{g}_1\mathrm{,}se{g}_2\mathrm{,}\dots \mathrm{,}se{g}_{w_{\text{q}}}\right\}$. $y_i\in \mathcal{Y}$ is the class label of $x_{seg}^i$. During the process of transformation, information of loc and t are discarded. The transformed spectrum data $x_{seg}^i$ is in the form of vector and represented by Eq.6. $x_{seg}^i=\left\{\left(x_{se{g}_1}^i\mathrm{,}{y}_1\right)\mathrm{,}\left(x_{se{g}_2}^i\mathrm{,}{y}_1\right)\mathrm{,}\cdots \mathrm{,}\left(x_{se{g}_{w_{\text{q}}}}^i\mathrm{,}{y}_1\right)\right\}$ (6)

$x_{se{g}_k}^i$ represents the transformed energy spectrum of the k-th segment of the i-th sample of the time-series detector response dataset. $x_{se{g}_k}^i=\left\{x_1^{\left(i\mathrm{,}k\right)}\mathrm{,}{x}_2^{\left(i\mathrm{,}k\right)}\mathrm{,}\cdots \mathrm{,}{x}_n^{\left(i\mathrm{,}k\right)}\right\}$ (7) i=1,2,⋅s,M and k=1,2,...,w_q. $x_n^{\left(i,k\right)}$ is the aggregated count of the n-th bin of $x_{se{g}_k}^i$.

The length of each $x_{se{g}_k}^i$, i.e. n, is the same because the same energy range is precomputed before the transformation process and a fixed number of bins is selected uniformly.

The value of n is determined based on the Expected value of the maximum and minimum energy values across all samples in T. The expected value of the maximum and minimum energy values are denoted as $\lceil \mathbb{E}\left[\max \left(e^i\right)\right]\rceil$ and $\lfloor \mathbb{E}\left[\min \left(e^i\right)\right]\rfloor$ i=1,2,⋅s,M. $n\in {\mathbb{Z}}^{+}$ and the value of n is calculated by Eq.8. $n=\lceil \mathbb{E}\left[\max \left(e^i\right)\right]\rceil -\lfloor \mathbb{E}\left[\min \left(e^i\right)\right]\rfloor$ (8) Furthermore, $x_n^{\left(i\mathrm{,}k\right)}$ indicates the photon count in $e_{{\text{seg}}_k}^i$ in the corresponding energy interval, and $x_n^{\left(i\mathrm{,}k\right)}$ is calculated by Eq. 9. Specifically, countif(A, B) is a function that searches the range A for items that match condition B and counts them. Additionally, α is applied to represent the energy range to fine-tune the transform accuracy of the energy spectrum. $\begin{array}{c}{x}_n^{\left(i\mathrm{,}k\right)}=\text{countif}(e_{{\text{seg}}_k}^i\mathrm{,}\lfloor \mathbb{E}\left[\min \left(e^i\right)\right]\rfloor +\left(n-1\right)\\ \times \alpha ⩽e_{{\text{seg}}_k}^i<\lfloor \mathbb{E}\left[\min \left(e^i\right)\right]\rfloor +n\times \alpha )\end{array}$ (9) Detection and identification of radioactive materials primarily hinge on nuclear radiation detectors, which capture gamma rays emitted during the decay process. The measurement of gamma ray energy is determined by registering the energy dispersed within the detector. The main mechanisms driving gamma energy spectrum measurements encompass three interactions between gamma rays and the detector medium: the photoelectric effect, the Compton effect, and pair production.

Low-energy gamma rays (0 – hundreds of keV) predominantly undergo the photoelectric effect, resulting in at least one distinct photoelectric peak. Medium-energy gamma rays (hundreds of keV – 3 MeV) primarily interact through the Compton effect. Conversely, high-energy gamma rays (5–10 MeV and beyond) are primarily subject to pair production. The photoelectric peak, when the energy of the incident gamma radiation is below 1.02 MeV, is often termed as the full energy peak. This peak is traditionally considered as the primary hallmark for identifying specific radioactive nuclides. The full energy peak arises from the sum of the photoelectric peak’s energy combined with energy from Compton electrons and photoelectrons stemming from Compton scattering interactions. In the spectrum of low-to-medium energy gamma rays, pair production is negligible. Instead, the energy spectrum is characterized by a Compton continuum and photoelectric peaks. When gamma rays possess intermediate energy, incident gamma photons experience multiple successive Compton scatterings. The energy from the recoil electrons, produced from these scatterings, is deposited in the detector. Notably, the cumulative energy of these recoil electrons can surpass the energy transfer’s upper limit in a single scattering event, filling regions between the Compton edge and photoelectric peaks [9].

From the prior discussion on the formation mechanism and measurement principle of the gamma-ray instrument spectrum, it is clear that gamma energy spectra contain both photoelectric peaks and the Compton continuum. Conventionally, the photoelectric peaks serve as the primary identifiers for radionuclides. Conversely, the Compton continuum, which often exhibits similar shapes across different contexts, is usually overlooked. However, relying solely on characteristic peaks for radioactive nuclide identification may fall short in complex background situations [31-33]. Drawing inspiration from Ref. [7], this subsection introduces peak-to-flat ratios and peak-to-peak ratios as descriptors for the spectral features.

Equation 7 defines the form of energy spectrum after binning. Here, $x_n^{(i\mathrm{,}k)}$ indicates the photon counts in the corresponding energy bins and can be calculated by Eq.9. Based on this, specific bins are selected according to the decay properties of the radionuclide material, which correspond to the area of theoretical Compton continuum, characteristic peaks, and auxiliary peaks, respectively. For $x_{{\text{seg}}_k}^i$, the area of the theoretical Compton continuum, characteristic peaks, and auxiliary peaks are represented as a_c, a_f, and a_u, respectively, and defined as Eq.10– Eq.12, respectively. i=1,2,⋅s,M and k=1,2,...,w_q. $a_{\text{c}}={\displaystyle \sum _{n^{\prime }\in \left[c_{\text{l}}\mathrm{,}{c}_{\text{r}}\right]}\left(x_{n^{\prime }}^{\left(i\mathrm{,}k\right)}\right)}$ (10) $a_{\text{f}}={\displaystyle \sum _{n^{\prime }\in \left[f_{\text{l}}\mathrm{,}{f}_{\text{r}}\right]}\left(x_{n^{\prime }}^{\left(i\mathrm{,}k\right)}\right)}$ (11) $a_{\text{u}}={\displaystyle \sum _{n^{\prime }\in \left[u_{\text{l}}\mathrm{,}{u}_{\text{r}}\right]}\left(x_{n^{\prime }}^{\left(i\mathrm{,}k\right)}\right)}$ (12) The boundaries of the Compton continuum are denoted by c_l and c_r, while the characteristic peaks are bounded by f_l and f_r, and the auxiliary peaks are bounded by a_l and a_r. These boundaries are selected based on the decay properties of the radionuclide material. The peak-to-flat ratio r₁ and peak-to-peak ratio r₂ are defined as Eq.13. $r_1=a_{\text{f}}/a_{\text{c}}\mathrm{,}{r}_2=a_{\text{u}}/a_{\text{f}}$ (13) Based on a macroscopic perspective, r₁ characterizes the capability to discern low-energy weak peaks amidst a complex background, while r₂ measures the likelihood of gamma rays experiencing multiple interactions within the detector, culminating in their contribution to the full-energy peaks.

Classification

In this subsection, considering the requirements for imbalanced multi-classification, we developed a detection and identification model using the weighted KNN architecture. By capitalizing on the inherent trait that energy spectra from identical radioactive materials exhibit minimal inter-class variability, the model significantly boosts the classification accuracy for underrepresented classes and improves the overall efficacy of the classifier.

Following the peak-ratio spectrum analysis, we derive a set of feature vectors comprised of aggregated counts, peak-to-flat ratios, and peak-to-peak ratios, denoted as T_app. $T_{\text{app}}=\left\{\left(c_{seg}^1\mathrm{,}{y}_1\right)\mathrm{,}\left(c_{seg}^2\mathrm{,}{y}_2\right)\mathrm{,}\cdots \mathrm{,}\left(c_{seg}^M\mathrm{,}{y}_M\right)\right\}$ $c_{seg}^i$ represents the set of feature vectors from $x_{seg}^i$, which is denoted by Eq.14. i=1,2,…,M and $seg=\left\{{\text{seg}}_1\mathrm{,}{\text{seg}}_2\mathrm{,}\dots \mathrm{,}{\text{seg}}_{w_q}\right\}$. $y_i\in \mathcal{Y}$ is the class label. $c_{seg}^i=\left\{\left(c_{se{g}_1}^i\mathrm{,}{y}_i\right)\mathrm{,}\left(c_{se{g}_2}^i\mathrm{,}{y}_i\right)\mathrm{,}\cdots \mathrm{,}\left(c_{se{g}_{w_{\text{q}}}}^i\mathrm{,}{y}_i\right)\right\}$ (14) Here, $c_{{\text{seg}}_k}^i$ represents the feature vector of the k-th segment of the i-th sample in the original time-series detector response dataset. $c_{{\text{seg}}_k}^i=\left[x_{{\text{seg}}_k}^i\mathrm{,}{r}_1^{\left(i\mathrm{,}k\right)}\mathrm{,}{r}_2^{\left(i\mathrm{,}k\right)}\right]$ (15) i=1,2,…, M and k=1,2,...,wq. $x_{{\text{seg}}_k}^i$, $r_1^{\left(i\mathrm{,}k\right)}$, and $r_2^{\left(i\mathrm{,}k\right)}$ denote the aggregated counts, peak-to-flat ratio, and peak-to-peak ratio of the k-th segment of the i-th sample in the original time-series detector response dataset, respectively.

The sample to be classified is represented as 𝑺₀, with its corresponding feature vector symbolized as 𝒄₀. 𝒄₀ and $c_{{\text{seg}}_k}^i$ are n-dimensional vectors, i.e., ${\text{c}}_0\in {ℝ}^n$ and ${\text{c}}_{se{g}_k}^i\in {ℝ}^n$.

The function f evaluated at the sample point $c_{{\text{seg}}_k}^i$ is yi, i.e., $y_i=f\left(c_{{\text{seg}}_k}^i\right)$. The vector of observations is defined as: $\text{y}=\left[\begin{array}{c}{y}_1\\ ⋮\\ y_M\end{array}\right]$ To construct a surrogate model $\widehat{f}$ of a function f, sample points $c_{{\text{seg}}_k}^i$ are acquired into the matrix. The dimension of the matrix $C$ is M × w_q × n. $C=\left[\begin{array}{c}{c}^{{\left(1,1\right)}^{\text{T}}}\\ ⋮\\ c^{{\left(M,w_{\text{q}}\right)}^{\text{T}}}\end{array}\right]$ Notation c^(I,k) signifies $c_{{\text{seg}}_k}^i$. Here, K denotes the number of nearest neighboring sample points. 𝒁 denotes the set of K sample points that are closest to 𝒄₀ in terms of distance, which is denoted by Eq.16. $\text{Z}\subseteq \text{C}$ and $\left|Z\right|=K$. $Z=\left\{c^{\left(i,k\right)}\in C:\text{rank}\left(d\left(c_0\mathrm{,}{c}^{\left(i,k\right)}\right)\right)\le K\right\}$ (16) Function $\text{rank}\left(d\left(c_0\mathrm{,}{c}^{\left(i,k\right)}\right)\right)$ represents the ranking of the distance $d\left(c_0\mathrm{,}{c}^{\left(i,k\right)}\right)$ in ascending order.

The surrogate model $\widehat{f}$ of function f is defined by Eq.17. $\widehat{f}\left(c_0\right)=\frac{{\displaystyle {\sum }_{{}^{c\left(i,k\right)}\in C}w}\left(c_0\mathrm{,}{c}^{\left(i,k\right)}\right)y_i}{{\displaystyle {\sum }_{c^{(i\mathrm{,}k)}\in C}w}\left(c_0\mathrm{,}{c}^{\left(i,k\right)}\right)}$ (17) w(·) is the inverse distance weight function, which is defined by Eq.18. $w{\left(c_0\mathrm{,}{c}^{\left(i,k\right)}\right)}_{c^{\left(i,k\right)}\in C}=\frac{1}{{\left[d\left(c_0\mathrm{,}{c}^{\left(i,k\right)}\right)\right]}^q}$ (18) where q denotes a normalization power and $d\left(c_0\mathrm{,}{c}^{\left(i,k\right)}\right)$ denotes the distance between the target point c₀ and sample point $c_{{\text{seg}}_k}^i$. The metric used to measure this distance is Lp-norm, which is defined by Eq.19. $d\left(c_0\mathrm{,}{c}^{\left(i,k\right)}\right)={\Vert c_0-c^{\left(i,k\right)}\Vert }_p={\left({\displaystyle \sum _{\beta =1}^m{\left|c_{\mathrm{0,}\beta }-c_{\beta }^{\left(i,k\right)}\right|}^p}\right)}^{1/p}$ (19) Here, β denotes the index of the dimension of the vector. $c_{\beta }^{\left(i,k\right)}$ denotes the β-th dimension of $c_{{\text{seg}}_k}^i$, and $c_{\mathrm{0,}\beta }$ denotes the β-th dimension of 𝒄₀. Specifically, L₁-norm (where p=1) represents the rectangular distance, while L₂-norm (where p=2) represents the Euclidean norm. Furthermore, $L_{\infty }$-norm (where $p\to \infty$) represents the maximum norm.

To limit the effect of farther samples points and also avoid divisions by zero, the distance function is implemented as follows. If ${\Vert c_0-c^{\left(i,k\right)}\Vert }_p=0$, then distance $d\left(c_0\mathrm{,}{c}^{\left(i,k\right)}\right)$ is set to ε. If $0<{\Vert c_0-c^{\left(i,k\right)}\Vert }_p\le R$, then distance $d\left(c_0\mathrm{,}{c}^{\left(i,k\right)}\right)$ is calculated by Eq.19. If ${\Vert {\text{c}}_0-{\text{c}}^{(i\mathrm{,}k)}\Vert }_p>R$, then distance $d\left(c_0\mathrm{,}{c}^{\left(i,k\right)}\right)=0$. Where ε is a small number, and R is the radius of the distance function $d(\cdot )$. Table 1 summarizes the overall flow of the algorithm.

Introduction of data source

The experimental data utilized in this study originated from a time-series detector response dataset, representing a NaI(Tl) detector’s movement within a simulated city block using the Monte Carlo method. This dataset was curated by J. M. Ghawaly Jr and his team at Oak Ridge National Laboratory (ORNL) [34]. Fig.2 offers a visual representation capturing the core features of the dataset. The model simulated seven interconnected city blocks in three dimensions, encompassing various buildings, sidewalks, roads, parking areas, and other urban elements. The Naturally Occurring Radioactive Materials (NORM) incorporated included ⁴⁰K, ²³²Th and its progeny, as well as ²³⁸U/²³⁵U and their respective offspring. The concentration of each component within the KUT (Potassium, Uranium, and Thorium) might vary depending on the specific material. Each block’s background radiation was individually computed. The radioactive materials were potentially concealed in 15 distinct spots. Each radioactive source could exist in one of two states: either unshielded or shielded by 1 cm of lead. A NaI(Tl) detector navigated through these city blocks without the interference of cars or other forms of clutter.

Fig. 2Full-screen View

(Color online) Schematic diagram of the fundamental characteristics of the dataset. This model consisted of seven modular city blocks, and the order of the blocks can be adjusted. Size of the model was 989–1047 m × 201 m × 158 m. For each component of the blocks, every NORM isotope in each material (asphalt, brick, granite, concrete, and soil) in its composition was modeled. These data form the background of the urban environment. A 2″×4″×16″ NaI(Tl) detector traversed the city block in the absence of cars or other forms of clutter. The velocity of the detector was a value in the range of 1–13.4 m/s and remains constant. The walls of the buildings in the model were 6 in (15.24 cm) thick [34].

The dataset comprises radioactive materials from two categories: special nuclear materials (SNMs) and common sources. The SNMs are represented by highly enriched uranium (HEU) and weapons-grade plutonium (WGPu), while the common sources are Technetium-99m (^99mTc), Iodine-131 (¹³¹I), and Cobalt-60 (⁶⁰Co). Both HEU and WGPu are characterized by energy spectra dominated by prompt fission neutrons and prompt gamma rays, which are emitted during fission. These gamma rays possess a broad energy range, spanning from several hundred keV up to multiple MeV. Conversely, ^99mTc releases gamma rays that predominantly linger around the 140-keV energy mark. ¹³¹I emits mainly beta particles accompanied by gamma rays; the beta particles peak at energies near 606 keV. The emitted gamma rays have varying energy levels, with the most notable peaks observed at 364 keV and 637 keV. Finally, ⁶⁰Co radiates gamma rays that prominently feature two energy peaks, one at 1.17 MeV and the other at 1.33 MeV [9].

Specifically, 9700 samples were labeled and listed in Table 2, of which 4900 were background samples without any radioactive materials, while the remaining 4800 samples contained radioactive materials.

Comparative experiments

In this subsection, to optimize and assess the model’s performance while ensuring its practical applicability, the time-series detector response dataset was partitioned into three distinct subsets: training (60%), validation (20%), and testing (20%). This division was subjected to ten-fold cross-validation. A stratified random split was adopted, guaranteeing a balanced representation of radioactive materials across all subsets. The model was implemented in Python, leveraging the capabilities of the PyTorch framework. For comparative analysis, the Weka machine learning toolkit was employed. All experiments were executed on a system furnished with an Intel Core i7 processor, 16GB of RAM, and an NVIDIA GeForce RTX 3070 graphics card.

To streamline the discussion, the proposed algorithm in the experimental results will be referred to as TPW. In the experiments detailed in this subsection, the values of K, p, and q were set to 5, 2, and 2, respectively. A comprehensive examination and discussion regarding the selection of these parameter values can be found in Sect. 3.5. The testing accuracy achieved was 99.1%, with an F1 score of 0.9868. The confusion matrix derived from the test data can be viewed in Fig. 3.

Fig. 3Full-screen View

Confusion matrix for the proposed algorithm’s test results. Each cell within the matrix’s core is normalized according to the total observations of the respective class, illustrating the proportion of correctly identified samples within the whole dataset. The column summary indicates the percentage of correct and incorrect classifications for each predicted class, scaled by the overall observations of that predicted class. Similarly, the row summary portrays the percentage of correct and incorrect classifications for each actual class, adjusted by the total observations of that specific class.

In summary, the TPW algorithm has shown promising results in both passive backgrounds and active scenarios. Analysis of the column and row summaries indicates that the most notable misclassifications occur between class 1(HEU), 5(^99mTc), and 6(HEU+^99mTc). The primary reason for this is the presence of radioactive material in the detection scene of class 6, which is also present in classes 1 and 5, causing ambiguity in the identification process.

Furthermore, to provide a comprehensive comparison, the standard KNN (KNN) [16, 35], support vector machine (SVM) [17, 36], Bayesian network (BayesNet) [18, 37], random tree (RandomTree) [19, 38], and the proposed algorithm (TPW) were applied for evaluation. The aforementioned methods were commonly utilized for radionuclide identification classification in recent years. Comparative experiments were conducted using the Weka machine learning toolkit [39] with a batch size of 100 and ten-fold cross-validation. The main parameters used for each method were as follows: For standard KNN, the number of neighbors was set to 5 with no distance weighting. For SVM, the Poly Kernel was used as the kernel function. The complexity and tolerance parameters were set to 1.0 and 0.001, respectively. For Bayesian network, the initial network structure used for learning was the Naive Bayes Network and the maximum number of parents that a node in the Bayes net can have was limited to 1. For Random Tree, the random number seed used for selecting attributes was 1, and the minimum total weight of instances in a leaf was set to 1.0. The maximum depth of the tree was unlimited.

Given the imbalanced nature of the dataset, relying solely on traditional classification accuracy can be misleading. For instance, a model might achieve high accuracy simply by categorizing all samples as the majority class, in this case, ’Background’. Hence, a variety of evaluation metrics were utilized in this subsection to provide a holistic view of model performance. Fig. 4 presents these metrics for different models. Comparing the TPW algorithm with four other methods (Standard KNN, Support Vector Machine, Bayesian Network, and Random Tree) across five distinct evaluation metrics, it became evident that TPW excels. Specifically, TPW consistently showcased superior accuracy, F1 score, MCC, ROC Area, and PRC Area when compared to its counterparts. This underlines TPW’s enhanced efficacy and reliability in tasks related to radionuclide identification.

Fig. 4Full-screen View

Multiple evaluation metrics across various models. The x-axis represents the evaluation metrics including accuracy (Acc), F1 measure (F-Measure), Matthews correlation coefficient (MCC), receiver operating characteristic area (ROC Area), and precision-recall curve area (PRC Area). The y-axis shows the performance values for each method under the corresponding metrics for every class of samples.

We further performed individual tests for each class of samples, contrasting the TPW algorithm’s performance with the four other methods using the F1 measure. Fig.5 showcases the classification results across different models for every sample class. Overall, the TPW algorithm emerged as the top performer amongst all the tested methods. In particular, it exhibited a commendable capacity to accurately classify samples from every class, underscoring its robustness and adaptability to various sample types.

Fig. 5Full-screen View

F1 measure for each class of samples across various models. The F1 measure values of different methods for each class of samples are plotted on the y-axis of the plot, while the x-axis of the plot indicates the corresponding class of samples as listed in Table 2.

Examining the results in depth, we note that the efficacy of different methods varies considerably across classes. In particular, for some classes, the TPW algorithm notably surpasses the F1 measures of its competitors, while for others, the performance differences are more nuanced. This indicates that the TPW algorithm is especially adept at processing certain sample types, although its relative advantage might be less distinct for other sample types. Notably, the F1 measure for the samples of class 0 (background), 3 (¹³¹I), and 4 (⁶⁰Co) is higher, whereas it is somewhat subdued for class 1 (HEU), 2 (WGPu), 5 (^99mTc), and 6 (^99mTc). By analyzing the peak energies of these radioactive materials, we discerned that their characteristic peaks are all below 200 keV. This suggests that their accurate detection and identification might be compromised by the Compton continuum. Nevertheless, the TPW algorithm excels over other models, surpassing them by at least 0.18% in multi-isotope scenarios like HEU+^99mTc, and showcases lower variability than other models when detecting radioactive materials.

Discussion on temporal energy window

In this subsection, we explored the effects of varying parameters associated with the temporal energy window, specifically focusing on the number and duration of these windows. The movement of the detector poses challenges in determining the ideal length for an energy window. A brief window, such as 1 s, might not capture sufficient relevant data, while an extended window, such as 20 s, could introduce substantial noise, potentially overshadowing crucial signals. Fig.6 depicts the energy spectra for a sample from class 4 (⁶⁰Co) at varying energy window durations. As the window length increased, there was a noticeable rise in the count of the energy spectrum. However, the count values across different energy points did not grow linearly with the expansion of the window length, possibly due to statistical fluctuations and other influencing factors [29, 30]. For a deeper insight into the nuances of the spectral lines, Fig.7 highlights the variations in the morphology and attributes of the energy spectra, as the window length transitions across five distinct durations.

Fig. 6Full-screen View

(Color online) Energy spectra with respect to different temporal energy window lengths. This figure illustrates the energy spectra of ⁶⁰Co for different temporal energy window lengths of 1, 3, 5, 10 and 20 s. The energy range shown is between 0–3000 keV. The count value displayed on the Z-axis shows an increasing trend with the varying window lengths.

Fig. 7Full-screen View

(Color online) Detailed energy spectra across varied window lengths. These five diagrams offer an in-depth examination of the energy spectra of ⁶⁰Co, captured at different temporal energy window durations, complementing Fig. 6. The inset images in the upper right magnify the spectrum within the highlighted red boxes. The irregular growth in count values at various energy addresses, as the window length expands, can be linked to statistical fluctuations and other potential influences.

Figure 8 illustrates the model’s performance changes in relation to varying window lengths and quantities. The data indicates that as the number of temporal energy windows increases, the model’s prediction accuracy also rises, particularly when the window quantity equals 5. With the increase in the number of temporal energy windows, the proposed algorithm captures a richer representation of the signal, enhancing the differentiation between signal and noise, and thus providing a more precise target prediction. Additionally, the model’s prediction accuracy trends upward with longer window durations. However, there is a significant decline in accuracy with excessively long windows. These experiments shed light on the influence of window lengths and quantities on classification efficacy, emphasizing the need for optimal temporal scale selection and feature extraction methods for accurate classification in the given task.

Fig. 8Full-screen View

Performance variations of models with different window lengths and quantities. The plot displays the F1 measure values for various parameter combinations, with x-axis representing the energy window lengths.

Discussion on peak-ratio spectrum analysis

In this subsection, we explored the influence of individual and combined features on the classification performance using ablation experiments. Table 3 contrasts the classification outcomes of the aggregated energy spectrum counts, peak-flat ratio, and peak-peak ratio features against those using the joint features, shedding light on their individual and combined impacts.

The F1 measure comparisons reveal that combined features outperform their individual counterparts. This superior performance of joint features arises from their capacity to seamlessly assimilate energy spectrum data from diverse viewpoints. By harmoniously harnessing the unique attributes of each feature, joint features amplify classification precision, overshadowing the results achieved by singular features. Additionally, joint features adeptly counteract challenges intrinsic to individual features, such as noise interference, data sparsity, or lack of comprehensive representation. Conversely, singular features often struggle to offer a holistic and resilient information foundation for classification [23]. Therefore, combined features furnish a more holistic and richer data representation, bolstering classification efficiency. In essence, the empirical findings underscore the merit of deploying joint features for more accurate radionuclide identification.

Discussion on classification model

In this subsection, we examined the proposed algorithm by assessing the impact of three factors: the distance metric, value of K, and distance weight. We experimented with various distance metrics, including Euclidean, City block, Chebyshev, Correlation, Spearman, Hamming, Jaccard, and Cosine. The value of K was varied between 1 and 20. For distance weighting, we considered three approaches: "Equal distance" (ED), which did not incorporate any weight; "Inverse distance" (ID), where the weight was based on the inverse of the distance to the data point; and "Inverse distance squared" (IDS), where the weight was determined by the inverse of the squared distance to the data point. The experimental results are listed in Table 4. During the experiments, K was set to 10 and all the samples were subjected to standardization.

Based on the observed F1 measure during validation and testing, it is evident that the Euclidean metric outperformed other distance metrics in classification accuracy. In terms of distance weighting, the IDS emerged as the superior performer. By significantly reducing the influence of far-off points on the classification decision, IDS led to more precise and dependable results. Collectively, the experimental data indicated that a combination of the Euclidean distance metric with the IDS weighting scheme stands as an optimal choice for the proposed algorithm when tackling classification tasks.

Selecting an appropriate value of K was crucial for optimal model performance. A very small K makes the model vulnerable to noise in the feature points, which can greatly influence classification outcomes. Conversely, an overly large K dilutes the specificity of the model as the neighborhood around the training instance becomes too expansive, increasing the likelihood of misclassifications [16, 35]. Thus, striking a balance between noise resistance and model precision by carefully adjusting the K value is imperative. We undertook a series of tests to discern the effect of different K values, ranging from 1 to 20, on the efficacy of the proposed algorithm. The outcomes of these tests are depicted in Fig. 9.

Fig. 9Full-screen View

Performance variations based on different K values. The figure showcases the F1 measure for the proposed algorithm for K values spanning from 1 to 20, plotted on the X-axis. The Y-axis indicates the F1 measure. Two distinct curves denote the outcomes from verification and test experiments, respectively.

From the results, it can be observed that the F1 score initially increases as the value of K increases from 1 to 5, reaching a peak value of 0.9868 at K=5. Then, F1 score slightly fluctuates and then starts to decline as K further increases. Additionally, the results suggest that the model has a high overall performance, with F1 scores consistently above 0.95 for all values of K. This indicates that the model is effective in accurately predicting the class labels of the input data. This pattern of results suggests that increasing the value of K can lead to better classification performance up to a certain point, beyond which overfitting may occur, resulting in a decline in performance.