Theoretical model of Bonner Spheres Spectrometer

The Bonner Spheres Spectrometer (BSS) is considered a typical neutron spectrum measurement instrument that offers a wide energy range, high sensitivity, isotropic response, consistent neutron flux, and ambient dose equivalent. The neutron is moderated by a polyethylene shell in the BSS. The center uses a thermal neutron detector, which can be one of the following three types: ³He, ¹⁰B and ⁶Li. The ³He neutron proportional detector with the highest thermal neutron cross section is selected as the detector. For unfolding the neutron spectrum, the response function for each sphere of BSS is required, and is commonly calculated by the Monte Carlo method. The measurement principle can be expressed as the first kind of Fredholm integral equation ^[21-23]: $\begin{array}{c}{N}_i+{\epsilon }_i=\underset{E_{\text{min}}}{\overset{E_{\text{max}}}{{\displaystyle \int }}}{R}_i\left(E_j\right)\Phi \left(E_j\right)\text{d}E,\\ i=1,2,\cdots ,n; j=1,2,\cdots ,m\end{array}$ (1) where, $N_i$ is the count rate of the ith Bonner Sphere; ${\epsilon }_i$ is the measurement error of ith Bonner Sphere. $E_{\text{max}}$ and $E_{\text{min}}$ are the highest and lowest energy of the neutron; $R_i\left(E_j\right)$ is the Bonner Sphere response function corresponding to the jth energy interval; $\Phi \left(E_j\right)$ is the fluence rate of neutrons in the jth energy interval; n is the number of Bonner Spheres; m is the number of energy intervals.

Since the $R_i\left(E_j\right)$ and $N_i$ are discretely distributed, Eq. (1) can be expressed in the discrete form of Eq. (2): $N_i+{\epsilon }_i={\displaystyle \sum _j{R}_{ij}{\phi }_j} i=1,2,\cdots ,n; j=1,2,\cdots ,m$ (2)

Eq. (2) can also be expressed in matrix form as follows: $N+\epsilon =R\phi ,$ (3) where, $N={(N_1,N_2,\dots ,N_n)}^T$, $\phi ={\left({\phi }_1,{\phi }_2,\dots ,{\phi }_m\right)}^T$, R is the response matrix of n×m.

Usually, N and R are known quantities, and the number of Bonner spheres n is often less than the number of the energy interval m. The Eq. (3) belongs to the underdetermined equations. In this work, the PSO-MLEM algorithm was used to solve the underdetermined problem in neutron spectrum unfolding.

Principle of PSO-MLEM algorithm

This paper combines the PSO algorithm with the MLEM algorithm, and the dynamic acceleration factor is used to balance the global and local search ability, which improves the convergence speed and calculation accuracy of the PSO-MLEM algorithm. The MLEM algorithm is used to limit the search direction of particles during each iteration of the PSO-MLEM algorithm. Figure 1(a) depicts the particle optimization process in the PSO-MLEM algorithm. The MLEM algorithm guides the flight direction of the particle swarm to ensure that the particles always fly to the solution space of the MLEM algorithm. The theoretical position calculated by the MLEM algorithm replaces the position determined by the particle inertia, which effectively reduces the possibility of search stagnation and missing the optimal solution. At the same time, the displacement of the particles is determined by both the PSO and MLEM algorithms during each iteration of the PSO-MLEM algorithm. As shown in the dynamic acceleration factor transformation Eq. (8), the PSO-MLEM algorithm adopts a smaller $c_1$ and a larger $c_2$ in the early stage, so that the particles can quickly narrow the range of the solution space with the help the MLEM algorithm. With the number of iterations increases, MLEM algorithm is easy to suffer the pitfalls of local optima. Therefore, the PSO-MLEM algorithm gradually increases $c_1$ and decreases $c_2$ in the iterative process, which strengthens the local search ability of the algorithm and helps MLEM algorithm to jump out of the pitfalls of local optima.

Fig. 1Full-screen View

Graphical depiction of PSO-MLEM algorithm. (a) Schematic diagram of PSO-MLEM particle flight; (b) Flow chart of PSO-MLEM algorithm.

The following 6 steps further explain the general procedure of the PSO-MLEM method for neutron spectrum unfolding:

1. Initialization of the particle swarm

