1 Introduction

The inverse transport problem, corresponding to the problem that acquires the unknown information of an object by analyzing the detected radiation fields, is always a strong challenge in nonproliferation, arms control verification, and international security [1, 2]. Thus, great efforts have been dedicated to investigate the inverse transport problem for various application circumstances since the 1960’s. It is found that the inverse transport problems can be solved through an optimization approach by searching the unknown system parameters that minimize a cost difference between the quantities of interest measured from a system and the quantities of interest calculated with a direct transport model using a set of postulated parameters [2].

Several optimization methods have been applied to solve the inverse transport problems, including the gradient-based techniques, e.g., Levenberg-Marquardt (LM) method [3, 4] and Schwinger method [5], and the derivative-free techniques, e.g., differential evolution (DE) method [6, 7] and mesh adaptive direct search (MADs) [8]. Based on gradient-based techniques, a few inverse radiation transport models have been developed to identify the unknown information of one-dimensional source/shield systems, e.g., the locations of interfaces between material layers, source composition, shield material identification, and material mass density, as well as the combinations of these unknown properties [3-6, 9-11]. It is found that the gradient-based inverse models can be applied to solve the inverse problems containing relatively few unknowns (approximately smaller than 4) with high success as long as the proper initial guesses for the unknown parameters are used. However, in the problems with several unknown parameters, the gradient-based inverse model is heavily dependent on the accuracy of the initial guess for the unknown parameters and can often fall into the local minima when random initial guesses are used (corresponding to the cases in which there were not enough prior information about the unknown parameter values) [10]. The requirement of relatively accurate initial guesses of unknowns can also restrict the application of such gradient-based inverse model to a certain extent.

On the other hand, with the development of the computer and optimization algorithm, several inverse radiation models based on derivative-free techniques have been developed to solve the inverse problem in both spherical and cylindrical geometries [6-8, 12]. It is found that the derivative-free inverse solvers can be applied for the source systems containing unknowns of more than 4 or having more complicated geometric properties with a much higher success rate than the ones based on gradient-based techniques. Furthermore, such derivative-free inverse calculations usually do not depend on the initial guesses of the unknowns [8, 13], indicating that they are more proper for the inverse problems with more unknowns or less prior information. However, in previous studies, it is also pointed out that further investigations of algorithm improvement, exploration of technical feasibility and validation, as well as experimental validation and testing for the derivative-free inverse radiation transport models are always needed.

Recently, based on an enhanced differential evolution algorithm with global and local neighborhoods (DEGL) [14, 15], which is suggested to perform better than the standard differential evolution algorithm adopted in previous inverse solver development [6-8, 12], an inverse radiation transport model (IRT-DEGL) is developed for estimating the unknown shielding layer thicknesses [16]. By conducting a series of experiments with a multi-shield point source ¹⁵²Eu and a 60% relative efficiency HPGe spectrometer, the numerical tests are carried out for IRT-DEGL model. It is found that in comparison with the traditional gamma-ray absorption method the IRT-DEGL model can provide a much more accurate shielding layer thicknesses and may have great potential for complicated inverse problems.

It is known that one of the major advantages for inverse solvers based on the differential evolution algorithm is that their inverse results in principle do not depend on the initial guesses of unknowns, thus the present IRT-DEGL model is extended to estimating the unknown layer thicknesses with random initial guesses of unknowns, and the influence on the iteration process and the calculation amount is discussed in detail. Furthermore, in real applications supporting nonproliferation, arms control verification, and international security, one often has to deal with the inverse problems needed to identify the material information or part of the geometric and material information that the traditional analytical algorithm and the basic iteration method can not provide these accuracy results easily. As an elementary attempt, the IRT-DEGL model is further attended to estimate the unknown material mass densities as well as the combination of unknown layer thickness and mass densities of multi-shielding layers. Using the detected gamma-ray spectra, the numerical tests are implemented and the main influence factors for calculation results are analyzed.

2 Theoretical framework of IRT-DEGL model

Inverse procedure of IRT-DEGL model starts with the construction of a direct transport model according to the real detection circumstance and usually containing several unknown parameters of source/shield system to be determined. In order to complicate the direct transport model, a set of vectors, vi (i = 1,..., P), where each vector contains a set of postulated values of the unknown parameters, are generated. The number of unknown parameters determines the dimension of each vector, ui, and P represents the total number of the vectors, i.e., the population size. The fitness of each population vector is determined by a cost function, f(u), between the measured gamma-ray spectrum and the calculated value using the parameters of the population vector, i.e.,

