2.1 Geometry of NaI(Tl) Scintillation Detector and Building of Theoretical Model

The structural model of the NaI(Tl) scintillation detector is shown in Fig.1, which is composed of a cylinder with 75 mm in diameter on the substrate surface and 75 mm in height, and a single layered aluminum shell within the cylinder that is 3 mm thick. The back-end of the cylinder has a photoconductive layer, which connects the photomultiplier tube with 75 mm in diameter and 20 mm in height, and other components. The detector is placed within a lead shield layer with an inner radius of 90 mm, an outer radius of 140 mm, and a height of 360 mm, respectively. A simulation model contains the NaI(Tl) crystal, a reflector, a parcel layer, an optical (SiO₂) and lead shield layer so as to close to the experimental conditions as much as possible. The crystal size is Φ75×75 mm, the aluminum shell thickness is 3 mm at the front facet and 2 mm on the lateral surface, and the MgO₂ reflective thickness is 0.5 mm, and its SiO₂ optical glass thickness is 2 mm on the back surface. In order to decompose the spectrum of consistency, the used resolution of the simulated detector is set to 8% and the detector resolution of the measured spectrum is about 8%. The simulated radioactive source is a point source with a distance of 50 mm from the detector. In the simulation calculations, in order to weigh the simulation time consuming and guarantee sufficient statistical accuracy, the incident photon number is set to 10⁸.

Fig.1Full-screen View

NaI(Tl) scintillation detector geometry model

The main purpose of the Monte-Carlo simulation is to produce a corresponding geometric model that is as close as possible to the response function of the actual detector. According to detector’s geometric model, after simulating the 661.6 keV of ¹³⁷Cs, the difference between the simulated spectrum and the measured ¹³⁷Cs radioactive source spectrum (the measured spectrum) is shown in Fig.2. The characteristic parameters of the full-energy peak, the Compton plateau and edge are all consistent. The full-energy peak (peak 1) is somewhat consistent but there are some gaps at the low-energy end. The Compton edge of the simulated spectrum (peak 2) is higher than the measured spectrum, and the backscattering peak (peak 3) is lower than the measured spectrum. On the one hand, the result is closer to an ideal state because the simulated photons in the transport processes are closer to practical effect of γ photons in the detector, and because they are similar to the γ-rays absorbed by the scintillator. On the other hand, the actual measurement is performed in the lead shield chamber. The measured spectrum also contains the X-ray peak (peak 4) released from the decay daughter ¹³⁷Ba of ¹³⁷Cs, which is unique for ¹³⁷Cs. The differences between the theoretical simulation and actual measurement environment and conditions are responsible for this. However, the theoretical simulation reflects the response characteristic of the detector, and does not affect the generation of the Monte-Carlo response-matrix. Therefore, the deviation is acceptable.

Fig.2Full-screen View

A simulated spectrum with the ¹³⁷Cs Source is contrasted with the measured spectrum

2.2 Energy Deposition and Response Function

When using the Monte-Carlo simulation and calculation, it is necessary to consider the energy deposition of each photon in the NaI(Tl) scintillation detector. The photon-electronic transport processes of energy deposition corresponding to γ-ray spectra track a large number of γ photons that are either absorbed by the detector or escape from the detector. If the MCNP software is used for the simulations, energy deposition calculations can be obtained after accumulating the energy deposition of each source’s particle event. Considering the resolution of the γ-ray spectrometer and Gaussian broadening of energy deposition, the relationship between Full Width at Half Maximum (FWHM) and γ-ray energy E_d can be written as the following formula:

a, b, and c are the calibration coefficients of resolution, a=0. 001, b=0.05086, c=0.030486. The energy is:

In Eq.(2), $\sigma =FWHM/2\sqrt{2\ln 2}=FWHM/2.355$, and x is obtained using the standard normal distribution sampling.

When calculating the detection efficiency of the γ–ray spectrometer, it is necessary to record the total number n of the full-spectra during the simulation process: including the number N of incident γ photons in the detector and the total number of photons to be recorded by the detector, and a ratio of full-energy peak counting n_p to the total number n of full-spectra is used to calculate the peak-to-total ratio (R[E]). The energy deposition spectrum approximates to the pulse response function of the γ-ray spectrometer. Simulating the response function using the Monte-Carlo method is actually the process of tracking particle trajectories and their interactions in the NaI(Tl) scintillation detector. With consideration to radioactive statistical fluctuations, when a photon’s tracking processes end, according to the variance that responds to deposition energy, a Gaussian distribution should be used to sample the deposited energy to obtain the pulse counting generated by photons interacting with each other in the NaI(Tl) crystal. Therefore, the total deposition energy of each incident photon in the NaI(Tl) crystal is equal to the sum of the deposited energy after the photons have interacted with the dielectric and deposited energy of each second photon [6], which is the photon’s response in the NaI(Tl) scintillation detector.

