Data mapping

In this subsection, the essential features of the different radionuclides are extracted from a novel perspective.

Gamma rays are the products of the de-excitation process of atomic excitation and manifest as short-wavelength electromagnetic radiation. The essence of gamma rays is a stream of particles, which are gamma photons. Gamma photons are uncharged particles, and their interaction with matter is a random event. Collecting gamma photons using a general signal acquisition device is complex, whereas collecting electrical signals is relatively easy; the amplitude of the electrical signal is proportional to the photon energy value. Gamma photons can be converted into electrical signals for signal processing to be collected by a nuclear radiation detector. The principle of this process is that photons emitted by the radioactive source interact with the atoms of the medium in the detector to produce charged particles. The detector collects the particles, converts them into electrical signals, and detects nuclear signals by measuring the electrical signals. The counts distributed with the energy value of the particles can be obtained by scaling the pulse amplitude by the energy; the energy spectrum is the curve of the distribution of the counts with the energy of the particles. As the fingerprints of radionuclides, the energy spectrum contains the distinguishable information of different radionuclides[2-4, 6].

For nuclear events, the count and time of events are random within a certain time interval. In radiation detection, the number of nuclear events measured over a certain period (e.g., detector counts) is also random. Because radioactive decay is a random process, each observation can be considered as a random experiment, and the count per unit time can be regarded as a random variable that obeys the Poisson or Gaussian distributions. For a spherical space with a single point source as the center of the sphere and a certain radius, the process of generating gamma photons by radioactive decay is random and continuous, and photons are uniformly emitted in all directions in space.

The gamma-ray energy spectrum 𝒔 is a broad stationary random vector in Euclidean space, $\text{s}=\left\{s_k\right\}\in H^l$. For the convenience of subsequent expressions, let k=0,1, ..., l-1, where sk is the k+1-th count of photons distributed over the energy value, $s_k\in \mathbb{N}$. As previously indicated, the photons generated by radioactive decay are uniformly emitted in all directions in spherical space with a single point source as the center of the sphere and a certain radius. Therefore, sk ideally positively correlates with the duration of the measurement, as formulated by Eq. (1). $s_k={\alpha }_k\cdot t$ (1) where t is the duration of the measurement and αk is a parameter over the k+1-th count value affected by the background noise of the environment, the measuring angle and distance, and the intensity of the radiation source.

Considering 𝒔 as a Markov chain or Markov process, which is a stochastic model describing a sequence of possible events in which the probability of each event depends only on the state attained in the previous event, sk is only related to s_k-1. In vector-based feature extraction methods, only part of the information of 𝒔 is considered, that is, the relative relationship between the corresponding counts of the front and rear energy addresses. This is not conducive to mining the mutual relationship between non-adjacent sk, which will cause difficulties in extracting discriminative and effective features.

Every radioactive decay produces photons of different energy values, and the photon counts of a certain energy obey the Poisson or Gaussian distributions. Thus, the total count of photons generated by radioactive decay per unit time also obeys the Poisson or Gaussian distributions. Let 𝒂 and 𝒃 be the gamma-ray energy spectra obtained by using different measurement durations in a measurement scenario with a fixed background noise, using a specified detector to measure the specified radioactive source at a fixed location and orientation. Therefore, as previously indicated, 𝒂 and 𝒃 are vectors from the same 𝑯, $\text{a}\in H^l$ and $\text{b}\in H^l$. The corresponding durations of the measurement are ta and tb, respectively. $a_k={\alpha }_k^{(a)}\cdot t_a$ (2) $b_k={\alpha }_k^{(b)}\cdot t_b$ (3) where the measuring conditions, that is, the background noise of the environment, measuring angle and distance, and intensity of the radiation source, are unchanged during the measurement, and we can assume that ${\alpha }_k^{(a)}={\alpha }_k^{(b)}$. Assuming that the measurement time is sufficient, $\frac{a_k}{b_k}$ converges in probability to $\frac{t_a}{t_b}$; that is, for ak and bk, $\forall \epsilon >0$, $P\left(\left|\frac{a_k}{b_k}-\frac{t_a}{t_b}\right|\ge \epsilon \right)=0$ (when photons are uniformly emitted in all directions).

To efficiently transfer discriminative information for identification from the gamma-ray energy spectra, mapping the gamma-ray energy spectra from a vector form in 𝑯 to a matrix form in Banach space 𝑩, the mapping f can be formulated by Eq. (4). $P=f(\text{s})$ (4) where $P={\left(p_{ij}\right)}_{m\times n}\in B^{m\times n}$, m × n=l. For the convenience of subsequent expressions, let i=0,1, ..., m-1 and j=0,1, ..., n-1, where pi j represents the element in the i+1-th row and j+1-th column in the matrix. The relationship between sk and pi j is represented by Eq. (5). $p_{ij}=p_{\lfloor k/n\rfloor \mathrm{,}\mathrm{mod}(k\mathrm{,}n)}=s_k$ (5) where k=0,1, ..., l-1.