where Rg is the measured data, Rg(ui) is the calculated value using postulated parameter set ui, and the sum over g runs over the full gamma-ray spectrum or the strong emission lines. In an inverse problem, one needs to seek to find the population vector with the global minimum, f(ui).

The IRT-DEGL model employs a generational process with mutation, crossover, and selection operations for optimization. The population vectors move from generation h to generation h+1 by first creating trial vectors with mutation operation. One trial vector vi (i = 1,..., P) is created for each member of the population according to

with the global neighborhood mutation, $g_i^{(h)}$, and local neighborhood mutation, $l_i^{(h)}$

$u_{\text{best}\mathrm{,}i}^{(h)}$ represents the best vector in the entire population at generation h and random integral numbers a, b ∈ [1, P] with a ≠ b ≠ i. $u_{\text{pbest}\mathrm{,}i}^{(h)}$ represents the best vector in the neighborhood of $u_i^{(h)}$ and random integral numbers c, d ∈ [i-k, i+k] with c ≠ d ≠ i, where the neighborhood of radius k is a nonzero integer from 0 to (P-1)/2. α and β are the scaling factors. In comparison with Ref. [16], the weight factor, ω, controlling the balance between the exploration and exploitation capability adopted here is calculated in a self-adaptive way [14],

where ${\omega }_{\text{best}\mathrm{,}i}^{(h)}$ is the weight factor associated with the best vector $u_{\text{best}\mathrm{,}i}^{(h)}$, and the value of ${{\omega }^{\prime }}_i^{(h)}$ is restricted in the range [0.05,0.95].

After P trial vectors are created, a crossover operation is used to create the child vectors vi (i = 1,..., P) with the following rule,

where rj is a random number in [0,1], and CR is a real number in [0,1] denoting the crossover probability and remains the constant throughout the optimization process. i and j denote the i’th child or trial or parent vector and j’th unknown parameter, respectively. In such a way, the i’th child vector could have characteristics of the i’th trial vector and the i’th parent.

After P child vector are created, the selection operation is implemented to choose the vectors in the next generation by a direct competition between the i’th parent vectors and the i’th child vectors,

that is, the i’th population vector of generation h+1 is chosen as the better fit (lower cost function value) between $y_i^{(h)}$ and $u_i^{(h)}$. Furthermore, according to the self-adaptation scheme, the weight factor, ω, in Eq. (4) for the next generation is updated following the selection operation in Eq. (6), that is, ${\omega }_i^{(h+1)}$ equals ${{\omega }^{\prime }}_i^{(h)}$ for $f(y_i^{(h)})\le f(u_i^{(h)})$ while ${\omega }_i^{(h+1)}$ equals ${\omega }_i^{(h)}$ for other cases.

The population of the first generation are created according to the initial guesses set by hand or randomly. Then the generational process proceeds until the converge criterion is achieved, e.g., either the cost function, f(u), of the best population vector is less than a set converge value or until the population has little diversity in the unknown parameter values. And finally, the best population vector is thought to be the inverse results.

3 Results and discussion

It is known that one of the advantages for the derivative-free inverse transport model is that its calculation results do not depend on the initial guesses of unknowns, therefore, one can, in principle, create the initial values randomly or just according to the few pieces of prior information. Although the proper initial guesses from either enough prior information or numerical results from the traditional analytical method may help the inverse calculation converge, sometimes one can not get enough information in one detection, or the source system is so complicated that the analytical method is not suit or can not be solved easily. Thus, for more general inverse cases, it is necessary to directly generate the initial values for unknown parameters randomly according to the prior information already obtained instead of setting an accurate initial guess. Therefore, in comparison with the inverse investigations in Ref. [16] where the relative proper initial guesses of the unknown layer thicknesses are input to create the population of the first generation, the IRT-DEGL model in this paper generates the first generation randomly just under the constriction of prior information.

