Generation of theoretical continuum

The measured gamma-ray spectrum is broadened owing to the statistical fluctuations of either light in the scintillation detector or electron-hole pairs in the semiconductor detector [18]. However, a spectrum without broadening can be synthesized using Monte Carlo (MC) simulations. When a simulated spectrum without broadening is generated, each peak assumes a single channel, thereby simplifying the removal of net peak counts. The remaining spectrum can then be manually broadened, and the theoretical continuum of the corresponding broadened spectrum can be obtained. The theoretical continuum is obtained by the procedure illustrated in Fig. 1 and is described below.

Fig. 1Full-screen View

Generation of theoretical continuum

Let $x=\left[x_1\mathrm{,}{x}_2\mathrm{,}\dots \mathrm{,}{x}_n\right]$ be a simulated gamma-ray spectrum and $x^{\text{'}}=\left[x_1^{\text{'}}\mathrm{,}{x}_2^{\text{'}}\mathrm{,}\dots \mathrm{,}{x}_n^{\text{'}}\right]$ be its corresponding spectrum without broadening, where n is the number of channels. If a peak exists at channel k in x′, the net peak counts are removed by replacing the counts of channel k with the average counts of its adjacent left and right channels, as follows: $\{\begin{array}{c}{x}_k^{\text{'}}=\frac{x_{k-1}^{\text{'}}+x_{k+1}^{\text{'}}}{2}1<k<n\\ x_k^{\text{'}}=x_2^{\text{'}}k=1\\ x_k^{\text{'}}=x_{n-1}^{\text{'}}k=n\end{array}$ (1) If x′ has multiple peaks, the above calculation is applied to each one. To obtain x^b according to the broadening function applied in the MC simulation, x′ is manually broadened. Consider a common Gaussian broadening function given by $E_{\text{p}}=E_{\text{d}}+\sigma X_{\text{f}}\mathrm{,}$ (2) where E_d and E_p are the deposited energies before and after broadening during simulation, respectively; σ is the standard deviation of E_d; and X_f is a Gaussian random number. Manual broadening proceeds as follows [19]. $x_i^b={\displaystyle \sum _{j=1}^n\frac{x_j^{\text{'}}}{\sqrt{2\pi }\sigma w}}{e}^{-\frac{{\left(E_i-E_j\right)}^2}{2{\sigma }^2}}\mathrm{,}$ (3) where $x_i^{\text{b}}$ represents the counts of channel i in the broadened spectrum x^b; Ej and Ej are the energy values of channels i and j, respectively; and w is the channel width expressed in terms of energy. The resulting x^b corresponds to the theoretical continuum of the gamma-ray spectrum x.

CNN for continuum estimation

CNNs are the most common architecture used in deep learning [20-24]. Compared to a fully connected network, a CNN can extract multiscale features over multiple convolutional layers and prevent overfitting through parameter sharing when treating high-dimensional data in computer vision and other areas [25-27]. Considering a spectrum with hundreds or thousands of channels as the input and a predicted continuum over the entire spectral range as the output, high-dimensional data are involved at both ends. Hence, a CNN is preferable to a fully connected network for extracting distinctive and stable shape features from the spectrum and relating them to a continuum. Additionally, the convolution operation is highly effective for handling data with localized correlations or local information, as has been extensively demonstrated in image processing, where 2D local correlations are prevalent. Therefore, we propose the incorporation of CNN to harness their inherent ability to exploit the 1D local correlations exhibited in the spectra.