3.1 Relationship between Response Function and Response-Matrix

In the theory of the continuous linear system [7], the principle of convolution methods is that the signal is decomposed into the sum of the impulse signals with the system response function h(t), which is used to solve the zero state response of the system to an arbitrary excitation signal.

The convolution formula is shown as follows:

where y(t) is the output function, h(t) is the response function, and x(t) is the input function. If any two of the three functions are known, a third one will be similarly sought; the formation of the γ-ray instrument response spectrum also can be seen as the form of convolution between the input function and the response function, which leads to the output function. The c energy spectrum form diagram is shown in Fig. 3.:

Fig.3.Full-screen View

The γ energy spectrum form diagram

The formation of instrument response spectrum can be represented as follows:

In Eq.(4), y(n) is the measured response spectrum, h(n) is the impulse response function, and x(n) is the incident spectrum. The decomposing process has turned into a process of solving the input x(n) through the measured instrument response spectrum y(n) and the predictable response function h(n) of the detector. Eq.(4) is expressed as a matrix or matrix equation as the following:

or

In Eq.(5), the vector y is the instrument response spectrum in channel n, R is a response matrix, and vector x is the original incident γ-ray spectrum in channel m, which is different from the instrument response spectrum. Each channel corresponds to one characteristic energy point of the incident γ-ray spectrum: the energy of incident γ-ray. Column vector $R_{\mathrm{...},i}$ is a response to the energy points in channel i, and is also a spectra line in channel n, $i=1\cdots m$.

Therefore, from a mathematical perspective, the spectra x of the incident γ-ray can be generated after decomposing the instrument response spectrum through an inversion decomposition on linear Eq.(5) or (6). The response matrix is composed of the response function, which includes various spectral characteristics, meaning that response function cannot be obtained through the transformation of a simple photoelectric peak function, and each γ-ray spectrum needs to generate the corresponding response function. Because the response functions may be obtained by several γ-rays with single energy, several single γ-rays with a single energy spectrum can be implemented by means of the Monte-Carlo simulation to obtain more energy-segment spectral lines [8].

3.2 Construction of Simulating Full Spectrum Response-Matrix

The measured spectrum was obtained in this experiment using a NaI(Tl) scintillation detector with a resolution of 8% and Φ75×75 mm in size; the distance from the front facet of the detector to the samples and to the source are both 50 mm. Considering that environmental impacts and statistical fluctuations occur in the actual measurements, and some deviations exist between the response matrix and the actual response, the system is not absolutely linear, so the response matrix must be constructed under the same environmental conditions (such as the size of the detector and the distance from the detector to source). Since the simulation on each response function is time consuming within the full spectrum (0~3 MeV), to save computation time, according to the energy interval of 100 keV, each corresponding response function is calculated one by one, and another response function is calculated using the interpolation algorithm between the two simulated response functions at an interval of 100 keV. For example, assuming that two response functions present [$R_1(E_1,e)$, $R_2(E_2,e)$], respectively, the corresponding energy values $e_i$ can be calculated using the interpolation algorithm in the response matrix $R_i$, as follows:

In general, it is assumed that the interpolation segment is divided into n parts. Each point of the response functions $R_i(E_i,e)$ in the interpolation can be classified and then calculated. All the response functions in a given area can be generated using the above algorithms. The whole response matrix R can be obtained by filling all the interval areas. The constructed part response matrix with the interpolation algorithm is as shown in Fig.4. Its interpolation diagram is as shown in Fig.4a, the thick and thin lines representing the simulated values and response functions of the interpolation, respectively, and a part of response matrix was obtained using the above algorithms, as shown in Fig.4b.

Fig. 4Full-screen View

(Color online). The constructed part response matrix with the interpolation algorithm