In the PSO-MLEM algorithm, the initial step is to define the number of particles, the number of particles is set to 600, and each particle is defined as a solution vector that must be optimized $\phi ={\left[{\phi }_1,{\phi }_2,\dots ,{\phi }_m\right]}^{ }$. The solution vector $\phi$ is the neutron spectrum to be optimized in the process of neutron spectrum unfolding, while the initial position and velocity of the particle are initialized by random numbers between 0 and 1.

2. Calculation of the cost function

The following cost function is specified to limit the search range of the PSO-MLEM algorithm: $\text{Cost}= \frac{{\left|\left|N_{\text{cal}}-N_{\text{obs}}\right|\right|}^2}{{\left|\left|N_{\text{obs}}\right|\right|}^2}+\frac{{\left|\left|\phi -{\phi }_{\text{ref}}\right|\right|}^2}{{\left|\left|{\phi }_{\text{ref}}\right|\right|}^2} ,$ (4) where, $N_{\text{obs}}$ is the count rate of Bonner spheres; $\phi$ is the solution of each particle in each iteration; ${\phi }_{\text{ref}}$ is the reference spectrum that contains the a priori information; $N_{\text{cal}}$ is calculated by using Eq. (5) as follows: $N_{\text{cal}}=R\phi .$ (5)

Through the reference spectrum ${\phi }_{\text{ref}}$ and the count rate $N_{\text{obs}}$, the above cost function limits the search direction of the particle. The first term of the Cost function is the distance between $N_{\text{cal}}$ and $N_{\text{obs}}$, and the second term is the distance between the updated spectrum $\phi$ and the reference spectrum ${\phi }_{\text{ref}}$. When the minimum of Cost function is obtained, the $\phi$ is the result of the neutron spectrum unfolding.

3. Update the local optima solution of a single particle (pbest)

The solution of each particle is calculated in the current iteration process. Then the cost value of the current solution is compared with the cost value of the historical local optima solution(pbest). If the current cost value is smaller than the historical cost value, the local optima solution is updated to the solution of the current particle.

4. Update the global optima solution of all particles (gbest)

The cost value of the local optima solution (pbest) is compared with the cost value of the global optima solution (gbest). If the cost value of the local optima solution is smaller than the cost value of the global optima solution, the global optima solution is updated to the local optima solution of the particle.

5. Update the velocity and solution vector of the particles

According to Eq. (6) and (7), the velocity $v_{ij}$ and solution ${\phi }_{ij}$ of the jth dimension of the ith particle are updated: $v_{ij}\left(t+1\right)=r_1\cdot \left(pbes{t}_i\left(t\right)-{\phi }_{ij}\left(t\right)\right)+r_2\cdot \left(gbest\left(t\right)-{\phi }_{ij}\left(t\right)\right),$ (6) ${\phi }_{ij}\left(t+1\right)=c1\cdot v_{ij}\left(t+1\right)+c2\cdot {\phi }_{ij}^{\text{MLEM}}\left(t+1\right),$ (7) where, $v_i$ is the velocity vector of the ith particle, ${\phi }_i$ is the solution vector of the ith particle, the $pbes{t}_i$ is the historical local optima solution of the ith particle, and $gbest$ is the global optima solution of all particles, $r$ is a random number between [0, 1], $t$ is the number of iterations, $c$ is the acceleration factor between [0, 1], as shown in Eq. (8): $\{\begin{array}{l}{c}_1=\frac{t}{T_{\text{max}}}\hfill \\ c_2=1-\left(\frac{t}{T_{\text{max}}}\right)\hfill \end{array}$ (8) where, $t$ is the current number of iterations, and $T_{\text{max}}$ is the maximum number of iterations. ${\phi }_{ij}^{\text{MLEM}}$ is calculated according to the MLEM algorithm using Eq. (9): ${\phi }_{ij}^{\text{MLEM}}\left(t+1\right)={\phi }_{ij}\left(t\right)·\frac{{\displaystyle {\sum }_{i=1}^n{R}_{ij}}\frac{N_i}{{\displaystyle {\sum }_{j=1}^m{R}_{ij}}{\phi }_j\left(t\right)}}{{\displaystyle {\sum }_{i=1}^n{R}_{ij}}}.$ (9)

6. Judging convergence

The search is terminated when the cost value of the global optima solution achieves the predefined iterative convergence condition or when the current number of iterations reaches the maximum number of iterations. If the above conditions are not met, steps 2)-6) are repeated.

The algorithm flow chart summarizes the preceding 6 steps as follows Fig.1(b):

Verification by Monte Carlo simulation