The proposed CNN is a modified version of ResNet-50 [28], the architecture of which is shown in Fig. 2. It consists of an input layer, multiple residual modules, a fully connected layer, and an output layer. The input spectrum is arranged in 1024 channels, a configuration commonly considered for a low-resolution detector. Spectra with different numbers of channels can be matched to the input via channel splitting or merging. Residual modules prevent the vanishing gradient problem during deep learning by establishing skip connections from the output of the front layer to the subsequent outputs across a convolutional layer [29, 30]. Two types of residual modules, denoted as R1 and R2, are used, and their architectures are shown in Fig. 2a and b, respectively. Module R1 is characterized by four parameters: R1($D_{\text{in}}\mathrm{,}{D}_{\text{out}}\mathrm{,}C{h}_{\text{in}}\mathrm{,}\text{ and }C{h}_{\text{out}}$), where D represents the data dimension, Ch represents the number of convolution kernels, and the subscripts in and out indicate the input and output, respectively. Module R2 is characterized by two parameters: R2(D and Ch), where D and Ch apply to both input and output data. In Fig. 2a and b, Conv1D(k, s) represents a 1D convolutional layer with a convolution kernel width k and stride s, and Ch₁, Ch₂, and Ch₃ represent the number of channels of the corresponding convolution kernels. In R1, $C{h}_1=C{h}_2=C{h}_{\text{out}}/4$ and $C{h}_3=C{h}_{\text{out}}$, and in R2, $C{h}_1=C{h}_2=Ch/4$ and $C{h}_3=Ch$. BN denotes the batch normalization applied to the batch training data, and ReLU denotes the rectified linear unit (ReLU) activation in each layer. ReLU activation is employed to introduce nonlinearity, thereby enhancing the mapping capability of the CNN and ensuring non-negativity in each channel of the continuum in the final layer. Each module R1 reduces the data dimensions by one-fourth and quadruples the number of convolution kernel channels (except for the first module), whereas R2 maintains the two parameters. Through four R1–R2 blocks, the spectral features are finally embedded into a 16 × 256 vector and mapped onto the continuum through the last fully connected layer. Because modules R1 and R2 have three convolutional layers, the entire CNN contains 51 layers, including the input and output layers, and more than 10⁵ parameters.

Fig. 2Full-screen View

Architectures of residual modules (a) R1, (b) R2, and (c) proposed CNN

Test setup

First, a laboratory test experiment was conducted using available equipment. The setup involved an in-house NaI-type WBC to measure the radioactivity from the human body. WBC is commonly used in occupational radiation monitoring in nuclear facilities (e.g., nuclear power plants) to determine the category and activity of radionuclides inside an exposed human body by detecting the emitted gamma-ray spectrum [31, 32].

This setup is suitable for evaluating the continuum estimation for the following reasons: First, owing to the limited energy resolution of the NaI detector, the peaks in its spectrum are strongly broadened, leading to a wide continuum that is more difficult to estimate accurately than narrow continuums. Second, the NaI detector used in the WBC is larger than that used in other devices. Specifically, the size of the WBC used in this study was 7.6 cm×12.7 cm×40.6 cm, which resulted in multiple Compton scattering events inside the detector and a much higher continuum within the peak region than that obtained from a small detector. Clearly, a higher continuum requires better estimation to obtain accurate peak net counts. Third, the human body with radionuclides is a large volumetric radioactive source, and the emitted gamma rays are considerably scattered before detection. Consequently, more continuum counts are recorded in the low- and middle-energy regions of the spectrum. Meanwhile, the detector in a WBC is usually well-shielded by thick stainless steel and lead, causing severe backscattering of gamma rays. Combining the two abovementioned scattering mechanisms, the continuum differs completely from that observed when measuring a simple point source without shielding. Fourth, radionuclides in the human body have low activity; thus, the peak-to-continuum count ratios are lower than those formed by a strong source. This increases the importance of accurate continuum estimation.

Dataset construction

The digital model of the setup described in Sect. 2.3 was constructed using a human body represented by a human phantom (see Fig. 3c). Accordingly, an MC simulation was conducted to generate a dataset for training and testing the CNN using the GEANT4 code [33-35].

Fig. 3Full-screen View

In-house (a) WBC, (b) human physical phantom, and (c) corresponding digital model constructed in GEANT4