As shown in the above figures, the current response functions are divided into four aspects: $(0,E_{\text{b}})$, $(E_{\text{b}},E_{\text{c}})$, $(E_{\text{c}},E_{\text{p}}-10\text{keV})$ and $(E_{\text{p}}-10\text{keV},E_{\text{p}}+10\text{keV})$, among them, $E_{\text{b}}=E_{\text{p}}/(1+2E_{\text{p}}/m_{\text{e}}{c}^2)$, $E_{\text{c}}=E_{\text{p}}-E_{\text{b}}$ and $E_{\text{p}}$ represent backscattering peak energy, the Compton edge and the photopeak, respectively. When energy is lower than $E=m_{\text{e}}{c}^2/2$, the Compton edge overlap with backscattering peak and the regions will turn into $(0,E_{\text{c}})$, $(E_{\text{c}},E_{\text{b}})$, $(E_{\text{b}},E_{\text{p}}-10\text{keV})$ and $(E_{\text{p}}-10\text{keV},E_{\text{p}}+10\text{keV})$, here $E_{\text{c}}<E_{\text{b}}$. Within the range of 200 keV to 300 keV, the Compton edge and the backscattering peak are aliasing. In order to achieve a better resolution and accuracy with the interpolation method, the response function interval obtained through the simulation should be appropriately reduced to 50 keV. When the energy is greater than 1.02 MeV, it is necessary to consider the electron pair effect, which results in the appearance of escape peaks on the γ energy spectrum. Therefore, when energy is greater than 1.02 MeV, it is necessary to interpolate the response functions separately and add escape peaks into accurate locations (E-511 keV and E-1022 keV).

Each spectral characteristic parameter of the γ-ray spectrum in the simulated response matrix is shown in Fig.5. The energy is between 100 keV~3000 keV, and it is clearly seen that γ photons with different energies produce the instrument response spectrum results in the NaI(Tl) scintillation detector. Abscissa corresponds to the channel of the instrument response spectrum (Energy), and the ordinate represents the energy of γ photons that are incident to the NaI(Tl) detector (Count). The full energy peak located near the diagonal line and its trend changes also can be clearly seen. The Compton platform is located at the bottom left triangle; single and double escape peaks are generated by γ photons when the energy is greater than 1.02 MeV, the effect of which is more obvious with the increase of γ photon energy. For example, there are backscattering peaks around the E-250 keV, and single and double escape peaks are located at E-511 keV and E-1022 keV, respectively, and the Compton edge and the full energy peak also can be clearly observed.

Fig.5Full-screen View

(Color online). The full-spectra response matrix simulated by the MCNP

4.1 Gold Decomposition Method

The analytical algorithm of the Gold analysis is described in detail in Ref. [9-14]. It is based on Van Cittert analytical algorithm [15-18]. By multiplying $R^T$ on both sides of the Eq. (5) one can obtain the Toeplitz linear equations:

in Eq. (8), $A=R^{T}R$, $z=R^{T}y$.

The results of the k+1 step iteration can be expressed as the following:

where $\mu$ is the relaxation factor, and the local variable relaxation factor is introduced:

Substituting the local variable relaxation factor ${\mu }_i$ into formula (9), one obtains:

namely:

L is the number of iterations. To further define the vector:

one obtains,

Iteration algorithms start with the initial solutions,

The algorithms are positively definite. If all quantities are positive values in Eq.(12), then all of the solutions are also positive values. As a result, the algorithms are suitable for all positive data, and the estimation of x vector is constrained by the sub space of positive solutions.

4.2 The Improved Boosted-Gold Decomposition Method

The result from the above proof shows that when the Gold analysis has converged to a stable value, no matter how the number of iterations is increased, the decomposed results will no longer change. If the width of the peaks continues to shrink, it should stop the iterations when the solutions reach a steady state, then, $x^{(L)}$ vector will change in some way, which is then taken as a new initial value; at that point, the Gold iterative solution of formula repeats (12). The nonlinear function enhancement (acceleration) should be used to change the particular solutions [16-19]. The power function has been shown to give better results; its enhanced algorithms of the improved iterative Gold are described as follows:

(1) According to equation (15), the initial value is set to $x^{(0)}$; (2) times R and iteration times L should be set repeatedly according to the needs; (3) repeat times r=1; (4) according to $k=0,1,\mathrm{...},L-1$, seek out $x^{(L)}$; (5) if r=R, stop calculating; otherwise, continue to the following operations: (i) accelerate enhanced operation, namely $x_i^{(0)}={\left|x_i^{(L)}\right|}^p$, $i=1,2,\mathrm{...},N$, p is accelerative index; (ii) r=r+1; (iii) repeat step (5).

5.1 Test of Decomposing the Simulated Multi-characteristic Energy γ-Ray Spectrum

According to Monte-Carlo numerical modeling and the detector’s geometric model, the mixed sources ¹³⁷Cs and ⁶⁰Co of the incident γ-ray are simulated and the activity ratio is set to 0.35:0.65. The number of emission photons is 10⁸, so the simulated mixed source spectrum is as shown in Fig.6, and the response matrix generated by the Monte-Carlo simulation is as shown in Fig.5. The decomposed results are shown in Fig.7.

Fig. 6Full-screen View