From the matrix perspective, there are more elements adjacent to pi j. When extracting features from gamma-ray energy spectra, not only is the relative relationship between the corresponding counts of the front and rear energy addresses considered, but also the relative relationship between the upper and lower counts and the diagonal counts, which can be easily used and is more conducive to mining the internal and mutual relationships between elements in 𝑷.

To apply 𝑷 as the input of a DCNN model and transfer image descriptors as features of the gamma-ray energy spectra, it is essential to map 𝑷 to the image form. Mapping g maps 𝑷 from the matrix form in 𝑩 to the image form in the RGB color space 𝑱, can be formulated by Eq. (6). $Q=g(\text{P})$ (6) where $Q={\left(q_{ij}\right)}_{m\times n}\in J^{m\times n}$, m × n=L. For the convenience of subsequent expressions, let i=0,1, ..., m-1, j=0,1, ..., n-1, where qi j represents the element in the i+1-th row and j+1-th column in the matrix. 𝑷 is normalized before the mapping process, which can be formulated using Eqs. (7). $q_{ij}=\frac{p_{ij}-p_{\min }}{p_{\max }-p_{\min }}$ (7) where $p_{\max }={{\displaystyle \max }}_{i\mathrm{,}j}\left\{p_{ij}\right\}$, $p_{\min }={{\displaystyle \min }}_{i\mathrm{,}j}\left\{p_{ij}\right\}$.

Eq. (6) maps the element values of 𝑸 onto the corresponding pixels of an image with specified colors. Each qi j corresponds to a rectangular area in the image, and the values of qi j are indices in the Parula color map [45] that determine the color of each patch. Eq. (6) maps the smallest value in 𝑸 to the first entry in the Parula color map and maps the largest value in 𝑸 to the last entry in the Parula color map. All intermediate values of 𝑸 are linearly scaled to the Parula color map in the ascending order. The relationship between the values of the elements in 𝑸 and the colors of the corresponding pixels in the Parula color map is shown in Fig. 2.

Fig. 2Full-screen View

(Color online) Relationship between the elements and the colormap. The Parula colormap is a three-column array with 64 rows, where each row in the array defines one color using an RGB triplet, i.e., contains the red, green, and blue intensities for a specific color. Each row in the matrix defines one color using an RGB triplet. The intensities are in the range [0,1], where a value of 0 indicates no color and a value of 1 indicates full intensity.

𝒂 and 𝒃 are two gamma-ray energy spectra obtained for different measurement durations under the same measurement conditions, and Qa=g(f(_a)) and Qb=g(f(b)) are the images of 𝒂 and 𝒃, respectively, under mapping f * g which are defined as previously mentioned. For $\forall q_{ij}^{(a)}\in Q_{\text{a}}$, $q_{ij}^{(a)}$ and t_a are related, as shown in Eq. (8). $q_{ij}^{(a)}=\frac{p_{ij}^{(a)}-p_{\min }^{(a)}}{p_{\max }^{(a)}-p_{\min }^{(a)}}=\frac{s_k^{(a)}-s_{\min }^{(a)}}{s_{\max }^{(a)}-s_{\min }^{(a)}}$ (8) where k=i × m + j and $\begin{array}{l}{s}_{\max }^{(a)}=\underset{k^{\prime }=\mathrm{0,1,}\dots \mathrm{,}l-1}{{\displaystyle \max }}\left\{s_{k^{\prime }}^{(a)}\right\}=\underset{k^{\prime }=\mathrm{0,1,}\dots \mathrm{,}l-1}{{\displaystyle \max }}\left\{{\alpha }_{k^{\prime }}^{(a)}\right\}\times t_a,\\ s_{\min }^{(a)}=\underset{k^{\prime }=\mathrm{0,1,}\dots \mathrm{,}l-1}{{\displaystyle \min }}\left\{s_{k^{\prime }}^{(a)}\right\}=\underset{k^{\prime }=\mathrm{0,1,}\dots \mathrm{,}l-1}{{\displaystyle \min }}\left\{{\alpha }_{k^{\prime }}^{(a)}\right\}\times t_a.\end{array}$ Obviously, $q_{ij}^{(a)}=\frac{{\alpha }_{i\times m+j}^{(a)}-\underset{k^{\prime }=\mathrm{0,1,}\dots \mathrm{,}l-1}{{\displaystyle \min }}\left\{{\alpha }_{k^{\prime }}^{(a)}\right\}}{\underset{k^{\prime }=\mathrm{0,1,}\dots \mathrm{,}l-1}{{\displaystyle \max }}\left\{{\alpha }_{k^{\prime }}^{(a)}\right\}-\underset{k^{\prime }=\mathrm{0,1,}\dots \mathrm{,}l-1}{{\displaystyle \min }}\left\{{\alpha }_{k^{\prime }}^{(a)}\right\}}$ Similarly, for $\forall q_{ij}^{(b)}\in Q_{\text{b}}$, we have $q_{ij}^{(b)}=\frac{{\alpha }_{i\times m+j}^{(b)}-\underset{k^{\prime }=\mathrm{0,1,}\dots \mathrm{,}l-1}{{\displaystyle \min }}\left\{{\alpha }_{k^{\prime }}^{(b)}\right\}}{\underset{k^{\prime }=\mathrm{0,1,}\dots \mathrm{,}l-1}{{\displaystyle \max }}\left\{{\alpha }_{k^{\prime }}^{(b)}\right\}-\underset{k^{\prime }=\mathrm{0,1,}\dots \mathrm{,}l-1}{{\displaystyle \min }}\left\{{\alpha }_{k^{\prime }}^{(b)}\right\}}$ This is because the measuring conditions, that is, the background noise of the environment, the measuring angle, and distance are unchanged between 𝑸a and 𝑸b. If 𝑸a and 𝑸b were from the same radioactive source, then ${\alpha }_k^{(a)}$= ${\alpha }_k^{(b)}$.

