2.1 Gaussian modeling

In order to fit the peak shape, many curves (e.g., Gaussian, Lorentzian, mixed Gaussian-Lorentzian, and modified Gaussian) have been proposed. Fits using curves of Gaussian, Lorentzian, and mixed Gaussian-Lorentzian (Voight) type were experimented with for comparison to a real chlorite, short-wavelength infrared (SWIR) spectrum obtained by Brown [10]. Miyamoto et al. found that a Gaussian function shows a better fit than a Lorentzian function for the decomposition of ordinary chondrite and eucrite spectra [11]. Furthermore, a modified Gaussian successfully modeled the electronic transition absorption in pyroxene, pyroxene mixtures, and olivine spectra [12].

Here, we use Gaussian fitting function to decompose gamma-ray spectra into individual full energy peaks, as described by Eq. (1):

where ${{\displaystyle A}}_i$, ${{\displaystyle B}}_i$, and ${{\displaystyle p}}_i$ are the peak height, the full width at half maximum (FWHM), and the center of the ith peak, respectively, and ${{\displaystyle x}}_j$ is the channel number [13,14].

2.2 Levenberg-Marquardt algorithm

The Levenberg-Marquardt algorithm is a method of finding the minimum value of a function even if the initial value is far from the final value [15]. It is a combination of a gradient descent algorithm and Gauss-Newton algorithm. It is more accurate than gradient descent algorithms and more robust than the Gauss-Newton algorithm and can effectively improve convergence performance. Therefore, the main advantage of the LM algorithm is the ability to make the values of model parameters range under strong control [8].

The LM Gaussian fitting algorithm is described as follows:

Where $\mu >0$, $I$ is a unit vector, $d$ is the mobile vector, and $a={{\displaystyle (A,B,p)}}^T$. We introduce a function ${{\displaystyle r}}_i(a)={{\displaystyle y}}_i-y({{\displaystyle x}}_i,a),{{\displaystyle }}^{}i=1,2,\cdots ,n$, and let $r(a)=({{\displaystyle r}}_1(a);{{\displaystyle r}}_2(a);\cdots ;{{\displaystyle r}}_n(a))$; $J$ is the Jacobi matrix of $r(a)$. ${{\displaystyle y}}_i$ and $y({{\displaystyle x}}_i,a)$ are, respectively, the original and fitted value, ${{\displaystyle \sigma }}_i$ is the standard deviation, and $n$ is the total number of channels. The purpose of the algorithm is to make ${{\displaystyle \chi }}^2(a)$ a minimum by searching for the optimal solution of $a$ with the given initial value of $a$, replacing $a$ into the Gaussian function after obtaining the optimal value, and then decomposing the overlapping peaks [16-21]. A flow chart of the LM algorithm is presented in Fig. 1.

Fig. 1.Full-screen View

Flow chart of Levenberg-Marquardt algorithm.

In this paper, a comparison with a genetic algorithm found in Matlab software was also made.

2.3 Data preprocessing

In gamma-ray spectrum analysis, a very large error may be introduced when the method applied to deduct the background is not accurate. Among background-subtraction methods, the statistics-sensitive, nonlinear iterative peak-clipping (SNIP) algorithm has been broadly applied in nuclear data calculations. Its main advantage is that it can compress the dynamic range of the channel counts and enhance small peaks. It is both swift and accurate. Therefore, in this paper we used the SNIP method to deduct the background.

We first transformed the original data $y(i)$, letting $Z(i)=\ln [1+\ln (1+\sqrt{y(i)+1})]$, where $i$ is the channel number. A natural logarithmic operation has advantages in processing high-count spectral peaks and a square-root operation has comparative superiority in strengthening the weak peaks. Second, in the wth iteration, we calculated the new value in channel $i$ as ${{\displaystyle Z}}_w(i)=\min \{{{\displaystyle Z}}_{w-1}(i),[{{\displaystyle Z}}_{w-1}(i-w)+{{\displaystyle Z}}_{w-1}(i+w)]/2\},{{\displaystyle }}^{}w=1,\cdots ,N$, where $N$ is the number of iterations. Finally, we obtained the background channel counts, $B(i)={{\displaystyle Z}}^{-1}(i)$[22]. Then, the net peak area can be easily determined.

3.1 In situ gamma-ray spectrometry

A NaI (Tl) scintillation spectrometer and a ¹³⁷Cs point source (3.7×10⁴ Bq) were used in our simulation experiment. The scintillation spectrometer is portable and of high precision for taking measurements in the field. The diameter of the NaI (Tl) detector is 3×3 in., and its energy resolution is 7.5% at 662 keV [23].

A schematic of the simulation experiment is shown in Fig. 2. The NaI (Tl) spectrometer was placed on the ground, and a ¹³⁷Cs point source was buried in the soil directly under the detector. The source-to-detector distance changed from 18 to 38 cm when the source was put at a designated depth and the hole then filled with soil. Spectra were recorded at every depth, with an acquisition time of 600 s. The simulation experiment at every depth was performed twice.

Fig. 2.Full-screen View

Schematic of ¹³⁷Cs point source simulation experiment.

3.2 Lack-of-fit test

