I. INTRODUCTION

Gamma energy spectrum analysis is a technique of using gamma spectrometer to measure γ-rays of radionuclides and determine directly categories and contents of particular radionuclides by processing the measured γ-ray spectra. Generally, the data processing includes smoothing, background deduction, peak searching, calculation of peak areas and energy calibration. Certainly, peak position and peak area are the key to qualitative analysis and quantitative analysis respectively. Traditional peak searching methods [1] includes IF function method, Gaussian product function method, derivative method, covariance method, and symmetrical zero area transformation method. Although they can be of high precision for strong radioactive γ-ray spectra and singlet, these peak searching methods are rough for weak radioactive γ-ray spectra and doublets. Especially for weak radioactive γ-ray spectra, IF function method, Gaussian product function method and covariance are not desirable [2], while derivative and symmetric zero area methods have limited effects.

Most of the traditional peak area methods, such as total peak area (TPA) method, Covell peak area method, Wasson peak area method, Sterlinski peak area method, Wasson-Sterlinski peak area method, Quittner peak area method, W-S-Q peak area method, and Q-S peak area method, request determination of the peak position and border, and can be only applied for singlet only [3].

The curve fitting method performs better [4, 5]. Its basic principle is as follows. A mathematical model, or a function, is established to describe the shape of total energy peak and baseline, then parameters of the function are determined from the experimental data, and finally the peak area is obtained by integrating the function. In practice, however, due to complexity of detectors and detection environment, measured γ-ray spectra are often too complicated to determine the function. In other words, application of the curve fitting method of peak area is extremely circumscribed.

In 1997, de-convolution [6] was proposed to process γ-ray spectra. For ill-posed problems, conventional de-convolution methods, such as maximum entropy method [7], is sensitive to the error of input data, namely, a small error can cause great oscillation. Fortunately, Tikhonov et al. studied ill-posed problems in 1977 on a strict mathematic basis by introducing the regularization theory and methods [8]. The idea is to transform ill-conditioned problems into well-posed problems, so as to ensure the solution acceptable and stable enough. Regularization falls into three categories generally: the least squares, smoothing methods, and iterative. The least-squares method is vulnerable to premature and the solutions have large error [9]. The degree of smoothing is not easy to control and smoothing itself has error [10]. In the iterative approach, one attempts to generate successive approximations which converge to an appropriate solution. Perhaps the most significant feature of iterative method, for the unfolding application, lies in the ability to incorporate readily most of the major physical implication of the description. Also, the necessity for an explicit calculation of the inverse matrix is avoided [11]. This is important, as the most response matrices are ill-conditioned, and the error in the elements of the inverse matrix can be prohibitive.

There are three commonly used regularization (iteration) de-convolution methods: Gold [12], Richardson-Lucy [13], and the maximum posteriori algorithms [14]. The direct demodulation method (DDM) [15, 16] was proposed by Li Tibei and Wu Mei. Its basic principle is to use some known physical conditions to control iterative process for solving modulation equation. The Richardson-Lucy iteration converges to the maximum likelihood solution. Compared with conventional de-convolution methods, DDM can be rarely restricted by shape of γ-ray spectra, radioactive intensity, and error of input, by taking full advantage of known information and utilizing nonlinear physical constraint conditions. DDM has succeeded in high-energy astronomy for reconstructing astrophysical images of low signal-to-noise ratio (SNR), low statistics and low-resolution [17-19]. But there are rarely reports about DDM for γ-ray spectra. In this paper, we propose to utilize DDM to reconstruct γ-ray spectra. Monte Carlo simulation is performed to calibrate response functions and establish response matrix. The purpose of this study is to improve the energy resolution of gamma spectrometry, and accuracy of γ-ray spectra analysis.

II. THE DIRECT DEMODULATION METHOD

The input and output of gamma spectrometer can be regarded as a complete system. The incident radiation can be imaged as a sum of δ-functions. So at the system output the spectrum represents a linear combination of δ-functions with various amplitudes located at various channels function, which are blurred by the system response functions. The purpose of DDM was to eliminate the influence of response function to the utmost. Ideally, a complete γ-ray spectra consisting of δ-functions can be obtained.

