Scintillation detector response

The energy resolution of the scintillation detector is related to Gaussian broadening and can be described by an unbroadening spectrum convoluted with a Gaussian response function. For discrete data, the convolution is given by $s(i)={\displaystyle \sum _{k=0}^{i}x}(k)h(i-k)=x\left(i\right)*h(i)i=\mathrm{0,1,}\cdots \mathrm{,}N+M-\mathrm{1,}$ (1) where $x(k)$ denotes the original input signal, $h(i)$ is the response function, $s(i)$ is the response signal. The same signal length as that of the input signal was obtained by discarding the boundary values.

In real measurements, noise is an inevitable influencing factor in gamma-ray detection systems. Therefore, the scintillation detector response spectrum can be expressed as a matrix equation, as shown in Eq. (2). $s=Hx+n\mathrm{,}$ (2) where s is the scintillation detector response-spectrum vector, x is the gamma-ray energy vector in the detector, n is the noise vector, and H is the Gaussian response matrix.

The Gaussian response matrix determines the extent of energy broadening, which can be written as $H={\left[\begin{array}{ccc}{G}_{11}& \cdots & G_{1m}\\ ⋮& \ddots & ⋮\\ G_{m1}& \cdots & G_{mm}\end{array}\right]}_{m\times m}\mathrm{,}$ (3) where m denotes the length of the response spectrum vector. The element Gij in the matrix H is expressed as a Gaussian function: $G_{ij}=\frac{1}{\sqrt{2\pi }\sigma }{\text{e}}^{-\frac{{\left(E_i-E_j\right)}^2}{2{\sigma }^2}}$ (4) It means that Eq.(4) represents the contribution of a gamma ray with energy Ej to the channel with energy Ei, and σ in Eq.(4) is represented by Eq.(5), $\sigma =\frac{1}{2\sqrt{2\ln 2}}FWHM=\frac{1}{2\sqrt{2\ln 2}}\left(a+b\sqrt{c{E}^2+E}\right)\mathrm{,}$ (5) where FWHM is full width at half maximum, it's a function of a, b, and c; E is the energy of the corresponding channel.

Therefore, the column vector in the matrix H can be regarded as a normalized Gaussian distribution centered on the corresponding diagonal element. Figure 1 shows an example of a Gaussian response matrix in three-dimensional coordinates with a=0.032, b=0.045, c=0.83.

Fig. 1Full-screen View

(Color online) The example of the Gaussian response matrix

Locally constrained regularization

The key to calibrating the spectral energy resolution is to solve the Gaussian response matrix that matches the response of the gamma-ray measurement system. The solution to the Gaussian response matrix is a typical discrete ill-posed problem, which implies that employing unconstrained least squares cannot yield satisfactory results. Hence, the classical method for solving this type of problem is regularization, which can find a reasonable approximate solution under given constraints. Equation (6) presents an example of a regularization objective function with a linear approximation constraint using weighted least squares (WLS). $F=\text{argmin}{\displaystyle \sum _{i=1}^m{\omega }_i}{\left(u_i-s_i\right)}^2+\lambda {\displaystyle \sum _{i=2}^{m-1}{\left({\text{Δ}}^2{s}_i\right)}^2}\mathrm{,}$ (6) where si and ui are the ith elements of the calculated and measured vectors, respectively, λ is a regularization parameter that balances both parts of Eq. (6), ωi is the weight of the ith channel, Δ² is the second-order difference operator, which can be expressed as a matrix D shown in Eq. (7). Similarly, quadratic or cubic approximations can be obtained by constructing the corresponding difference matrix [35]. $D={\left[\begin{array}{cccccccc}1& -2& 1& 0& \cdots & 0& 0& 0\\ 0& 1& -2& 1& \cdots & 0& 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ 0& 0& 0& \cdots & 1& -2& 1& 0\\ 0& 0& 0& \cdots & 0& 1& -2& 1\end{array}\right]}_{m-2\times m}$ (7) Therefore, the objective function can be expressed as a matrix equation, as shown in Eq. (8), $F(s)={\left(u-s\right)}^{\text{T}}{W}_1\left(u-s\right)+\lambda s^{\text{T}}{D}^{\text{T}}Ds\mathrm{,}$ (8) where u is the measured vector, s is the calculated vector, W₁ is the diagonal weight matrix.