Nine common radionuclides used in routine internal exposure monitoring in nuclear power plants were selected to simulate the spectra, as detailed in Table 1. Per simulation run, a random selection of one to five radionuclides was made from a given set. The activity of each radionuclide was then randomly assigned a value ranging from hundreds to thousands of becquerels, based on realistic activity levels. Simulated source particles were then included with energies equal to the gamma-ray theoretical energies and quantities equal to the radionuclide activities multiplied by the gamma-ray branching ratios for 5 min, which is a typical acquisition time in practice. In addition, the source particles were uniformly distributed in the phantom, which is consistent with real conditions. The real spectrum-broadening function in Eq. (4) was used during the simulation, and the gamma-ray energies shifted slightly from their theoretical values (as presented in Eq. (5)) to mimic the temperature shifts that occur in NaI detectors. $\sigma =-3.46+0.98\sqrt{E_{\text{d}}}$ (4) where σ and E_d are both expressed in keV according to the definitions provided in Eq. (2). $E_{\gamma }^{\text{s}}=E_{\gamma }+0.022E_{\gamma }\xi \mathrm{,}$ (5) where Eγ and $E_{\gamma }^{\text{s}}$ are the gamma-ray energies before and after shifting, respectively; ξ is a random number between -1 and 1; and the value of 0.022 is also determined through measurements.

We obtained a training set with 10⁶ spectrum samples, 20% of which were reserved as test samples. The corresponding continuums were generated using the method described in Sect. 2.1.

Evaluation measures

Instead of directly assessing the error in the estimated continuum counts, the proposed method was evaluated more intuitively by comparing the deduced radionuclide activity values with the theoretical values defined during the simulation. The activity relative error (AcE) and mean AcE (MAcE) obtained from the test set were used for evaluation.

Consider a peak region that includes n channels: The activities of the radionuclides were determined as follows: $A=\frac{S}{\epsilon T\eta }=\frac{{\displaystyle {\sum }_{i=1}^n\left(Y_i-C_i\right)}}{\epsilon T\eta },$ (6) where S is the sum of the peak net counts; Yi and Ci are the total counts and estimated continuum counts in channel i, respectively; n is the number of channels in the peak region; ε is the photoelectric efficiency determined by the MC simulation; T is the acquisition time (5 min for a test spectrum); and η is the branch ratio. For a multiplet, S is determined per peak using nonlinear least-squares fitting, similar to the fitting method for continuum estimation. When multiple peaks are involved for one radionuclide, its activity is given by the weighted average of the activity across the peaks, as demonstrated below. $\overline{A}=\frac{{\displaystyle {\sum }_{p=1}^P\frac{A_p}{{\sigma }_{A_p}^2}}}{{\displaystyle {\sum }_{p=1}^P\frac{1}{{\sigma }_{Ap}^2}}}$ (7) where Ap and ${\sigma }_{A_p}^2$ are the activity and its uncertainty estimated based on peak p, respectively, and P is the total number of peaks of this radionuclide.

${\sigma }_{A_p}^2$ is determined by the error propagation based on Eq. (7). Because of the challenge of accurately evaluating the relative error of Ci, the Poisson distribution was utilized as an approximation to simplify the distribution of Ci. Consequently, ${\sigma }_{A_p}^2$ can be expressed as follows: ${\sigma }_{A_p}^2=\frac{{\displaystyle {\sum }_{i=1}^n\left(Y_i+C_i\right)}}{\epsilon T\eta }\mathrm{.}$ (8) Based on the activity of each radionuclide, the AcE_j and MAcE are defined as follows: ${\text{AcE}}_j=\frac{{\widehat{A}}_j-A_j}{A_j}\mathrm{,}$ (9) $\text{MAcE}={\displaystyle \sum _{j=1}^m\frac{\left|{\widehat{A}}_j-A_j\right|}{m{A}_j}}\mathrm{,}$ (10) where ${\widehat{A}}_j$ and Aj are the estimated and theoretical activity values of radionuclide j, respectively; and m is the number of radionuclides in the test spectrum samples.

Comparison methods

To demonstrate the high performance of the proposed method, we compared it with three existing continuum estimation methods—the fitting, step-type, and peak erosion methods—on the same training and test sets.

The fitting and step-type methods are only applicable to the peak regions. Thus, performance evaluation was limited to the peak regions of each spectrum. The peak region was defined as the range from left to right of the peak centroid with a full width at half maximum of 1.5. Overlapping peak regions formed by adjacent peaks were treated as a single region.

