The study of a neutron spectrum unfolding method based onparticle swarm optimization combined with maximum likelihood expectation maximization

Hong-Fei Xiao ·Qing-Xian Zhang ·He-Yi Tan ·Bin Shi ·Jun Chen ·Zhi-Qiang Cheng ·Jian Zhang ·Rui Yang

Abstract The neutron spectrum unfolding by Bonner sphere spectrometer (BSS) is considered a complex multidimensional model,which requires complex mathematical methods to solve the first kind of Fredholm integral equation.In order to solve the problem of the maximum likelihood expectation maximization (MLEM) algorithm which is easy to suffer the pitfalls of local optima and the particle swarm optimization (PSO) algorithm which is easy to get unreasonable flight direction and step length of particles, which leads to the invalid iteration and affect efficiency and accuracy, an improved PSO-MLEM algorithm, combined of PSO and MLEM algorithm, is proposed for neutron spectrum unfolding.The dynamic acceleration factor is used to balance the ability of global and local search, and improves the convergence speed and accuracy of the algorithm.Firstly, the Monte Carlo method was used to simulated the BSS to obtain the response function and count rates of BSS.In the simulation of count rate, four reference spectra from the IAEA Technical Report Series No.403 were used as input parameters of the Monte Carlo method.The PSO-MLEM algorithm was used to unfold the neutron spectrum of the simulated data and was verified by the difference of the unfolded spectrum to the reference spectrum.Finally, the 252Cf neutron source was measured by BSS, and the PSO-MLEM algorithm was used to unfold the experimental neutron spectrum.Compared with maximum entropy deconvolution (MAXED), PSO and MLEM algorithm, the PSO-MLEM algorithm has fewer parameters and automatically adjusts the dynamic acceleration factor to solve the problem of local optima.The convergence speed of the PSO-MLEM algorithm is 1.4 times and 3.1 times that of the MLEM and PSO algorithms.Compared with PSO, MLEM and MAXED, the correlation coefficients of PSO-MLEM algorithm are increased by 33.1%, 33.5% and 1.9%, and the relative mean errors are decreased by 98.2%, 97.8% and 67.4%.

Keywords Particle swarm optimization ·Maximum likelihood expectation maximization ·Neutron spectrum unfolding ·Bonner spheres spectrometer ·Monte Carlo method

1 Introduction

The precise measurement of the neutron spectrum distribution in the radiation field is a basic process in nuclear physics research, radiation protection and monitoring, reactor safety and control, particle accelerator and other nuclear technology fields.During the process of neutron spectrum measurement, it is necessary to unfold the neutron spectrum.Because the neutron spectrum unfolding procedure is an underdetermined problem with numerous viable solutions, a specialized optimization algorithm must be used to find the optimal solution [1-4].

The neutron spectrum unfolding is a process of solving underdetermined equations, and at the earliest, the least square method is used to unfold the neutron spectrum.However, the unfolded spectrum obtained through multiple iterations is quite different from the reference spectrum.However, by adding the prior information, the results could be improved to a certain extent [5].Engl et al.used a regularization method to unfold the neutron spectrum, which constrained the solution by adding the prior information to improve the stability and continuity of the unfolded spectrum [6].Reginatto et al.used the NE213 scintillation detector to measure the neutron spectrum.Moreover, the authors combined the Bayesian analysis method with the maximum entropy principle to unfold the neutron spectrum and to achieving relatively high accuracy [7, 8].Artificial intelligence algorithms were frequently utilized in the neutron spectrum unfolding [9].Freeman et al.based on BSS used a genetic algorithm to unfold the neutron spectrum.This method used a computer-simulated "chromosome" to encode the neutron spectrum, and then, the individuals with the best fitness was chosen as the final result.Compared with other neutron spectrum unfolding methods, the extracted results show that the genetic algorithm can obtain more accurate unfolded spectrum in the absence of prior information [10].Subsequently, Sharghi Ido et al.developed two different artificial neural network methods to unfold the neutron spectrum(multilayer and single-layer).With this method, the human nervous system was simulated and the neural network model was trained with known inputs and outputs to adjust its internal parameters to appropriate values.Finally, Am-Be and252Cf neutron sources were measured, and the unfolded spectrum have a good agreement [11].M.Shahmohammadi Beni et al.proposed the maximum likelihood expectation maximization algorithm to unfold the neutron spectrum and unfolded the Am-Be neutron spectrum in the absence of prior information and got good results [3].The MLEM algorithm has strong convergence ability and high stability in the early stage of iteration, but the later iteration is easy to suffer the pitfalls of local optima.