By solving the matrix equation, the optimal solution vector can be obtained when the error between the calculated and measured vectors is minimized.

For spectral resolution calibration, it is difficult to determine an accurate Gaussian response matrix using only the WLS method, particularly in the case of a poor SNR of the gamma spectrum. Fortunately, incident gamma-ray energies are useful prior knowledge in gamma spectral analysis. The energies of neutron-induced or natural gamma rays are related to the inherent properties of nuclei. Based on this premise, the characteristic peaks of gamma rays in the energy spectrum can be used to set additional constraints to determine more accurate energy-resolution parameters.

In the proposed locally constrained regularization method, the objective function can provide additional weights to the peak regions so that the determined energy resolution is more suitable for broadening the characteristic peaks. The new objective function is expressed by Eq. (9): $F=\text{argmin}{\displaystyle \sum _{i=1}^m{\omega }_{1i}}{\left(u_i-s_i\right)}^2+\lambda {\displaystyle \sum _{j=1}^n{\displaystyle \sum _{i=p_{j1}}^{p_{jk}}{\left[{\omega }_{2ji}\left(u_{ji}-s_{ji}\right)\right]}^2}}\mathrm{,}$ (9) where λ is the regularization parameter, n is the number of selected peak regions, pjk is the number of channels in the jth peak region, $w_{1i}$ is the weight of the ith channel, $w_{2ji}$ is the additional weight of the ith channel in the jth peak region.

The weight of each channel can be designed as a diagonal matrix W₁, as expressed in Eq. (10), which is determined by the relative count of the corresponding channel [36]. $W_1={\left[\begin{array}{ccccc}{u}_1^{-1}& & & & \\ & u_2^{-1}& & & \\ & & \ddots & & \\ & & & u_{m-1}^{-1}& \\ & & & & u_m^{-1}\end{array}\right]}_{m\times m}$ (10) The additional weights can also be expressed as matrix W₂ shown in Eq. (11). The number of columns is equal to the spectral length m, and the number of rows is determined by the number of selected peak regions n. $W_2={\left[\begin{array}{cccccc}0& \cdots & \stackrel{FWHM}{\overbrace{w_{1p_1}\cdots Peak\cdots w_{1p_k}}}& \cdots & \cdots & 0\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ 0& \cdots & \cdots & \underset{FWHM}{\underbrace{w_{n{p}_1}\cdots Peak\cdots w_{n{p}_k}}}& \cdots & 0\end{array}\right]}_{n\times m}$ (11) Only the channels within the FWHM of the selected peak region were assigned additional weights. A nonlinear weighting scheme based on a sigmoid function was designed as shown in Eq. (12): $w_{ip}=\text{Sigmoid}\left(d_{ip}\right)=\frac{1}{1+{\text{e}}^{2\left(d_{ip}-FWHM/2\right)/\left(FWHM/6\right)}}\mathrm{,}$ (12) where wip is the weight of the pth channel in the ith peak region and dip is the distance between the pth channel and the peak position of the ith peak region.

In the proposed weighting scheme, the weight of the channel is determined by the distance from the peak position and gradually decreases with an increase in the distance. The channel of the peak position had the maximum weight, and the boundaries of the FWHM had the minimum weight. The weight distributions in the peak regions are shown in Fig. 2(a). Channels within the FWHM of the spectral peak were assigned additional weights using this weighting scheme, as shown in Fig. 2(b)). Similarly, the new objective function of the locally constrained regularization can be expressed as a matrix equation: $F(s)={\left(u-s\right)}^{\text{T}}{W}_{\text{1}}\left(u-s\right)+\lambda {\left(u-s\right)}^{\text{T}}{W}_2^{\text{T}}{W}_2\left(u-s\right)\mathrm{.}$ (13) Finally, the objective function is converted into a function of a Gaussian response matrix. $F(H)={\left(u-Hx\right)}^{\text{T}}{W}_1\left(u-Hx\right)+\lambda {\left(u-Hx\right)}^{\text{T}}{W}_2^{\text{T}}{W}_2\left(u-Hx\right)$ (14)

Fig. 2Full-screen View

(Color online) (a)The weight distribution of the peak region; (b)The additional weight region of the peak region

Energy resolution calibration

