Extraction of full energy peak of ¹³⁷Cs from in situ NaI (Tl) gamma-ray spectrum

Extraction of full energy peak of 137 Cs from in situ NaI (Tl) gamma-ray spectrum 1 first-author 1 corresp Corresponding author. E-mail address: npwang@cugb.edu.cn Corresponding author. E-mail address: npwang@cugb.edu.cn npwang@cugb.edu.cn 1 School of Geophysics and Information Technology, China University of Geosciences, Beijing 100083, China. School of Geophysics and Information Technology, China University of Geosciences, Beijing 100083, China. A Levenberg-Marquardt Gaussian fitting algorithm has been used for analyzing the overlap of three peaks (the 583-keV peak of 208 Tl, the 609-keV peak of 214 Bi, and the 662-keV peak of 137 Cs) using an in situ NaI (Tl) scintillation spectrometer. The algorithm, in addition, was compared with a genetic algorithm used for multiple deconvolution. The three fitted peak areas (583, 609 and 662 keV) were calculated from the measured gamma-ray spectra obtained from a simulation experiment in which a 137 Cs source was buried at different soil depths (from 18 to 38 cm). The application of the Levenberg-Marquardt algorithm yielded similar results compared to the genetic algorithm. A lack-of-fit test showed that the fitting is good when the instrumental noise levels were estimated from replicated analyses. The relative fitting error of the total net area and the residual standard deviation were within 5% and 0.04, respectively, and the goodness of the fitting was better than 0.98. While the methods used in this paper give high performance, the results may lead to incorrect estimation when the signal-to-noise ratio is smaller than −30 dB. This study is useful for the determination of radioactive specific activity of 137 Cs by in situ spectrometry. NaI (Tl) detector Gamma-ray spectrum Levenberg-Marquardt algorithm Overlapping-peak decomposition Curve fitting 1 INTRODUCTION 1 Radioactive contamination has drawn unprecedented research attention with the continuous development of the nuclear power industry. The artificial radionuclide 137 Cs is one of the sources of contamination found in environmental materials derived from nuclear accidents, such as that which occurred in Chernobyl (Ukraine, 1986) and Fukushima (Japan, 2011), both of which caused worldwide radioactive contamination [ 1 1 , 2 2 ]. In particular, 137 Cs has a long half-life period, and thus could be retained in the environment for a long time, contaminating surface water, soils, foodstuffs, etc.. Thus, determinations regarding contamination by 137 Cs have been an important focus of research. Gamma-ray spectrometry is a highly efficient and accurate method of measuring environmental radioactivity. NaI (Tl) scintillation spectrometers play an important role in measuring 137 Cs. However, the 662-keV peak of 137 Cs normally overlaps with the 583-keV peak of 208 Tl in the 232 Th series and the 609-keV peak of 214 Bi in the 238 U series in spectra obtained with a NaI (Tl) detector, because of the lower-energy resolution of NaI (Tl) detectors and a high Compton-scattering background. In this case, more complex and accurate techniques are required to replace the simple summation method to extract the single-peak area of 137 Cs from the overlap of three peaks. So far, techniques for fitting peaks have been studied by many workers, and there are differences in the calculation results and the reliability of application effects. Yamada et al . used the least-squares algorithm [ 3 3 ] and maximum likelihood estimation algorithm [ 4 4 ] to analyze 134 Cs and 137 Cs, respectively. Wang et al . also showed the decomposition of pulse waveform using a weighted least-squares algorithm [ 5 5 ]. All of these workers achieved good results, but these methods still need further research when the peaks over lap on a large scale. Guillot applied the filtering technique to extract the full energy peak of the aviation gamma-ray spectrum [ 6 6 ], but the pure spectrum profile method could not separate the overlapping peaks. A genetic algorithm approach for multiple deconvolution was also used to decompose the overlapping peaks [ 7 7 ], and while the result was good, the approach was more costly in computer time than traditional techniques. Parente et al . showed that the Levenberg-Marquardt algorithm provides good results for solving overlapping-peak problems because of its global rapid convergence feature [ 8 8 ]. Singer also used the same technique to fit a smooth polynomial continuum and three Gaussian absorptions to an olivine spectrum [ 9 9 ]. Brown used it to fit hyperspectral data [ 10 10 ]. All of the aforementioned methods separate a wave into two or more waves, and they all stress the importance of a correct model on the accuracy and stability of the curve fitting. Therefore, in this paper, in order to respond to such a situation in the use of a NaI (Tl) scintillation spectrometer, a Levenberg-Marquardt (LM) optimization algorithm with a Gaussian model is proposed for data obtained in the field, and is demonstrated experimentally for the determination of the single-peak area of 137 Cs from the overlapping of three peaks. The LM algorithm is compared with a genetic algorithm (GA). The performances of the two methods are also assessed. 2 BACKGROUND 1 2.1 Gaussian modeling 2 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 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 11 ]. Furthermore, a modified Gaussian successfully modeled the electronic transition absorption in pyroxene, pyroxene mixtures, and olivine spectra [ 12 12 ]. Here, we use Gaussian fitting function to decompose gamma-ray spectra into individual full energy peaks, as described by Eq. (1): <math display="block" id="M1"><mrow><msub><mstyle mathsize="140%" displaystyle="true"><mi>y</mi></mstyle><mi>i</mi></msub><mo stretchy="false">(</mo><msub><mstyle mathsize="140%" displaystyle="true"><mi>x</mi></mstyle><mi>j</mi></msub><mo stretchy="false">)</mo><mo>=</mo><msub><mstyle mathsize="140%" displaystyle="true"><mi>A</mi></mstyle><mi>i</mi></msub><mi>exp</mi><mo stretchy="false">[</mo><mo>−</mo><mn>4</mn><mo stretchy="false">(</mo><mi>ln</mi><mn>2</mn><mo stretchy="false">)</mo><mo>∗</mo><msup><mstyle mathsize="140%" displaystyle="true"><mrow><mo stretchy="false">(</mo><mfrac><mrow><msub><mstyle mathsize="140%" displaystyle="true"><mi>x</mi></mstyle><mi>j</mi></msub><mo>−</mo><msub><mstyle mathsize="140%" displaystyle="true"><mi>p</mi></mstyle><mi>i</mi></msub></mrow><mrow><msub><mstyle mathsize="140%" displaystyle="true"><mi>B</mi></mstyle><mi>i</mi></msub></mrow></mfrac><mo stretchy="false">)</mo></mrow></mstyle><mn>2</mn></msup><mo stretchy="false">]</mo><mo>,</mo></mrow></math> where <math id="M2"><mrow><msub><mstyle mathsize="140%" displaystyle="true"><mi>A</mi></mstyle><mi>i</mi></msub></mrow></math> , <math id="M3"><mrow><msub><mstyle mathsize="140%" displaystyle="true"><mi>B</mi></mstyle><mi>i</mi></msub></mrow></math> , and <math id="M4"><mrow><msub><mstyle mathsize="140%" displaystyle="true"><mi>p</mi></mstyle><mi>i</mi></msub></mrow></math> are the peak height, the full width at half maximum (FWHM), and the center of the i th peak, respectively, and <math id="M5"><mrow><msub><mstyle mathsize="140%" displaystyle="true"><mi>x</mi></mstyle><mi>j</mi></msub></mrow></math> is the channel number [ 13 13 , 14 14 ]. 2.2 Levenberg-Marquardt algorithm 2 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 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 8 ]. The LM Gaussian fitting algorithm is described as follows: <math display="block" id="M6"><mrow><mo stretchy="false">(</mo><msup><mstyle mathsize="140%" displaystyle="true"><mrow><mi>J</mi><mo stretchy="false">(</mo><mi>a</mi><mo stretchy="false">)</mo></mrow></mstyle><mi>T</mi></msup><mi>J</mi><mo stretchy="false">(</mo><mi>a</mi><mo stretchy="false">)</mo><mo>+</mo><mi>μ</mi><mi>I</mi><mo stretchy="false">)</mo><mi>d</mi><mo>=</mo><mo>−</mo><msup><mstyle mathsize="140%" displaystyle="true"><mrow><mi>J</mi><mo stretchy="false">(</mo><mi>a</mi><mo stretchy="false">)</mo></mrow></mstyle><mi>T</mi></msup><mi>r</mi><mo stretchy="false">(</mo><mi>a</mi><mo stretchy="false">)</mo><mo>,</mo></mrow></math> <math display="block" id="M7"><mrow><msup><mi>χ</mi><mn>2</mn></msup><mo stretchy="false">(</mo><mi>a</mi><mo stretchy="false">)</mo><mo>=</mo><msup><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><mrow><mo stretchy="false">[</mo><mfrac><mrow><msub><mi>y</mi><mi>i</mi></msub><mo>−</mo><mi>y</mi><mo stretchy="false">(</mo><msub><mi>x</mi><mi>i</mi></msub><mo>,</mo><mi>a</mi><mo stretchy="false">)</mo></mrow><mrow><msub><mi>σ</mi><mi>i</mi></msub></mrow></mfrac><mo stretchy="false">]</mo></mrow></mstyle></mrow><mn>2</mn></msup><mo>,</mo></mrow></math> Where <math id="M8"><mrow><mi>μ</mi><mo>></mo><mn>0</mn></mrow></math> , <math id="M9"><mi>I</mi></math> is a unit vector, <math id="M10"><mi>d</mi></math> is the mobile vector, and <math id="M11"><mrow><mi>a</mi><mo>=</mo><msup><mstyle mathsize="140%" displaystyle="true"><mrow><mo stretchy="false">(</mo><mi>A</mi><mo>,</mo><mi>B</mi><mo>,</mo><mi>p</mi><mo stretchy="false">)</mo></mrow></mstyle><mi>T</mi></msup></mrow></math> . We introduce a function <math id="M12"><mrow><msub><mstyle mathsize="140%" displaystyle="true"><mi>r</mi></mstyle><mi>i</mi></msub><mo stretchy="false">(</mo><mi>a</mi><mo stretchy="false">)</mo><mo>=</mo><msub><mstyle mathsize="140%" displaystyle="true"><mi>y</mi></mstyle><mi>i</mi></msub><mo>−</mo><mi>y</mi><mo stretchy="false">(</mo><msub><mstyle mathsize="140%" displaystyle="true"><mi>x</mi></mstyle><mi>i</mi></msub><mo>,</mo><mi>a</mi><mo stretchy="false">)</mo><mo>,</mo><msup><mstyle mathsize="140%" displaystyle="true"><mrow/></mstyle><mrow/></msup><mi>i</mi><mo>=</mo><mn>1</mn><mo>,</mo><mn>2</mn><mo>,</mo><mo>⋯</mo><mo>,</mo><mi>n</mi></mrow></math> , and let <math id="M13"><mrow><mi>r</mi><mo stretchy="false">(</mo><mi>a</mi><mo stretchy="false">)</mo><mo>=</mo><mo stretchy="false">(</mo><msub><mstyle mathsize="140%" displaystyle="true"><mi>r</mi></mstyle><mn>1</mn></msub><mo stretchy="false">(</mo><mi>a</mi><mo stretchy="false">)</mo><mo>;</mo><msub><mstyle mathsize="140%" displaystyle="true"><mi>r</mi></mstyle><mn>2</mn></msub><mo stretchy="false">(</mo><mi>a</mi><mo stretchy="false">)</mo><mo>;</mo><mo>⋯</mo><mo>;</mo><msub><mstyle mathsize="140%" displaystyle="true"><mi>r</mi></mstyle><mi>n</mi></msub><mo stretchy="false">(</mo><mi>a</mi><mo stretchy="false">)</mo><mo stretchy="false">)</mo></mrow></math> ; <math id="M14"><mi>J</mi></math> is the Jacobi matrix of <math id="M15"><mrow><mi>r</mi><mo stretchy="false">(</mo><mi>a</mi><mo stretchy="false">)</mo></mrow></math> . <math id="M16"><mrow><msub><mstyle mathsize="140%" displaystyle="true"><mi>y</mi></mstyle><mi>i</mi></msub></mrow></math> and <math id="M17"><mrow><mi>y</mi><mo stretchy="false">(</mo><msub><mstyle mathsize="140%" displaystyle="true"><mi>x</mi></mstyle><mi>i</mi></msub><mo>,</mo><mi>a</mi><mo stretchy="false">)</mo></mrow></math> are, respectively, the original and fitted value, <math id="M18"><mrow><msub><mstyle mathsize="140%" displaystyle="true"><mi>σ</mi></mstyle><mi>i</mi></msub></mrow></math> is the standard deviation, and <math id="M19"><mi>n</mi></math> is the total number of channels. The purpose of the algorithm is to make <math id="M20"><mrow><msup><mstyle mathsize="140%" displaystyle="true"><mi>χ</mi></mstyle><mn>2</mn></msup><mo stretchy="false">(</mo><mi>a</mi><mo stretchy="false">)</mo></mrow></math> a minimum by searching for the optimal solution of <math id="M21"><mi>a</mi></math> with the given initial value of <math id="M22"><mi>a</mi></math> , replacing <math id="M23"><mi>a</mi></math> into the Gaussian function after obtaining the optimal value, and then decomposing the overlapping peaks [ 16 16 - 21 21 16-21 ]. A flow chart of the LM algorithm is presented in Fig. 1 Fig. 1 . In this paper, a comparison with a genetic algorithm found in Matlab software was also made. 2.3 Data preprocessing 2 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 <math id="M24"><mrow><mi>y</mi><mo stretchy="false">(</mo><mi>i</mi><mo stretchy="false">)</mo></mrow></math> , letting <math id="M25"><mrow><mi>Z</mi><mo stretchy="false">(</mo><mi>i</mi><mo stretchy="false">)</mo><mo>=</mo><mi>ln</mi><mo stretchy="false">[</mo><mn>1</mn><mo>+</mo><mi>ln</mi><mo stretchy="false">(</mo><mn>1</mn><mo>+</mo><msqrt><mrow><mi>y</mi><mo stretchy="false">(</mo><mi>i</mi><mo stretchy="false">)</mo><mo>+</mo><mn>1</mn></mrow></msqrt><mo stretchy="false">)</mo><mo stretchy="false">]</mo></mrow></math> , where <math id="M26"><mi>i</mi></math> 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 w th iteration, we calculated the new value in channel <math id="M27"><mi>i</mi></math> as <math id="M28"><mrow><msub><mstyle mathsize="140%" displaystyle="true"><mi>Z</mi></mstyle><mi>w</mi></msub><mo stretchy="false">(</mo><mi>i</mi><mo stretchy="false">)</mo><mo>=</mo><mi>min</mi><mo>{</mo><msub><mstyle mathsize="140%" displaystyle="true"><mi>Z</mi></mstyle><mrow><mi>w</mi><mo>−</mo><mn>1</mn></mrow></msub><mo stretchy="false">(</mo><mi>i</mi><mo stretchy="false">)</mo><mo>,</mo><mtext>  </mtext><mo stretchy="false">[</mo><msub><mstyle mathsize="140%" displaystyle="true"><mi>Z</mi></mstyle><mrow><mi>w</mi><mo>−</mo><mn>1</mn></mrow></msub><mo stretchy="false">(</mo><mi>i</mi><mo>−</mo><mi>w</mi><mo stretchy="false">)</mo><mo>+</mo><msub><mstyle mathsize="140%" displaystyle="true"><mi>Z</mi></mstyle><mrow><mi>w</mi><mo>−</mo><mn>1</mn></mrow></msub><mo stretchy="false">(</mo><mi>i</mi><mo>+</mo><mi>w</mi><mo stretchy="false">)</mo><mo stretchy="false">]</mo><mo>/</mo><mn>2</mn><mo>}</mo><mo>,</mo><msup><mstyle mathsize="140%" displaystyle="true"><mrow/></mstyle><mrow/></msup><mi>w</mi><mo>=</mo><mn>1</mn><mo>,</mo><mo>⋯</mo><mo>,</mo><mi>N</mi></mrow></math> , where <math id="M29"><mi>N</mi></math> is the number of iterations. Finally, we obtained the background channel counts, <math id="M30"><mrow><mi>B</mi><mo stretchy="false">(</mo><mi>i</mi><mo stretchy="false">)</mo><mo>=</mo><msup><mstyle mathsize="140%" displaystyle="true"><mi>Z</mi></mstyle><mrow><mo>−</mo><mn>1</mn></mrow></msup><mo stretchy="false">(</mo><mi>i</mi><mo stretchy="false">)</mo></mrow></math> [ 22 22 ]. Then, the net peak area can be easily determined. 3 METHODS 1 3.1 In situ gamma-ray spectrometry 2 A NaI (Tl) scintillation spectrometer and a 137 Cs point source (3.7×10 4 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 23 ]. A schematic of the simulation experiment is shown in Fig. 2 Fig. 2 . The NaI (Tl) spectrometer was placed on the ground, and a 137 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. 3.2 Lack-of-fit test 2 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: <math display="block" id="M31"><mtable columnalign="left"><mtr><mtd><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mi>m</mi></munderover><mrow><msubsup><mstyle mathsize="140%" displaystyle="true"><mi>ε</mi></mstyle><mrow><mi>i</mi><mi>j</mi></mrow><mn>2</mn></msubsup></mrow></mstyle></mrow></mstyle><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mi>m</mi></munderover><mrow><msup><mstyle mathsize="140%" displaystyle="true"><mrow><mo stretchy="false">(</mo><msub><mstyle mathsize="140%" displaystyle="true"><mi>y</mi></mstyle><mrow><mi>i</mi><mi>j</mi></mrow></msub><mo>−</mo><msub><mstyle mathsize="140%" displaystyle="true"><mover accent="true"><mi>y</mi><mo>^</mo></mover></mstyle><mi>i</mi></msub><mo stretchy="false">)</mo></mrow></mstyle><mn>2</mn></msup></mrow></mstyle></mrow></mstyle></mtd></mtr><mtr><mtd><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mi>m</mi></munderover><mrow><msup><mstyle mathsize="140%" displaystyle="true"><mrow><mo stretchy="false">(</mo><msub><mstyle mathsize="140%" displaystyle="true"><mi>y</mi></mstyle><mrow><mi>i</mi><mi>j</mi></mrow></msub><mo>−</mo><msub><mstyle mathsize="140%" displaystyle="true"><mover accent="true"><mi>y</mi><mo>¯</mo></mover></mstyle><mrow><mi>i</mi><mo>•</mo></mrow></msub><mo stretchy="false">)</mo></mrow></mstyle><mn>2</mn></msup></mrow></mstyle></mrow></mstyle><mo>+</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><mrow><mi>m</mi><msup><mstyle mathsize="140%" displaystyle="true"><mrow><mo stretchy="false">(</mo><msub><mstyle mathsize="140%" displaystyle="true"><mover accent="true"><mi>y</mi><mo>¯</mo></mover></mstyle><mrow><mi>i</mi><mo>•</mo></mrow></msub><mo>−</mo><msub><mstyle mathsize="140%" displaystyle="true"><mover accent="true"><mi>y</mi><mo>^</mo></mover></mstyle><mi>i</mi></msub><mo stretchy="false">)</mo></mrow></mstyle><mn>2</mn></msup></mrow></mstyle><mo>,</mo><msup><mstyle mathsize="140%" displaystyle="true"><mrow/></mstyle><mrow/></msup><mi>i</mi><mo>=</mo><mn>1</mn><mo>,</mo><mo>⋯</mo><mo>,</mo><mi>n</mi><mo>;</mo><msup><mstyle mathsize="140%" displaystyle="true"><mrow/></mstyle><mrow/></msup><mi>j</mi><mo>=</mo><mn>1</mn><mo>,</mo><mo>⋯</mo><mo>,</mo><mi>m</mi><mo>,</mo></mtd></mtr></mtable></math> where <math id="M32"><mi>n</mi></math> is the total channel number, <math id="M33"><mi>m</mi></math> is the number of times the experiment is repeated, <math id="M34"><mrow><msub><mstyle mathsize="140%" displaystyle="true"><mi>y</mi></mstyle><mrow><mi>i</mi><mi>j</mi></mrow></msub></mrow></math> and <math id="M35"><mrow><msub><mstyle mathsize="140%" displaystyle="true"><mover accent="true"><mi>y</mi><mo>^</mo></mover></mstyle><mi>i</mi></msub></mrow></math> are the observed value and fitted value, respectively, and <math id="M36"><mrow><msub><mstyle mathsize="140%" displaystyle="true"><mover accent="true"><mi>y</mi><mo>¯</mo></mover></mstyle><mrow><mi>i</mi><mo>•</mo></mrow></msub></mrow></math> is the average of observed y values in channel <math id="M37"><mi>i</mi></math> . The first component is called pure error. The second component is called lack of fit [ 24 24 ]. It then follows that, statistically, <math display="block" id="M38"><mrow><mi>F</mi><mo>=</mo><mfrac><mrow><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></msubsup><mrow><msup><mstyle mathsize="140%" displaystyle="true"><mrow><mi>m</mi><mo stretchy="false">(</mo><msub><mstyle mathsize="140%" displaystyle="true"><mover accent="true"><mi>y</mi><mo>¯</mo></mover></mstyle><mrow><mi>i</mi><mo>•</mo></mrow></msub><mo>−</mo><msub><mstyle mathsize="140%" displaystyle="true"><mover accent="true"><mi>y</mi><mo>^</mo></mover></mstyle><mi>i</mi></msub><mo stretchy="false">)</mo></mrow></mstyle><mn>2</mn></msup><mo>/</mo><mo stretchy="false">(</mo><mi>n</mi><mo>−</mo><mi>p</mi><mo stretchy="false">)</mo></mrow></mstyle></mrow><mrow><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></msubsup><mrow><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mi>m</mi></msubsup><mrow><msup><mstyle mathsize="140%" displaystyle="true"><mrow><mo stretchy="false">(</mo><msub><mstyle mathsize="140%" displaystyle="true"><mi>y</mi></mstyle><mrow><mi>i</mi><mi>j</mi></mrow></msub><mo>−</mo><msub><mstyle mathsize="140%" displaystyle="true"><mover accent="true"><mi>y</mi><mo>¯</mo></mover></mstyle><mrow><mi>i</mi><mo>•</mo></mrow></msub><mo stretchy="false">)</mo></mrow></mstyle><mn>2</mn></msup><mo>/</mo><mo stretchy="false">(</mo><mi>N</mi><mo>−</mo><mi>n</mi><mo stretchy="false">)</mo></mrow></mstyle></mrow></mstyle></mrow></mfrac><mo>,</mo><msup><mstyle mathsize="140%" displaystyle="true"><mrow/></mstyle><mrow/></msup><mi>N</mi><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><mi>m</mi></mstyle><mo>,</mo></mrow></math> where <math id="M39"><mi>N</mi></math> is the total number of observations and <math id="M40"><mi>p</mi></math> is the number of parameters in the model. The degrees of freedom of pure error and lack of fit are <math id="M41"><mrow><mi>N</mi><mo>−</mo><mi>n</mi></mrow></math> and <math id="M42"><mrow><mi>n</mi><mo>−</mo><mi>p</mi></mrow></math> , 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 <math id="M43"><mrow><msub><mstyle mathsize="140%" displaystyle="true"><mi>d</mi></mstyle><mn>1</mn></msub><mo>=</mo><mi>n</mi><mo>−</mo><mi>p</mi></mrow></math> and <math id="M44"><mrow><msub><mstyle mathsize="140%" displaystyle="true"><mi>d</mi></mstyle><mn>2</mn></msub><mo>=</mo><mi>N</mi><mo>−</mo><mi>n</mi></mrow></math> . 3.3 Signal-to-noise ratio 2 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 <math display="block" id="M45"><mrow><mtext>SNR(dB</mtext><mo stretchy="false">)</mo><mo>=</mo><mn>20</mn><mo>×</mo><msub><mstyle mathsize="140%" displaystyle="true"><mrow><mi>log</mi></mrow></mstyle><mrow><mn>10</mn></mrow></msub><mo stretchy="false">(</mo><mfrac><mrow><msub><mstyle mathsize="140%" displaystyle="true"><mi>A</mi></mstyle><mrow><mtext>peak</mtext></mrow></msub></mrow><mrow><msub><mstyle mathsize="140%" displaystyle="true"><mi>A</mi></mstyle><mrow><mtext>background</mtext></mrow></msub></mrow></mfrac><mo stretchy="false">)</mo><mo>,</mo></mrow></math> where <math id="M46"><mrow><msub><mstyle mathsize="140%" displaystyle="true"><mi>A</mi></mstyle><mrow><mtext>peak</mtext></mrow></msub></mrow></math> and <math id="M47"><mrow><msub><mstyle mathsize="140%" displaystyle="true"><mi>A</mi></mstyle><mrow><mtext>background</mtext></mrow></msub></mrow></math> are the amplitude of the 137 Cs peak and background, respectively. The smaller the SNR, the lower the amplitude of the 137 Cs peak. When the value of SNR is too small, smaller than even −30 dB, the amplitude of the 137 Cs peak is smaller than 10 counts. It is difficult to say that the peak of 137 Cs was the source for this; assuming so may lead to incorrect estimation. 3.4 Evaluation methods 2 The relative fitting error of the total net area <math id="M48"><mi>δ</mi></math> , residual standard deviation <math id="M49"><mi>S</mi></math> , and fitting goodness <math id="M50"><mrow><msup><mstyle mathsize="140%" displaystyle="true"><mi>R</mi></mstyle><mn>2</mn></msup></mrow></math> were also used to evaluate the curve-fitting effect. Their mathematical expressions are <math display="block" id="M51"><mrow><mi>δ</mi><mo>=</mo><mrow><mo>|</mo><mrow><mfrac><mrow><msub><mstyle mathsize="140%" displaystyle="true"><mrow><mi>A</mi><mn>2</mn></mrow></mstyle><mrow><mi>T</mi><mi>h</mi><mo>+</mo><mi>R</mi><mi>a</mi><mo>+</mo><mi>C</mi><mi>s</mi></mrow></msub></mrow><mrow><msub><mstyle mathsize="140%" displaystyle="true"><mrow><mi>A</mi><mn>1</mn></mrow></mstyle><mrow><mi>T</mi><mi>h</mi><mo>+</mo><mi>R</mi><mi>a</mi><mo>+</mo><mi>C</mi><mi>s</mi></mrow></msub></mrow></mfrac><mo>−</mo><mn>1</mn></mrow><mo>|</mo></mrow><mo>×</mo><mn>100</mn><mi>%</mi><mi>，</mi></mrow></math> <math display="block" id="M52"><mrow><mi>S</mi><mo>=</mo><msqrt><mrow><mfrac><mrow><mstyle displaystyle="true"><mo>∑</mo><mrow><msup><mstyle mathsize="140%" displaystyle="true"><mrow><mo stretchy="false">(</mo><mi>y</mi><mo>−</mo><mover accent="true"><mi>y</mi><mo>^</mo></mover><mo stretchy="false">)</mo></mrow></mstyle><mn>2</mn></msup></mrow></mstyle></mrow><mrow><mi>n</mi><mo>−</mo><mn>2</mn></mrow></mfrac></mrow></msqrt><mi>，</mi></mrow></math> <math display="block" id="M53"><mrow><msup><mstyle mathsize="140%" displaystyle="true"><mi>R</mi></mstyle><mn>2</mn></msup><mo>=</mo><mn>1</mn><mo>−</mo><mfrac><mrow><mstyle displaystyle="true"><mo>∑</mo><mrow><msup><mstyle mathsize="140%" displaystyle="true"><mrow><mo stretchy="false">(</mo><mi>y</mi><mo>−</mo><mover accent="true"><mi>y</mi><mo>^</mo></mover><mo stretchy="false">)</mo></mrow></mstyle><mn>2</mn></msup></mrow></mstyle></mrow><mrow><mstyle displaystyle="true"><mo>∑</mo><mrow><msup><mstyle mathsize="140%" displaystyle="true"><mrow><mo stretchy="false">(</mo><mi>y</mi><mo>−</mo><mover accent="true"><mi>y</mi><mo>¯</mo></mover><mo stretchy="false">)</mo></mrow></mstyle><mn>2</mn></msup></mrow></mstyle></mrow></mfrac><mi>，</mi></mrow></math> where <math id="M54"><mrow><msub><mstyle mathsize="140%" displaystyle="true"><mrow><mi>A</mi><mn>1</mn></mrow></mstyle><mrow><mi>T</mi><mi>h</mi><mo>+</mo><mi>R</mi><mi>a</mi><mo>+</mo><mi>C</mi><mi>s</mi></mrow></msub></mrow></math> and <math id="M55"><mrow><msub><mstyle mathsize="140%" displaystyle="true"><mrow><mi>A</mi><mn>2</mn></mrow></mstyle><mrow><mi>T</mi><mi>h</mi><mo>+</mo><mi>R</mi><mi>a</mi><mo>+</mo><mi>C</mi><mi>s</mi></mrow></msub></mrow></math> are the total net area of the three full energy peaks at 583, 609, and 662 keV before fitting and after fitting, respectively, <math id="M56"><mi>y</mi></math> is the original value, <math id="M57"><mover accent="true"><mi>y</mi><mo>^</mo></mover></math> is the fitting value, and <math id="M58"><mover accent="true"><mi>y</mi><mo>¯</mo></mover></math> is the average of the original value. The closer <math id="M59"><mi>δ</mi></math> and <math id="M60"><mi>S</mi></math> are to 0 and the closer <math id="M61"><mrow><msup><mstyle mathsize="140%" displaystyle="true"><mi>R</mi></mstyle><mn>2</mn></msup></mrow></math> is to 1, the better the fitting curve. 4 RESULTS AND DISCUSSIONS 1 4.1 137 Cs source experiment 2 Based on the data from the 137 Cs point source simulation experiment ( Fig. 2 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 Fig. 3 and Table 1 Table 1 . Figure 3 Figure 3 and Table 1 Table 1 show that the net area of the 137 Cs peak (662 keV) became smaller as the source-to-detector distance increased, while the full energy peak areas of 208 Tl (583 keV) and 214 Bi (609 keV) emerged. The decomposed net full energy peak areas of 208 Tl and 214 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. 1 2 Algorithm 1 2 Distance of source to detector (cm) 3 1 583 keV ( 208 Tl) 3 1 609 keV ( 214 Bi) 3 1 662 keV ( 137 Cs) 1 1 Peakheight(cps) 1 1 FWHM 1 1 Net area (cps) 1 1 Peak height (cps) 1 1 FWHM 1 1 Net area (cps) 1 1 Peak height (cps) 1 1 FWHM 1 1 Net area (cps) 1 6 L-M 1 1 No source 1 1 0.099 1 1 11.821 1 1 1.137 1 1 0.071 1 1 10.747 1 1 0.811 1 1 0.016 1 1 3.568 1 1 0.060 1 1 18 1 1 0.070 1 1 9.336 1 1 1.016 1 1 0.175 1 1 13.544 1 1 2.527 1 1 1.033 1 1 14.010 1 1 15.398 1 1 24 1 1 0.099 1 1 9.612 1 1 1.138 1 1 0.063 1 1 10.763 1 1 0.716 1 1 0.347 1 1 14.791 1 1 5.468 1 1 30 1 1 0.087 1 1 10.167 1 1 1.042 1 1 0.073 1 1 11.232 1 1 0.874 1 1 0.117 1 1 10.769 1 1 1.343 1 1 34 1 1 0.097 1 1 9.854 1 1 0.989 1 1 0.088 1 1 9.607 1 1 0.896 1 1 0.052 1 1 11.358 1 1 0.627 1 1 38 1 1 0.081 1 1 11.458 1 1 1.018 1 1 0.056 1 1 11.769 1 1 0.707 1 1 0.011 1 1 6.243 1 1 0.072 1 6 GA 1 1 No source 1 1 0.099 1 1 11.822 1 1 1.137 1 1 0.071 1 1 10.747 1 1 0.811 1 1 0.016 1 1 3.566 1 1 0.060 1 1 18 1 1 0.069 1 1 9.650 1 1 1.000 1 1 0.175 1 1 13.551 1 1 2.519 1 1 1.033 1 1 13.998 1 1 15.392 1 1 24 1 1 0.099 1 1 9.619 1 1 1.139 1 1 0.063 1 1 10.783 1 1 0.718 1 1 0.347 1 1 14.790 1 1 5.467 1 1 30 1 1 0.087 1 1 10.172 1 1 1.043 1 1 0.073 1 1 11.246 1 1 0.875 1 1 0.117 1 1 10.764 1 1 1.343 1 1 34 1 1 0.097 1 1 9.860 1 1 0.989 1 1 0.088 1 1 9.617 1 1 0.896 1 1 0.052 1 1 11.350 1 1 0.627 1 1 38 1 1 0.081 1 1 11.464 1 1 1.018 1 1 0.056 1 1 11.786 1 1 0.707 1 1 0.011 1 1 6.214 1 1 0.072 4.2 Discussions 2 Table 2 Table 2 presents the results of evaluating the data in Table 1 Table 1 for the decomposition and fitting of 137 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 1 =62 and d 2 =71, F 0.05 ( d 1 , d 2 )=1.498, and all of the F values obtained and presented in Table 2 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 137 Cs peak is becoming smaller and smaller. The relative error of total net area <math id="M62"><mi>δ</mi></math> and residual standard deviation <math id="M63"><mi>S</mi></math> were no more than 5% and 0.04, respectively, and the fitting goodness <math id="M64"><mrow><msup><mstyle mathsize="140%" displaystyle="true"><mi>R</mi></mstyle><mn>2</mn></msup></mrow></math> 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 137 Cs when the SNR value of in Table 1 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 137 Cs peak of the two distances may lead to incorrect estimation in Fig. 2 Fig. 2 and Table 1 Table 1 . 1 1 Algorithm 1 1 Distance of source to detector (cm) 1 1 F 1 1 <math id="M65"><mrow><mi>δ</mi><mo>/</mo><mi>%</mi></mrow></math>   1 1 <math id="M66"><mrow><mi>S</mi><mo>/</mo><mo>×</mo><msup><mstyle mathsize="140%" displaystyle="true"><mrow><mn>10</mn></mrow></mstyle><mrow><mo>−</mo><mn>2</mn></mrow></msup></mrow></math>   1 1 <math id="M67"><mrow><msup><mstyle mathsize="140%" displaystyle="true"><mi>R</mi></mstyle><mn>2</mn></msup><mo>/</mo><mi>%</mi></mrow></math>   1 1 SNR (dB) 1 6 LM 1 1 No source 1 1 0.072 1 1 4.69 1 1 0.49 1 1 98.06 1 1 −30.9 1 1 18 1 1 0.646 1 1 3.51 1 1 3.93 1 1 98.55 1 1 7.3 1 1 24 1 1 0.229 1 1 3.63 1 1 1.30 1 1 98.52 1 1 −0.7 1 1 30 1 1 0.070 1 1 4.63 1 1 0.40 1 1 98.77 1 1 −10.2 1 1 34 1 1 0.101 1 1 1.46 1 1 0.52 1 1 98.24 1 1 −15.8 1 1 38 1 1 0.033 1 1 2.11 1 1 0.33 1 1 98.99 1 1 −30.3 1 6 GA 1 1 No source 1 1 0.072 1 1 4.70 1 1 0.49 1 1 98.06 1 1 −30.9 1 1 18 1 1 0.644 1 1 3.35 1 1 3.93 1 1 98.56 1 1 7.3 1 1 24 1 1 0.230 1 1 3.65 1 1 1.30 1 1 98.52 1 1 −0.7 1 1 30 1 1 0.070 1 1 4.65 1 1 0.40 1 1 98.77 1 1 −10.2 1 1 34 1 1 0.101 1 1 1.48 1 1 0.52 1 1 98.24 1 1 −15.8 1 1 38 1 1 0.033 1 1 2.16 1 1 0.33 1 1 98.99 1 1 −30.3 The comparison between the LM algorithm and genetic algorithm is presented in Table 1 Table 1 and Table 2 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. 5 CONCLUSIONS 1 A 137 Cs point source simulation experiment was undertaken and the results were reported in this paper. The net area of 137 Cs was extracted by use of a Levenberg-Marquardt Gaussian fitting algorithm, and the performance of the algorithm was compared with that of a genetic algorithm. After analyzing the data, the following conclusions have been drawn: 1. The extraction of the net area of 137 Cs is feasible using a Levenberg-Marquardt Gaussian fitting algorithm for in situ NaI (Tl) gamma-ray spectrum data. 2. The results of the application of our approach based on a Levenberg-Marquardt algorithm yielded similar results compared to other studies based on a genetic algorithm. The Levenberg-Marquardt algorithm gives better performance than the genetic algorithm for a large quantity of data. 3. The method used in this paper gives high performance, but the results may lead to incorrect estimation when the signal-to-noise ratio is smaller than −30 dB. REFERENCES Steinhauser G , Brandl A , Johnson T E . Comparison of the Chernobyl and Fukushima nuclear accidents: A review of the environmental impacts . Sci Total Environ , 2014 , 470-471 : 800 - 817 . DOI: 10.1016/j.scitotenv.2013.10.029 http://doi.org/10.1016/j.scitotenv.2013.10.029 Comparison of the Chernobyl and Fukushima nuclear accidents: A review of the environmental impacts Saegusa J , Kikuta Y , Akino H . Observation of gamma-rays from fallout collected at Ibaraki, Japan, during the Fukushima nuclear accident . Appl Radiat Isot , 2013 , 77 : 56 - 60 . DOI: 10.1016/j.apradiso.2013.02.018 http://doi.org/10.1016/j.apradiso.2013.02.018 Observation of gamma-rays from fallout collected at Ibaraki, Japan, during the Fukushima nuclear accident Yamada T , Ishizu H , Kawada Y . A simple method for activity determination of 134Cs and 137Cs in foodstuffs using NaI(Tl) scintillation spectrometer . Appl Radiat Isot , 2013 , 81 : 353 - 355 . DOI: 10.1016/j.apradiso.2013.03.048 http://doi.org/10.1016/j.apradiso.2013.03.048 A simple method for activity determination of 134Cs and 137Cs in foodstuffs using NaI(Tl) scintillation spectrometer Yamada T , Takano N . Doublet peak area determination in NaI(Tl) scintillation spectrometry using maximum likelihood estimation . Appl Radiat Isot , 2014 , 87 : 407 - 409 . DOI: 10.1016/j.apradiso.2013.11.069 http://doi.org/10.1016/j.apradiso.2013.11.069 Doublet peak area determination in NaI(Tl) scintillation spectrometry using maximum likelihood estimation Wang L , Xu L S , Feng S T , et al . Mult-Gaussian fitting for pulse waveform using Weighted Least Squares and multi-criteria decision making method . Computers in biology and medicine , 2013 , 43 : 1661 - 1672 . DOI: 10.1016/j.compbiomed.2013.08.004 http://doi.org/10.1016/j.compbiomed.2013.08.004 Mult-Gaussian fitting for pulse waveform using Weighted Least Squares and multi-criteria decision making method Guillot L . Extraction of full absorption peaks in airborne gamma-spectrometry by filtering techniques coupled with a study of the derivatives . Comparison with the window method. J Environ Radioactiv , 2001 , 53 : 381 - 398 . DOI: 10.1016/S0265-931X(00)00144-2 http://doi.org/10.1016/S0265-931X(00)00144-2 Extraction of full absorption peaks in airborne gamma-spectrometry by filtering techniques coupled with a study of the derivatives Garcia-Talavera M , Ulicny B . A genetic algorithm approach for multiplet deconvolution in γ-ray spectra . Nucl Instrum Methods Phys Res Sect A , 2003 , 3 : 585 - 594 . DOI: 10.1016/S0168-9002(03)02052-7 http://doi.org/10.1016/S0168-9002(03)02052-7 A genetic algorithm approach for multiplet deconvolution in γ-ray spectra Parente M , Makarewicz H D , Bishop J L . Decomposition of mineral absorption bands using nonlinear least squares curve fitting: Application to Martian meteorites and CRISM data . Planet Space Sci , 2011 , 59 : 423 - 442 . DOI: 10.1016/j.pss.2011.01.009 http://doi.org/10.1016/j.pss.2011.01.009 Decomposition of mineral absorption bands using nonlinear least squares curve fitting: Application to Martian meteorites and CRISM data Singer R B . Near-infrared spectral reflectance of mineral mixtures: Systematic combinations of pyroxenes, olivine, and iron oxides . J Geophys Res , 1981 , 86 : 7967 - 7982 . DOI: 10.1029/JB086iB09p07967 http://doi.org/10.1029/JB086iB09p07967 Near-infrared spectral reflectance of mineral mixtures: Systematic combinations of pyroxenes, olivine, and iron oxides Brown A J . Spectral curve fitting for automatic hyperspectral data analysis . IEEE T Geosci Remote , 2006 , 44 : 1601 - 1608 . DOI: 10.1109/TGRS.2006.870435 http://doi.org/10.1109/TGRS.2006.870435 Spectral curve fitting for automatic hyperspectral data analysis Miyamoto M , Arai T , Komatsu M , et al . Evaluation of a curve-fitting method for diffuse reflectance spectra in the UV-Visible-NIR wavelength region . Polar Sci , 2009 , 3 : 110 - 116 . DOI: 10.1016/j.polar.2009.06.004 http://doi.org/10.1016/j.polar.2009.06.004 Evaluation of a curve-fitting method for diffuse reflectance spectra in the UV-Visible-NIR wavelength region Sunshine J M , Pieters C M , Pratt S F . Deconvolution of mineral absorption bands: An improved approach . J Geophys Res , 1990 , 95 : 6955 - 6966 . DOI: 10.1029/JB095iB05p06955 http://doi.org/10.1029/JB095iB05p06955 Deconvolution of mineral absorption bands: An improved approach Fu C , Wang N P . Gaussian fitting in gamma-ray spectral decomposition . Nucl Sci Tech , 2010 , 21 : 214 - 217 . Gaussian fitting in gamma-ray spectral decomposition Huang H Q , Shao Q , Ding W C , et al . Simulation of gamma spectrum-shifting based on the parameter adjustment of Gaussian function space . Nucl Sci Tech , 2013 , 24 : 060405 . DOI: 10.13538/j.1001-8042/nst.2013.06.009 http://doi.org/10.13538/j.1001-8042/nst.2013.06.009 Simulation of gamma spectrum-shifting based on the parameter adjustment of Gaussian function space Gill P R , Murray W , Wright M H . The Levenberg-Marquardt method . Practical Optimization, Academic Press , London , 1981 . The Levenberg-Marquardt method Dan H , Yamashita N , Fukushima M . Convergence Properties of the Inexact Levenberg-Marquardt Method under Local Error Bound Conditions . Optim Methods Softw , 2002 , 17 : 605 - 626 . DOI: 10.1080/1055678021000049345 http://doi.org/10.1080/1055678021000049345 Convergence Properties of the Inexact Levenberg-Marquardt Method under Local Error Bound Conditions Behling R , Fischer A , Herrich M , et al . A Levenberg-Marquardt method with approximate projections . Comput Optim Appl , 2013 , 59 : 5 - 26 . DOI: 10.1007/s10589-013-9573-4 http://doi.org/10.1007/s10589-013-9573-4 A Levenberg-Marquardt method with approximate projections Ling C , Wang G F , He H J . A new Levenberg-Marquardt type algorithm for solving nonsmooth constrained equations . Appl Math Comput , 2014 , 229 : 107 - 122 . DOI: 10.1016/j.amc.2013.12.015 http://doi.org/10.1016/j.amc.2013.12.015 A new Levenberg-Marquardt type algorithm for solving nonsmooth constrained equations He Y D , Ma C F , Fan B . A corrected Levenberg-Marquardt algorithm with a nonmonotone line search for the system of nonlinear equations . Appl Math Comput , 2015 , 260 : 159 - 169 . DOI: 10.1016/j.amc.2015.03.076 http://doi.org/10.1016/j.amc.2015.03.076 A corrected Levenberg-Marquardt algorithm with a nonmonotone line search for the system of nonlinear equations Xiao Y . A higher-order Levenberg-Marquardt method for nonlinear equations . Appl Math Comput , 2013 , 219 : 10682 - 10694 . DOI: 10.1016/j.amc.2013.04.033 http://doi.org/10.1016/j.amc.2013.04.033 A higher-order Levenberg-Marquardt method for nonlinear equations Amini K , Rostami F . A modified two steps Levenberg-Marquardt method for nonlinear equations . J Comput Appl Math , 2015 , 288 : 341 - 350 . DOI: 10.1016/j.cam.2015.04.040 http://doi.org/10.1016/j.cam.2015.04.040 A modified two steps Levenberg-Marquardt method for nonlinear equations Ryan C G , Clayton E , Griffin W L , et al . SNIP, a statistics-sensitive background treatment for the quantitative analysis of PIXE spectra in geoscience applications . Nucl Instrum Methods Phys Res Sect B , 1988 , 34 : 396 - 402 . DOI: 10.1016/0168-583X(88)90063-8 http://doi.org/10.1016/0168-583X(88)90063-8 SNIP, a statistics-sensitive background treatment for the quantitative analysis of PIXE spectra in geoscience applications Wang N P , Xiao L , Li C P , et al . Determination of Radioactivity Level of 238U, 232Th and 40K in Surface Medium in Zhuhai City by in-situ Gamma-ray Spectrometry . J Nucl Sci Technol , 2005 , 42 : 888 - 896 . DOI: 10.1080/18811248.2005.9711040 http://doi.org/10.1080/18811248.2005.9711040 Determination of Radioactivity Level of 238U, 232Th and 40K in Surface Medium in Zhuhai City by in-situ Gamma-ray Spectrometry Neter J , Kutner M H , Nachtsheim C J , et al . Applied Linear Statistical Models . Chicago (USA) : McGraw-Hill/Irwin , 1996 , 121 - 122 . Applied Linear Statistical Models 10.1007/s41365-016-0090-x This work was supported by the National Natural Science Foundation of China (No. 41474107) . 2015-07-20 2015-11-15 2015-11-22 2016-07-08 2016-08-20 2021-04-11 © China Science Publishing & Media Ltd. (Science Press), Shanghai Institute of Applied Physics, the Chinese Academy of Sciences, ChineseNuclear Society and Springer Nature Singapore Pte Ltd. 2016 This is a post-peer-review, pre-copyedit version of an article published in Nuclear Science and Techniques. The final authenticated version is available online at: http://dx.doi.org/10.1007/s41365-016-0090-x http://dx.doi.org/10.1007/s41365-016-0090-x . 2016 Nuclear Science and Techniques 04 27