The PSO algorithm was developed by Eberhart and Kennedy in 1995 and is based on bird foraging behavior [12].More specifically, the PSO begins with a randomly generated initial population (known as a particle swarm) that is randomly distributed as a candidate solution in the search space.The PSO contains the velocity and position for each particle(solution vector).The speed of each particle must be dynamically modified based on its own flight experience, as well as the flight experiences of the surrounding particles.On top of that, the position of each particle must be updated by using the cost function [13].PSO algorithm has the advantage of global convergence, strong randomness of search, and adjustable local search ability and global search ability in the whole iteration process.H.Shahabinejad et al.developed a novel neutron spectrum unfolding code (SDPSO) based on particle swarm optimization algorithm and performed neutron spectrum unfolding for Am-Be neutron sources [14].The disadvantage is that the highly random search process is easy to get unreasonable flight direction and step length of particles, which leads to the invalid iteration and affect efficiency and accuracy.In this regard, many scholars have studied traditional PSO from the perspective of parameter setting, convergence and combining with other algorithms to improve algorithms.Some authors use the inertia weights to balance the global search capability and the local search capability [15-17].Some authors propose combination algorithms, which are based on the global search advantage of PSO algorithm, such as the combination of PSO and genetic algorithm [18], PSO and k-clustering algorithm [19].

In this paper, an improved PSO-MLEM algorithm, combined of PSO algorithm and MLEM algorithm, is proposed for neutron spectrum unfolding.The PSO-MLEM algorithm has the advantages of PSO and MLEM, which can avoid the pitfalls of local optima through the adjustment of dynamic acceleration factors, and improve the convergence speed and accuracy.

2 Method andmaterials

2.1 Theoretical model of Bonner spheres spectrometer

The Bonner sphere 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:3He,10B and6Li.The3He 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 [20-22]:

whereNiis the count rate of theith Bonner sphere;εiis the measurement error ofith Bonner sphere;EmaxandEminare the highest and lowest energy of the neutron;is the Bonner sphere response function corresponding to thejth energy interval;s the fluence rate of neutrons in thejth energy interval;nis the number of Bonner spheres; andmis the number of energy intervals.

whereN=(N1,N2,…,Nn)T,φ=(φ1,φ2,…,φm)TandRis the response matrix ofn×m.

Usually,NandRare known quantities, and the number of Bonner spheresnis often less than the number of the energy intervalm.Equation (3) belongs to the underdetermined equations.In this work, the PSO-MLEM algorithm was used to solve the underdetermined problem in neutron spectrum unfolding.

2.2 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 PSOMLEM algorithm.The MLEM algorithm is used to limit the search direction of particles during each iteration of the PSO-MLEM algorithm.Figure 1a 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 smallerc1and a largerc2in 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 PSOMLEM algorithm gradually increasesc1and decreasesc2in 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.

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φ=[φ1,φ2,…,φm].The solution vectorφ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:

whereNobsis the count rate of Bonner spheres;φis the solution of each particle in each iteration;φrefis the reference spectrum that contains the a priori information; andNcalis calculated by using Eq.(5) as follows:

Through the reference spectrumφrefand the count rateNobs, the above cost function limits the search direction of the particle.The first term of the cost function is the distance betweenNcalandNobs, and the second term is the distance between the updated spectrumφand the reference spectrumφref.When the minimum of Cost function is obtained, theφ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)