A strong correlation exists between the elements in the Gaussian response matrix, which is determined by a broadening mechanism. Thus, the objective function can be regarded as a function of resolution calibration parameters a, b, and c. Calibrating the energy resolution can be considered a high-dimensional constrained optimization problem, which can be expressed as Eq.(15), where N is the number of measured spectra: $\text{min}{\displaystyle \sum _{i=1}^N{f}_i}\left(a\mathrm{,}b\mathrm{,}c\right)s\mathrm{.}t\mathrm{.}a\mathrm{,}b\mathrm{,}c\ge 0$ (15) The swarm intelligence optimization algorithm is an iterative method that updates feasible solutions based on their fitness. It has the advantage that it does not depend on the initial guess and does not require a derivative of the objective function. Therefore, we adopt the particle swarm optimization (PSO) to achieve the global optimal solution of the objective function.

PSO is a classical optimization algorithm inspired by bird predation behaviors. PSO initializes a group of particles as potential solutions, and the fitness of the solution is determined by the objective function. In the kth iteration, the ith particle updates its velocity $V_i^k$ and position $X_i^k$ by tracking the individual optimum position $P_i^k$ and global optimum position $P_g^k$ to search for the optimal solution, as shown in Eq. (16) and Eq. (17) [37]. $V_i^{k+1}=\omega V_i^k+c_1{r}_1\left(P_i^k-X_i^k\right)+c_2{r}_2\left(P_g^k-X_i^k\right)\mathrm{,}$ (16) $X_i^{k+1}=X_i^k+V_i^{k+1}\mathrm{,}$ (17) where ω is the inertia weight, c₁ and c₂ are the acceleration factors, r₁ and r₂ are random numbers distributed between [0, 1].

Another key parameter in the objective function is the incident gamma ray energy vector x. As mentioned above, Monte Carlo simulation is a useful method for obtaining the response of a gamma-ray measurement system, which can solve neutron, photon, electron, or coupled neutron/photon/electron transport problems in a three-dimensional complex geometry. The gamma-ray energy distribution in the scintillation detector can be effectively obtained by establishing a simulation model consistent with the actual measurement system. Then, utilizing the simulated gamma ray energy distribution, the error between the measured and calibrated spectra was computed based on the established objective function under the preset constraints. Normalization is required in the calculation to ensure that the measured and calibrated spectra are of the same order of magnitude. This error can be used as a fitness evaluation criterion to update the particle swarm using the PSO algorithm. A flowchart of the energy-resolution calibration using the locally constrained regularization method is shown in Fig. 3.

Fig. 3Full-screen View

(Color online) The flowchart of energy resolution calibration by locally constrained regularization method

Experiments on simulated spectra

To verify the performance of the proposed locally constrained regularization method, a Monte Carlo simulation model was constructed based on an actual NaI detector structure, as shown in Fig. 4. The diameter and length of the NaI scintillation crystals were 3.18 cm, and 10.16 cm. The NaI crystal was coated with a MgO reflecting layer with a thickness of 0.05 cm. The outermost layer of the detector was an aluminum shell with a thickness of 0.1cm. One end of the detector was a light guide made of SiO₂, which is 0.2 cm thick. A gamma source that can emit gamma rays with a specified energy distribution was set in front of the detector. As neutrons are transported in a medium, they gradually thermalize because they collide numerous times with the nuclei of the medium. Thermal neutrons are then captured by the target nucleus to release gamma photons. The energies of the neutron-induced gamma rays determined by the type of the elements can be regarded as the elemental “fingerprint.” In the simulation experiments, the thermal neutron-captured gamma rays of four common elements, H, Si, Al, and Ca, were used as the incident gamma-ray energy distributions [38], as shown in Fig. 5. The gamma source of the simulation model emitted gamma rays that followed the above four distributions, and the simulated gamma rays with energies ranging from 0.088 MeV to 10 MeV were recorded in 1024 channels as the experimental test set. Random noise was added to each simulated spectrum to illustrate the applicability of the proposed method.

Fig. 4Full-screen View

(Color online) The model for the simulated spectra

Fig. 5Full-screen View

(Color online) The capture gamma-ray energy distributions of H, Si, Al, and Ca elements