Let us suppose that the intensity distribution function of source is f(x). From a signal processing point of view, measurement process of spectrum can be regarded as the modulation process of spectrum of radioactive source. The measured γ-ray y(j) is the modulated result. The modulation equation can be expressed as

where h (j, x) represents the system response function.

The matrix equation can be expressed as

where H is system response matrix. For γ-ray spectra, H is the lower triangular Toeplitz matrix:

Equation (2) means that spectrum of radioactive source F is modulated by gamma spectrometer, and the measured spectrum Y is obtained. The reconstruction process is an inverse process: demodulate the spectrum Y to reconstruct spectrum of radioactive source F.

For gamma spectrometry, y(j) and f(x) are both nonnegative discrete points. It should be emphasized that this is vitally important for DDM. Besides, h(j, x) is discretized in the numerical calculation. Hence the discrete form of Eq. (1) can be expressed as

where h(k, j) is the element at the k^th row and j^th column of response matrix H, f(j) is the j^th element of F and y(k) is the k^th element of Y. However, Eq. (4) is an ill-posed problem. So DDM is applied for Eq. (4).

Richardson-Lucy iteration algorithm (R-L) is based on Bayesian statistical theory. It can be expressed as

where d⁽ⁿ⁾(k) = h(k,j)f⁽ⁿ⁾(j). However, Eq. (5) is apt to "premature". So an enhanced Richardson-Lucy iterative algorithm [8] is applied in this paper.

The nonlinear physical constraints, including upper limit and lower limit, can be expressed as [20]

Obviously, each element of F is less than the upper limit. So the lower limit only needs to be considered. For γ-ray spectra, background is regarded as the lower limit b(j).

Equations (5) and (6) consist a complete representation of DDM. Referring to Ref. [8], exponential factor p was introduced to DDM. The principles of enhanced DDM can be expressed as follows:

Step 1: for n=0, set the iterative initial value f (0) = [1, 1,, 1]^T;

Step 2: setting the number of iterations L and cycle number R;

Step 3: initializing the cycle number r=1;

Step 4: according to Eq. (4) to calculate and seek the solution f^(L);

Step 5: introducing parameter p, and f⁽⁰⁾(i) = [f^(L)(i)]^p, p (1, 2), i=0, 1,, N-1;

Step 6: introducing physical constraints;

Step 7: if r=R, then stop calculating. Otherwise, r=r+1, continue to do Step 4.

III. RESULTS AND EVALUATION

The difficulty evaluation in separation of doublets of IEC standard [21] is given in Table 1. Referring to it, the Gaussian function was used to simulate nine full-energy peaks. Their positions, heights and areas are listed in Table 2. The synthetic spectrum is shown in Fig. 1(a). From Table 2, the biggest ratio in height is 10:1, and the minimum interval of peak positions is less than 1/3 FWHM. All cases except for area ratio 1:100 in Table 1 are considered. The spectral line was reconstructed by applying DDM with 1000 (L=50, R=20, p=1.8), 5000 (L=100, R=50, p=1.8), 10000 (L=200, R=50, p=1.8), and 50000 (L=500, R=100, p=1.8) iterations, respectively. Then, peak areas were calculated by accumulating counts for "isolated" singlet of reconstructed spectrum. The result of 50000 iterations is shown in Fig. 1(b).

Fig. 1.Full-screen View

(Color online) synthetic spectrum (a) and reconstructed spectrum using DDM algorithm after 50000 iterations (b).

Also, width of the reconstructed peak decreased with increasing number of iterations. Theoretically, a complete spectrum consisting of δ-function (single pulse) can be obtained after a large number of iterations. And corresponding counts of pulse are peak areas. But a bigger iteration means more time and memory consumed. Thus, from the perspective of efficiency, iteration is stopped when doublets are completely separated. The peak areas are not the corresponding counts of pulse, but the sum of corresponding counts of a few "isolated" points. The results show that DDM can effectively decompose the peaks overlapped by just 1/3 FWHM. It should be emphasized that this is an extremely difficult task.

Table 3 lists the peak positions and areas of the DDM-reconstructed spectrum. For area ratio of 1:1 and peak position interval of about 1/3 FWHM, the error of peak position is 2 channels; while it is 1 channel for area ratio of 1:10 and peak position interval of about FWHM. The peak area error is less than 4% for all cases. This indicates that positions and areas can be precisely obtained with DDM.