A lack-of-fit test is a useful method in statistics, often used to determine whether the proposed model fits well. We partition the error into two components:

where $n$ is the total channel number, $m$ is the number of times the experiment is repeated, ${{\displaystyle y}}_{ij}$ and ${{\displaystyle \widehat{y}}}_i$ are the observed value and fitted value, respectively, and ${{\displaystyle \overline{y}}}_{i•}$ is the average of observed y values in channel $i$. The first component is called pure error. The second component is called lack of fit [24].

It then follows that, statistically,

where $N$ is the total number of observations and $p$ is the number of parameters in the model. The degrees of freedom of pure error and lack of fit are $N-n$ and $n-p$, respectively. We reject the model if the value of the F statistic is larger than the critical F value. Otherwise, we accept the model. We find the critical value in the critical value table of F distribution with the desired confidence level, and the degrees of freedom ${{\displaystyle d}}_1=n-p$ and ${{\displaystyle d}}_2=N-n$.

3.3 Signal-to-noise ratio

The signal to noise ratio (SNR) is the ratio of the power of signal and power of noise, so it is the square of the amplitude of signal and noise. Its units generally are dB, and it is expressed in this paper as

where ${{\displaystyle A}}_{\text{peak}}$ and ${{\displaystyle A}}_{\text{background}}$ are the amplitude of the ¹³⁷Cs peak and background, respectively. The smaller the SNR, the lower the amplitude of the ¹³⁷Cs peak. When the value of SNR is too small, smaller than even −30 dB, the amplitude of the ¹³⁷Cs peak is smaller than 10 counts. It is difficult to say that the peak of ¹³⁷Cs was the source for this; assuming so may lead to incorrect estimation.

3.4 Evaluation methods

The relative fitting error of the total net area $\delta$, residual standard deviation $S$, and fitting goodness ${{\displaystyle R}}^2$ were also used to evaluate the curve-fitting effect. Their mathematical expressions are

where ${{\displaystyle A1}}_{Th+Ra+Cs}$ and ${{\displaystyle A2}}_{Th+Ra+Cs}$ are the total net area of the three full energy peaks at 583, 609, and 662 keV before fitting and after fitting, respectively, $y$ is the original value, $\widehat{y}$ is the fitting value, and $\overline{y}$ is the average of the original value. The closer $\delta$ and $S$ are to 0 and the closer ${{\displaystyle R}}^2$ is to 1, the better the fitting curve.

4.1 ¹³⁷Cs source experiment

Based on the data from the ¹³⁷Cs point source simulation experiment (Fig. 2), the decomposed and fitted curves, and the results for three overlapping peaks as determined by the LM algorithm, are shown in Fig. 3 and Table 1. Figure 3 and Table 1 show that the net area of the ¹³⁷Cs peak (662 keV) became smaller as the source-to-detector distance increased, while the full energy peak areas of ²⁰⁸Tl (583 keV) and ²¹⁴Bi (609 keV) emerged. The decomposed net full energy peak areas of ²⁰⁸Tl and ²¹⁴Bi were larger, fluctuating within statistics, when the source-to-detector distance was larger than 24 cm. The average value of the fitting net peak area of 583 and 609 keV is 1.028 counts per second (cps) and 0.798 cps, respectively, but the net area of 662 keV becomes smaller and smaller.

Fig. 3.Full-screen View

Decomposed and fitted curves for the overlapping peaks at different source-to-detector distances (a, 18 cm; b, 24 cm; c, 30 cm; d, 34 cm; e, 38 cm; f, no source).

4.2 Discussions

Table 2 presents the results of evaluating the data in Table 1 for the decomposition and fitting of ¹³⁷Cs source simulation spectrum data. In a lack-of-fit test, the instrumental noise levels were estimated from replicated analyses. The degrees of freedom d₁=62 and d₂=71, F_0.05(d₁,d₂)=1.498, and all of the F values obtained and presented in Table 2 were smaller than 1.498. Therefore, the two models proposed in the paper are feasible. The SNR values were even negative and showed a decreasing trend when the source-to-detector distance increased. The SNR is smaller than −10dB when the source-to-detector distance is larger than 30 cm. This also indicates that the amplitude of the ¹³⁷Cs peak is becoming smaller and smaller. The relative error of total net area $\delta$ and residual standard deviation $S$ were no more than 5% and 0.04, respectively, and the fitting goodness ${{\displaystyle R}}^2$ is larger than 0.98. These results show that the methods used in the paper give high performance. However, it is difficult to say that the peak at 662 keV was a contribution of ¹³⁷Cs when the SNR value of in Table 1 is smaller than −30 dB. It may stem from the statistics fluctuation from the background for an amplitude of only approximately 6 counts. This situation is suited for a source-to-detector distance of 38 cm and no source. The determination of the ¹³⁷Cs peak of the two distances may lead to incorrect estimation in Fig. 2 and Table 1.

The comparison between the LM algorithm and genetic algorithm is presented in Table 1 and Table 2. The results of the application of our approach based on the LM algorithm yielded similar results compared to other studies based on the genetic algorithm. The results were in good agreement, but the genetic algorithm is costlier in terms of computer time. Thus, the LM algorithm gives better performance than the genetic algorithm based on these results.