Based on the aforementioned analysis, the effects of different measurement durations on 𝑸a and 𝑸b can be ignored under the same measurement conditions. More explicitly, the mappings of the gamma-ray energy spectra in the RGB color space from the same radionuclide would present nearly identical images in an ideal situation, which can reduce the intraclass differences caused by different measuring durations. Even in a real measuring situation, the information processing of mapping the gamma-ray energy spectra from a Euclidean space to a Banach space and then to an RGB color space remains to have a strong applicability and can reduce intraclass differences, thereby helping to extract the essential features from the gamma-ray energy spectra of different radionuclides.

Feature transferring

In this subsection, a DCNN model was constructed and trained on the ImageNet dataset. The mapped gamma-ray energy spectrum images were applied as inputs to the DCNN model, and the corresponding activation vectors of fully connected layers and activations from convolution layers were transferred as descriptors of images to construct a new classification model for radionuclide identification.

VGG is a widely used convolutional neural network (CNN) model proposed by Karen Simonyan and Andrew Zisserman at the University of Oxford[46]. The VGG has various configurations; instance, VGG-11, VGG-16, VGG-19, etc.. Of all the configurations, VGG-16 was identified as the best-performing model on the ImageNet dataset. The basic building block of VGG can be summarized as a stack of multiple (usually one, two, or three) convolution layers with a filter size of 3 × 3, one stride, and one padding, followed by a max-pooling layer of size 2 × 2. Different configurations of this stack were repeated in the network to achieve various depths. The number associated with each configuration is the number of layers with the weight parameters. The convolution stacks are followed by three fully connected layers, two with a size of 4096, and the last one with a size of 1000. The last layer is the output layer with a Softmax activation. A size of 1000 refers to the total number of possible classes in ImageNet. In the method proposed in this study, the structure of the VGG consists of U convolution layers, V max pooling layers, and N fully connected layers. The convolution layers use convolution kernels to convolve the input image, and the results of the convolution constitute the feature maps of the input image; thus, the local features of the image are extracted after the convolution layers. The max pooling layers are arranged after the convolution layers, and the maximum values of the corresponding positions in the feature maps are calculated. This can reduce the dimension of the extracted feature information to make the feature maps smaller, simplify the computing complexity of the network, and avoid overfitting. The U convolution layers can be divided into M groups using the max pooling layers as separations. The fully connected layers map the distributed feature representation learned from the convolution layers and max pooling layers to the sample marker space, the essence of which is to perform a weighted sum of the features to integrate the local features and output them as a value to reduce the influence of the local feature position on the classification.

The convolution layers and max pooling layers have parameters such as the height of the convolution kernel, width of the convolution kernel, number of input channels, number of output channels (number of convolution kernels), padding, and stride. The padding parameter refers to the boundary padding of the original matrix with ′0′ before the convolution; thus, the convolution kernel can extend to the pseudo-pixels beyond the edge when scanning the input image, thereby avoiding the loss of edge information. The stride parameter refers to the length of each movement of the convolution kernel, and the stride size affects the efficiency of the model. Each convolution kernel strictly has a bias parameter. To simplify the calculation, the bias is omitted in the following calculation.

For the convolution layers, the height and width of the convolution kernel are represented as $h_{\text{k}}^{\text{conv}}$ and $w_{\text{k}}^{\text{conv}}$, respectively; the number of input and output channels are represented as $c_{\text{i}}^{\text{conv}}$ and $c_{\text{o}}^{\text{conv}}$, respectively; and the parameters of the padding and stride are represented as pad^conv and str^conv. The number of neurons in the fully connected layer is represented by $c_{\text{o}}^{fc}$.