(Color online). The simulated spectrum of ¹³⁷Cs/⁶⁰Co mixed sources

Fig. 7Full-screen View

(Color online). The decomposed results of the ¹³⁷Cs/⁶⁰Co mixed spectrum with Gold and Boosted-Gold method

Table 1 lists the inversion-decomposing results of the mixed-source response spectrum for the radioactive ¹³⁷Cs and ⁶⁰Co with different activity ratios, performed through the Monte-Carlo simulation. The deviations of the decomposition results with this method are small, so the deviations are smaller, allowing the results to be more stable. This method is also used to accurately obtain the characteristic energy of nuclides. According to the information in Table 1, both of the activity ratios in the Monte-Carlo simulation and the decomposing mixed sources for ¹³⁷Cs and ⁶⁰Co show that one of positive features: ‘Features of maintaining content equal’ (the total count of decomposed spectra is roughly equal to the number of photons that come from un-decomposed photon sources that are incident to the crystal detector). This is because the whole process, from the interaction of γ photons simulated by the response matrix in the detector to the formation of electrical signals and radionuclide activity, can be directly solved out in theory. Therefore, under ideal conditions, the inversion-decomposing results that come from using the Gold algorithms are an activity of the nuclide.

5.2 Test of Decomposition to Synthetic Overlapping Peak Spectrum

The purpose of testing ¹³⁷Cs synthetic spectrum is to validate the decomposition capability of analytical methods for overlapped peaks. The synthetic spectrum is obtained using different amplitude attenuation and displacement with ¹³⁷Cs spectra lines, and then adding Gaussian noise equivalent to one percent of a maximum, shown in Fig.8. In the figure, spectral lines of channel address have been extracted at the 429 ch, 450 ch and 600 ch in the response matrix, the amplitudes of which represent 1, 0.5 and 0.5 respectively. The response matrix made by the Monte-Carlo simulation and interpolation is shown in Fig.5, and decomposed results are shown in Fig.9.

Fig.8Full-screen View

(Color online). The synthetic overlapping peak spectrum

Fig. 9Full-screen View

(Color online). The decomposed results of the synthetic spectrum with Gold and Boosted-Gold method

The resolution of the simulated detector is 8%, which corresponds approximately to the 50 channels of the FWHM; Fig.9 shows that the Gold and the Boosted-Gold methods for decomposing the spectrum have an obvious effect, and using the improved Boosted-Gold method to decompose overlapped peaks whose FWHM is less than one can obtain accurate results. Further, under ideal conditions (such as instance in which the simulated spectra and measured spectra are close), this method is used to decompose complex spectrum obtained by the lower resolution detector.

5.3 Tests of Decomposing ²³⁸U Series, ²³²Th Series and ²³⁸U-²³²Th Mixed Spectra

The purpose of testing ²³⁸U series, ²³²Th series and ²³⁸U-²³²Th mixed spectra is to validate whether the method is feasible in practical applications. The sample spectrum is measured using the NaI(Tl) detector with Φ75×75 mm in size and energy calibration is that E=2.9577 keV/ch, and the resolution of the detector is about 8%. The sample sources used in the experiment adopted the ²³⁸U-²²⁶Ra equilibrium source with a mass fraction of 1.86×10^-4 and a specific activity of 2308.26 Bq/kg, and the ²³²Th equilibrium source with the mass fraction of 1.57×10^-4 and specific activity of 638.99 Bq/kg, and the ²³⁸U -²³²Th mixed source with a mass fraction of 2.20×10^-4 and a specific activity of 1754.72 Bq/kg; additionally, their physical conditions are solid and shapes are globular. The measured time is 20000 s. For a natural γ-ray spectrum analysis, it is common to choose energy characteristic peak at 0.583, 0.609, 1.120, 1.76 and 2.62 MeV energies photopeak as counting windows, and their energy widths expand as the window of uranium and thorium, then using the inversion-decomposition method to calculate the activity or content of the energy window nuclides. The experimental result of the environmental natural background is as shown in Fig.10, and the measured time is 33890 s. The response matrix was generated by the Monte-Carlo simulation and interpolation, and was constructed in accordance with the method of the above sections 2 and 3. The decomposed results of the measured spectra for the ²³⁸U-²²⁶Ra equilibrium source and ²³⁸U -²³²Th mixed source with response matrix are shown in Fig.11 and Fig.12.

Fig.10Full-screen View

(Color online). The background spectrum(2.9577keV/ch)

Fig. 11Full-screen View

The decomposed results of the measured ²³⁸U-²²⁶Ra equilibrium source with response matrix

Fig. 12Full-screen View