For the regularization method, the regularization parameter is the key factor affecting performance. The main function of the regularization parameter is to balance the importance of the sum of the weighted squares of the errors and constraints. Therefore, mathematical experiments were conducted to determine the most appropriate regularization parameters. In the experiments, the weight of each channel in the matrix W₁ was determined by the corresponding intensity expressed in Eq. (10), and part of the characteristic gamma-ray energies of the four elemental energy distributions whose intensities are larger than 10% are selected to construct the regularization weight matrix W₂. Regularization parameters 10, 20, 30, 50, 70, and 100 were used to establish the objective function, and the PSO algorithm was adopted to search for the optimal resolution calibration parameters with different regularization parameters. The average RMSEs of the spectra of the H, Si, Al, and Ca elements between the calibrated and simulated spectra with different regularization parameters are shown in Fig. 6. It can be found that the calibrated spectra have the minimum RMSE when the regularization parameter is about 30. These conclusions can guide the application of the proposed method.

Fig. 6Full-screen View

(Color online) The average RMSEs between the original spectra and the calibrated spectra with different regularization parameters

A comparison between the simulated and calibrated spectra is shown in Fig. 7. Even for spectra with noise, the resolution calibration parameters determined by the proposed locally constrained regularization method still have a high consistency with the simulated spectra. The correlation coefficients between all the simulated spectra and the corresponding calibrated spectra were greater than 0.99, which indicates the applicability of the proposed method.

Fig. 7Full-screen View

(Color online) The resolution extraction results of the locally constrained regularization method for different elemental spectra

The performances of the locally constrained regularization and traditional weighted least squares were also compared to illustrate their applicability. Resolution calibration results for the Al spectra obtained using the two methods are shown in Fig. 8. This is shown in Fig. 8(a), the resolution calibration parameters determined using the locally constrained regularization method cause each characteristic peak of the calibration spectrum to match the corresponding characteristic peak of the simulation spectrum. In contrast, although the spectrum calibrated by the weighted least squares is well consistent with the simulated spectrum in the energy range of 3.5 MeV–6 MeV, the resolution error of the characteristic peaks in the energy range of 0.5 MeV–3.5 MeV and 6 MeV–8 MeV is relatively large because of the lack of sufficient constraints, as shown in Fig. 8(b).

Fig. 8Full-screen View

(Color online) The resolution calibration results of the Al spectrum using the locally constrained regularization method and weight least squares

The optimal resolution calibration parameters of the two methods were determined by the PSO algorithm using 1000 iterations. The iterative curves of the locally constrained regularization method and the weighted least-squares method are shown in Fig. 9. We found that the RMSEs of the two methods gradually decreased with the iterations, indicating that the fitness of the objective function was enhanced. Although the iterative curve of the weighted least squares method has a faster convergence rate during the iteration, the final RMSE is higher than that of the proposed method, which means that the optimization is trapped in a locally optimal solution.

Fig. 9Full-screen View

(Color online) Comparison of the iterative curves of the locally constrained regularization method and weight least squares

Experiments on measured spectra

Pulsed neutron logging is a powerful method for detecting and analyze reservoir elemental compositions and has been widely used in oilfield exploration and development. The logging instrument was a typical neutron-induced gamma-ray measurement system. It usually uses a D-T neutron source to excite the nuclei of formation elements to release gamma rays and uses scintillation detectors to obtain the corresponding gamma energy spectra. Therefore, test pit experiments for pulsed neutron logging were performed to verify the performance of the proposed locally constrained regularization method for spectral resolution calibration. The experimental pulsed neutron logging instrument was set in different artificial pits with known compositions and geometric structures. As shown in Fig. 10 (a), several test pits are located at the experimental site. The configurations of the pure rock and simple substance test pits are shown in Fig. 10 (b). The depth of these test pits was 6050 mm and they were divided into multiple intervals. The intervals are filled with different materials, including aluminum (Al), magnesium (Mg), iron (Fe), titanium (Ti), graphite (C), limestone (LS), and granite (GS). The instrument was inserted into a borehole with a diameter of 200 mm, which was filled with fresh water, to measure neutron-induced gamma-rays. In the experiments, four test pits were selected as the measurement objects, three of which were pure rock or simple substance test pits (LS, Al, and C), and the fourth was a water tank. As described in the Methodology section, a Monte Carlo simulation model was constructed based on the actual experimental instrument structure and test pit configurations to simulate the energy deposition distribution in the scintillation crystal, as shown in Fig. 11. The diameter and length of the instrument were SI41mmand 2358 mm. The D-T neutron source consisted of an accelerating system and a neutron-producing target. The instrument has three scintillation detectors with different source-detector spacings (29 cm, 59 cm and 91 cm), and there is a tungsten shield between the source and the nearby detector. The far and near detectors use LaBr₃ crystals with crystal sizes of 27.5 mm×152 mm and 27.5 mm×38 mm (diameter×length), respectively. The extra detector adopts a NaI crystal of size 27 mm×152 mm. The extra detector has the maximum spacing; therefore, it is mainly used to collect natural gamma-rays rather than neutron-induced gamma-rays. In the experiments, considering the influence of radioactive statistical fluctuations on gamma-ray spectral measurements, the energy spectrum collected by the near detector was selected to test the proposed method.