The size of the input image is represented as $c_{\text{i}}^{\text{im}}\times h_{\text{i}}\times w_{\text{i}}$, and the size of the output image is represented as $c_{\text{o}}^{\text{im}}\times h_{\text{o}}\times w_{\text{o}}$, where $c_{\text{i}}^{\text{im}}$, h_i, w_i, and $c_{\text{o}}^{\text{im}}$, h_o, w_o represent the channel, height, and width of the input and output images, respectively. The height and width of the output image of the convolution and max pooling layers can be formulated as Eq. (9), and Eq. (10). $h_{\text{o}}=\lfloor h_{\text{i}}-h_{\text{k}}+2\times \text{pad}\rfloor /\text{str}+\mathrm{1,}$ (9) $w_{\text{o}}=\lfloor w_{\text{i}}-w_{\text{k}}+2\times \text{pad}\rfloor /\text{str}+\mathrm{1,}$ (10) where the Eq. (9), and Eq. (10) can be generalized to the max-pooling layers to calculate the height and width of the convolution kernels. $c_{\text{i}}^{\text{conv}}$ of a convolution layer depends on the channel of the input image, and $c_{\text{o}}^{\text{conv}}$ of a convolution layer depends on the number of convolution kernels.

VGG is training on the ImageNet dataset to obtain the classification model PM. ImageNet is an image database organized according to the WordNet hierarchy, in which hundreds or thousands of images depict each node of the hierarchy. The dataset has been instrumental in advancing computer vision and deep learning research[47]. Let 𝑺 be a set of the gamma-ray energy spectra in Euclidean space, $\text{S}\subseteq {\text{H}}^l$, 𝒔 be the gamma-ray energy spectrum, $\text{s}\in \text{S}$. 𝒔 is converted into 𝑸 using Eqs. (4), and Eq. (6), and 𝑸 is applied as an input to the PM to transfer features. Generally, PM will only have one final output, and the corresponding activation vectors of the fully connected layers and activations are transferred from the convolution layers in the process as the image descriptors of 𝑸.

The convolution kernel of PM corresponds to a receptive field, and a small part of the image (local receptive area) is used as the input of the lowest convolution layer, which makes each neuron output by the convolution layer only experience the local image area and does not need to experience the global image. This operation is equivalent to passing through a digital filter to obtain the most salient features of the observed data. In the fully connected layers, different local features from the convolution layers are synthesized through the weight matrix to form a representation of the global information. Therefore, the transferred activation vectors of the fully connected layers are regarded as global feature representations for the gamma-ray energy spectrum images, and the transferred activation maps from the convolution layers are regarded as local features describing particular gamma-ray energy spectrum image regions.

Specifically, the activated maps transferred from the convolution layers, that is, the set of feature maps transferred from the i-th group of convolution layers of PM, is denoted as 𝒄i, $i\in [\mathrm{1,2,}\cdots \mathrm{,}M]$, which is essentially a set of multiple square matrices. 𝒄i contains d⁽ⁱ⁾ feature maps, where d⁽ⁱ⁾ depends on the number of convolution kernels in one of the convolution layers in the i-th group of convolution layers. The size of the feature map is $h_{\text{o}}\times w_{\text{o}}$, where h_o and w_o can be calculated using Eq. (9), and Eq. (10). 𝒄i can be written as: ${\text{c}}_i=\left\{{\text{c}}_1^{(i)}\mathrm{,}{\text{c}}_2^{(i)}\mathrm{,}\dots \mathrm{,}{\text{c}}_{d^{(i)}}^{(i)}\right\}\mathrm{.}$ To convert the set of feature maps into a set of feature vectors, which can be further used for classifier training, the conversion operation can be formulated as follows: $\begin{array}{l}GAP\left({\text{c}}_i\right)=\left(\text{ave}\left({\text{c}}_1^{(i)}\right)\mathrm{,}\text{ave}\left({\text{c}}_2^{(i)}\right)\mathrm{,}\dots \mathrm{,}\text{ave}\left({\text{c}}_{d^{(i)}}^{(i)}\right)\right)\mathrm{,}\hfill \end{array}$ where $GAP\left({\text{c}}_i\right)$ is the global average pooling on the set of feature maps, and $\text{ave}\left(c^{(i)}\right)$ is the average operation on all elements of the feature maps, where $\text{ave}\left({\text{c}}^{(i)}\right)=\text{ave}\left(\left[\begin{array}{ccc}{c}_{\mathrm{1,1}}^{(i)}& \cdots & c_{\mathrm{1,}{w}_{\text{o}}}^{(i)}\\ ⋮& \ddots & ⋮\\ c_{h_{\text{o}}\mathrm{,1}}^{(i)}& \cdots & c_{h_{\text{o}}\mathrm{,}{w}_{\text{o}}}^{(i)}\end{array}\right]\right)\mathrm{.}$ The transferred local features of 𝑸 are expressed by Eqs. (11). $F_L^{(Q)}=\left\{GAP\left({\text{c}}_i\right)\right\}\mathrm{,}i\in [\mathrm{1,2,}\cdots \mathrm{,}M]$ (11) where $\text{GAP}\left({\text{c}}_i\right)$ is a vector of length d⁽ⁱ⁾ and $F_L^{(Q)}$ is a set of multiple $\text{GAP}\left({\text{c}}_i\right)$, $i\in [\mathrm{1,2,}\cdots \mathrm{,}M]$, which are essentially the local image descriptors of Q transferred from PM.

