The study of a neutron spectrum unfolding method based on Particle Swarm Optimization combined with Maximum Likelihood Expectation Maximization

The study of a neutron spectrum unfolding method based on Particle Swarm Optimization combined with Maximum Likelihood Expectation Maximization 1 first-author https://orcid.org/0000-0003-1391-8243 1 corresp * shinecore@163.com * shinecore@163.com shinecore@163.com 1 https://orcid.org/0000-0001-5949-1500 2 2 1 1 1 Applied Nuclear Technology in Geosciences Key Laboratory of Sichuan Province, Chengdu University of Technology, Chengdu 610059, China Applied Nuclear Technology in Geosciences Key Laboratory of Sichuan Province, Chengdu University of Technology, Chengdu 610059, China National Key Laboratory for Metrology and Calibration Techniques, China Institute of Atomic Energy, Beijing 102413, China National Key Laboratory for Metrology and Calibration Techniques, China Institute of Atomic Energy, Beijing 102413, China The neutron spectrum unfolding by Bonner Spheres 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 252 Cf 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 increase by 33.1%, 33.5% and 1.9%, and the relative mean errors decreased by 98.2%, 97.8% and 67.4%. Particle swarm optimization Maximum Likelihood Expectation Maximization Neutron spectrum unfolding Bonner Spheres Spectrometer Monte Carlo method Introduction 1 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 1 - 4 4 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 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 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 7 - 8 8 7-8 ] . Artificial intelligence algorithms were frequently utilized in the neutron spectrum unfolding [ 9 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 10 ] . Subsequently, Sharghi Ido et al. developed two different artificial neural network methods to unfold the neutron spectrum (multi-layer 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 and 252 Cf neutron sources were measured, and the unfolded spectrum have a good agreement [ 11 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 [ 12 12 ] . 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 [ 13 13 ] . 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 [ 14 14 ] . 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 [ 15 15 ] . 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 [ 16 16 - 18 18 16-18 ] . 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 [ 19 19 ] , PSO and k-clustering algorithm [ 20 20 ] . 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. Method and Materials 1 Theoretical model of Bonner Spheres Spectrometer 2 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: 3 He, 10 B and 6 Li. The 3 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 21 - 23 23 21-23 ] : <math display="block" id="M1"><mtable><mtr><mtd><msub><mi>N</mi><mi>i</mi></msub><mo>+</mo><msub><mi>ε</mi><mi>i</mi></msub><mo>=</mo><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∫</mo></mstyle><mrow><msub><mi>E</mi><mrow><mtext>min</mtext></mrow></msub></mrow><mrow><msub><mi>E</mi><mrow><mtext>max</mtext></mrow></msub></mrow></munderover><msub><mi>R</mi><mi>i</mi></msub><mrow><mo>(</mo><mrow><msub><mi>E</mi><mi>j</mi></msub></mrow><mo>)</mo></mrow><mi>Φ</mi><mrow><mo>(</mo><mrow><msub><mi>E</mi><mi>j</mi></msub></mrow><mo>)</mo></mrow><mtext>d</mtext><mi>E</mi><mo>,</mo></mtd></mtr><mtr><mtd><mi>i</mi><mo>=</mo><mn>1</mn><mo>,</mo><mn>2</mn><mo>,</mo><mo>⋯</mo><mo>,</mo><mi>n</mi><mo>;</mo><mo> </mo><mi>j</mi><mo>=</mo><mn>1</mn><mo>,</mo><mn>2</mn><mo>,</mo><mo>⋯</mo><mo>,</mo><mi>m</mi></mtd></mtr></mtable></math> where, <math id="M2"><mrow><msub><mi>N</mi><mi>i</mi></msub></mrow></math> is the count rate of the i th Bonner Sphere; <math id="M3"><mrow><msub><mi>ε</mi><mi>i</mi></msub></mrow></math> is the measurement error of i th Bonner Sphere. <math id="M4"><mrow><msub><mi>E</mi><mrow><mtext>max</mtext></mrow></msub></mrow></math> and <math id="M5"><mrow><msub><mi>E</mi><mrow><mtext>min</mtext></mrow></msub></mrow></math> are the highest and lowest energy of the neutron; <math id="M6"><mrow><mtext> </mtext><msub><mi>R</mi><mi>i</mi></msub><mrow><mo>(</mo><mrow><msub><mi>E</mi><mi>j</mi></msub></mrow><mo>)</mo></mrow></mrow></math> is the Bonner Sphere response function corresponding to the j th energy interval; <math id="M7"><mrow><mi>Φ</mi><mrow><mo>(</mo><mrow><msub><mi>E</mi><mi>j</mi></msub></mrow><mo>)</mo></mrow></mrow></math> is the fluence rate of neutrons in the j th energy interval; n is the number of Bonner Spheres; m is the number of energy intervals. Since the <math id="M8"><mrow><msub><mi>R</mi><mi>i</mi></msub><mrow><mo>(</mo><mrow><msub><mi>E</mi><mi>j</mi></msub></mrow><mo>)</mo></mrow></mrow></math> and <math id="M9"><mrow><msub><mi>N</mi><mi>i</mi></msub></mrow></math> are discretely distributed, Eq. (1) can be expressed in the discrete form of Eq. (2): <math display="block" id="M10"><mrow><msub><mi>N</mi><mi>i</mi></msub><mo>+</mo><msub><mi>ε</mi><mi>i</mi></msub><mo>=</mo><mstyle displaystyle="true"><munder><mo>∑</mo><mi>j</mi></munder><mrow><msub><mi>R</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub><msub><mi>φ</mi><mi>j</mi></msub></mrow></mstyle><mo> </mo><mo> </mo><mo> </mo><mi>i</mi><mo>=</mo><mn>1</mn><mo>,</mo><mn>2</mn><mo>,</mo><mo>⋯</mo><mo>,</mo><mi>n</mi><mo>;</mo><mo> </mo><mo> </mo><mi>j</mi><mo>=</mo><mn>1</mn><mo>,</mo><mn>2</mn><mo>,</mo><mo>⋯</mo><mo>,</mo><mi>m</mi></mrow></math> Eq. (2) can also be expressed in matrix form as follows: <math display="block" id="M11"><mrow><mi>N</mi><mo>+</mo><mi>ε</mi><mo>=</mo><mi>R</mi><mi>φ</mi><mo>,</mo></mrow></math> where, <math id="M12"><mrow><mi>N</mi><mo>=</mo><msup><mrow><mo stretchy="false">(</mo><msub><mi>N</mi><mn>1</mn></msub><mo>,</mo><msub><mi>N</mi><mn>2</mn></msub><mo>,</mo><mo>…</mo><mo>,</mo><msub><mi>N</mi><mi>n</mi></msub><mo stretchy="false">)</mo></mrow><mi>T</mi></msup></mrow></math> , <math id="M13"><mrow><mi>φ</mi><mo>=</mo><msup><mrow><mrow><mo>(</mo><mrow><msub><mi>φ</mi><mn>1</mn></msub><mo>,</mo><msub><mi>φ</mi><mn>2</mn></msub><mo>,</mo><mo>…</mo><mo>,</mo><msub><mi>φ</mi><mi>m</mi></msub></mrow><mo>)</mo></mrow></mrow><mi>T</mi></msup></mrow></math> , 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 2 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) 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 <math id="M14"><mrow><msub><mi>c</mi><mn>1</mn></msub></mrow></math> and a larger <math id="M15"><mrow><msub><mi>c</mi><mn>2</mn></msub></mrow></math> 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 <math id="M16"><mrow><msub><mi>c</mi><mn>1</mn></msub></mrow></math> and decreases <math id="M17"><mrow><msub><mi>c</mi><mn>2</mn></msub></mrow></math> 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. 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 <math id="M18"><mrow><mi>φ</mi><mo>=</mo><msup><mrow><mrow><mo>[</mo><mrow><msub><mi>φ</mi><mn>1</mn></msub><mo>,</mo><msub><mi>φ</mi><mn>2</mn></msub><mo>,</mo><mo>…</mo><mo>,</mo><msub><mi>φ</mi><mi>m</mi></msub></mrow><mo>]</mo></mrow></mrow><mo> </mo></msup></mrow></math> . The solution vector <math id="M19"><mi>φ</mi></math> 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: <math display="block" id="M20"><mrow><mtext>Cost</mtext><mo>=</mo><mo> </mo><mfrac><mrow><msup><mrow><mrow><mo>|</mo><mrow><mrow><mo>|</mo><mrow><msub><mi>N</mi><mrow><mtext>cal</mtext></mrow></msub><mo>−</mo><msub><mi>N</mi><mrow><mtext>obs</mtext></mrow></msub></mrow><mo>|</mo></mrow></mrow><mo>|</mo></mrow></mrow><mn>2</mn></msup></mrow><mrow><msup><mrow><mrow><mo>|</mo><mrow><mrow><mo>|</mo><mrow><msub><mi>N</mi><mrow><mtext>obs</mtext></mrow></msub></mrow><mo>|</mo></mrow></mrow><mo>|</mo></mrow></mrow><mn>2</mn></msup></mrow></mfrac><mo>+</mo><mfrac><mrow><msup><mrow><mrow><mo>|</mo><mrow><mrow><mo>|</mo><mrow><mi>φ</mi><mo>−</mo><msub><mi>φ</mi><mrow><mtext>ref</mtext></mrow></msub></mrow><mo>|</mo></mrow></mrow><mo>|</mo></mrow></mrow><mn>2</mn></msup></mrow><mrow><msup><mrow><mrow><mo>|</mo><mrow><mrow><mo>|</mo><mrow><msub><mi>φ</mi><mrow><mtext>ref</mtext></mrow></msub></mrow><mo>|</mo></mrow></mrow><mo>|</mo></mrow></mrow><mn>2</mn></msup></mrow></mfrac><mo> </mo><mo>,</mo></mrow></math> where, <math id="M21"><mrow><msub><mi>N</mi><mrow><mtext>obs</mtext></mrow></msub></mrow></math> is the count rate of Bonner spheres; <math id="M22"><mi>φ</mi></math> is the solution of each particle in each iteration; <math id="M23"><mrow><msub><mi>φ</mi><mrow><mtext>ref</mtext></mrow></msub></mrow></math> is the reference spectrum that contains the a priori information; <math id="M24"><mrow><msub><mi>N</mi><mrow><mtext>cal</mtext></mrow></msub></mrow></math> is calculated by using Eq. (5) as follows: <math display="block" id="M25"><mrow><msub><mi>N</mi><mrow><mtext>cal</mtext></mrow></msub><mo>=</mo><mi>R</mi><mi>φ</mi><mo>.</mo></mrow></math> Through the reference spectrum <math id="M26"><mrow><msub><mi>φ</mi><mrow><mtext>ref</mtext></mrow></msub></mrow></math> and the count rate <math id="M27"><mrow><msub><mi>N</mi><mrow><mtext>obs</mtext></mrow></msub></mrow></math> , the above cost function limits the search direction of the particle. The first term of the Cost function is the distance between <math id="M28"><mrow><msub><mi>N</mi><mrow><mtext>cal</mtext></mrow></msub></mrow></math> and <math id="M29"><mrow><msub><mi>N</mi><mrow><mtext>obs</mtext></mrow></msub></mrow></math> , and the second term is the distance between the updated spectrum <math id="M30"><mi>φ</mi></math> and the reference spectrum <math id="M31"><mrow><msub><mi>φ</mi><mrow><mtext>ref</mtext></mrow></msub></mrow></math> . When the minimum of Cost function is obtained, the <math id="M32"><mi>φ</mi></math> 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 <math id="M33"><mrow><msub><mi>v</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub></mrow></math> and solution <math id="M34"><mrow><msub><mi>φ</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub></mrow></math> of the j th dimension of the i th particle are updated: <math display="block" id="M35"><mrow><msub><mi>v</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub><mrow><mo>(</mo><mrow><mi>t</mi><mo>+</mo><mn>1</mn></mrow><mo>)</mo></mrow><mo>=</mo><msub><mi>r</mi><mn>1</mn></msub><mo>⋅</mo><mrow><mo>(</mo><mrow><mi>p</mi><mi>b</mi><mi>e</mi><mi>s</mi><msub><mi>t</mi><mi>i</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>−</mo><msub><mi>φ</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow><mo>)</mo></mrow><mo>+</mo><msub><mi>r</mi><mn>2</mn></msub><mo>⋅</mo><mrow><mo>(</mo><mrow><mi>g</mi><mi>b</mi><mi>e</mi><mi>s</mi><mi>t</mi><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>−</mo><msub><mi>φ</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> <math display="block" id="M36"><mrow><msub><mi>φ</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub><mrow><mo>(</mo><mrow><mi>t</mi><mo>+</mo><mn>1</mn></mrow><mo>)</mo></mrow><mo>=</mo><mi>c</mi><mn>1</mn><mo>⋅</mo><msub><mi>v</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub><mrow><mo>(</mo><mrow><mi>t</mi><mo>+</mo><mn>1</mn></mrow><mo>)</mo></mrow><mo>+</mo><mi>c</mi><mn>2</mn><mo>⋅</mo><mo> </mo><msubsup><mi>φ</mi><mrow><mi>i</mi><mi>j</mi></mrow><mrow><mtext>MLEM</mtext></mrow></msubsup><mrow><mo>(</mo><mrow><mi>t</mi><mo>+</mo><mn>1</mn></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> where, <math id="M37"><mrow><msub><mi>v</mi><mi>i</mi></msub></mrow></math> is the velocity vector of the i th particle, <math id="M38"><mrow><msub><mi>φ</mi><mi>i</mi></msub></mrow></math> is the solution vector of the i th particle, the <math id="M39"><mrow><mi>p</mi><mi>b</mi><mi>e</mi><mi>s</mi><msub><mi>t</mi><mi>i</mi></msub></mrow></math> is the historical local optima solution of the i th particle, and <math id="M40"><mrow><mi>g</mi><mi>b</mi><mi>e</mi><mi>s</mi><mi>t</mi></mrow></math> is the global optima solution of all particles, <math id="M41"><mi>r</mi></math> is a random number between [0, 1], <math id="M42"><mi>t</mi></math> is the number of iterations, <math id="M43"><mi>c</mi></math> is the acceleration factor between [0, 1], as shown in Eq. (8): <math display="block" id="M44"><mrow><mrow><mo>{</mo><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow><msub><mi>c</mi><mn>1</mn></msub><mo>=</mo><mfrac><mi>t</mi><mrow><msub><mi>T</mi><mrow><mtext>max</mtext></mrow></msub></mrow></mfrac></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow><msub><mi>c</mi><mn>2</mn></msub><mo>=</mo><mn>1</mn><mo>−</mo><mrow><mo>(</mo><mrow><mfrac><mi>t</mi><mrow><msub><mi>T</mi><mrow><mtext>max</mtext></mrow></msub></mrow></mfrac></mrow><mo>)</mo></mrow></mrow></mtd></mtr></mtable></mrow></mrow></mrow></math> where, <math id="M45"><mi>t</mi></math> is the current number of iterations, and <math id="M46"><mrow><msub><mi>T</mi><mrow><mtext>max</mtext></mrow></msub></mrow></math> is the maximum number of iterations. <math id="M47"><mrow><msubsup><mi>φ</mi><mrow><mi>i</mi><mi>j</mi></mrow><mrow><mtext>MLEM</mtext></mrow></msubsup></mrow></math> is calculated according to the MLEM algorithm using Eq. (9): <math display="block" id="M48"><mrow><msubsup><mi>φ</mi><mrow><mi>i</mi><mi>j</mi></mrow><mrow><mtext>MLEM</mtext></mrow></msubsup><mrow><mo>(</mo><mrow><mi>t</mi><mo>+</mo><mn>1</mn></mrow><mo>)</mo></mrow><mo>=</mo><msub><mi>φ</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>·</mo><mfrac><mrow><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></msubsup><mrow><msub><mi>R</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub></mrow></mstyle><mfrac><mrow><msub><mi>N</mi><mi>i</mi></msub></mrow><mrow><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mi>m</mi></msubsup><mrow><msub><mi>R</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub></mrow></mstyle><msub><mi>φ</mi><mi>j</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow></mrow></mfrac></mrow><mrow><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></msubsup><mrow><msub><mi>R</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub></mrow></mstyle></mrow></mfrac><mo>.</mo></mrow></math> 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) Fig.1(b) : Verification by Monte Carlo simulation 2 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 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 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. Bonner Spheres Spectrometer response function simulation 3 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 26 ] . Each sphere has a thermal neutron detector in the center. The center detector chooses the SP90 3 He proportional counter manufactured by Centronics of the UK. The Bonner sphere model is shown in Fig.2(a) Fig.2(a) . 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 3 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 27 ] : <math display="block" id="M49"><mrow><msub><mi>R</mi><mtext>d</mtext></msub><mrow><mo>(</mo><mrow><msub><mi>E</mi><mtext>n</mtext></msub></mrow><mo>)</mo></mrow><mo>=</mo><mfrac><mrow><msub><mi>M</mi><mtext>d</mtext></msub></mrow><mrow><mi>Φ</mi><mrow><mo>(</mo><mrow><msub><mi>E</mi><mtext>n</mtext></msub></mrow><mo>)</mo></mrow></mrow></mfrac><mo>,</mo></mrow></math> where, <math id="M50"><mrow><msub><mi>R</mi><mtext>d</mtext></msub><mrow><mo>(</mo><mi>E</mi><mo>)</mo></mrow></mrow></math> is the energy response (cm2) when the incident neutron energy is E, <math id="M51"><mrow><msub><mi>M</mi><mtext>d</mtext></msub></mrow></math> represents the count rate of the Bonner sphere with diameter d , and <math id="M52"><mrow><mi>Φ</mi><mrow><mo>(</mo><mi>E</mi><mo>)</mo></mrow></mrow></math> is the neutron fluence (cm – 2 ) when the incident neutron energy is E . Figure 2(b) Figure 2(b) depicts the response functions of 15 Bonner spheres as simulated by the Monte Carlo method. As Fig. 2(b) 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. Validation of the unfolded spectrum 3 As shown in Fig.3(a1) 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) 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) 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. Figure 3(b1) - (e1) 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) Figure 3(b2) - (e2) is the ratio relation between <math id="M77"><mrow><msub><mi>N</mi><mrow><mtext>def</mtext></mrow></msub><mo>/</mo><msub><mi>N</mi><mrow><mtext>obs</mtext></mrow></msub></mrow></math> , <math id="M78"><mrow><msub><mi>N</mi><mrow><mtext>PSO</mtext><mo>_</mo><mtext>MLEM</mtext></mrow></msub><mo>/</mo><msub><mi>N</mi><mrow><mtext>obs</mtext></mrow></msub></mrow></math> , <math id="M79"><mrow><msub><mi>N</mi><mrow><mtext>PSO</mtext></mrow></msub><mo>/</mo><msub><mi>N</mi><mrow><mtext>obs</mtext></mrow></msub></mrow></math> , <math id="M80"><mrow><msub><mi>N</mi><mrow><mtext>MLEM</mtext></mrow></msub><mo>/</mo><msub><mi>N</mi><mrow><mtext>obs</mtext></mrow></msub></mrow></math> and <math id="M81"><mrow><msub><mi>N</mi><mrow><mtext>MAXED</mtext></mrow></msub><mo>/</mo><msub><mi>N</mi><mrow><mtext>obs</mtext></mrow></msub></mrow></math> . <math id="M82"><mrow><msub><mi>N</mi><mrow><mtext>obs</mtext></mrow></msub></mrow></math> are the count rates of 15 Bonner spheres, <math id="M83"><mrow><msub><mi>N</mi><mrow><mtext>def</mtext></mrow></msub><mrow><mo>(</mo><mrow><msub><mi>N</mi><mrow><mtext>def</mtext></mrow></msub><mo>=</mo><mi>R</mi><msub><mi>φ</mi><mrow><mtext>input</mtext></mrow></msub></mrow><mo>)</mo></mrow></mrow></math> is the result of the response function times the input spectrum. Besides, <math id="M84"><mrow><msub><mi>N</mi><mrow><mtext>PSO</mtext><mo>_</mo><mtext>MLEM</mtext></mrow></msub><mrow><mo>(</mo><mrow><msub><mi>N</mi><mrow><mtext>PSO</mtext><mo>_</mo><mtext>MLEM</mtext></mrow></msub><mo>=</mo><mi>R</mi><msub><mi>φ</mi><mrow><mtext>PSO</mtext><mo>_</mo><mtext>MLEM</mtext></mrow></msub></mrow><mo>)</mo></mrow></mrow></math> is the result of the response function times the unfolded spectrum of the PSO-MLEM algorithm, <math id="M85"><mrow><msub><mi>N</mi><mrow><mtext>PSO</mtext></mrow></msub><mrow><mo>(</mo><mrow><msub><mi>N</mi><mrow><mtext>PSO</mtext></mrow></msub><mo>=</mo><mi>R</mi><msub><mi>φ</mi><mrow><mtext>PSO</mtext></mrow></msub></mrow><mo>)</mo></mrow></mrow></math> is the result of the response function times the unfolded spectrum of the PSO algorithm, <math id="M86"><mrow><msub><mi>N</mi><mrow><mtext>MLEM</mtext></mrow></msub><mrow><mo>(</mo><mrow><msub><mi>N</mi><mrow><mtext>MLEM</mtext></mrow></msub><mo>=</mo><mi>R</mi><msub><mi>φ</mi><mrow><mtext>MLEM</mtext></mrow></msub></mrow><mo>)</mo></mrow></mrow></math> is the result of the response function times the unfolded spectrum of the MLEM algorithm, and <math id="M87"><mrow><msub><mi>N</mi><mrow><mtext>MAXED</mtext></mrow></msub><mrow><mo>(</mo><mrow><msub><mi>N</mi><mrow><mtext>MAXED</mtext></mrow></msub><mo>=</mo><mi>R</mi><msub><mi>φ</mi><mrow><mtext>MAXED</mtext></mrow></msub></mrow><mo>)</mo></mrow></mrow></math> 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 <math id="M88"><mrow><msup><mi>φ</mi><mrow><mtext>cal</mtext></mrow></msup></mrow></math> and the input spectrum <math id="M89"><mrow><msup><mi>φ</mi><mrow><mtext>input</mtext></mrow></msup></mrow></math> , and its value ranges from +1 and –1. The r can calculate as following equation: <math display="block" id="M90"><mrow><msub><mi>r</mi><mrow><msup><mi>φ</mi><mrow><mtext>cal</mtext></mrow></msup><msup><mi>φ</mi><mrow><mtext>input</mtext></mrow></msup></mrow></msub><mo>=</mo><mfrac><mrow><mi>m</mi><mo>∑</mo><msubsup><mi>φ</mi><mi>i</mi><mrow><mtext>cal</mtext></mrow></msubsup><msubsup><mi>φ</mi><mi>i</mi><mrow><mtext>input</mtext></mrow></msubsup><mo>−</mo><mo>∑</mo><msubsup><mi>φ</mi><mi>i</mi><mrow><mtext>cal</mtext></mrow></msubsup><mo>∑</mo><msubsup><mi>φ</mi><mi>i</mi><mrow><mtext>input</mtext></mrow></msubsup></mrow><mrow><msqrt><mrow><mi>m</mi><mo>∑</mo><msup><mrow><mrow><mo>(</mo><mrow><msubsup><mi>φ</mi><mi>i</mi><mrow><mtext>cal</mtext></mrow></msubsup></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup><mo>−</mo><msup><mrow><mrow><mo>(</mo><mrow><mo>∑</mo><msubsup><mi>φ</mi><mi>i</mi><mrow><mtext>cal</mtext></mrow></msubsup></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow></msqrt><msqrt><mrow><mi>m</mi><mo>∑</mo><msup><mrow><mrow><mo>(</mo><mrow><msubsup><mi>φ</mi><mi>i</mi><mrow><mtext>input</mtext></mrow></msubsup></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup><mo>−</mo><msup><mrow><mrow><mo>(</mo><mrow><mo>∑</mo><msubsup><mi>φ</mi><mi>i</mi><mrow><mtext>input</mtext></mrow></msubsup></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow></msqrt></mrow></mfrac></mrow></math> 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: <math display="block" id="M91"><mrow><mi>R</mi><mi>M</mi><msub><mi>E</mi><mrow><msup><mi>φ</mi><mrow><mtext>cal</mtext></mrow></msup><msup><mi>φ</mi><mrow><mtext>input</mtext></mrow></msup></mrow></msub><mo>=</mo><mfrac><mn>1</mn><mi>m</mi></mfrac><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>m</mi></munderover><mrow><mrow><mo>|</mo> <mrow><mfrac><mrow><msubsup><mi>φ</mi><mi>i</mi><mrow><mtext>cal</mtext></mrow></msubsup><mo>−</mo><msubsup><mi>φ</mi><mi>i</mi><mrow><mtext>input</mtext></mrow></msubsup></mrow><mrow><msubsup><mi>φ</mi><mi>i</mi><mrow><mtext>input</mtext></mrow></msubsup></mrow></mfrac></mrow> <mo>|</mo></mrow></mrow></mstyle><mo>,</mo></mrow></math> where, <math id="M92"><mi>m</mi></math> is the sample size, represents the number of energy intervals, and <math id="M93"><mrow><msubsup><mi>φ</mi><mi>i</mi><mrow><mtext>cal</mtext></mrow></msubsup></mrow></math> and <math id="M94"><mrow><msubsup><mi>φ</mi><mi>i</mi><mrow><mtext>input</mtext></mrow></msubsup></mrow></math> respectively represent the normalized value of the i th channel corresponding to the unfolded spectra and the input spectrum. For Eq. (11), when <math id="M95"><mi>r</mi></math> =0, the two variables are uncorrelated, when <math id="M96"><mi>r</mi></math> > 0, there is a positive correlation, when <math id="M97"><mi>r</mi></math> < 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 ( <math id="M98"><mrow><msup><mi>φ</mi><mrow><mtext>cal</mtext></mrow></msup><mo>=</mo><msup><mi>φ</mi><mrow><mtext>input</mtext></mrow></msup></mrow></math> ), <math id="M99"><mrow><mi>R</mi><mi>M</mi><mi>E</mi></mrow></math> is 0. As can be seen from Table 1 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 <math id="M100"><mrow><mi>R</mi><mi>M</mi><mi>E</mi></mrow></math> 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%. 1 1 Unfolding method 1 1 r and RME 4 1 Input spectrum 1 1   1 1   1 1 Input-Spec1 1 1 Input-Spec2 1 1 Input-Spec3 1 1 Input-Spec4 1 1 PSO-MLEM 1 1 <math id="M101"><mrow><msub><mi>r</mi><mrow><msup><mi>φ</mi><mrow><mtext>cal</mtext></mrow></msup><msup><mi>φ</mi><mrow><mtext>input</mtext></mrow></msup></mrow></msub></mrow></math> 1 1 1.0000 1 1 0.9989 1 1 0.9940 1 1 0.9991 1 1   1 1 <math id="M102"><mrow><mi>R</mi><mi>M</mi><msub><mi>E</mi><mrow><msup><mi>φ</mi><mrow><mtext>cal</mtext></mrow></msup><msup><mi>φ</mi><mrow><mtext>input</mtext></mrow></msup></mrow></msub></mrow></math> 1 1 0.0132 1 1 0.0203 1 1 0.0664 1 1 0.0195 1 1 PSO 1 1 <math id="M103"><mrow><msub><mi>r</mi><mrow><msup><mi>φ</mi><mrow><mtext>cal</mtext></mrow></msup><msup><mi>φ</mi><mrow><mtext>input</mtext></mrow></msup></mrow></msub></mrow></math> 1 1 0.9862 1 1 0.8123 1 1 0.7467 1 1 0.8815 1 1   1 1 <math id="M104"><mrow><mi>R</mi><mi>M</mi><msub><mi>E</mi><mrow><msup><mi>φ</mi><mrow><mtext>cal</mtext></mrow></msup><msup><mi>φ</mi><mrow><mtext>input</mtext></mrow></msup></mrow></msub></mrow></math> 1 1 1.6720 1 1 2.1514 1 1 3.7010 1 1 0.2062 1 1 MLEM 1 1 <math id="M105"><mrow><msub><mi>r</mi><mrow><msup><mi>φ</mi><mrow><mtext>cal</mtext></mrow></msup><msup><mi>φ</mi><mrow><mtext>input</mtext></mrow></msup></mrow></msub></mrow></math> 1 1 0.9877 1 1 0.9481 1 1 0.7447 1 1 0.8655 1 1   1 1 <math id="M106"><mrow><mi>R</mi><mi>M</mi><msub><mi>E</mi><mrow><msup><mi>φ</mi><mrow><mtext>cal</mtext></mrow></msup><msup><mi>φ</mi><mrow><mtext>input</mtext></mrow></msup></mrow></msub></mrow></math> 1 1 0.6627 1 1 1.0882 1 1 3.0150 1 1 0.1978 1 1 MAXED 1 1 <math id="M107"><mrow><msub><mi>r</mi><mrow><msup><mi>φ</mi><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>c</mi><mi>a</mi><mi>l</mi></mstyle></mrow></msup><msup><mi>φ</mi><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>i</mi><mi>n</mi><mi>p</mi><mi>u</mi><mi>t</mi></mstyle></mrow></msup></mrow></msub></mrow></math> 1 1 0.9917 1 1 0.9037 1 1 0.9753 1 1 0.9914 1 1   1 1 <math id="M108"><mrow><mi>R</mi><mi>M</mi><msub><mi>E</mi><mrow><msup><mi>φ</mi><mrow><mtext>cal</mtext></mrow></msup><msup><mi>φ</mi><mrow><mtext>input</mtext></mrow></msup></mrow></msub></mrow></math> 1 1 0.1305 1 1 0.3528 1 1 0.2035 1 1 0.0719 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 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 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. Experiment 1 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 27 ] . As Fig.5(a) Fig.5(a) , the Bonner Spheres Spectrometer consists of one SP9 3 He proportional counter and ten polyethylene spheres with diameters of 1, 2.5, 3, 3.5, 4, 5, 6, 8, 10, 12 inches. The 3 He 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 –10 and 18.84 MeV, which is divided into 207 energy groups. The response function of BSS is shown in Fig.5(b) Fig.5(b) . The BSS of the China Institute of Atomic Energy is used to measure the reference 252 Cf neutron source, and the count rates of BSS <math id="M109"><mrow><mi>N</mi><mo>=</mo><msup><mrow><mrow><mo>(</mo><mrow><msub><mi>N</mi><mn>1</mn></msub><mo>,</mo><msub><mi>N</mi><mn>2</mn></msub><mo>,</mo><mo>…</mo><mo>,</mo><msub><mi>N</mi><mrow><mn>11</mn></mrow></msub></mrow><mo>)</mo></mrow></mrow><mi>T</mi></msup></mrow></math> is obtained and used to unfold the neutron spectrum by the PSO-MLEM algorithm. The reference spectrum and unfolded spectrum by the PSO-MLEM algorithm are shown in in Fig. 6 Fig. 6 . The Table 2 Table 2 is 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 <math id="M113"><mrow><msup><mi>φ</mi><mrow><mtext>cal</mtext></mrow></msup></mrow></math> is more than the 0.99. The relative mean error of the unfolded spectrum is 0.1718%. 1 1   2 1 Correlation coefficient ( r ) and relative mean error ( RME ) 1 1 252 Cf neutron source 1 1 <math id="M114"><mrow><msub><mi>r</mi><mrow><msup><mi>φ</mi><mrow><mtext>cal</mtext></mrow></msup><msup><mi>φ</mi><mrow><mtext>ref</mtext></mrow></msup></mrow></msub></mrow></math> 1 1 <math id="M115"><mrow><mi>M</mi><mi>R</mi><msub><mi>E</mi><mrow><msup><mi>φ</mi><mrow><mtext>cal</mtext></mrow></msup><msup><mi>φ</mi><mrow><mtext>ref</mtext></mrow></msup></mrow></msub></mrow></math> 1 1   1 1 0.9980 1 1 0.1718 Conclusion 1 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 fluctuate 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 of 252 Cf neutron source, and the result shows that the PSO-MLEM algorithm is competent for neutron spectrum unfolding. References G. Medkour Ishak-Boushaki , M. Allab , Comparison between standard unfolding and Bayesian methods in Bonner spheres neutron spectrometry . Nucl. Instrum. Meth. A. 689 , 35 - 39 ( 2012 ). doi: 10.1016/j.nima.2012.06.016 http://doi.org/10.1016/j.nima.2012.06.016 Comparison between standard unfolding and Bayesian methods in Bonner spheres neutron spectrometry H. Song , Y. Yuan , T. Peng , et al. , Optimization study on neutron spectrum unfolding based on the least-squares method . Nucl. Sci. Tech. 29 , 118 ( 2018 ). doi: 10.1007/s41365-018-0454-5 http://doi.org/10.1007/s41365-018-0454-5 Optimization study on neutron spectrum unfolding based on the least-squares method M. Shahmohammadi Beni , D. Krstic , D. Nikezic , Studies on unfolding energy spectra of neutrons using maximum-likelihood expectation–maximization method . Nucl. Sci. Tech. 30 , 134 ( 2019 ). doi: 10.1007/s41365-019-0662-7 http://doi.org/10.1007/s41365-019-0662-7 Studies on unfolding energy spectra of neutrons using maximum-likelihood expectation–maximization method M. Reginatto ,  P. Goldhagen ,  MAXED, a computer code for maximum entropy deconvolution of multisphere neutron spectrometer data .  Health Phys.  77 ,  579 - 583 ( 1999 ).  doi: 10.1097/00004032-199911000-00012 http://doi.org/10.1097/00004032-199911000-00012  MAXED, a computer code for maximum entropy deconvolution of multisphere neutron spectrometer data F.G. Perey ,  Least-squares dosimetry unfolding: The program STAY'SL , ( Oak Ridge National Laboratory/TM-6062 ,  1977 ). doi: 10.2172/5443582 http://doi.org/10.2172/5443582 Least-squares dosimetry unfolding: The program STAY'SL H.W. Engl , Discrepancy principles for Tikhonov regularization of ill-posed problems leading to optimal convergence rates . J. Optimiz. Theory App. 52 , 209 - 215 ( 1987 ). doi: 10.1007/BF00941281 http://doi.org/10.1007/BF00941281 Discrepancy principles for Tikhonov regularization of ill-posed problems leading to optimal convergence rates M. Reginatto ,  P. Goldhagen ,  S. Neumann ,  Spectrum unfolding, sensitivity analysis and propagation of uncertainties with the maximum entropy deconvolution code MAXED .  Nucl. Instrum. Meth. A.  476 , 242 - 246 ( 2002 ).  doi: 10.1016/S0168-9002(01)01439-5 http://doi.org/10.1016/S0168-9002(01)01439-5  Spectrum unfolding, sensitivity analysis and propagation of uncertainties with the maximum entropy deconvolution code MAXED M. Reginatto , A. Zimbal , Bayesian and maximum entropy methods for fusion diagnostic measurements with compact neutron spectrometers . Rev. Sci. Instrum. 79 , 23505 ( 2008 ). doi: 10.1063/1.2841695 http://doi.org/10.1063/1.2841695 Bayesian and maximum entropy methods for fusion diagnostic measurements with compact neutron spectrometers J. Wang , Y. Zhou , Z. Guo , et al. , Neutron spectrum unfolding using three artificial intelligence optimization methods . Appl. Radiat. Isot. 147 , 136 - 143 ( 2019 ). doi: 10.1016/j.apradiso.2019.03.009 http://doi.org/10.1016/j.apradiso.2019.03.009 Neutron spectrum unfolding using three artificial intelligence optimization methods D.W. Freeman ,  D.R. Edwards ,  A.E. Bolon ,  Genetic algorithms - a new technique for solving the neutron spectrum unfolding problem .  Nucl. Instrum. Meth. A.  425 ,  549 - 576 ( 1999 ).  doi: 10.1016/S0168-9002(98)01427-2 http://doi.org/10.1016/S0168-9002(98)01427-2  Genetic algorithms - a new technique for solving the neutron spectrum unfolding problem A. Sharghi Ido , M.R. Bonyadi , G.R. Etaati , et al. ,  Unfolding the neutron spectrum of a NE213 scintillator using artificial neural networks .  Appl. Radiat. Isotopes.  67 , 1912 - 1918 ( 2009 ).  doi: 10.1016/j.apradiso.2009.05.020 http://doi.org/10.1016/j.apradiso.2009.05.020  Unfolding the neutron spectrum of a NE213 scintillator using artificial neural networks R. Eberhart , J. Kennedy , A new optimizer using particle swarm theory . Proceedings of the sixth international symposium on micro machine and human science . IEEE 39 - 43 ( 1995 ). A new optimizer using particle swarm theory H. Wang , W. Wang , Z. Wu , Particle swarm optimization with adaptive mutation for multimodal optimization . Appl. Math. Comput. 221 , 296 - 305 ( 2013 ). doi: 10.1016/j.amc.2013.06.074 http://doi.org/10.1016/j.amc.2013.06.074 Particle swarm optimization with adaptive mutation for multimodal optimization H. Shahabinejad , M. Sohrabpour , A novel neutron energy spectrum unfolding code using particle swarm optimization . Radiat. Phys. Chem. 136 , 9 - 16 ( 2017 ). doi: 10.1016/j.radphyschem.2017.03.033 http://doi.org/10.1016/j.radphyschem.2017.03.033 A novel neutron energy spectrum unfolding code using particle swarm optimization Y. Shi , R.C. Eberhart , Parameter selection in particle swarm optimization . International conference on evolutionary programming, Springer , 591 – 600 ( 1998 ). Parameter selection in particle swarm optimization L. Zhang , Y. Tang , C. Hua , et al. , A new particle swarm optimization algorithm with adaptive inertia weight based on Bayesian techniques . Appl. Soft Comput. 28 , 138 - 149 ( 2015 ). doi: 10.1016/j.asoc.2014.11.018 http://doi.org/10.1016/j.asoc.2014.11.018 A new particle swarm optimization algorithm with adaptive inertia weight based on Bayesian techniques I.K. Gupta , A. Choubey , S. Choubey , Particle swarm optimization with selective multiple inertia weights . 2017 8th International Conference on Computing, Communication and Networking Technologies (ICCCNT) . IEEE Computer Society pp. 1 - 6 ( 2017 ). doi: 10.1109/icccnt.2017.8204132 http://doi.org/10.1109/icccnt.2017.8204132 Particle swarm optimization with selective multiple inertia weights M.A. Memon , M.D. Siddique , S. Mekhilef et al. , Asynchronous Particle Swarm Optimization-Genetic Algorithm (APSO-GA) Based Selective Harmonic Elimination in a Cascaded H-Bridge Multilevel Inverter . IEEE T. Ind. Electron. 69 , 1477 - 1487 ( 2022 ). doi: 10.1109/TIE.2021.3060645 http://doi.org/10.1109/TIE.2021.3060645 Asynchronous Particle Swarm Optimization-Genetic Algorithm (APSO-GA) Based Selective Harmonic Elimination in a Cascaded H-Bridge Multilevel Inverter J. Li , L. Zhang , X. Gu , A composite particle swarm optimization algorithm for hospital equipment management risk control optimization and prediction . J. Env. Public Health 2022 , 1 - 9 ( 2022 ). doi: 10.1155/2022/5268887 http://doi.org/10.1155/2022/5268887 A composite particle swarm optimization algorithm for hospital equipment management risk control optimization and prediction T. Belgya , L. Szentmiklósi , Monte-Carlo calculated detector response functions to unfold radiative neutron capture spectra . Nucl. Instrum. Meth. A. 991 , 165018 ( 2021 ). doi: 10.1016/j.nima.2021.165018 http://doi.org/10.1016/j.nima.2021.165018 Monte-Carlo calculated detector response functions to unfold radiative neutron capture spectra V. Suman , P.K. Sarkar , Neutron spectrum unfolding using genetic algorithm in a Monte Carlo simulation . Nucl. Instrum. Meth. A. 737 , 76 - 86 ( 2014 ). doi: 10.1016/j.nima.2013.11.012 http://doi.org/10.1016/j.nima.2013.11.012 Neutron spectrum unfolding using genetic algorithm in a Monte Carlo simulation S. Maeda , H. Tomita , J. Kawarabayashi et al. , Fundamental study on neutron spectrum unfolding using maximum entropy and maximum likelihood method . Prog. Nucl. Sci. Technol. 1 , 233 - 236 ( 2011 ). doi: 10.15669/pnst.1.233 http://doi.org/10.15669/pnst.1.233 Fundamental study on neutron spectrum unfolding using maximum entropy and maximum likelihood method International Atomic Energy Agency , Compendium of neutron spectra and detector responses for radiation protection purposes . Technical Reports Series No.403 ( 2001 ). Compendium of neutron spectra and detector responses for radiation protection purposes R. Li , J.B. Yang , X.G. Tuo et al. , Unfolding neutron spectra from water-pumping-injection multilayered concentric sphere neutron spectrometer using self-adaptive differential evolution algorithm . Nucl. Sci. Tech. 32 , 26 ( 2021 ). doi: 10.1007/s41365-021-00864-4 http://doi.org/10.1007/s41365-021-00864-4 Unfolding neutron spectra from water-pumping-injection multilayered concentric sphere neutron spectrometer using self-adaptive differential evolution algorithm B. Wiegel , A.V. Alevra , B. Siebert , Calculations of the response functions of Bonner spheres with a spherical 3 He proportional counter using a realistic detector model . SCAN-9503183 , ( 1994 ). Calculations of the response functions of Bonner spheres with a spherical 3He proportional counter using a realistic detector model X. Jiang , J. Cao , C. Jiang , et al. , Geant4 simulation of multi-sphere spectrometer response function and the detection of 241 Am–Be neutron spectrum . Nucl. Sci. Tech. 28 , 174 ( 2017 ). doi: 10.1007/s41365-017-0328-2 http://doi.org/10.1007/s41365-017-0328-2 Geant4 simulation of multi-sphere spectrometer response function and the detection of 241Am–Be neutron spectrum J. Chen , C.J. LI , W. Li et al. , Neutron spectra measurement of IHNI-I BNCT beam with multi- sphere spectrometer . Atomic Energy Science and Technology 48 , 2127 - 2132 ( 2014 ). doi: 10.7538/yzk.2014.48.11.2127 http://doi.org/10.7538/yzk.2014.48.11.2127 ( in Chinese ) Neutron spectra measurement of IHNI-I BNCT beam with multi- sphere spectrometer All 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 and Jian Zhang. The first draft of the manuscript was written by Hong-Fei Xiao, Qing-Xian Zhang and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript. 10.1007/s41365-023-01200-8 This work was supported by the National Natural science Foundation of China (No. 42127807), the Sichuan Science and Technology Program (No. 2020YJ0334), and the Sichuan Science and Technology Breeding Program (No. 2022041). 2022-10-18 2023-02-18 2023-02-25 2023-04-26 2023-06-07 © China Science Publishing & Media Ltd. (Science Press), Shanghai Institute of Applied Physics, the Chinese Academy of Sciences, ChineseNuclear Society and Springer Nature Singapore Pte Ltd. 2023 This is a post-peer-review, pre-copyedit version of an article published in Nuclear Science and Techniques. The final authenticated version is available online at: http://dx.doi.org/10.1007/s41365-023-01200-8 http://dx.doi.org/10.1007/s41365-023-01200-8 . Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law. 2023 38408226 The data that support the findings of this study are openly available in Science Data Bank at https://doi.org/10.57760/sciencedb.07741 and https://cstr.cn/31253.11.sciencedb.07741. Data Availability Statement Nuclear Science and Techniques 4 34