Fig.1 Graphical depiction of PSO-MLEM algorithm.a Schematic diagram of PSO-MLEM particle flight; b flowchart of PSO-MLEM algorithm

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 velocityvijand solutionφijof thejth dimension of theith particle are updated:

whereviis the velocity vector of theith particle,φiis the solution vector of theith particle, thepbestiis the historical local optima solution of theith particle,gbestis the global optima solution of all particles,ris a random number between [0, 1],tis the number of iterations, andcis the acceleration factor between [0, 1], as shown in Eq.(8):

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 flowchart summarizes the preceding 6 steps in Fig.1b.

2.3 Verification by Monte Carlo simulation

In order to verify the accuracy and convergence of the PSOMLEM 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 [23] 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 [24].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 [25].Each sphere has a thermal neutron detector in the center.The center detector chooses the SP903He proportional counter manufactured by Centronics of the UK.The Bonner sphere model is shown in Fig.2a.

The response functions of Bonner sphere 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 and3He 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 ofEis defined as follows [26]:

whereRd(E) is the energy response (cm2) when the incident neutron energy isE,Mdrepresents the count rate of the Bonner sphere with diameterdand Φ(E) is the neutron fluence (cm-2) when the incident neutron energy isE.Figure 2b depicts the response functions of 15 Bonner spheres as simulated by the Monte Carlo method.

As shown in Fig.2b, 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 toward higher energies.

2.3.2 Validation oftheunfolded spectrum

As shown in Fig.3a1, the Input-Spec1 is from the CERNCEC 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 Fig.3a1, 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 PSOMLEM, PSO, MLEM and MAXED were used to unfold the neutron spectrum.

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

Fig.3 (Color online) 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: ratio of NPSO_MLEM ,NPSO , NMLEM , NMAXED and Ndef to Nobs ); c unfolded spectrum of the Spec2 (c1: Comparison of PSO-MLEM, PSO, MLEM and MAXED unfolded spectra with the input spectrum; c2: ratio of NPSO_MLEM ,NPSO , NMLEM , NMAXED and Ndef to Nobs ); d unfolded spectrum of the Spec3 (d1: comparison of PSO-MLEM, PSO, MLEM and MAXED unfolded spectra with the input spectrum; d2: ratio of NPSO_MLEM ,NPSO , NMLEM , NMAXED and Ndef to Nobs ); e unfolded spectrum of the Spec4 (e1: Comparison of PSO-MLEM, PSO, MLEM and MAXED unfolded spectra with the input spectrum; e2: ratio of NPSO_MLEM ,NPSO , NMLEM , NMAXED and Ndef to Nobs)

Fig.3 (continued)

Figure 3(b1)-(e1) depicts 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) presents the ratio relation betweenNdef∕Nobs,NPSO_MLEM∕Nobs,NPSO∕Nobs,NMLEM∕NobsandNMAXED∕Nobs.Nobsare the count rates of 15 Bonner spheres, andNdef(Ndef=Rφinput) is the result of the response function times the input spectrum.Besides,NPSO_MLEM(NPSO_MLEM=RφPSO_MLEM) is the result of the response function times the unfolded spectrum of the PSO-MLEM algorithm,NPSO(NPSO=RφPSO) is the result of the response function times the unfolded spectrum of the PSO algorithm,NMLEM(NMLEM=RφMLEM) is the result of the response function times the unfolded spectrum of the MLEM algorithm, andNMAXED(NMAXED=RφMAXED) 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φcaland the input spectrumφinput, and its value ranges from + 1 to -1.Thercan be calculated using the following equation:

The relative mean error (RME) was also used to calculate the difference of the unfolded spectrum from the input spectrum.TheRMEcan be calculated using the following equation:

wheremis the sample size, represents the number of energy intervalsi, andand, respectively, represent the normalized value of theith channel corresponding to the unfolded spectra and the input spectrum.For Eq.(11), whenr=0, the two variables are uncorrelated, whenr＞0, there is a positive correlation, whenr＜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 (φcal=φinput),RMEis 0.