The transferred activation vectors of the fully connected layers, that is, the feature vectors transferred from the j-th fully connected layer of $\text{PM}$, are denoted as 𝒇j, $j\in [\mathrm{1,2,}\cdots \mathrm{,}N]$, whose length is the same as the number of neurons in the j-th fully connected layer $c_{\text{o}}^{fc}$. The transferred global features of 𝑸 are expressed by Eqs. (12). $F_G^{(Q)}=\left\{{\text{f}}_j\right\}\mathrm{,}j\in [\mathrm{1,2,}\cdots \mathrm{,}N]\mathrm{.}$ (12) where $F_G^{(Q)}$ is a set of multiple 𝒇j, which are the global image descriptors of 𝑸 transferred from PM. Table 1 summarizes the overall flow of the algorithm.

Illustration

In this subsection, we present a simple example of the proposed method. For illustration, we selected two synthetic gamma-ray energy spectra 𝒔₁ and 𝒔₂ of ⁶⁰Co. The essence of 𝒔₁ and 𝒔₂ are two vectors in Euclidean space 𝑯l, where l = 4096. The simulation settings of 𝒔₁ and 𝒔₂ are essentially the same, and the only difference is the total number of simulated particles, which is equivalent to the difference in the measurement duration. The number of simulated particles correlates with the duration of the measurement because the process of generating gamma photons by radioactive decay is random and continuous, and photons are uniformly emitted in all directions in space. The simulated particle numbers corresponding to 𝒔₁ and 𝒔₂ are represented as t₁ and t₂, where t₁ is five hundred thousand and t₂ is two hundred and fifty thousand. To illustrate the difference between the two spectra more vividly, we plot 𝒔₁ and 𝒔₂ in Fig. 3, where the blue spectrum is the description of 𝒔₁, and the red spectrum is the description of 𝒔₂. Figure. 3 clearly shows differences between 𝒔₁ and 𝒔₂ due to the duration of measurement, where the counts distributed over the energy addresses of the entire spectrum is significantly more prominent in 𝒔₁ than 𝒔₂.

Fig. 3Full-screen View

(Color online) Two original spectra of ⁶⁰Co. The simulation settings of 𝒔₁ and 𝒔₂ are essentially the same, and the only difference is the total number of simulated particles t₁ and t₂, which is equivalent to the difference in the measurement duration.

To extract more discriminative image descriptors from the gamma-ray energy spectra for identification, 𝒔₁ and 𝒔₂ are mapped to the matrix form using Eq. (4). The mapped matrices are represented as 𝑷₁ and 𝑷₂, which are in the Banach space $B^{m\times n}$, where m = n = 64. To reduce the discrepancy between 𝑷₁ and 𝑷₂ caused by the difference in the measurement duration, and to facilitate the matrices as the inputs of the DCNN for the transfer of image descriptors, 𝑷₁ and 𝑷₂ are mapped to the image form using Eq. (6), and the results are represented by 𝑸₁ and 𝑸₂. The gamma-ray energy spectra of ¹³⁷Cs and ¹⁵²Eu were randomly selected for comparison, and their corresponding Q_Cs and Q_Eu values were obtained using Eq. (4). The intraclass similarity and interclass distinction of the different gamma-ray energy spectra are shown in Fig. 4 vividly. Specifically, 𝑸₁ and 𝑸₂ exhibited nearly identical images in Fig. 4(a) and (b), which indicates that the difference between 𝒔₁ and 𝒔₂ in Fig. 3 is significantly reduced when 𝒔₁ and 𝒔₂ are mapped to RGB color space. By comparing Fig. 4(c), (d) with Fig. 4(a), (b), (c), (d) clearly have completely different characteristics than Fig. 4(a), (b). Namely, the mapping of the gamma-ray energy spectra of the randomly selected ¹³⁷Cs and ¹⁵²Eu in the RGB color space exhibits completely different characteristics from those of ⁶⁰Co. Through the aforementioned comparison and analysis in Sect. 2, the effects of different measuring durations on 𝑸₁ and 𝑸₂ can apparently be ignored under the same measuring conditions. Therefore, the information processing of mapping the gamma-ray energy spectra from Euclidean space to Banach space and then to the RGB color space remains to have a strong applicability and can reduce intraclass differences, thereby helping to extract the essential features from the gamma-ray energy spectra of different radionuclides. A further quantitative analysis is presented in Sect. 3.

Fig. 4Full-screen View

(Color online) Gamma-ray energy spectrum images. (a) and (b) are the images of 𝑸₁ and 𝑸₂, (c) and (d) are the images of Q_Cs and Q_Eu.