In order to verify the accuracy and convergence of the PSO-MLEM algorithm, the Monte Carlo method was used to simulate the response function of the Bonner Spheres Spectrometer. And the count rates of Bonner spheres were simulated while four reference spectra from the IAEA Technical Report Series No. 403^[24] were used as the input spectra of the Monte Carlo method. The MXD_FC33 program in the UMG unfolding software package V3.3, which was based on the method of maximum entropy deconvolution (MAXED), was used to unfold the neutron spectrum^[25]. The reference spectrum, adding a random error, was set as the prior information of the PSO-MLEM algorithm and MAXED. Combined with the response function and the count rates of Bonner spheres obtained by simulation, the neutron spectrum was unfolded using PSO-MLEM, PSO, MLEM and MAXED, respectively.

2.3.1

Bonner Spheres Spectrometer response function simulation

The Bonner spheres model refers to the parameters used by the Physikalisch Technische Bundesanstalt (PTB). The Bonner spheres spectrometer consists of 15 polyethylene moderated spheres of various diameters, with diameters of 0, 2.5, 3, 3.5, 4, 4.5, 5, 6, 7, 8, 9.5, 10, 12, 15, 18 inches. The response function has an energy range of 1×10^–9~398 MeV, which is divided into 60 energy groups ^[26]. Each sphere has a thermal neutron detector in the center. The center detector chooses the SP90 ³He proportional counter manufactured by Centronics of the UK. The Bonner sphere model is shown in Fig.2(a).

Fig. 2Full-screen View

(Color online) Bonner sphere model and Responses functions. (a) Bonner sphere model and ³He proportional counter; (b) Responses functions of Bonner spheres spectrometer simulated by Monte Carlo method

The response functions of Bonner spheres spectrometer refer to the counts caused by unit incident neutrons emitted from parallel source terms. The number of reactions between the neutrons entering the detector and ³He must be recorded when the Monte Carlo method was used to simulate the response function. The response functions of Bonner spheres spectrometer to a single-energy parallel neutron beam with a diameter of d and incident energy of E is defined as follows^[27]: $R_{\text{d}}\left(E_{\text{n}}\right)=\frac{M_{\text{d}}}{\Phi \left(E_{\text{n}}\right)},$ (10) where, $R_{\text{d}}\left(E\right)$ is the energy response (cm2) when the incident neutron energy is E, $M_{\text{d}}$ represents the count rate of the Bonner sphere with diameter d, and $\Phi \left(E\right)$ is the neutron fluence (cm^–2) when the incident neutron energy is E. Figure 2(b) depicts the response functions of 15 Bonner spheres as simulated by the Monte Carlo method.

As Fig. 2(b), the Bonner sphere with a small diameter has a high response in the low energy range, and the Bonner sphere with a large diameter has a high response in the high energy range. In addition, as the diameter of the Bonner sphere increases, the maximum response gets narrower and moves towards higher energies.

2.3.2

Validation of the unfolded spectrum

As shown in Fig.3(a1), the Input-Spec1 is from the CERN-CEC high energy reference field facility, 40 cm thick Fe shield in the IAEA Technical Report Series No. 403, the Input-Spec2 is from Booster-Synchrotron of the Stanford Synchrotron Radiation Laboratory (SSRL) SPEAR in the IAEA Technical Report Series No. 403, the Input-Spec3 is from CERN Pu-Be calibration spectra in the IAEA Technical Report Series No. 403 and the Input-Spec4 is from Final Focus Test Beam Facility (FFTB) in the IAEA Technical Report Series No. 403. In the Fig.3(a1), the Ref-Spec1, Ref-Spec2, Ref-Spec3 and Ref-Spec4, which are the spectra from Input-Spec1 to Input-Spec2 adding the ±10% random error, are the prior information for the PSO-MLEM algorithm and MAXED. Figure 3(a2) shows the count rates of Bonner spheres spectrometer obtained when the Input-Spec1 to Input-Spec4 were used for the input spectrum of simulation. With the response function and the count rates of Bonner spheres spectrometer obtained by simulation, the PSO-MLEM, PSO, MLEM and MAXED were used to unfold the neutron spectrum.

Fig. 3Full-screen View