Following the typical experiment conducted in Ref. [16], where two shielding layers of copper and aluminium with thicknesses of 0.308 cm and 0.57 cm, respectively, locates between a point source of ¹⁵²Eu and a 60% relative efficiency HPGe spectrometer, the gamma-ray spectrum detected for the inverse study is shown in Fig. 1. Considering the cases where the shielding layers are encapsulated as a whole, it is not easy to distinguish the layer thicknesses. The thicknesses of copper layer and aluminium layer are both assumed to be unknown. The densities of the copper layer and aluminium layer are 8.96 g/cm³ and 2.7 g/cm³, respectively, and the material in the layer is assumed to be homogeneous. The two layers are close to each other and the copper layer is on the side close to the HPGe detector. The distance between the radiation source and the HPGe detector and the distance between the HPGe detector and one surface of copper layer are already known. Referring to Ref. [16], the 7 strongest emission lines of ¹⁵²Eu (122, 344, 779, 964, 1086, 1112, and 1408 keV) of detected gamma-ray spectrum are used for inverse calculation, and the simulated gamma-ray spectrum during the inverse calculation is provided by a Monte Carlo code MCNP5 [17]. For all inverse investigations in the present work, the statistical errors for both experimental and simulated full-energy peak areas for these emission lines are smaller than 1%.

Figure 1:Full-screen View

(Color online) Gamma-ray spectrum detected for the inverse study.

The IRT-DEGL calculation converges after 27 generation evolution and the final optimal thicknesses of the Cu layer (0.319 cm) and Al layer (0.551 cm) are obtained. The normalized counts calculated with the inverse result for the 7 emission lines are shown in comparison with the experimental values in table 1 as well as the corresponding cost difference, f(ui). The calculated and experimental normalized counts are the same up to 0.01 with the relativistic difference smaller than 2.126%, which indicates that adopting the final inverse results of the strength of the 7 emission lines of the simulated gamma-ray spectrum can reproduce the detected values well.

In Table 2, the thicknesses of the Cu and Al layers calculated by the IRT-DEGL model with and without inputting initial guesses are shown in comparison with the actual values. One can see that both IRT-DEGL calculations converge quickly and provide the inverse results within the relative difference from the actual value smaller than 6% and 4%, respectively. The major reason for the deviation between two results comes from the convergence criterion that the change of cost function, f(ui), is smaller than 0.1% for 5 generations. When an inverse calculation converges to the area around global minimum and the cost function, f(u), varies slightly in this area, such criterion can stop the inverse calculation in a few generations and provide a value close to the global minimum. Therefore, using such convergence criterion, both the IRT-DEGL calculations with and without the initial guesses could not provide the correct global minimum value, but converge to the global minimum area along the different pathes and finally obtain the different results around the global minimum point. Although this criterion is rough and can only provide an approximate global minimum value, for most applications these inverse result with such degree of accuracy is acceptable. Of course, if one needs to obtain the inverse result closer to the global minimum, the more precision converge criterion should be to adopt, e.g., using criterion where the cost function, f(u_opt.), of the best population vector is less than a small converge value, or allowing the cost function, f(u), to change a little for large generation steps, or using the combination of them.

Furthermore, in Table 2 one can also see that with proper initial guesses, i.e. 0.508 cm for copper layer (corresponding to the relative difference of 64.9%) and 0.700 cm for aluminium layer (corresponding to the relative difference of 22.8%), the IRT-DEGL calculation converges after 19 generation evolutions, while for the case without an initial guess, where the initial population is generated randomly with the prior information introduced above, the inverse calculation converges after 27 generation evolutions and provides relatively accurate results as the inverse values obtained with proper initial guesses. It indicates that the IRT-DEGL model can be applied for the inverse problems with and without enough prior information, and the computing time of the inverse calculation can be saved greatly when the proper initial guesses are used.

In addition to identifing the unknown shielding layer thicknesses, the IRT-DEGL model in the present work is further applied for the inverse problems of estimating the unknown matter mass densities of shielding layers. The detected gamma-ray spectrum for identifying unknown shielding layer thicknesses are used for the present inverse study of the densities. Assuming that the geometric parameters of the source system are already known, three test cases for mass densities of shielding layers, including only copper density (8.96 g/cm³) or aluminium density (2.7 g/cm³) unknowns, as well as both copper and aluminium densities unknown, are considered.