In the feature transfer phase, the structure of VGG and the corresponding feature sizes are listed in Table 2, which consists of five groups of convolution layers and two fully connected layers. The layer number only considers the number of convolution layers, and the last layer of each group is the max-pooling layer. The number of channels in the feature size is not changed in the convolution layers, but changes in the max pooling layer owing to the size of the stride parameter (the number of output channels in different groups can be calculated by Eq. (9), and Eq. (10)). For fully connected layers, the number of output channels is equal to the number of neurons in the fully connected layers. The partial parameters of the convolution, max pooling, and fully connected layers are shown in Table 3.

The VGG framework was trained on the ImageNet dataset to obtain the classification model PM. 𝑸₁ and 𝑸₂ are applied as input images for the PM to transfer features, and the activation vectors of the fully connected layers and activations from the convolution layers in the process are transferred as the image descriptors of 𝑸a and 𝑸b. Image descriptors transferred from the fifth group of convolution layers 𝒄₅ and the first fully connected layer 𝒇₁ were selected as illustrations. 𝒄₅ and 𝒇₁ from the synthetic and measured datasets are shown in a low-dimensional space through the t-SNE[48] in Fig. 7. Each color in the figure represents one type of radionuclide, which intuitively reflects the difference in the expressive and discriminative abilities between the different features. The transferred image descriptors have a strong intra-class similarity and inter-class differentiation, and a further analysis of the feature performance is presented in the next section.

Fig. 7Full-screen View

(Color online) Visualization of various features by t-SNE. Each color in the figure represents one type of radionuclide, which intuitively reflects the difference in the expressive ability and discrimination ability between different features.

Data preparation

The production of the synthetic dataset consisted of the following three steps: (1) Background data acquisition. A self-made 3-inch NaI detector was used to measure the ambient background. The measurement was conducted in a laboratory environment without the presence of a separate radioactive source. The detector was then placed at a fixed position for 12 h. Two measurements were performed, one with lead bricks placed around the detector and the other without; two sets of background data were obtained. (2) Single-nuclide energy-spectrum acquisition. Based on the Geant4 platform, the transport process of the radioactive gamma-ray particles of 26 single nuclides was simulated using the Monte Carlo method. We constructed simulated scenarios in Geant4 containing only different types of single radioactive point sources and a 3-inch NaI detector. The positions and relative distances of the radioactive point sources and detector were fixed. A total of one million photons were simulated, and their trajectories were recorded by the detector and converted into the gamma-ray energy spectrum. (3) Data synthesis. To simulate real measurements as closely as possible, two sets of background data and 26 sets of synthetic spectra were linearly superimposed with a random signal-to-noise ratio (SNR). $\text{SNR}=N_{\text{nc}}/N_{\text{bg}}$, where N_nc is the sum of the photon counts emitted by the radioactive point source and N_bg is the sum of the photon counts of the background. In the linear superposition process, the SNR value is a random number ranging between 0.3-1. A total of 2080 synthetic gamma-ray energy spectra of 26 radionuclides were obtained and named as dataset 1. Table 4 demonstrates a list of common radionuclides, which includes a total of 26 radionuclides in the following four categories: SNM, industrial, medical, and NORM.

The measured gamma-ray energy spectra were obtained from radioactive sources in a laboratory environment with lead brick shielding. The measurements took advantage of a cadmium zinc telluride (CZT) gamma-ray spectrometer from Kromek. The spectrometer has 4096 measuring channels; the measurable energy ranges between 25 keV-3.0 MeV, and the electronic noise is lower than 10 keV. Two types of V radiation sources, ¹³⁷Cs and ⁶⁰Co, and one type of IV radiation source, ¹⁵²Eu were used in the measurement. The spectrometer was carried on the Turtlebot robot to quantitatively control the measuring distance and reduce radiation damage to the experimental operators. As shown in Table 5, seven groups of samples were established; the measuring object, measuring distance, and measuring duration were varied in the process. A total of 150 measured gamma-ray energy spectra of three single radiation sources and a total of 200 measured gamma-ray energy spectra of four mixed radiation sources were obtained and named dataset 2. Figure 5 presents the spectrum of ⁶⁰Co in different sample sets. Figure 6 presents the gamma-ray energy spectrum images in different sample sets.

Fig. 5Full-screen View

Examples of synthetic and measured spectra. (a) displays a ⁶⁰Co synthetic spectrum, and (b) displays a real measured ⁶⁰Co spectrum.

Fig. 6Full-screen View

(Color online) Examples of synthetic and measured gamma-ray energy spectrum images. Group(a) are synthetic gamma-ray energy spectrum images of 26 radionuclides, group(b) are measured gamma-ray energy spectrum images of 7 single and mixed radionuclides.

Feature performance comparison and analysis