Unfolded spectra of Monte Carlo simulation data. (a) Input spectrum and count rates of Bonner spheres (a1: The solid line is the input spectrum of Monte Carlo method, and the dashed line is the prior information with a random error; a2: The count rates of 15 Bonner spheres); (b) Unfolded spectrum of the Spec1 (b1: Comparison of PSO-MLEM, PSO, MLEM and MAXED unfolded spectra with the input spectrum; b2: The ratio of $N_{PSO\_MLEM}$, $N_{PSO}$, $N_{MLEM}$, $N_{MAXED}$ and $N_{def}$ to $N_{obs}$); (c) Unfolded spectrum of the Spec2 (c1: Comparison of PSO-MLEM, PSO, MLEM and MAXED unfolded spectra with the input spectrum; c2: The ratio of $N_{PSO\_MLEM}$, $N_{PSO}$, $N_{MLEM}$, $N_{MAXED}$ and $N_{def}$ to $N_{obs}$); (d) Unfolded spectrum of the Spec3 (d1: Comparison of PSO-MLEM, PSO, MLEM and MAXED unfolded spectra with the input spectrum; d2: The ratio of $N_{PSO\_MLEM}$, $N_{PSO}$, $N_{MLEM}$, $N_{MAXED}$ and $N_{def}$ to $N_{obs}$); (e) Unfolded spectrum of the Spec4 (e1: Comparison of PSO-MLEM, PSO, MLEM and MAXED unfolded spectra with the input spectrum; e2: The ratio of $N_{PSO\_MLEM}$, $N_{PSO}$, $N_{MLEM}$, $N_{MAXED}$ and $N_{def}$ to $N_{obs}$)

Figure 3(b1) - (e1) depict the results of the unfolded spectra by the PSO-MLEM, PSO, MLEM and MAXED with the prior information (Ref-Spec1 to Ref-Spec4). Figure 3(b2) - (e2) is the ratio relation between $N_{\text{def}}/N_{\text{obs}}$, $N_{\text{PSO}\_\text{MLEM}}/N_{\text{obs}}$, $N_{\text{PSO}}/N_{\text{obs}}$, $N_{\text{MLEM}}/N_{\text{obs}}$ and $N_{\text{MAXED}}/N_{\text{obs}}$. $N_{\text{obs}}$ are the count rates of 15 Bonner spheres, $N_{\text{def}}\left(N_{\text{def}}=R{\phi }_{\text{input}}\right)$ is the result of the response function times the input spectrum. Besides, $N_{\text{PSO}\_\text{MLEM}}\left(N_{\text{PSO}\_\text{MLEM}}=R{\phi }_{\text{PSO}\_\text{MLEM}}\right)$ is the result of the response function times the unfolded spectrum of the PSO-MLEM algorithm, $N_{\text{PSO}}\left(N_{\text{PSO}}=R{\phi }_{\text{PSO}}\right)$ is the result of the response function times the unfolded spectrum of the PSO algorithm, $N_{\text{MLEM}}\left(N_{\text{MLEM}}=R{\phi }_{\text{MLEM}}\right)$ is the result of the response function times the unfolded spectrum of the MLEM algorithm, and $N_{\text{MAXED}}\left(N_{\text{MAXED}}=R{\phi }_{\text{MAXED}}\right)$ is the result of the response function times the unfolded spectrum of the MAXED. The closer the ratio is to 1, the smaller the difference between the unfolded spectrum and the input spectrum is proved.

The Pearson correlation coefficient was utilized for the evaluation of the correlation between the unfolded neutron spectrum ${\phi }^{\text{cal}}$ and the input spectrum ${\phi }^{\text{input}}$, and its value ranges from +1 and –1. The r can calculate as following equation: $r_{{\phi }^{\text{cal}}{\phi }^{\text{input}}}=\frac{m\sum {\phi }_i^{\text{cal}}{\phi }_i^{\text{input}}-\sum {\phi }_i^{\text{cal}}\sum {\phi }_i^{\text{input}}}{\sqrt{m\sum {\left({\phi }_i^{\text{cal}}\right)}^2-{\left(\sum {\phi }_i^{\text{cal}}\right)}^2}\sqrt{m\sum {\left({\phi }_i^{\text{input}}\right)}^2-{\left(\sum {\phi }_i^{\text{input}}\right)}^2}}$ (11)