IV. DISCUSSION

Actually, a measured γ-ray spectrum usually contains noise and background. Therefore, it is necessary to study the influence of noise and background on DDM.

A. Noise

A Gaussian peak is added with random (Fig. 2(a)) and gauss (Fig. 2(c)) noises, in amplitudes of 3% and 1% of the Gaussian peak height, respectively. The results are shown in Figs. 2(b) and 2(d), respectively. One sees that DDM has a strong inhibitory effect on noise. Because the least-squares, smoothing methods, and iterative can solve ill-conditioned problem, and the Richardson-Lucy iteration algorithm is adopted, noises can be intangibly inhibited in the process of iteration. If the noise is not serious, a spectrum can be reconstructed without de-noising.

Fig. 2.Full-screen View

(Color online) A Gaussian peak and its DDM reconstructed spectrum. (a) and (b) peak with random noise and the DDM result; (c) and (d), peak with gauss noise and the DDM result.

B. Background

A spectrum is added with straight-line, oblique-line, and step backgrounds, respectively. The DDM-reconstructed results are shown in Fig. 3. The spectra were not reconstructed correctly. This indicates that background affected the accuracy of reconstructed spectra, but researches show that if the response matrix can be obtained under the same condition of background, or the background can be removed from the measured spectrum and response matrix by applying the same way, the lines can be reconstructed accurately.

Fig. 3.Full-screen View

(Color online) Spectra with backgrounds in (a) straight line, (b) oblique-line and (c) steps, and the reconstructed spectra.

V. CALIBRATION OF NAI GAMMA SPECTROMETER BASED ON DDM

DDM was used to calibrate the field NaI gamma spectrometer. Generally, spectrum stripping was used. The principles of calibration can be seen in Ref. [25]. IED-3000A NaI gamma spectrometer (104#) was used. The dimensions of NaI crystal were 75 mm×75 mm, with an energy resolution of 10.3% (661 keV). The measuring time was 3600 s. Substance contents of quasi-saturation model were listed in Table 4, which were measured by Analysis and Test Center of China Geological Survey.

The measured spectrum of hybrid model was shown in Fig. 4(a). SNIP [26] was applied to deduct background. DDM was used to reconstruct the spectrum. The result after 50000 iterations is shown in Fig. 4(b). Calibration coefficients and net peak areas of 1.46 MeV (⁴⁰K), 1.76 MeV (²¹⁴Bi of U series), and 2.62 MeV (²⁰⁸Tl of Th series) are listed in Table 4.

Fig. 4.Full-screen View

(Color online) The measured spectrum of the U model (a) and the reconstructed spectrum (b).

Table 5 shows the elemental analysis results of a γ-ray analysis standard by applying DDM and spectrum stripping method (SSM), It can be seen that the DDM results are of smaller deviations from the given value of U, Th and K, being respectively -1.02%, -4.19% and -2.41%, while the SSM results deviated by -8.06%, 10.44% and 8.84%, respectively.

VI. CONCLUSION

Due to response function and iterative algorithm, DDM can be rarely restricted by shape of γ-ray spectra, radioactive intensity. Besides, DDM can significantly pseudo-improve energy resolution, precisely search peaks, calculate areas, and improve the performance of gamma spectrometer without hardware cost. It is worth mentioning that response matrix is the key of DDM. However, in practical, it is tedious and difficult to use the standard radioactive sources to establish response matrix. Hence GEANT 4 was proposed to calibrate response function. It is convenient, safe, and effective to employ Monte Carlo method to establish response matrix.

On the other hand, in theory, DDM can pseudo-improve energy resolution to infinity after a large number of iterations. But the bigger iterations, the more time and memory will be consumed. Thus, from the perspective of efficiency, we should stop iteration when doublets were completely separated.

It should be emphasized that, theoretically, DDM is not only appropriate for γ-ray spectra but also for any other spectrum as long as response matrix can be established correctly. In fact, DDM is not only just applied to NaI (Tl) detector but also to any other detector. It is time-consuming to establish response matrix and run DDM, especially for the complex spectrum. Thus, DDM is usually regarded as a means for offline analysis. In this respect, how to improve the running speed will be focused during the next study.