As given in Table 1, 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 errorRMEis 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 fluctuates more strongly, and 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%.

Fig.4 (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

The time complexity of the algorithm is represented by BigOnotation, and the time complexity of the PSOMLEM algorithm isO(T*N), which is linearly related to the number of iterationsTand populationNof the algorithm.Figure 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 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, and 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.

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

Fig.5 Experimental Bonner sphere model and response functions.a Photograph of Bonner spheres spectrometer; b response functions of Bonner spheres spectrometer

3 Experiment

The Bonner spheres spectrometer designed by the China Institute of Atomic Energy was used for experiment to further evaluate the viability of the algorithm [27].As shown in Fig.5a, the Bonner spheres spectrometer consists of one SP93He proportional counter and ten polyethylene spheres with diameters of 1, 2.5, 3, 3.5, 4, 5, 6, 8, 10, 12 inches.The3He neutron proportional counter as the central detector is wrapped in 10 polyethylene spheres.The neutron energy range of the response function is between 9.441 × 10-10and 18.84 MeV, which is divided into 207 energy groups.The response function of BSS is shown in Fig.5b.

The BSS of the China Institute of Atomic Energy is used to measure the reference252Cf neutron source, and the count rates of BSSN=(N1,N2,…,N11)Tare obtained and used to unfold the neutron spectrum by the PSO-MLEM algorithm.The reference spectrum and unfolded spectrum by the PSOMLEM algorithm are shown in in Fig.6.

Table 2 presents the correlation coefficient (r) and relative mean error (RME) of unfolded spectrum by PSO-MLEM algorithm to the reference spectrum.The error analysis of the unfolded spectrum shows that the correlation coefficient of the unfolded spectrumφcalis more than the 0.99.The relative mean error of the unfolded spectrum is 0.1718%.

4 Conclusion

Fig.6 (Color online).252Cf neutron source unfolded spectrum (a: Comparison of PSO-MLEM unfolded spectrum with reference spectrum; b:Ratio of Ncal、Ndef and Nobs of each detector)

Table 2 Correlation coefficient (r) and relative mean error (RME) of unfolded spectrum

In this work, to solve the problem of the MLEM algorithm which is easy to suffer the pitfalls of local optima and the problem of the PSO algorithm which is easy to get unreasonable flight direction and step length of particles, which leads to the invalid iteration and affect efficiency and accuracy, an improved PSO-MLEM algorithm was proposed for neutron spectrum unfolding.Based on the traditional PSO algorithm, PSO and MLEM algorithm are mixed, and the convergence speed and accuracy of the algorithm are improved by dynamic acceleration factor, and the algorithm will not suffer the pitfalls of the local optima.

The Monte Carlo method was employed to simulate the response functions of the BBS, and four reference spectra from the IAEA Technical Report Series No.403 were used as the input spectra.The effects caused by prior information with random error and fluctuated more strongly of spectrum on the unfolded spectrum were investigated.The results show that the PSO-MLEM algorithm is not sensitive to the fluctuation of prior information, and the unfolded spectrum has a strong correlation with the input spectrum of Monte Carlo method.In addition, the PSO-MLEM algorithm has strong robustness to the neutron spectrum with large fluctuations.Besides, PSO-MLEM algorithm is used to unfold the experimental data of252Cf neutron source, and the result shows that the PSO-MLEM algorithm is competent for neutron spectrum unfolding.

Author contributionsAll authors contributed to the study conception and design.Material preparation, data collection and analysis were performed by Hong-Fei Xiao, Qing-Xian Zhang, He-Yi Tan, Bin Shi,Jun Chen, Zhi-Qiang Cheng, Jian Zhang and Rui Yang.The first draft of the manuscript was written by Hong-Fei Xiao and Qing-Xian Zhang,and all authors commented on previous versions of the manuscript.All authors read and approved the final manuscript.

Data AvailabilityThe data that support the findings of this study are openly available in Science Data Bank at https:// doi.org/ 10.57760/ scien cedb.07741 and https:// cstr.cn/ 31253.11.scien cedb.07741.