The relative mean error (RME) was also used to calculate the difference of the unfolded spectrum from the input spectrum. The RME can calculate as following equation: $RM{E}_{{\phi }^{\text{cal}}{\phi }^{\text{input}}}=\frac{1}{m}{\displaystyle \sum _{i=1}^m\left|\frac{{\phi }_i^{\text{cal}}-{\phi }_i^{\text{input}}}{{\phi }_i^{\text{input}}}\right|},$ (12) where, $m$ is the sample size, represents the number of energy intervals, and ${\phi }_i^{\text{cal}}$ and ${\phi }_i^{\text{input}}$ respectively represent the normalized value of the ith channel corresponding to the unfolded spectra and the input spectrum. For Eq. (11), when $r$=0, the two variables are uncorrelated, when $r$ >0, there is a positive correlation, when $r$ <0, there is a negative correlation, and the degree of correlation is as follows: very strong correlation (0.8–1.0), strong correlation (0.6–0.8), moderate correlation (0.4–0.6), weak correlation (0.2–0.4), very weak correlation or no correlation (0.0–0.2). For equation Eq. (12), under ideal conditions (${\phi }^{\text{cal}}={\phi }^{\text{input}}$), $RME$ is 0.

As can be seen from Table 1 that even if the prior information with a random error from the input spectra of the Monte Carlo method is used in the PSO-MLEM algorithm, the correlation coefficients of the PSO-MLEM unfolded spectrum of the four neutron sources are all above the value of 0.99 and the relative mean error $RME$ are all less than 0.07%, indicating that the unfolded spectrum by the PSO-MLEM has a strong correlation with the input spectrum of the Monte Carlo method. The correlation coefficient of the PSO and MLEM algorithm is significantly lower than the PSO-MLEM algorithm. The Input-Spec3 fluctuate more strongly, the results of PSO and MLEM algorithms are lower accuracy. The correlation coefficient of Spec3 is only 0.74 in PSO and MLEM algorithms, which is lower than the PSO-MLEM algorithm (0.99), and the relative mean error is also more than the PSO-MLEM algorithm. The result shows that the PSO-MLEM algorithm has better robustness when the neutron spectrum is complex. At the same time, compared with the unfolded spectrum by MAXED, the correlation coefficient of PSO-MLEM is increased by 1.9%, and the relative mean error is decreased by 67.4%.

Table 1Full-screen View

Correlation coefficient (r) and relative mean error (RME) of PSO-MLEM, PSO, MLEM and MAXED unfolded spectra

Unfolding method	r and RME	Input spectrum
		Input-Spec1	Input-Spec2	Input-Spec3	Input-Spec4
PSO-MLEM	$r_{{\phi }^{\text{cal}}{\phi }^{\text{input}}}$	1.0000	0.9989	0.9940	0.9991
	$RM{E}_{{\phi }^{\text{cal}}{\phi }^{\text{input}}}$	0.0132	0.0203	0.0664	0.0195
PSO	$r_{{\phi }^{\text{cal}}{\phi }^{\text{input}}}$	0.9862	0.8123	0.7467	0.8815
	$RM{E}_{{\phi }^{\text{cal}}{\phi }^{\text{input}}}$	1.6720	2.1514	3.7010	0.2062
MLEM	$r_{{\phi }^{\text{cal}}{\phi }^{\text{input}}}$	0.9877	0.9481	0.7447	0.8655
	$RM{E}_{{\phi }^{\text{cal}}{\phi }^{\text{input}}}$	0.6627	1.0882	3.0150	0.1978
MAXED	$r_{{\phi }^{cal}{\phi }^{input}}$	0.9917	0.9037	0.9753	0.9914
	$RM{E}_{{\phi }^{\text{cal}}{\phi }^{\text{input}}}$	0.1305	0.3528	0.2035	0.0719

Show more

The time complexity of the algorithm is represented by Big O notation, and the time complexity of the PSO-MLEM algorithm is O(T*N), which is linearly related to the number of iterations T and population N of the algorithm. Fig. 4 shows a plot between the minimum cost of Eq. (5) and the number of iterations of the PSO-MLEM, PSO and MLEM algorithms, which indicate the convergence of algorithm. From the Fig.4, the minimum cost of the PSO-MLEM algorithm is less than that of the PSO and MLEM algorithms when the convergence condition is achieved. And when the convergence condition is achieved, the iterations number of the PSO-MLEM algorithm is 500, the iteration number of the MLEM algorithm is 700, the iterations number of the PSO algorithm is 1550. The convergence speed of the PSO-MLEM algorithm is 1.4 times and 3.1 times higher than that of MLEM and PSO algorithms.

Fig. 4Full-screen View

(Color online) A plot of the algorithm convergence between the minimum cost and the number of iterations of the PSO-MLEM, PSO and MLEM algorithms.