Details of the comparison methods can be found in corresponding studies, and we provide brief descriptions for convenience.

Validation experiment

After testing, a validation experiment was conducted by measuring a human physical phantom using WBC (Fig. 3). The in-house phantom contained uniformly distributed ¹³⁴Cs, ¹³⁷Cs, ⁵⁷Co, and ⁶⁰Co with the known activity listed in Table 2. We performed 100 repeated measurements, and the average activity of each radionuclide estimated by the proposed method was validated by comparing it with the true value and the results of the three comparison methods. The required detection efficiency for each gamma ray was determined in the same manner as in the test step.

Test

The AcE distribution and MAcE of each method are shown in Fig. 4 and Table 3. The AcE of the proposed method is within ±3% for all test samples, leading to the smallest MAcE of 1.5%. The fitting method provides a relatively small AcE of -6% to 10% under most conditions but shows some outliers up to -40%–60%, resulting in an MAcE of 5.5%. The step-type method provided an AcE within ±8%, achieving the second-best results among all the methods with an MAcE of 3.1%. The peak erosion method had the worst performance, with its AcE ranging from -20% to 90%, resulting in the largest MAcE of 18.2%.

Fig. 4Full-screen View

AcE distribution of proposed (CNN), fitting (Fit), step-type (Step), and peak erosion (Erosion) methods

Typical test scenarios

During testing, the performance of the compared methods was evaluated for three typical scenarios (see Fig. 5).

Fig. 5Full-screen View

Estimation results of each method on singlet of ¹³⁷Cs ((a) full estimation, (b) magnified view, (c) theoretical vs. estimated counts), on singlet of ⁵⁷Co ((d) full estimation, (e) magnified view, (f) theoretical vs. estimated counts), and on multiplet of ¹³⁴Cs and ¹³⁷Cs ((g) full estimation, (h) magnified view, (i) theoretical vs. estimated counts)

Validation

The results obtained from the measurements listed in Table 4 show the largest AcE of -5.1%, -55.1%, -8.3%, and ~99.1% for the proposed, fitting, step-type, and peak erosion methods, respectively. Similar to testing results, the proposed method provides the best estimation, whereas the step-type method yields the second-best results with AcE of less than 10% for the four radionuclides, and the fitting method provides the third-best results with good estimation for radionuclides ¹³⁴Cs, ¹³⁷Cs, and ⁶⁰Co but poor estimation for ⁵⁷Co. Additionally, the peak erosion method again yields the worst results with its high AcE of 99.1% for ⁵⁷Co.

Test results

The AcE distributions and MAcE values reported in Sect. 3.2 can be explained by the limitations of existing methods and the advantages of the proposed method.

The fitting method involves multiparameter optimization, which is highly nonlinear, nonconvex, and sensitive to parameter initialization. Hence, this method can easily fail for a wide and complex continuum when performed automatically without manual adjustments, resulting in considerable errors as shown in Fig. 4. In addition, the estimation is determined by a fitting function that limits the representable continuums. Overall, the fitting method is unstable for continuum estimation in low-resolution gamma-ray spectra because its unpredictable results may be suitable under simple conditions but unacceptable for high and complex continuums.

The step-type method can describe continuum counts within a peak region formed by multiple Compton scattering events of the concerned gamma ray but cannot determine the background counts. The detailed theory of this method can be found in existing literature [40, 41]. This method can only provide a continuum decline from left to right across the peak region and correct results when the background counts are negligible. However, it deviates when the continuum shape changes significantly from the ideal step curve for low-energy singlets or complex multiplets.

The peak erosion method considers the convexity of the peak structure and generates a relatively flat curve across the entire spectrum. However, when the continuum is convex, a large estimation error is observed. Moreover, the generated curve exhibits a random shape and cannot represent the details of a real continuum. Consequently, this method exhibits the worst performance in most scenarios.

