Abstract
NaI(Tl) detectors are frequently operated under unstable temperature conditions when used in an open environment. Temperature changes would result in a peak shift and spectral distortion during measurement. Two easy-to-implement methodologies are proposed to stabilize the measured spectrum without the necessity of adjusting the gain, which are a correction algorithm for temperature-caused peak-shift based on multiple characteristic peak area weighting factors and an interpolation correction algorithm based on multicharacteristic peak sequence. Both of them can be used when the relative channel displacement of characteristic peaks in the spectrum due to temperature changes is not constant. Experimental data obtained under controlled temperature conditions in the laboratory were adopted to correct a spectrum, with joint consideration of some known characteristic peaks, such as 40K, U (214Bi), or Th (208Tl) peaks. Through constructing a reversible temperature coefficient matrix, one can easily obtain the coefficients of the n-th polynomial describing the influence of temperature on peak position, which presents their nonlinear mathematical relationship. Then, corrections of these two effects can also be easily calculated. Comparing the experimental results, peak positions before and after correction, it is proved that the interpolation correction algorithm based on multicharacteristic peak sequence has better correction accuracy, but the temperature-caused peak shift correction algorithm based on the multicharacteristic peak area weighting factor has a shorter calibration time.
1. Introduction
NaI(Tl) scintillation detectors are still widely used in many applications nowadays [1–3]. In most applications, they are used in a laboratory or in an open environment for short-time measurements, where temperature changes are pretty inappreciable [4]. However, when used in an open outdoor environment for a long time, detectors are likely to encounter a wide range of temperature variations [4, 5]. In this case, their performance will be obviously affected either by the temperature effect on the detector crystal itself or the electronic components [6, 7]. Temperature changes may ultimately cause gain instabilities and probably lead to a peak shift and then spectral distortion [8]. There are many reasons for spectrum peak shifts, such as the temperature effect of the photomultiplier tube [4, 7, 9, 10] and the crystal’s light output [11–13], the decay time constant of the scintillator [7, 11, 14], and the temperature characteristics of other associated electronic components [8]. It is difficult to distinguish the primary cause for spectrum peak shift only from the simple temperature effects between crystals and electronics.
There are already several kinds of available spectrum stabilization methods in the literature [4]. The first category is to stabilize the spectrum by adjusting the gain; Shepard et al. achieved it by applying an electronic reference pulse to generate a known equivalent energy in the spectrum with an easily recognized peak position [15]. Pausch et al. attached an external radioactive source to the detector [5, 15], utilized some radioactive isotopes from the natural background, and analyzed the temperature dependence of a light pulse’s decay time using LEDs as reference light sources [5, 16]. Casanovas et al. used some radioactive isotopes from an inner contaminant [4]. Zeng et al. presented a frequency spectrum analysis method for the natural gamma-ray background spectrum by calculating the overall spectral drift of the natural gamma-ray spectrum [17].
The previously described methods are all based on an automatic gain adjusting system, which is not suitable for systems only including analogue gain control. The second category is to stabilize the spectrum without adjusting the gain. Csurgai et al. elaborated on a simple method for energy recalibration of scintillation detectors at different temperatures [18], which can be used in software processing. Bu et al. adopted some simple offline methods to correct the drift of a gamma spectrum [19]. Casanovas et al. proposed some software algorithms [4] to correct spectrum drift by adopting experimental data obtained under controlled temperatures in the laboratory, such as using one known peak (K peak for the NaI(Tl) detector) in the spectrum or using an external source. Under normal circumstances, the spectrum is usually corrected in the channel domain. Wang et al. presented a correction method for gamma spectra based on the system transformation theory of random signals. Using this method, the theoretical deposition energy spectrum (i.e., the corrected spectrum) can be derived from the corresponding measured spectrum without any mathematical approximation [20]. Mitra et al. proposed another method based on the assumption that the spectrum obtained at measurement temperature represents the same statistical distribution as that at reference temperature but with different energy scales [21]. Rezaei Moghaddam et al. suggested a “waiting time” before long-time measurements to avoid the channel shifting effects, which finally improved the detector’s resolution [22].
All of these software algorithms are based on the assumption that the relative channel displacement due to temperature changes is approximately the same for all characteristic peaks. In this study, we describe two easy-to-implement methodologies to correct the peak shift of the spectrum without continuously adjusting the gain.
2. Nonlinear Correction Methods of Temperature-Caused Peak Shift
In the past, for most of the correction methods, it has been assumed that the drift direction of each peak and the relative channel displacement due to temperature changes are approximately the same. But the actual situation may be different. When this assumption is not true, those previously studied methods will meet challenges. Channel shifts caused by temperature changes may invalidate the energy calibration and lead to the wrong identification of radionuclides.
2.1. Model Establishment of Temperature and Peak Position
It is assumed here that for a fixed bias and gain, the channel shift of peak position only depends on the temperature.
For the detector, at the constant reference temperature T0 (e.g., T0 = 25°C), the peak position is also constant. Setting the sampling time and then getting the background counts at different preset temperatures Tj, , in the approximate temperature range of 0°C to 50°C, we need to find the peak positions , , of the spectrum after deducting the background at the corresponding temperature Tj,. The relationship between temperature Tj, and peak position can be expressed by binomial expansion as in equation (1).
In equation (1), , and ,, …… are coefficients of the temperature and peak position’s nth polynomial. Equation (2) is the matrix representation of equation (1).
After the experiments are carried out at different temperatures from to , equation (2) can be extended to equation (3) by using nth order polynomial fitting.
The mathematical relationship between temperature and peak position may be nonlinear. After obtaining the experimental data of temperature and peak position, the correction coefficients of nth order polynomial can be quickly obtained by matrix operations. It is necessary to ensure that the constructed temperature coefficient matrix T, which is shown in equation (4), is a square matrix. Therefore, when the n value of the nth order polynomial is chosen to be m − 1, the temperature matrix T becomes an m-order square matrix. When the temperature coefficient matrix T is a reversible nonsingular matrix, the correction coefficients , , …, of the nth order polynomial, shown in vector A, are determined by the following equation (5).
When the temperature coefficient matrix is an irreversible singular matrix, its Moore–Penrose generalized inverse matrix , which is shown in equations (6) and (7), can be used instead. is just an intermediate process matrix.
The unified expression of the correction coefficients vector A of the nth order polynomial is described by equation (8).
2.2. Calculation Method of Temperature-Caused Peak-Drift Correction Amount
Suppose is a peak position at the standard temperature , and is a peak position at temperature . Then, the peak position at temperature needs to be corrected to a new peak position at standard temperature . Equation (9) can be easily deduced from equation (1) replacing by in it.
So, the temperature-caused peak-shift correction amount from to can be calculated by using equation (1) minus equation (9), which is listed in equation (10) and in matrix form in equation (11):
Therefore, the temperature-caused peak shift can be corrected by applying equation (12), in which is the peak position before correction and is the peak position after correction.
2.3. Temperature-Caused Peak-Shift Correction Algorithm Based on the Multicharacteristic Peak Area Weighting Factor
According to the correction algorithm described in chapters 2.1 and 2.2, one appointed characteristic peak can be precisely adjusted to the desired characteristic peak position at a specified temperature. In some cases, the relative channel displacement caused by temperature changes is not the same for every characteristic peak. When the temperature-caused peak-shift correction amount is calculated by using a single characteristic peak, such as the K peak of 40K, a 1.46 MeV gamma-ray, the K peak will be corrected perfectly. But it cannot guarantee that other characteristic peaks such as U peak (214Bi, 1.76 MeV gamma-ray) or Th peak (208Tl, 2.62 MeV gamma-ray) can be synchronously corrected to the corresponding appointed peak position. That means when various radionuclides, such as K, U, and Th elements, and their contents in a mineral sample need to be determined, using a single characteristic peak (e.g., K peak) for correction will sacrifice the correction effect of other characteristic peaks (e.g., Th or U peak), thus finally affecting the measurement accuracy. Therefore, the correction algorithm based on the construction of the reversible temperature coefficient matrix in chapters 2.1 and 2.2 needs to be improved, and then, temperature-caused peak-shift correction algorithm based on a multicharacteristic peak area weighting factor is proposed.
According to the correction algorithm mentioned previously, the reversible temperature coefficient matrix is constructed, and correction coefficients of the temperature of K peak, U peak, and Th peak position models are calculated, respectively. Then, correction amounts , , and are obtained, respectively. The total correction amount of temperature-caused peak-shift correction algorithm based on the multicharacteristic peak area weighting factor can be described by equations (13)–(16).
In equation (13), , , and are characteristic peak area weighting factors, and the sum of the three is equal to 1. , , and are the peak areas of K, U, and Th peaks at standard temperature. The method chosen to solve the area of each characteristic peak should be the same one. This method also allows two characteristic peaks to be selected for correction. At this time, only two corresponding characteristic peaks related parameters are retained in the calculation equation (13). When there is only one characteristic peak chosen for correction, temperature-caused peak-shift correction algorithm based on the multicharacteristic peak area weighting factor just simplifies the method described in chapters 2.1 and 2.2. When = = , equation (13) can be transformed to equations (17)–(19).
From equations (17)–(19), the total correction amount is equal to the correction amount of every single characteristic peak, shown in equation (20).
When the K peak position corresponding to any temperature is corrected to the K peak position corresponding to the standard temperature , the total correction amount , which can be positive or negative, needs to be subtracted, as shown in equation (21). The similar correction formulas for characteristic peaks of U and Th are shown in equations (22) and (23).
2.4. Interpolation Correction Algorithm Based on Multicharacteristic Peak Sequence
Applying an interpolation correction algorithm, the correction amount corresponding to a single characteristic peak correction ( can be , , and ) needs to be calculated first, which has been described in Chapter 2.2. Secondly, at least two characteristic peaks need to be chosen. Using the first chosen characteristic peak to finish the first correction step is to apply the correction amount . Then, an interpolation offset , which can be obtained from equation (24), is used in the second correction step between the first and second appointed characteristic peaks.
In equation (24), the correction amount is used in the second correction step with another chosen characteristic peak. If there are other characteristic peaks existing between the first and second chosen ones in this method, the area calculation accuracy of these characteristic peaks, which are not chosen as corrected characteristic peaks, may be affected to some degree. Therefore, it is better that no other characteristic peaks exist between the first and the second chosen characteristic peaks for correction. The amount can be positive or negative. When the offset is positive, it means to insert channels between the two appointed characteristic peaks. Conversely, when the offset is negative, it means to delete channels between the two appointed characteristic peaks.
In order to reduce the potential adverse effects on peak area calculations, the interpolation interval selected between the two characteristic peak positions cannot overlap with the calculation interval of characteristic peak areas. The count value of each channel inserted is set to be a constant or follow another specified algorithm. In this way, nonimportant peak information is deleted, and important characteristic peak information is retained, which can ensure the accuracy of subsequent quantitative analysis.
3. Experiment and Discussion
In the process of natural gamma-ray measurement, the radioactive series are usually considered to be balanced. Therefore, the activity (or content) of the parents (K, U, or Th) or other offspring can be calculated by measuring the activity (content) of any of the offspring. In this section, the 1.46 MeV gamma-ray of 40K, the 1.76 MeV gamma-ray of 214Bi, and the 2.62 MeV gamma-ray of 208Tl are chosen to represent the typical K peak, U peak, and Th peak.
3.1. Experimental Setup
In this experiment, the two detectors used are both Φ3’ × 3’ NaI(Tl) scintillation detectors with 137Cs (0.662 MeV) gamma-ray. The IED-3000B digital gamma spectrometers with a full width at half maximum (FWHM) of 8.1% (@661 keV, Cs-137) are made in China and are identified as IED_309 and IED_310, respectively. The range of spectra is from 30 keV to 3000 keV. To validate these methods, two sets of 8 spectra for each detector are collected in the approximate temperature range of 0°C to 42°C, which is the manufacturer’s recommended operating range. A refrigerator and an oven are used to control the temperature changes, and all of the temperatures were measured by using a Pt100 temperature probe with an AA tolerance level, which is about ±0.1°C, and temperature values are eventually kept as integers. The experimental data are obtained using a radioactive ore sample as the source. Each spectrum is collected after thermal stability has been achieved (at least 1 h of constant temperature).
It is known that temperature changes can affect the detector’s performance either in the crystal itself or the electronics, which finally can lead to instabilities and result in a peak shift and spectral distortion during measurement. In this study, it is assumed that the peak position drift is only related to temperature. For each detector, two sets of spectra are collected at 8 different temperatures; the first set of spectra is collected using an ore sample of known composition, and the second set measures the background. After obtaining all the gamma spectra of these two detectors, the shift of K and Th peak positions with temperature are analyzed.
3.2. Nonlinear Relationship of Temperature and Characteristic Peak Position
After performing peak searching (e.g., K peak and Th peak) using the same method in the spectra obtained by detector IED_310, it is obvious that the relationship between temperature and corresponding characteristic peak position is nonlinear. Peak positions of K and Th characteristic peaks at different temperatures for IED_310 are recorded in Table 1. Compared with the peak position under standard temperature, K peak and Th peak positions have obviously drifted, which are shown in Figure 1(a). The two curves of the K and Th peak positions at different temperatures are nonlinear.
(a)
(b)
Relative deviations () of the uncorrected peak positions of K or Th characteristic peak positions at other temperatures compared to their reference peak positions at the standard temperature are shown in Figure 1(b).
3.3. Validation of Temperature-Caused Peak-Shift Correction Algorithm by Using a Single Characteristic Peak
In order to better display and compare the peak positions before and after correction, only counts of some continuous channels (from channel 300 to channel 700, covering three characteristic peaks of interest) are selected for display at different temperatures in Figures 2–5.
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
During the experiment, 26°C is preset as the standard temperature. Compared with the corresponding characteristic peak positions of K and Th under the standard temperature, the characteristic peak positions of K and Th peaks at other temperatures have been dramatically shifted.
According to the general fitting rule, after 8 spectra for each detector have been collected in the approximate temperature range of 0°C to 42°C, the 7th degree polynomial is selected to model the nonlinear relationship of temperature and peak position, and then 8th order temperature coefficient square matrix is constructed, which can achieve the best single peak position correction effect. Two correction amounts and are calculated, which are shown in Table 2. Then, equations (11) and (12) are used to correct the characteristic peak position to the one under standard temperature.
Figure 2(a) shows the uncorrected spectrum of an ore sample measured by the NaI(Tl) detector at 26°C (IED_310). Figure 2(b) shows the raw measured spectra as a function of temperature for the NaI(Tl) detector IED_310 (from channel 300 to channel 700). In Figures 2–5, mixed peaks are composed of various energies of uranium series and thorium series, with energies between 2.0–2.3 MeV. Figure 3(a) shows the measured spectra (IED_310) are corrected using the correction amount obtained from the K peak. Under the circumstances, the corrected K peak fits well, but the corrected Th peak may not meet the requirements. Figure 3(b) shows the measured spectra (IED_310) are corrected using the correction amount obtained from the Th peak. In this case, the corrected Th peak fits well, but the corrected K peak may not meet the requirements. In Figures 2–5, the horizontal axis represents the channel number, the vertical axis means the temperature at acquisition, and the colour scale indicates the counting rate, which is red for higher values and blue for lower ones. The red, green, and blue lines on the waterfall plots (bottom) indicate which spectra (top) are represented by way of example. Figure 3 proves that the method based on the assumption that the relative channel displacement due to temperature changes is approximately the same for all channels cannot obtain a satisfied correction effect.
3.4. Validation of the Temperature-Caused Peak-Shift Correction Algorithm Based on the Multicharacteristic Peak Area Weighting Factor
K and Th characteristic peaks are chosen to correct the temperature drift. The correction amount shown in Table 1 is calculated by using the transformation of equations (13)–(16), which takes no account of the U element. Then, we use equations (21) and (23) to correct the characteristic peak position under nonstandard temperature to the one under standard temperature.
Figure 4(a) shows that when adopting temperature-caused peak-shift correction algorithm based on the multicharacteristic peak area weighting factor, the peak positions of the K peak and Th peak have been improved to some extent, especially for the K peak. That is, because the single peak area of the K characteristic peak is larger, and the counting rate is higher. When calculating the total correction amount , the area weighting factor of the K characteristic peak is larger than that of the Th characteristic peak.
3.5. Validation of the Interpolation Correction Algorithm Based on Multicharacteristic Peak Sequence
K and Th characteristic peaks are selected for correction; the correction amount and is calculated first by the method described in chapters 2.1 and 2.2. The is adopted to finish the first correction step, and then the correction amount is obtained from equation (25), which is the interpolation amount between the K and Th characteristic peaks. This interpolation amount can be piecewise linear or constant.
When is positive, it first finds the midpoint position between K and Th characteristic peaks, inserts channels before the midpoint, and inserts the other channels after it. The count value of each channel inserted is equal to the count value at the midpoint channel. When is negative, we need to find the midpoint position between K and Th characteristic peaks, delete channels before the midpoint, and delete the other channels after it. In general, inserting or deleting channels means that signal quantification must then be based on peak height fitting rather than peak integration. In order to improve the peak integration accuracy, the interpolation position should avoid the peak area calculation intervals of characteristic peak areas used for correction.
This method is verified by two detectors with a satisfactory correction effect. For the spectra of IED_310, , , and are calculated respectively, which are shown in Table 1. Also, the corrected spectra are shown in Figure 4(b).
For the spectra of IED_309, , and are shown in Table 3, and the raw measured spectra and the corrected spectra are shown in Figure 5. If U peaks need to be corrected together with the K and Th peaks, the calculated correction amount such as , as its interpolation amount will be needed for the next correction step.
4. Conclusions
In this study, several methods are presented to correct the peak shift without requiring gain adjustment, which are helpful to stabilize the gamma-ray spectrum obtained under unstable temperature conditions. When used in an open environment, all of these methods are applied by software after spectrum acquisition instead of adjusting the gain during the measurements. It is acceptable for detectors that relative deviation (RD%) is less than 3% in absolute value in the temperature range of 5°C and 50°C, or a smaller relative deviation (RD%) is below 2% in absolute value in the normal operating range of 0°C to 40°C. When the relative channel displacement due to temperature changes is not the same for all channels in the spectrum, the average relative deviation (RD%) of K and Th peak position’s channel has been still up to 3.20% after peak position correction if standard peak correction methods based on fitting a single characteristic peak for reference and a polynomial for the nonlinear shift of peak position with temperature is used.
The first method uses multiple characteristic peaks for temperature-caused peak-shift correction by calculating the total correction amount based on multicharacteristic peak area weighting factor, so in the calculation process of correction amount; priority is given to the characteristic peaks with larger peak area and higher count rate. In multiple characteristic peaks selection, the peaks with similar peak areas and count rates are advised to be selected, which will help to improve the comprehensive correction effect. Using a total correction amount to correct the whole spectrum at one time cannot ensure that each characteristic peak can be corrected perfectly. However, a comprehensive multipeak correction effect can be obtained by using it without multiple correction iterations. Furthermore, its average relative deviation (RD%) of the corrected peak position in the channel can bee decreased to 1.93% after correction.
The second method uses an interpolation correction algorithm based on multicharacteristic peak sequence to correct the spectrum. First, a selected characteristic peak is used to finish the first correction step. Then, it calculates the interpolation amount and finds the interpolation position. In order to improve measurement accuracy, the interpolation position should avoid the peak area calculation intervals of characteristic peaks areas used for correction. In this paper, the midpoint position between two characteristic peak is chosen as the interpolation position. Through multiple iterative interpolation correction, its average relative deviation (RD%) has been dramatically decreased, which falls within the range of 0∼0.21% after peak position correction, and the absolute deviation is less than 1 channel.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by the National Key Research and Development Program of China (Grant no. 2017YFC0602105), the Sichuan Science and Technology Program (Grant no. 2021JDRC0112), and the Natural Science Foundation of Sichuan Province of China (Grant no. 2022NSFSC0280).