The IRT-DEGL calculations are performed for the three inverse cases with the detected full energy peak areas of the 7 strongest emission lines of ¹⁵²Eu (122, 344, 779, 964, 1086, 1112, and 1408 keV). The inverse results are shown in Table 3 with the actual values for comparison. It can be seen that for three test cases the IRT-DEGL calculations can provide the relatively accurate results with the differences from the actual values smaller than 9%. The deviation between the copper and aluminium densities calculated in different cases mainly comes from the definition of convergence criterion as discussed above, while the deviation between the inverse results and actual values is mainly induced by statistical error, extracting error for the experimental full energy peak area, theoretical simulation of the HPGe detector and detection process, and the lack of enough spectrum characteristics to be analyzed in inverse calculation. Of course, it should be pointed out that the discussion about the mass density of the shielding material only having one component is insignificant, however, the density is related to the material composition and therefore, the inverse investigation for the density of the shielding material containing several components can help to determine the weight fractions of each material component.

Besides the cases where only geometric information or material information is unknown, in the applications that support nonproliferation, arms control verification, and international security, sometimes one has to deal with the inverse problem that part of the geometric and material information are needed to be identified simultaneously. Again, adopting the detected gamma-ray spectrum shown in Fig. 1 and assuming the cases where the thicknesses of both aluminium and copper layers as well as the copper density are unknown, the inverse calculations are performed with IRT-DEGL model.

The optimal thicknesses for the aluminium layer and copper layer as well as the optimal copper density estimated by the IRT-DEGL model are shown as a function of generation in Fig. 2. The actual values are represented as dotted lines in each panel. Following the definition in Ref. [16], the 0 generation is used to denote the initial guesses for the unknown parameters, and since the IRT-DEGL calculations in the present work generate the first generation randomly instead of inputting the initial values, the inverse results in each panel of Fig. 2 start from 1 generation.

Figure 2:Full-screen View

(Color online) The optimal aluminium layer thickness (a), optimal copper layer thickness (b), optimal copper density (c), estimated by the IRT-DEGL model as a function of generation. The dotted line denotes actual thicknesses of copper and aluminium layers in panel (a) and (b), and the actual copper density in panel (c).

In the figure, one can see that the estimated aluminium and copper layer thicknesses and copper density approach the actual values gradually. After 10 generations, all calculated values are located in the area around the actual values with relative difference from actual values smaller than 10%, e.g., 8.42% for aluminium layer thickness, 4.33% for copper layer thickness and 6.46% for copper density, indicating the IRT-DEGL model can search the correct value area rapidly. For the applications that the rough values are required, the inverse calculation can be stopped manually and provide the estimating values. The whole inverse calculation stops after 67 generations of evolution, because the value of the cost function, f(u_opt.), is less than the converged value. As shown in Table 4, the final optimal thicknesses of the aluminium layer and copper layer are 0.314 cm and 0.551 cm corresponding to the relative difference from the actual value 1.95% and 3.33%, and the final optimal copper density is 9.090 g/cm³ corresponding to the relative difference of 1.45%.

4 Summary and conclusion

In summary, the derivative-free inverse radiation transport model based on the differential evolution algorithm with global and local neighborhoods (IRT-DEGL) developed recently has been further extended for the inverse problems of identifying the unknown thicknesses with random initial guesses and material mass densities of multi-shielding layers as well as their combination. Using the detected gamma-ray spectrum, the illustration calculations are implemented. It is found that when using either random or proper initial guesses the IRT-DEGL model can provide the same shielding layer thicknesses, but the computing time of the inverse calculation can be saved greatly with proper initial values. For identifying the unknown mass densities and the combination of unknown layer thickness and mass densities of multi-shielding layers, the IRT-DEGL model can provide the relatively accurate values within 9% and 5%, respectively, of their actual values.

Furthermore, as discussed above, the IRT-DEGL model solves unknown parameters of the source system by iteratively adjusting the hypothetical direct transport model to match the measured spectrum, thus the variance for solving different inverse problems is only correlated with the construction of the direct transport model while the framework and calculation process of the solver are identical, and therefore, the framework of IRT-DEGL model is general in its applications and thus it in principle can be easily extended for the complicated source systems. However, it should be also pointed out that the application of the IRT-DEGL model so far is limited to the simple source system, e.g., point source with one or two shielding layers. When it is applied to the complicated source systems, e.g., multi-shielded volume source or irregular source systems, more numerical problems, including the construction of the direct transport model and the method to reduce the numerical deviation will have to be dealed with, and therefore, a great deal of work will need to be done in the future.