Research have shown that image descriptors transferred from DCNNs can provide a reliable performance for image classification problems. However, choosing different features for a specific classification domain remains worth discussing. For the domain of the radioactive gamma-ray energy spectral classification presented in this study, and considering the limitation of the computing resources and the scale of parameters, VGG-16 was chosen as the DCNN framework owing to the characteristic that VGG-16 applies a significantly small 3 × 3 receptive field (filters) throughout the entire network with the stride of 1 pixel. A combination of multiple 3 × 3 filters and nonlinear activation layers can replace a receptive area of a larger size, which makes the decision functions more discriminative to the characteristics of the spectrum. This imparts the ability of the network to converge faster[46]. In addition, the consistent use of 3 × 3 convolutions across the network makes the network significantly simple, elegant, and conveniently transfers image descriptors.

In this subsection, 28 classification methods based on the Weka machine-learning toolkit [49] were applied in a series of experiments to verify whether image descriptors transferred from a DCNN model can be used to construct a classification model with a strong discrimination and advanced accuracy. Owing to the large number of transferred features (five sets of local features and two sets of global features for a single energy spectrum), in an actual training process, using only a certain set of features is sufficient to build a satisfactory classifier. Therefore, in the experiments described in this subsection, single sets of local and global features and combinations of high-level features were used to conduct the experiments. Multiple groups of transferred image descriptors from the synthetic and measured datasets were applied in the training and testing processes of the 28 classification models. All classification models were trained with the default parameters and settings specified in the toolkit, and the experiments applied the ten-fold cross-validation method to avoid the imbalance caused by random data segmentation. The percentage of misclassified cases for transferred features on the synthetic dataset and measured dataset are listed in Table 6 and Table 7. The results demonstrate that although the gamma-ray energy spectrum images are completely unfamiliar to the DCNN model and have not been used in the pre-training process, the transferred image descriptors achieved good classification results.

For local image descriptors, a higher group of convolution layers can provide more discriminative descriptors than a lower group of convolution layers. Part of the multigroup union image descriptors can provide better resolution, and different arrangements of image descriptors from the same groups also have an impact on the results. The best average classification effect of the local image descriptors on the synthetic dataset was the multi-group union image descriptors composed of c₃, c₄, and c₅, which was 93.08%. Local image descriptors achieved the best average classification effect on the measured dataset c₅, which was 93.51%. The best classification effects of the local image descriptor on the synthetic and measured datasets were 100.00% and 99.71%, respectively.

For the global image descriptors, 𝒇₁, which is transferred from the first fully connected layer, has strong semantic information, achieving an average accuracy of 92.76% and 94.86% on the synthetic dataset and measured dataset, respectively. The best classification effects of the global image descriptors on the synthetic and measured datasets were both 100.00%. In contrast, the classification model trained on the global image descriptors transferred from the second fully connected layer achieved a poor classification performance. The aforementioned experiments have preliminarily proved that image descriptors transferred from a DCNN model can be applied as the image descriptors of radionuclide gamma-ray energy spectrum images for classification; however, further comparisons of the features of gamma-ray energy spectra and different image descriptors of the gamma-ray energy spectrum images are needed.

Statistical comparison

Characteristic peaks are the most important features of traditional radionuclide identification methods. Four groups of characteristic peaks were extracted from the gamma-ray energy spectra by changing the peak properties, that is, distance, prominence, width, and threshold. Scale-invariant feature transform (SIFT)[50] and a histogram of oriented gradients (HOG)[51] are two classical feature extraction algorithms used in computer vision. Two groups of features extracted from the gamma-ray energy spectrum images were obtained using the HOG and SIFT algorithms. Fig. 7 demonstrates various features of the synthetic and measured data sets in a low-dimensional space through the t-SNE[48], which including characteristic peaks, features extracted by HOG and SIFT, present local image descriptors 𝒄₅ and global image descriptors 𝒇₁. As shown in Fig. 7, the distribution of the characteristic peaks is chaotic, and no apparent clustering center is observed for the different radionuclides. While HOG and SIFT features can form apparent clustering centers in certain radionuclides, the features of these radionuclides exhibit a significant crossover. Simultaneously, the transferred local and global image descriptors were significantly enhanced in feature discrimination, and the aggregation between similar feature points was also stronger.

The aforementioned features from one- and two-dimensional spectral images were applied in 21 classification methods using the Weka machine learning toolkit[49]. All experiments applied the ten-fold cross-validation method to avoid the imbalance caused by random data segmentation. The proportions of misclassified samples in the experiments are listed in Table 8 and Table 9. These results provide a basis for evaluating the performance of the transferred features. However, from a statistical point of view, these results do not provide strong support for identifying a group of features. For a statistical comparison of features, the non-parametric test method recommended by Demšar[52] was considered for comparing multiple features using several classification methods. First, we tested whether there were significant differences among the seven groups of features. The Friedman test was then applied to compare several features of the multiple classification methods.