Unlike existing methods, the proposed method estimates the continuum using global spectrum features extracted by a CNN, which provides small-scale count variance and other statistical characteristics, as well as large-scale count correlation and shape characteristics over the entire spectrum. Thus, it outperforms the fitting and step-type methods, which use limited local counts within peak regions, and the peak erosion method, which uses counts on each side of the concerned channel within 1.5 times the full width at half maximum. In fact, over a spectrum, the continuum of a local region is highly related to the counts in other regions; however, this relationship is too complex and nonlinear to be modeled by conventional methods. In contrast, a high-performance CNN is suitable for complex nonlinear mapping. By linking the primary spectrum to its continuum via multiple convolutional layers, the photoelectric peak, Compton scattering content, background radiation, backscattering counts, and other components are integrated by the CNN for prediction, establishing an end-to-end continuum estimation without any explicit regression or additional data processing. Therefore, the proposed method is convenient and accurate.

Analysis of selected scenarios

The limitations of the existing methods and the advantages of the proposed method can be further demonstrated by considering the scenarios detailed in Sect. 3.2.

4.2.2

Singlet of ⁵⁷Co on high-energy background

The simulated spectrum for the scenario reported in Sect. 3.2.2 without ⁵⁷Co is shown in Fig. 6a and Fig. 6b. The counts from the scattered gamma rays of 661.7 keV, 1099.2 keV, and 1292.6 keV form the background of the peak of ⁵⁷Co at 85.3-158.7 keV. When added to the original step-type continuum formed by multiple Compton scatterings of a 122.1 keV gamma ray, the additional background counts substantially change the final continuum shape. In Fig. 6b, the background curve fluctuates at 85.3-120 keV and drops sharply afterward. Thus, the continuum of ⁵⁷Co initially experiences a slower decline than expected, followed by a rapid decline. This indicates poor performance of the step-type method, as shown in Fig. 5f, which first underestimates and then overestimates the continuum.

Fig. 6Full-screen View

Simulated spectrum of Fig. 5d without ⁵⁷Co: ((a) full and (b) magnified views) of Fig. 5g with broadening and ((c) full and (d) magnified views) without broadening

Owing to least-squares optimization, the fitting method can adjust its shape more flexibly than the step-type method and thus displays a better result under this scenario when the correct fitting is obtained. However, a large error in the peak erosion method is evident.

Existing methods fail to predict the background of ⁵⁷Co contributed by ¹³⁷Cs and ⁵⁹Fe, thereby limiting their performance when multiple radionuclides are involved. In contrast, the proposed method bridges the counts in different energy regions through deep learning by relying on training samples and can implicitly estimate the background curve based on the peaks of ¹³⁷Cs and ⁵⁹Fe, as well as other spectrum characteristics, resulting in the best estimation.

Validation results

Figure 7 shows the continuum estimated by each method for the singlet of ⁵⁷Co (102.1-151.9 keV), multiplet of ¹³⁴Cs and ¹³⁷Cs (502.6-873.3 keV), and multiplet of ⁶⁰Co (1080.4-1430.7 keV). The theoretical continuum is not shown in Fig. 7 because it is not visible in the measured spectrum, as explained in Sect. 2.1.

Fig. 7Full-screen View

Estimation results of each method on measured spectrum. (a) Full estimation and magnified views for (b) multiplet of ¹³⁴Cs and ¹³⁷Cs, (c) singlet of ⁵⁷Co, and (d) multiplet of ⁶⁰Co

For the singlet, the performances of the proposed, step-type, and peak erosion methods are similar to those reported in Sect. 3.2.2. However, the fitting method performs significantly worse (Fig. 7c), possibly due to incorrect fitting, leading to a 55.1% underestimation of the activity of ⁵⁷Co. Moreover, the AcE of the proposed method (Table 4) is slightly higher than its MAcE (Table 3) because of the additional error induced by the difference between the simulated detection efficiency and the true value, as in the step-type method.

The results for the two multiplets agree with those reported in Sect. 3.2.3. The overestimation of the step-type method near the minima between the two overlapping peaks (approximately 725 keV and 1250 keV) in Fig. 7b and d and the negative estimation of the fitting method above 1350 keV in Fig. 7d are also observed.