Fig. 10Full-screen View

(Color online) The test pit experimental site (a) and the configurations of pure rock or substance test pits (b)

Fig. 11Full-screen View

(Color online) The test pit experimental Monte Carlo model

Based on the mechanism of nuclear reactions between fast neutrons and nuclei, two types of interactions can produce secondary gamma-rays: inelastic collisions and capture reactions. The experimental instrument utilized a special measurement sequence to measure and capture inelastic secondary gamma rays. As shown in Fig. 12, the D-T neutron source of the instrument periodically emits neutrons with energies of 14 MeV. The complete measurement period of 20ms consists of 200 subperiods of 100 μs, where the neutron source emitted neutrons during the first 40 μs of each subperiod, whereas the neutron source did not work for the next 60 μs. In the measurement, the scintillation detector recorded the gamma ray energy spectra during the G2 and G3 time gates and then accumulated the spectra of all sub-periods to obtain the final energy spectra. In general, inelastic collision and capture reactions occur during the neutron burst time (G2), whereas only the capture reaction can occur during the G3 time gate, when the neutron source stops working. Therefore, the spectrum obtained in time gate G2 is called the total spectrum, which contains inelastic secondary gamma-rays and captures secondary gamma-rays, whereas the spectrum obtained in time gate G3 is called the capture spectrum, which only contains captured secondary gamma-rays. The proposed locally constrained regularization method was used to calibrate the energy resolution using two types of spectra. All the total and capture spectra measured in the four test pits were used to establish the objective function. After optimization, the optimal resolution calibration parameters were a= 0.0479, b= 0.0146, and c= 8.3718. Obviously, the determined resolution calibration parameters apply to all the spectra measured by the instrument.

Fig. 12Full-screen View

(Color online) The measurement sequence of the experimental instrument

The calibration results for the entire spectrum are shown in Fig. 13. It can be seen that each total spectrum has an apparent characteristic peak with an energy of 2.23 MeV, which is produced by the capture reaction of hydrogen in the borehole fluid. In addition, some characteristic peaks of other elements were observed in the total spectra, such as the inelastic characteristic peaks of oxygen (6.13 MeV) and carbon (4.43 MeV) in the water tank and graphite test pit. All the characteristic peaks, including the corresponding escape peaks, can be used to set additional constraints as useful prior knowledge.

Fig. 13Full-screen View

(Color online) The comparison of the experimental spectra and the calibrated spectra of different test pits (total spectrum)

The selected characteristic peaks of each total spectrum, RMSE, and correlation coefficient (R-Square) between the normalized experimental spectrum and normalized calibrated spectrum are shown in Table 1. In the experiment, distinct peaks were selected for each spectrum to set the constraint conditions. We found that there were good correlations between the experimental and calibrated spectra from the comparison results. All correlation coefficients were greater than 0.99, and the average RMSE was only 4.22e-4 after resolution calibration using the locally constrained regularization method.

Comparisons of the normalized experimental capture spectra with the normalized calibrated capture spectra are shown in Fig. 14, and the capture spectral calibration results for the four test pits are listed in Table 2. The characteristic peak of hydrogen was also observed in all four measured spectra. In addition, the characteristic peaks of calcium (6.42 MeV) and aluminum (3.03 MeV and 4.69 MeV) were selected to set constraints for the spectra of the corresponding test pits. From the comparison results, it can be observed that the correlations between the calibrated and measured spectra are sufficient. All correlation coefficients were greater than 0.99, and the average RMSE was 4.03e-4. The calibration results of the capture spectra also proved the performance of the locally constrained regularization method in spectral-resolution calibration.

Fig. 14Full-screen View

(Color online) The comparison of the experimental spectra and the calibrated spectra of different test pits (capture spectrum)