Assuming that Num types of classification methods and k groups of features are involved in this experiment, the implementation of the Friedman test assigns a value $r_i^j$ to each group of features, which indicates the rank of the feature j on the i-th classification method. The average rank R of feature j was computed using Eq. (13). $R_j=\frac{1}{\text{Num}}{\displaystyle \sum _i}{r}_i^j\mathrm{.}$ (13) The null hypothesis was established, which states that all features have the same statistical performance. The Friedman statistic is asymptotically χ² distributed with k-1 degrees of freedom and is computed using Eq. (14). ${\chi }_F^2=\frac{12\text{Num}}{k\left(k+1\right)}\left({\displaystyle \sum _j}{R}_j^2-\frac{k{\left(k+1\right)}^2}{4}\right)$ (14) The null hypothesis is rejected if the size of ${\chi }_F^2$ exceeds the critical value, which indicates that there is a statistically significant difference between the classifiers. Conversely, the null hypothesis is accepted if the size of ${\chi }_F^2$ does not exceed the critical value. A post-hoc test was applied to determine the nature of the differences and compare the relative performances of different features when the null hypothesis was rejected.

In our experiments, the significance level α value was set to 0.05, the corresponding critical value of dataset 1 was calculated as 2.2541 while ${\chi }_F^2$= 44.8442, and the corresponding critical value of dataset 2 was calculated as 2.2541 while ${\chi }_F^2$= 2.6218. The result of the Friedman test indicated that the null hypothesis was rejected; therefore, further Holm tests [52] were required to compare the performance of the different features.

The test statistics of the Holm test for comparing features in a pair-wise manner are formulated in Eq. (15). $\text{z}=\left(R_i-R_j\right)/\sqrt{\frac{k\left(k+1\right)}{6\text{Num}}}$ (15) Using the z value obtained from Eq. (15), the corresponding probability p can be determined from the normal distribution table and is compared with the appropriate α to determine whether the hypothesis is rejected, that is, whether the proposed model is better than a certain classifier. Let p₁, p₂, ..., denote the ordered corresponding probability p-values; thus, p from all the features have the following relationship: p₁ ≤ p₂ ≤ ... ≤ p_k-1. The step-down procedure of the Holm test begins with the most significant p-value and compares each pi with the adjusted α, which is calculated using α/(k-i).

If the adjusted α corresponding to the feature, α/(k-1), is higher than p₁, the corresponding null hypothesis is rejected and indicates that both features have the same performance, and subsequently proceeds to compare p₂ with the corresponding adjusted α, α/(k-2). If the second null hypothesis remains rejected, the experiment continues to test the third null hypothesis, and so on. Once a certain null hypothesis is accepted, all remaining null hypotheses are preserved.

In our case, 𝒇₁ was chosen as the representative of the transferred image descriptors and compared with other groups of features. For data set 1, $R_{f_1}$ = 1.1818, $R_{\text{peaks}-\text{1}}$ = 2.0909, $R_{\text{peaks}-\text{2}}$ = 3.3636, $R_{\text{peaks}-\text{3}}$ = 3.9091, $R_{\text{peaks}-\text{4}}$ = 4.7727, R_{H O G} = 5.5455, R_{S I F T} = 7.0000. For data set 2, $R_{f_1}$ = 1.0000, $R_{\text{peaks}-\text{1}}$ = 3.5455, $R_{\text{peaks}-\text{2}}$ = 3.2273, $R_{\text{peaks}-\text{3}}$ = 4.5000, $R_{\text{peaks}-\text{4}}$ = 4.1818, R_{H O G} = 5.3636, R_{S I F T} = 5.2273. α = 0.05, Num = 11, and k = 7, and the standard error was $SE=\sqrt{\frac{7\times (7+1)}{6\times 11}}$=0.9211. The results of the Holm test are listed in Table 10.

For Dataset 1, the Holm procedure rejects the first, second, third, fourth, and fifth hypotheses because the corresponding p values are smaller than the adjusted α. Thus, the final hypothesis cannot be rejected. This indicates that the transferred image descriptors perform significantly better than the characteristic peaks searched by prominence, width, and threshold, and features extracted by HOG and SIFT at the significance level of α=0.05. The transferred image descriptors were not significantly better than the characteristic peaks searched for by distance. For Dataset 2, the Holm procedure rejects all the hypotheses, which indicates that the transferred image descriptors perform significantly better than the characteristic peaks and features extracted by HOG and SIFT at the significance level α = 0.05. The filtering basis of searching the characteristic peaks by using the distance removes peaks with smaller distances by defining the minimum horizontal distance between adjacent peaks until all remaining peaks satisfy the distance condition. This peak-finding method is effective for the synthesized gamma-ray energy spectra using the NaI detector and has a larger classification error in the measured gamma-ray energy spectra using the CZT detector. Because the resolution of the NaI detector is not excellent, and the synthetic energy spectra comes from an ideal environmental simulation, the simulated energy spectrum has fewer interference peaks, while the measured gamma-ray energy spectra fluctuate significantly.

The results of the aforementioned comparative experiments prove that image descriptors transferred from DCNNs are better than the characteristic peaks extracted from the gamma-ray energy spectra and the HOG and SIFT shallow features extracted from the gamma-ray energy spectrum images. These transferred image descriptors are essential, discriminative, and can provide a reliable performance for gamma-ray energy spectrum image classification problems.