(Color online). The decomposed results of the measured ²³⁸U-²³²Th mixed source with response matrix

For a characteristic peak analysis of the γ-ray spectrum, as shown in Fig.11, the decomposition results show that the characteristic peaks are located at 205 and 590 ch respectively. The linear calibration of energy is: E=2.9577 keV/ch, and the corresponding energy to be calculated represent 0.607 and 1.745 MeV respectively. According to a decay photon number and the characteristic energy of the γ-ray spectrum, the ²³⁸U series γ-ray spectra consist of the 0.609 and 1.764 MeV characteristic peaks generated by ²¹⁴Bi. The corresponding errors of characteristic energy are -0.32% and -0.85% respectively, and the decomposition of the instrument response spectrum has achieved ideal effects.

Given Fig.11, it is particularly worth mentioning that the corresponding energies of characteristic peaks in the low energy segment, located at 59, 82, 98 and 119 ch, represent, respectively, 0.179 MeV (0.1862 MeV characteristic peak generated by ²²⁶Ra), 0.242 MeV (0.242 MeV characteristic peak generated by ²¹⁴Pb), 0.292 MeV (0.295 MeV characteristic peak generated by ²¹⁴Pb), 0.351 MeV (0.352 MeV characteristic peak generated by ²¹⁴Pb). It is feasible to use the Monte-Carlo response matrix to decompose ²³⁸U-²²⁶Ra equilibrium source in a low energy segment; the effects are significant.

The expected results of the ²³²Th equilibrium source are also obtained in the same way. Its spectral decomposition results are accurate and acceptable.

The ²³⁸U-²³²Th mixed source spectrum are measured using the NaI(Tl) detector with Φ75×75 mm; energy calibration was E=2.9577 keV/ch. Since the 609 keV characteristic peaks generated by ²¹⁴Bi of the ²³⁸U series and the 583 keV characteristic peaks generated by ²⁰⁸Tl of the ²³²Th series in energy emerged with overlapping spectral peaks, the results of using the Gold and the Boosted Gold method for decomposing mixed spectrum of ²³⁸U series and ²³²Th series are shown in Fig.12.

From Fig.12, the characteristic peaks are located at 198 ch and 205 ch, respectively, and the corresponding energy respectively represents 0.581 MeV (0.583 MeV characteristic peak generated by ²⁰⁸Tl), and 0.607 MeV (0.609 MeV characteristic peak generated by ²¹⁴Bi). Therefore, ideal effects have been achieved regarding the decomposition of the overlapped peaks.

A measured spectrum of the ²³²Th equilibrium source with 40×40 cm detector is shown in Fig.13. The energy calibration is E=3.17 keV/ch, and the resolution of the detector is about 7.5%, the specific activity is 139.65 Bq/kg; their physical conditions are solid and shapes are globular. The measured time is 200 s. The response matrix was generated using the Monte-Carlo re-simulation, and it was constructed in accordance with the method shown in the above Sect. 2 and 3. The results of using the Gold and Boosted Gold methods for decomposing ²³²Th equilibrium source is shown in Fig.13. The decomposition results show that the characteristic peaks are located at 183 and 820 ch, respectively, and the corresponding energy to be calculated represents 0.581 and 2.599 MeV, respectively. According to a decay photon number and the characteristic energy of γ-ray spectrum, the ²³²Th series γ-ray spectra consist of the 0.583 and 2.62 MeV characteristic peaks generated by ²⁰⁸Tl, the corresponding errors of characteristic energy are -0.34% and -0.80% respectively.

Fig. 13Full-screen View

(Color online). The decomposed results of the measured ²³²Th equilibrium source with response matrix (40×40 cm)

From Fig.13, it is particularly worth mentioning that the corresponding energies of characteristic peaks in the whole energy segment, located at 71, 103, 183, 285, 302 and 820 ch, represent, respectively, 0.225 MeV (0.239 MeV characteristic peak generated by ²¹²Pb), 0.327 MeV (0.338 MeV characteristic peak generated by ²²⁸Ac), 0.581 MeV (0.583 MeV characteristic peak generated by ²⁰⁸Tl), 0.903 MeV (0.908 MeV characteristic peak generated by ²²⁸Ac), 0.957 MeV (0.960 MeV characteristic peak generated by ²²⁸Ac), 2.599 MeV (2.62 MeV characteristic peak generated by ²⁰⁸Tl). Similarly, the decomposition of the instrument response spectrum has achieved satisfactory effects.

In summary, in the decomposition tests, the decomposed spectrum can better distinguish characteristic peaks with similar energy, and provide a theoretical basis for studying the analytical algorithm of energy spectrum as well.