Numerical solution of a parameter estimation problem arising in Prompt-Gamma Neutron Activation Analysis

Abstract

We study a parameter estimation problem which arises from measurements with QUANTOM, a measurement facility utilizing Prompt-Gamma Neutron Activation Analysis (PGNAA) to analyze elemental masses in 200-l nuclear waste drums. The presence of neutron flux measurements in QUANTOM offers additional information for elemental composition reconstruction, motivating an alternative to the standard iterative approach. Therefore, the measurement is modelled using approximate models for neutron and photon transport. We introduce a mathematical formulation for solving the resulting parameter estimation problem of QUANTOM measurements and employ an adjoint-based method to determine the elemental composition. Uniqueness of solutions in this context is studied, and the results offer novel insights for PGNAA applications.

KEYWORDS:

• Parameter estimation problem
• neutron transport
• uniqueness of solutions
• PDE-constrained optimization
• Prompt-Gamma Neutron Activation Analysis (PGNAA)

MATHEMATICS SUBJECT CLASSIFICATIONS:

• 49M40
• 82D75
• 35R30

1. Introduction

For the disposal of low level (LLW) and intermediate level (ILW) radioactive waste in Germany a final deep geological repository called Konrad is currently expected to go into operation around the year 2030. In the process of disposal of nuclear wastes in Konrad one of the major challenges are the compliance with waste acceptance criteria in order to assure a safe and sustainable storage of the waste, and especially the proof of compliance. Radiological waste acceptance criteria can be verified using passive measurements of the gamma emitting radionuclide inventory of waste drums [1,2] and correlation of the remaining radionuclides. For the waste acceptance criteria regarding the material (elemental) composition no metrological solution is applied as a standard so far requiring a recourse on existing documentation or visual inspection by opening waste drums. To overcome this situation, a measurement device called QUANTOM [3,4] has been developed that is capable of performing a non-destructive analysis of elemental masses of radioactive wastes in 200-l drums. The parameter estimation problem arising in the evaluation of the spatially resolved measurement data from QUANTOM is investigated in this work.

The determination of the elemental composition of radioactive waste drums with QUANTOM is based on Prompt-Gamma Neutron Activation Analysis (PGNAA). Free neutrons are produced by a built-in neutron generator, which is in fact a small linear accelerator where Deuterium is accelerated to a Deuterium target. The resulting neutrons are moderated to thermal energies by a massive surrounding graphite layer and interact with the material. When a single neutron is captured by an atom, it becomes energetically excited. In the process of de-excitation gamma quants at element-specific energies are emitted [5]. This is the measurement signal resulting in a gamma spectrum ranging between 100 keV up to 12 MeV at max. Therefore, peaks in such a PGNAA spectrum can be correlated with an element, even specific isotopes. For the determination of the elemental composition from the gamma spectra the gamma and neutron transport in the complete measurement facility need to be considered, as shown in [6].

Because of the large volume samples under consideration in QUANTOM (typically 200-l-drums in Germany), a spatially resolved measurement is carried out corresponding to a complete surface scan. This is achieved by lifting the drum into 4 different positions and rotating the drum so that it can be measured from 6 different angles for each position. With two gamma detectors, a total of 48 gamma spectra is recorded. The gamma spectra are evaluated jointly based on a model of evaluation which accounts for the different measurement positions and angles [4]. In order to quantitatively evaluate the neutron flux in the measurement chamber, additional neutron detectors are included in QUANTOM. Gamma spectra are recorded using high purity Germanium detectors (HPGe), while the total neutron flux is measured using 3He counter tubes, as shown in Figure 1.

Figure 1. Cut-view of the QUANTOM measurement device (taken from [4]).

The standard use of PGNAA comprises the analysis of small samples at research reactors [5]. In the presence of a strong neutron flux and a small sample size, a constant neutron spectrum can be assumed in the sample. Then the model of evaluation reduces to a simple proportionality of element mass $m_l$ (index l indicating a specific element) and net photopeak count rate $P_{E_{\gamma ,l}}$ which is obtained by evaluating the gamma spectrum. The proportionality constant contains the so-called photopeak-efficiency ${ϵ}_{E_{\gamma ,l}}$ (fraction of gammas which are detected at full energy), ${\sigma }_{E_{\gamma ,l}}$, the so-called partial (n,γ)-cross section, total neutron flux in the sample Φ and the constant fraction of Avogadro constant $N_A$ and molar mass $M_l$: (1) $P_{E_{\gamma ,l}}=\frac{N_A}{M_l}\cdot m_l\cdot {ϵ}_{E_{\gamma ,l}}\cdot {\sigma }_{E_{\gamma ,l}}\cdot \mathrm{\Phi }.$(1) Hence, element masses can be obtained by inverting the proportionality. The computation of neutron flux and the photopeak-efficiency will be discussed in detail in Section 2.

In more complex cases, where the dependence of the neutron flux on the elemental composition can not be neglected, PGNAA evaluations are solved using an iterative formulation. Hereby, the neutron flux for a given elemental composition is computed and, using Equation (1(1) $P_{E_{\gamma ,l}}=\frac{N_A}{M_l}\cdot m_l\cdot {ϵ}_{E_{\gamma ,l}}\cdot {\sigma }_{E_{\gamma ,l}}\cdot \mathrm{\Phi }.$(1) ), the updated elemental composition is determined. The element masses in the $(n+1)$-th iteration are thus computed as follows: (2) $(m_l{)}_{n+1}=\frac{P_{E_{\gamma ,l}}{M}_l}{N_A\cdot {ϵ}_{E_{\gamma ,l}}\cdot {\sigma }_{E_{\gamma ,l}}\cdot \mathrm{\Phi }\left[\right(m_1{)}_n,\dots ,\left(m_L{)}_n\right)]},$(2) where L is the total number of elements. The iterative process terminates when the relative change of compositions between two iterations becomes neglectable. In [7,8] this approach has been studied for the evaluation of one single PGNAA measurement using a HPGe detector. Based on a relative formulation the authors were able to show the existence of a solution under mild assumptions. Nevertheless, the uniqueness of a solution cannot be guaranteed in this setting, although, in accompanying experimental studies ambiguous solutions did not occur. The setup of our paper differs in several aspects from the one considered in [7,8]. The sample volume in QUANTOM is larger than considered in standard literature before, hence spatial distribution of neutron flux and photopeak efficiencies need to be carefully represented in the simulation. Measurements for the neutron flux in QUANTOM present additional information regarding the reconstruction of elemental compositions. Furthermore, the problem formulation of [8], where relative compositions are determined, cannot be applied to the more complex measurement task of QUANTOM. The model of evaluation in this work therefore aims at the determination of absolute quantities.

Hence, we present an alternative to the standard iterative approach which is motivated by the spatially resolved measurement and additional measurement equipment for the neutron flux. We consider a mathematical formulation of the evaluation of QUANTOM measurements as a parameter estimation problem. Based on this formulation, we make use of gradient based methods to solve for the elemental composition. We only consider the parameter estimation problem with respect to neutron transport model, the dependence of photon transport to element composition is much weaker and can therefore be neglected.

This problem is studied in a simplified setting which makes it possible to also address the question of uniqueness of a solution. The material is assumed to be homogeneous, so only one measurement position of the HPGe detector is considered. Both types of detectors, for gamma spectra and total neutron flux, are considered in this setting. Elemental compositions are represented by the vector ρ of all corresponding elemental densities ${\rho }_l$, which are parameters of the neutron transport equation and therefore a natural candidate for the reconstruction. In this setting we are able to prove uniqueness of solutions to the parameter estimation problem. To our knowledge, this is the first mathematically rigorous result for PGNAA applications. Throughout this work we restrict the numerical setting to two space dimensions. This is sufficient to study and present all major effects arising in the problem setup. The problem setup is in parts comparable to Diffuse Optical Tomography (DOT), as a similar governing equation is used as an approximate model for particle transport. The connection between the two problems are discussed in the conclusive part of this paper.

This paper is structured as follows. In Section 2 we introduce the SP₁-approximation of the neutron transport equation and model functions for the detector response which, combined, represent the forward model. Section 3 comprises a formulation of the parameter estimation problem for QUANTOM in the form of a PDE-constrained optimization problem and a solution method via a formal Lagrange method. The key results of this work are contained in Section 4, where we prove uniqueness of solutions to the parameter estimation problem. In the following Section 5, we describe the computational implementation of the parameter estimation problem, based on an already existing neutron transport solver. Numerical results which support our findings are presented in Section 6, and we conclude the paper in Section 7 by discussing both the generalizability of our results and remaining open questions.

2. Forward model

As introduced in Section 1, two different kinds of detectors are used in QUANTOM to detect the total neutron flux and the energy resolved emitted gammas, respectively. For both detector types a model for the counting rates is required. Figure 2 shows the abstract measurement setup. Neutrons are constantly emitted by a neutron source. As a result, a steady state neutron flux field is formed in the graphite and the sample, while the graphite moderator provides a thermal neutron spectrum both in the graphite and the sample. It is important to note here, that the neutron flux field depends on the elemental composition of the sample, since the neutrons directly interact with its atoms. Both the sample and graphite are activated by the free neutrons, resulting in element-specific gamma emissions. The detectors are placed in proximity to the sample in order to collect the related measurement signal. Activation of the graphite can be neglected since scattering is dominant in the moderator. Furthermore, detectors can be collimated to the sample to further decrease this effect.

Figure 2. Abstract sketch of the measurement setup, including the neutron source, sample and both types of detectors in a surrounding graphite moderator.

In order to evaluate the counting rates of the detectors, it is necessary to introduce models for the gamma and photon transport first. The transport of gammas and neutrons can be described by the linear Boltzmann equation [9]. For both types of transport we make use of approximate models. For the neutron transport equation, the SP₁-equation, also known as neutron diffusion equation, is used as an approximation. The applicability of the SP₁-equation for this setup was shown in a previous publication [6]. Note that in this work we apply the one-group SP₁-equation, since we can assume a thermal neutron flux. The equation then reads (3) $-\mathrm{\nabla }\cdot \left(\frac{1}{3{\mathrm{\Sigma }}_t\left(x\right)}\mathrm{\nabla \varphi }\right(x\left)\right)+{\mathrm{\Sigma }}_a\left(x\right)\varphi \left(x\right)=Q\left(x\right).$(3) with the Marshak boundary conditions (4) $\frac{1}{2}\varphi \left(x\right)+\frac{1}{3{\mathrm{\Sigma }}_t\left(x\right)}\frac{\mathrm{\partial }}{\mathrm{\partial }n}\varphi \left(x\right)=0$(4) where $\varphi \left(x\right)$ is the neutron flux and $Q\left(x\right)$ the source term while ${\mathrm{\Sigma }}_t\left(x\right)$ and ${\mathrm{\Sigma }}_a\left(x\right)$ are the total and absorption cross sections of the neutrons. It is important to note that the total and absorption cross sections within the sample can be expressed using the elemental densities ${\rho }_l$ as well as the total and absorption cross sections of the elements, ${\sigma }_{\mathrm{tl}}$ and ${\sigma }_{\mathrm{al}}$, while these cross sections are independent of ρ anywhere else, leading to the representation: (5) $\begin{array}{rl}{\mathrm{\Sigma }}_t\left(x\right)& =\{\begin{array}{ll}{\displaystyle \sum _l{\sigma }_{\mathrm{tl}}{\rho }_l}& x\in V_{\mathrm{sam}}\\ M_t\left(x\right)& \mathrm{else}\end{array}\\ {\mathrm{\Sigma }}_a\left(x\right)& =\{\begin{array}{ll}{\displaystyle \sum _l{\sigma }_{\mathrm{al}}{\rho }_l}& x\in V_{\mathrm{sam}}\\ M_a\left(x\right)& \mathrm{else}\end{array}\end{array}$(5) Here $M_t\left(x\right)$ and $M_a\left(x\right)$ are piecewise constant functions that represent the cross sections in the actual facility. The dependence of the cross sections on ρ in the sample leads to a non-linearity of this parameter estimation problem.

According to Equation (1(1) $P_{E_{\gamma ,l}}=\frac{N_A}{M_l}\cdot m_l\cdot {ϵ}_{E_{\gamma ,l}}\cdot {\sigma }_{E_{\gamma ,l}}\cdot \mathrm{\Phi }.$(1) ), photon transport is relevant for the evaluation of a single peak as it determines the so-called photo-peak efficiency ${ϵ}_{E_{\gamma ,l}}$. This quantity represents the fraction of gammas which are detected at full energy. Consequently, scattered photons do not contribute to ${ϵ}_{E_{\gamma ,l}}$, allowing relatively simple models to calculate photo-peak efficiencies, rather than employing full photon transport simulations. The main effects which have to be considered are absorption, the geometry of the measurement setup and detector physics. Photopeak efficiency is a parameter dependent on gamma energy. Therefore, evaluating a measurement requires either point-wise values for all gamma energies considered or determining a (gamma-)energy-dependent calibration curve for a given measurement facility or setup, accounting for all these effects. This calibration curve provides photopeak efficiencies for any given gamma-energy and can be obtained experimentally, by measuring at given energies using calibration emitters and then interpolating between measured values. Alternatively, calibration curves can be determined using either analytical models or Monte Carlo simulations. As prompt gamma energies range between 50 keV up to 12 MeV, experimental calibration is of limited use, making modeling or simulation-based calibrations preferable. Optimal performance is typically achieved using Monte Carlo simulations with MCNP [10] to account for the intricate physics within a HPGe detector. Further details on the calibration of gamma detectors can be found in [11].

In principle, the photopeak efficiency also depends on the composition of the sample, which changes in every iteration of the evaluation (Equation (2(2) $(m_l{)}_{n+1}=\frac{P_{E_{\gamma ,l}}{M}_l}{N_A\cdot {ϵ}_{E_{\gamma ,l}}\cdot {\sigma }_{E_{\gamma ,l}}\cdot \mathrm{\Phi }\left[\right(m_1{)}_n,\dots ,\left(m_L{)}_n\right)]},$(2) )) (as the current result of the measurement, not the physical composition!). However, the tabulated values in [12] suggest that the mass attenuation coefficients for most of the materials relevant for the application are similar. Therefore, the sample's impact on the photo-peak efficiency is primarily governed by its average density, a known parameter that remains constant during the evaluation. Thus, we assume photo-peak efficiencies to be constant throughout the iterative evaluation.

In this work, we employ a simplistic model for the photopeak efficiencies in order to simplify notations and calculations, but all results remain valid when using experimentally obtained calibration curves, more sophisticated analytical models, or calibration curves based on extensive Monte Carlo simulations. The photopeak efficiency is given by (6) ${ϵ}_{E_{\gamma ,l}}\left(x\right)={ϵ}_{\mathrm{geom},E_{\gamma ,l}}\left(x\right)\cdot {ϵ}_{\mathrm{att},E_{\gamma ,l}}\left(x\right)\cdot {ϵ}_{\det ,E_{\gamma ,l}},$(6) where ${ϵ}_{\mathrm{att},E_{\gamma ,l}}\left(x\right)$ is the photon attenuation, ${ϵ}_{\mathrm{geom},E_{\gamma ,l}}\left(x\right)$ is the geometric efficiency and ${ϵ}_{\det ,E_{\gamma ,l}}$ is the detector efficiency. It should be noted that the photopeak efficiency is dependent of the location x where the photon is emitted. Since there exists no closed analytical expression for the full photopeak efficiency, we use a model that captures its main features. The geometric efficiency in this setup is given by (7) ${ϵ}_{\mathrm{geom},E_{\gamma ,l}}\left(x\right)=\frac{1}{\pi }\arctan \left(\frac{r_{\det }\left(x\right)}{d_{\det }\left(x\right)}\right),$(7) where $r_{\det }\left(x\right)$ is the radius of the detector and $d_{\det }\left(x\right)$ the distance to the center of the detector. The photon attenuation is given by (8) ${ϵ}_{\mathrm{att},E_{\gamma ,l}}\left(x\right)={\mathrm{e}}^{-\mu \cdot d\left(x\right)},$(8) where µ is the attenuation coefficient and d the traversed distance. Hereby the attenuation coefficient is mainly determined by the average density of the material, which can be easily determined in practice. The detector efficiency ${ϵ}_{\det ,E_{\gamma ,l}}$ can be interpreted as the probability for a photon entering the detector to contribute to the measured signal and has to be determined by photon transport simulations. For a given energy and detector, it is a constant. For simplicity it is set to one in the following.

Using both the neutron flux $\varphi \left(x\right)$ and the photopeak efficiency ${ϵ}_{E_{\gamma ,l}}$, we can define the counting rates of the detectors. The neutron detectors measure neutron counting rates given by (9) $R_n\left[\varphi \right]={\int }_{V_{\det ,n}}{\sigma }_d\cdot \varphi \left(x\right)\mathrm{d}x,$(9) where ${\sigma }_d$ is the so-called detector cross section. The integral is carried out over the volume of the nth detector $V_{\det ,n}$.

In contrast to the neutron detectors, the gamma detectors do resolve the energy of the measured photons. The resulting gamma spectrum can then be analyzed for element-specific photopeaks. Assuming a homogeneous sample, the counting rate for photons forming a photopeak at energy $E_{\gamma }$ specific for element l is given by (10) $P_{E_{\gamma ,l}}[\rho ,\varphi ]={\int }_{V_{\mathrm{sam}}}{ϵ}_{E_{\gamma ,l}}{\rho }_l{\sigma }_{E_{\gamma ,l}}\varphi \left(x\right)\mathrm{d}x,$(10) where ${\sigma }_{E_{\gamma ,l}}$ is the so-called partial (n,γ)-cross-section and ${\rho }_l$ is the density of the element l within the sample. Equation (10(10) $P_{E_{\gamma ,l}}[\rho ,\varphi ]={\int }_{V_{\mathrm{sam}}}{ϵ}_{E_{\gamma ,l}}{\rho }_l{\sigma }_{E_{\gamma ,l}}\varphi \left(x\right)\mathrm{d}x,$(10) ) is a generalization of Equation (1(1) $P_{E_{\gamma ,l}}=\frac{N_A}{M_l}\cdot m_l\cdot {ϵ}_{E_{\gamma ,l}}\cdot {\sigma }_{E_{\gamma ,l}}\cdot \mathrm{\Phi }.$(1) ), accounting for the additional spatial variation of the parameters ${ϵ}_{E_{\gamma ,l}}$ and ϕ. The integral in Equation (10(10) $P_{E_{\gamma ,l}}[\rho ,\varphi ]={\int }_{V_{\mathrm{sam}}}{ϵ}_{E_{\gamma ,l}}{\rho }_l{\sigma }_{E_{\gamma ,l}}\varphi \left(x\right)\mathrm{d}x,$(10) ) is carried out over the volume where the gamma photons are produced, the sample volume $V_{\mathrm{sam}}$.

3. Parameter estimation problem

The parameter estimation problem to be considered in this work is the following: given measurements from the neutron and gamma detectors, determine the elemental densities ${\rho }_l$. These are themselves parameters of the neutron transport equation. In the following subsections we present a rigorous formulation and a corresponding solution method.

3.1. Formulation as PDE-constrained optimization problem

The parameter estimation problem can be formulated as a PDE-constrained optimization problem with SP₁-equation as a constraint. Therefore, it can be written as an optimal control problem [13] with the vector of the material densities ρ as the control variable and the neutron flux ϕ as the state variable.

The forward model provides a map from parameter to data space. Given the measurements for the peak counting rates $\overline{P_{E_{\gamma ,l}}}$ and the neutron counting rates $\overline{R_n}$, an objective function including peak counting rates for all elements and the counting rates for all neutron detectors can be defined : (11) $J_{\det }[\rho ,\varphi ]=\sum _{l=1}^L\frac{w_{\gamma }}{{\overline{P_{E_{\gamma ,l}}}}^2}\left(P_{E_{\gamma ,l}}\right[\rho ,\varphi ]-\overline{P_{E_{\gamma ,l}}}{)}^2+\sum _{n=1}^N\frac{w_d}{{\overline{R_n}}^2}(R_n\left[\varphi \right]-\overline{R_n}{)}^2.$(11) Here, $w_{\gamma }$ and $w_d$ are weights that determine the contribution of the different detectors to the objective function, while L is the number of elements within the sample and N the number of neutron detectors. The weight of the neutron detectors is set to $w_d=\frac{1}{2}$, while the weight of the gamma detectors is set to $w_{\gamma }=\frac{3}{4}$ to ensure an equal weight of both terms despite the higher number of neutron detectors. In addition, we include a term punishing deviations from the known total mass $\overline{M}$: (12) $J_m\left[\rho \right]=\frac{w_m}{{\overline{M}}^2}(c_1{\rho }_1+c_2{\rho }_2-\overline{M}{)}^2$(12) Here $c_1$ and $c_2$ are material-specific constants and $w_m=\frac{1}{2}$. Since the total mass $\overline{M}$ can be measured with high accuracy, we assume in the following that the error on the total mass is negligible. This term effectively acts as a regularization. Its effect is investigated in Section 6.2.

The complete objective function is then given by (13) $J[\rho ,\varphi ]=J_{\det }[\rho ,\varphi ]+J_m\left[\rho \right].$(13) The parameter estimation problem can then be formulated as minimizing the objective function with the SP₁-equation as a constraint, so we want to solve the problem (14) $\begin{array}{rl}\underset{\rho }{min}J[\rho ,\varphi ]\mathrm{subjectto}D\left(x\right)\varphi \left(x\right)& =Q\left(x\right)\mathrm{in}D\\ \mathrm{and}\frac{1}{2}\varphi \left(x\right)+\frac{1}{3{\mathrm{\Sigma }}_t\left(x\right)}\frac{\mathrm{\partial }}{\mathrm{\partial }n}\varphi \left(x\right)& =0\mathrm{on}\mathrm{\partial D}.\end{array}$(14) The differential operator in Equation (14(14) $\begin{array}{rl}\underset{\rho }{min}J[\rho ,\varphi ]\mathrm{subjectto}D\left(x\right)\varphi \left(x\right)& =Q\left(x\right)\mathrm{in}D\\ \mathrm{and}\frac{1}{2}\varphi \left(x\right)+\frac{1}{3{\mathrm{\Sigma }}_t\left(x\right)}\frac{\mathrm{\partial }}{\mathrm{\partial }n}\varphi \left(x\right)& =0\mathrm{on}\mathrm{\partial D}.\end{array}$(14) ) is defined by (15) $D\left(x\right)\varphi \left(x\right)=-\mathrm{\nabla }\cdot \left(\frac{1}{3{\mathrm{\Sigma }}_t\left(x\right)}\mathrm{\nabla \varphi }\right(x\left)\right)+{\mathrm{\Sigma }}_a\left(x\right)\varphi \left(x\right),$(15) while D is the domain with corresponding boundary $\mathrm{\partial D}$. Outside the sample, ${\mathrm{\Sigma }}_t$ and ${\mathrm{\Sigma }}_a$ are independent of ρ. Thus, the differential operator D depends on ρ on a part of the domain.

3.2. Solution via formal Lagrange method

In order to solve the problem posed in Equation (14(14) $\begin{array}{rl}\underset{\rho }{min}J[\rho ,\varphi ]\mathrm{subjectto}D\left(x\right)\varphi \left(x\right)& =Q\left(x\right)\mathrm{in}D\\ \mathrm{and}\frac{1}{2}\varphi \left(x\right)+\frac{1}{3{\mathrm{\Sigma }}_t\left(x\right)}\frac{\mathrm{\partial }}{\mathrm{\partial }n}\varphi \left(x\right)& =0\mathrm{on}\mathrm{\partial D}.\end{array}$(14) ), we use a formal Lagrange method, from which we derive an expression for the gradient of the objective function, which in turn can be used for a gradient based optimization scheme. In the following, D is the computational domain with the corresponding boundary $\mathrm{\partial D}$.

The Lagrangian for this problem is given by the objective function with a Lagrange multiplier ${\varphi }^{\star }$ and the SP₁-equation as a constraint: (16) $L[\rho ,\varphi ,{\varphi }^{\star }]=J[\rho ,\varphi ]-{\int }_D{\varphi }^{\star }\left(x\right)\left(D\right(x\left)\varphi \right(x)-Q(x\left)\right)\mathrm{d}x.$(16) By formal integration by parts and using Equation (4(4) $\frac{1}{2}\varphi \left(x\right)+\frac{1}{3{\mathrm{\Sigma }}_t\left(x\right)}\frac{\mathrm{\partial }}{\mathrm{\partial }n}\varphi \left(x\right)=0$(4) ), we obtain (17) $\begin{array}{rl}L[\rho ,\varphi ,{\varphi }^{\star }]& =J[\rho ,\varphi ]-{\int }_D\varphi \left(x\right)\left(D\right(x\left){\varphi }^{\star }\right(x\left)\right)\mathrm{d}x-{\int }_{D}Q\left(x\right){\varphi }^{\star }\left(x\right)\mathrm{d}x\\ & +{\int }_{\mathrm{\partial }D}\frac{1}{2}\varphi \left(x\right){\varphi }^{\star }\left(x\right)+\frac{1}{3{\mathrm{\Sigma }}_t\left(x\right)}\varphi \left(x\right)\frac{\mathrm{\partial }}{\mathrm{\partial }n}{\varphi }^{\star }\left(x\right)\mathrm{d}s.\end{array}$(17) Since the operator $D\left(x\right)$ is self-adjoint, i.e. $D\left(x\right)=D^{\star }\left(x\right)$, we can identify ${\varphi }^{\star }$ with the adjoint flux. The boundary terms in Equation (17(17) $\begin{array}{rl}L[\rho ,\varphi ,{\varphi }^{\star }]& =J[\rho ,\varphi ]-{\int }_D\varphi \left(x\right)\left(D\right(x\left){\varphi }^{\star }\right(x\left)\right)\mathrm{d}x-{\int }_{D}Q\left(x\right){\varphi }^{\star }\left(x\right)\mathrm{d}x\\ & +{\int }_{\mathrm{\partial }D}\frac{1}{2}\varphi \left(x\right){\varphi }^{\star }\left(x\right)+\frac{1}{3{\mathrm{\Sigma }}_t\left(x\right)}\varphi \left(x\right)\frac{\mathrm{\partial }}{\mathrm{\partial }n}{\varphi }^{\star }\left(x\right)\mathrm{d}s.\end{array}$(17) ) need to vanish. This is only the case if the relation (18) $\frac{1}{2}{\varphi }^{\star }\left(x\right)+\frac{1}{3{\mathrm{\Sigma }}_t\left(x\right)}\frac{\mathrm{\partial }}{\mathrm{\partial }n}{\varphi }^{\star }\left(x\right)=0$(18) is fulfilled on the boundary. The derivatives of the Lagrangian can be used within an iterative scheme to reconstruct the elemental densities ρ. Differentiation of Equation (16(16) $L[\rho ,\varphi ,{\varphi }^{\star }]=J[\rho ,\varphi ]-{\int }_D{\varphi }^{\star }\left(x\right)\left(D\right(x\left)\varphi \right(x)-Q(x\left)\right)\mathrm{d}x.$(16) ) with respect to the adjoint flux ${\varphi }^{\star }$ yields the weak form of the forward SP₁-equation, while the differentiation with respect to the flux ϕ yields the weak form of the adjoint SP₁-equation (19) $D^{\star }\left(x\right){\varphi }^{\star }\left(x\right)=Q^{\star }\left(x\right)$(19) with the adjoint source (20) $\begin{array}{rl}& Q^{\star }=\frac{\mathrm{\delta J}}{\mathrm{\delta \varphi }}[\rho ,\varphi ]\\ & =\sum _{l=1}^L\frac{2w_{\gamma }}{{\overline{P_{E_{\gamma ,l}}}}^2}\left(P_{E_{\gamma ,l}}\right[\rho ,\varphi ]-\overline{P_{E_{\gamma ,l}}}){\int }_{V_{\mathrm{sam}}}{ϵ}_{E_{\gamma ,l}}\left(x\right){\rho }_l{\sigma }_{E_{\gamma ,l}}\mathrm{d}x\\ & +\sum _{n=1}^N\frac{2w_d}{{\overline{R_n}}^2}\left(R_n\right[\varphi ]-\overline{R_n}){\int }_{V_{\det ,n}}{\sigma }_d\mathrm{d}x.\end{array}$(20) and the boundary condition defined in Equation (18(18) $\frac{1}{2}{\varphi }^{\star }\left(x\right)+\frac{1}{3{\mathrm{\Sigma }}_t\left(x\right)}\frac{\mathrm{\partial }}{\mathrm{\partial }n}{\varphi }^{\star }\left(x\right)=0$(18) ). The derivative of Equation (16(16) $L[\rho ,\varphi ,{\varphi }^{\star }]=J[\rho ,\varphi ]-{\int }_D{\varphi }^{\star }\left(x\right)\left(D\right(x\left)\varphi \right(x)-Q(x\left)\right)\mathrm{d}x.$(16) ) with respect to the density of the lth element is given by (21) $\begin{array}{rl}\frac{\mathrm{\delta L}}{\delta {\rho }_l}[\rho ,\varphi ,{\varphi }^{\star }]& =\frac{{\sigma }_{t,l}}{\sum _j{\sigma }_{t,j}{\rho }_j}{\int }_{V_{\mathrm{sam}}}\mathrm{\nabla \varphi }\cdot \mathrm{\nabla }{\varphi }^{\star }\mathrm{d}x-{\sigma }_{a,l}{\int }_{V_{\mathrm{sam}}}\varphi {\varphi }^{\star }\mathrm{d}x\\ & +\frac{2w_{\gamma }}{{\overline{P_{E_{\gamma ,l}}}}^2}({\int }_{V_{\mathrm{sam}}}{ϵ}_{E_{\gamma ,l}}{\sigma }_{E_{\gamma ,l}}{\rho }_l\varphi \mathrm{d}x-\overline{P_{E_{\gamma ,l}}}){\int }_{V_{\mathrm{sam}}}{ϵ}_{E_{\gamma ,l}}{\sigma }_{E_{\gamma ,l}}\varphi \mathrm{d}x\\ & +\frac{2w_m{c}_l}{{\overline{M}}^2}(c_1{\rho }_1+c_2{\rho }_2-\overline{M}).\end{array}$(21) Using the results for ϕ and ${\varphi }^{\star }$, the gradient of the Lagrangian with respect to the different elemental density ${\rho }_l$ is computed and can be used in a gradient-based optimization scheme.

4. Uniqueness

In this chapter, we investigate the uniqueness of the solution of the unperturbed parameter estimation problem, i.e. without noise in the measurement data. First, we introduce a change of variables to simplify the discussion. We then investigate the objective function to motivate our proof for the uniqueness result that we finally present.

4.1. Change of variables

Within this chapter, we introduce a change of variables [14] to simplify the following discussion: (22) $\psi =\mathrm{\gamma \varphi }\mathrm{with}\gamma =\frac{1}{\sqrt{3{\mathrm{\Sigma }}_t}}.$(22) By substituting Equation (22(22) $\psi =\mathrm{\gamma \varphi }\mathrm{with}\gamma =\frac{1}{\sqrt{3{\mathrm{\Sigma }}_t}}.$(22) ) into (3(3) $-\mathrm{\nabla }\cdot \left(\frac{1}{3{\mathrm{\Sigma }}_t\left(x\right)}\mathrm{\nabla \varphi }\right(x\left)\right)+{\mathrm{\Sigma }}_a\left(x\right)\varphi \left(x\right)=Q\left(x\right).$(3) ) we obtain the Helmholtz-like equation (23) ${\mathrm{\nabla }}^2\psi +\mathrm{\eta \psi }=\frac{Q}{\gamma }\mathrm{with}\eta =\frac{{\mathrm{\nabla }}^2\gamma }{\gamma }+\frac{{\mathrm{\Sigma }}_a}{{\gamma }^2}.$(23) Since the cross sections as defined in Equation (5(5) $\begin{array}{rl}{\mathrm{\Sigma }}_t\left(x\right)& =\{\begin{array}{ll}{\displaystyle \sum _l{\sigma }_{\mathrm{tl}}{\rho }_l}& x\in V_{\mathrm{sam}}\\ M_t\left(x\right)& \mathrm{else}\end{array}\\ {\mathrm{\Sigma }}_a\left(x\right)& =\{\begin{array}{ll}{\displaystyle \sum _l{\sigma }_{\mathrm{al}}{\rho }_l}& x\in V_{\mathrm{sam}}\\ M_a\left(x\right)& \mathrm{else}\end{array}\end{array}$(5) ) are dependent on ρ only within the sample volume $V_{\mathrm{sam}}$, both the source term and the boundary conditions are independent of ρ. The parameter η in Equation (23(23) ${\mathrm{\nabla }}^2\psi +\mathrm{\eta \psi }=\frac{Q}{\gamma }\mathrm{with}\eta =\frac{{\mathrm{\nabla }}^2\gamma }{\gamma }+\frac{{\mathrm{\Sigma }}_a}{{\gamma }^2}.$(23) ) is only varied within $V_{\mathrm{sam}}$ and effectively reduces there to (24) $\eta =\frac{{\mathrm{\Sigma }}_a}{{\gamma }^2}=3{\mathrm{\Sigma }}_a{\mathrm{\Sigma }}_t$(24) due to the homogeneity of the sample. For the investigation of uniqueness, the term punishing deviations from the total mass M is not needed, so we set $w_m$ in Equation (12(12) $J_m\left[\rho \right]=\frac{w_m}{{\overline{M}}^2}(c_1{\rho }_1+c_2{\rho }_2-\overline{M}{)}^2$(12) ) to zero. Using the definition of ψ in Equation (22(22) $\psi =\mathrm{\gamma \varphi }\mathrm{with}\gamma =\frac{1}{\sqrt{3{\mathrm{\Sigma }}_t}}.$(22) ), we then obtain a transformed cost functional (25) $J[\rho ,\psi ]=\sum _{l=1}^L\frac{w_{\gamma }}{{\overline{P_{E_{\gamma ,l}}}}^2}\left(P_{E_{\gamma ,l}}\right[\rho ,\psi ]-\overline{P_{E_{\gamma ,l}}}{)}^2+\sum _{n=1}^N\frac{w_d}{{\overline{R_n}}^2}(R_n\left[\psi \right]-\overline{R_n}{)}^2.$(25) with neutron counting rates given by (26) $R_n\left[\psi \right]={\int }_{V_{\det ,n}}{\sigma }_d\frac{1}{\gamma }\psi \left(x\right)\mathrm{d}x,$(26) and gamma counting rates given by (27) $P_{E_{\gamma ,l}}[\rho ,\psi ]={\int }_{V_{\mathrm{sam}}}{ϵ}_{E_{\gamma ,l}}\left(x\right){\rho }_l{\sigma }_{E_{\gamma ,l}}\frac{1}{\gamma }\psi \left(x\right)\mathrm{d}x.$(27)

4.2. Objective function

Figure 3 shows an exemplary objective function where the cross sections of the elements are set to ${\sigma }_{t1}={\sigma }_{t2}=1$, ${\sigma }_{a1}=0.1$ and ${\sigma }_{a2}=0.5$. The different terms in the objective function defined in Equation (25(25) $J[\rho ,\psi ]=\sum _{l=1}^L\frac{w_{\gamma }}{{\overline{P_{E_{\gamma ,l}}}}^2}\left(P_{E_{\gamma ,l}}\right[\rho ,\psi ]-\overline{P_{E_{\gamma ,l}}}{)}^2+\sum _{n=1}^N\frac{w_d}{{\overline{R_n}}^2}(R_n\left[\psi \right]-\overline{R_n}{)}^2.$(25) ) are displayed separately. Figure 3 (a) and (b) show the objective function including only the first and the second photopeak, respectively. It can be seen that there exists no unique minimum in both figures, but instead a large ‘valley’ and thus many admissible values for ρ.

Figure 3. Subfigures (a)–(c) depict the different terms of the objective function after decomposition, subfigure (d) depicts the full objective function. The analytically computed minimum of the objective function of the neutron detectors is green. A logarithmic scale is used for all plots. (a) Objective function including the first Photopeak. (b) Objective function including the second Photopeak. (c) Objective function including all neutron detectors and (d) Full objective Function.

Since the neutron counting rates, as can be seen from Equations (23(23) ${\mathrm{\nabla }}^2\psi +\mathrm{\eta \psi }=\frac{Q}{\gamma }\mathrm{with}\eta =\frac{{\mathrm{\nabla }}^2\gamma }{\gamma }+\frac{{\mathrm{\Sigma }}_a}{{\gamma }^2}.$(23) ) and (26(26) $R_n\left[\psi \right]={\int }_{V_{\det ,n}}{\sigma }_d\frac{1}{\gamma }\psi \left(x\right)\mathrm{d}x,$(26) ), only depend on the parameter η, they can only be used to reconstruct this single parameter, independent of the number of detectors. Thus additional detectors can possibly reduce the measurement error, but cannot reconstruct additional parameters. Therefore, in Figure 3(c) the sum over the partial objective functions for all detectors is shown. It should be noted that all combinations of ${\mathrm{\Sigma }}_t$ and ${\mathrm{\Sigma }}_a$ that lead to the same parameter η give rise to the same solution ψ, since both the boundary conditions and the source term are independent of ρ. Thus the neutron detectors alone are not sufficient to uniquely reconstruct the material composition of the sample, although they might reconstruct the parameter η. In fact, for this specific problem one can derive an analytic expression that allows to compute the density of one material for a given density of the other material that leads to a desired parameter η. By inserting Equation (5(5) $\begin{array}{rl}{\mathrm{\Sigma }}_t\left(x\right)& =\{\begin{array}{ll}{\displaystyle \sum _l{\sigma }_{\mathrm{tl}}{\rho }_l}& x\in V_{\mathrm{sam}}\\ M_t\left(x\right)& \mathrm{else}\end{array}\\ {\mathrm{\Sigma }}_a\left(x\right)& =\{\begin{array}{ll}{\displaystyle \sum _l{\sigma }_{\mathrm{al}}{\rho }_l}& x\in V_{\mathrm{sam}}\\ M_a\left(x\right)& \mathrm{else}\end{array}\end{array}$(5) ) into (24(24) $\eta =\frac{{\mathrm{\Sigma }}_a}{{\gamma }^2}=3{\mathrm{\Sigma }}_a{\mathrm{\Sigma }}_t$(24) ) for the case $x\in V_{\mathrm{sam}}$, we obtain a relation between ${\rho }_2$ and ${\rho }_1$ (28) ${\rho }_2=\frac{-({\sigma }_{a1}{\sigma }_{t2}+{\sigma }_{a2}{\sigma }_{t1}){\rho }_1+\sqrt{({\sigma }_{a1}{\sigma }_{t2}+{\sigma }_{a2}{\sigma }_{t1}{)}^2{\rho }_1^2-4({\sigma }_{a1}{\sigma }_{t1}{\rho }_1^2-\eta /3){\sigma }_{a2}{\sigma }_{t2}}}{2{\sigma }_{a2}{\sigma }_{t2}}$(28) If the value of η belonging to the minimum is used, Equation (28(28) ${\rho }_2=\frac{-({\sigma }_{a1}{\sigma }_{t2}+{\sigma }_{a2}{\sigma }_{t1}){\rho }_1+\sqrt{({\sigma }_{a1}{\sigma }_{t2}+{\sigma }_{a2}{\sigma }_{t1}{)}^2{\rho }_1^2-4({\sigma }_{a1}{\sigma }_{t1}{\rho }_1^2-\eta /3){\sigma }_{a2}{\sigma }_{t2}}}{2{\sigma }_{a2}{\sigma }_{t2}}$(28) ) can be used to compute the expected position of the minimal values. It is depicted in Figure 3(c) as the green curve. It can be seen that the green curve lies within the valley of the numerical minimal values.

From Figure 3(a)–(c) it is apparent that neither the measurements of the different gamma lines nor the neutron flux measurements alone have a unique minimum. In combination however, we obtain an objective function with a minimum contained within a small region, well suitable for a reconstruction, as depicted in Figure 3(d).

4.3. Uniqueness of the parameter estimation problem

In the following, we prove our main result regarding the uniqueness of a solution to the considered parameter estimation problem. Throughout, we implicitly assume existence of a solution.

Theorem 4.1

Let $J[\rho ,\psi ]$ be the objective function as defined in Equation (25(25) $J[\rho ,\psi ]=\sum _{l=1}^L\frac{w_{\gamma }}{{\overline{P_{E_{\gamma ,l}}}}^2}\left(P_{E_{\gamma ,l}}\right[\rho ,\psi ]-\overline{P_{E_{\gamma ,l}}}{)}^2+\sum _{n=1}^N\frac{w_d}{{\overline{R_n}}^2}(R_n\left[\psi \right]-\overline{R_n}{)}^2.$(25) ), ψ being a solution of Equation (23(23) ${\mathrm{\nabla }}^2\psi +\mathrm{\eta \psi }=\frac{Q}{\gamma }\mathrm{with}\eta =\frac{{\mathrm{\nabla }}^2\gamma }{\gamma }+\frac{{\mathrm{\Sigma }}_a}{{\gamma }^2}.$(23) ), whereas $\rho ,{\mathrm{\Sigma }}_t,{\mathrm{\Sigma }}_a>0$. Then J has a unique minimum.

The proof of Theorem 4.1 relies on the following two lemmas, which are consequences of the weak and strong maximum principle, respectively.

Lemma 4.2

Let D be a bounded domain with smooth boundary $\mathrm{\partial D}$. Let u be a weakly differentiable function satisfying a mild regularity assumption (cf. Appendix A for details) and (29) $-\mathrm{\Delta u}\left(x\right)+\eta \left(x\right)u\left(x\right)=Q\left(x\right)\mathrm{in}D$(29) with bounded, non-negative η and Q that are strictly positive on a subset of D that is not a null set. Then u is strictly positive in D.

Lemma 4.3

Let D be a bounded domain with smooth boundary $\mathrm{\partial D}$. Then the problem (30) $\begin{array}{rl}-\mathrm{\Delta u}\left(x\right)+\alpha \left(x\right)u\left(x\right)& =Q\left(x\right)\mathrm{in}D\\ u\left(x\right)+\text{β}\left(x\right)\frac{\mathrm{\partial }u}{\mathrm{\partial }n}\left(x\right)& =0\mathrm{on}\mathrm{\partial D}\end{array}$(30) has at most one solution u if $\alpha \left(x\right),\text{β}\left(x\right)>0$.

The proofs of both lemmas are carried out in Appendices A and B. We can now present the proof of the main result.

Proof.

In the following, the different terms in the objective function (Equation (25(25) $J[\rho ,\psi ]=\sum _{l=1}^L\frac{w_{\gamma }}{{\overline{P_{E_{\gamma ,l}}}}^2}\left(P_{E_{\gamma ,l}}\right[\rho ,\psi ]-\overline{P_{E_{\gamma ,l}}}{)}^2+\sum _{n=1}^N\frac{w_d}{{\overline{R_n}}^2}(R_n\left[\psi \right]-\overline{R_n}{)}^2.$(25) )) are investigated separately. We define $J_{P,\min }$ as the set of values ρ that minimizes the first term in Equation (25(25) $J[\rho ,\psi ]=\sum _{l=1}^L\frac{w_{\gamma }}{{\overline{P_{E_{\gamma ,l}}}}^2}\left(P_{E_{\gamma ,l}}\right[\rho ,\psi ]-\overline{P_{E_{\gamma ,l}}}{)}^2+\sum _{n=1}^N\frac{w_d}{{\overline{R_n}}^2}(R_n\left[\psi \right]-\overline{R_n}{)}^2.$(25) ) (corresponding to the gamma detectors) and $J_{R,\min }$ as the set of values ρ the minimizes the second term in Equation (25(25) $J[\rho ,\psi ]=\sum _{l=1}^L\frac{w_{\gamma }}{{\overline{P_{E_{\gamma ,l}}}}^2}\left(P_{E_{\gamma ,l}}\right[\rho ,\psi ]-\overline{P_{E_{\gamma ,l}}}{)}^2+\sum _{n=1}^N\frac{w_d}{{\overline{R_n}}^2}(R_n\left[\psi \right]-\overline{R_n}{)}^2.$(25) ) (corresponding to the neutron detectors). Since both terms are positive, the intersection of $J_{P,\min }$ and $J_{R,\min }$ is then the set of values that minimizes the full objective function. The proof works as follows: showing that there exists a unique solution ${\psi }_{\min }$ that minimizes $J_{R,\min }$ using monotonicity arguments, we find explicit expressions for the elemental densities ${\rho }_l$ which are again unique.

Since the ψ is determined by the parameter η, the neutron counting rates defined in Equation (26(26) $R_n\left[\psi \right]={\int }_{V_{\det ,n}}{\sigma }_d\frac{1}{\gamma }\psi \left(x\right)\mathrm{d}x,$(26) ) can be seen as implicitly dependent of the parameter η. In the following we show that ψ and thus the counting rates are strictly monotone in η. We assume in the following that we have two parameters with $\tilde{\eta }>\eta$, and corresponding $\tilde{\psi }$ and ψ, such that Equation (23(23) ${\mathrm{\nabla }}^2\psi +\mathrm{\eta \psi }=\frac{Q}{\gamma }\mathrm{with}\eta =\frac{{\mathrm{\nabla }}^2\gamma }{\gamma }+\frac{{\mathrm{\Sigma }}_a}{{\gamma }^2}.$(23) ) is satisfied: (31) $\begin{array}{rl}-\mathrm{\nabla }\cdot \mathrm{\nabla \psi }+\mathrm{\eta \psi }& =\frac{Q}{\gamma }\\ -\mathrm{\nabla }\cdot \mathrm{\nabla }\tilde{\psi }+\tilde{\eta }\tilde{\psi }& =\frac{Q}{\gamma }.\end{array}$(31) By subtracting Equation  (31(31) $\begin{array}{rl}-\mathrm{\nabla }\cdot \mathrm{\nabla \psi }+\mathrm{\eta \psi }& =\frac{Q}{\gamma }\\ -\mathrm{\nabla }\cdot \mathrm{\nabla }\tilde{\psi }+\tilde{\eta }\tilde{\psi }& =\frac{Q}{\gamma }.\end{array}$(31) ) and defining (32) $\lambda =\psi -\tilde{\psi }\mathrm{and}\delta =\eta -\tilde{\eta }$(32) we obtain (33) $-\mathrm{\nabla }\cdot \mathrm{\nabla \lambda }+\mathrm{\eta \psi }-\tilde{\eta }\tilde{\psi }=0,$(33) which is equivalent to (34) $-\mathrm{\nabla }\cdot \mathrm{\nabla \lambda }+\tilde{\eta }\lambda =-\mathrm{\delta \psi }.$(34) By applying Lemma 4.2, we obtain the strict positivity of ψ and since $\delta <0$, the strict positivity of λ follows directly with the same lemma. Since λ is strictly positive, we have $\tilde{\psi }<\psi$ for $\tilde{\eta }>\eta$. The counting rate $R_n\left[\psi \right]$ defined in Equation (26(26) $R_n\left[\psi \right]={\int }_{V_{\det ,n}}{\sigma }_d\frac{1}{\gamma }\psi \left(x\right)\mathrm{d}x,$(26) ) is a positive linear functional of ψ, thus we also obtain monotonicity of $R_n\left[\psi \right]$. This ensures that there is just one ${\eta }_{\min }$ minimizing $J_n$.

With Lemma 4.3 it follows directly that there is only one ${\psi }_{\min }$ that minimizes the second term in Equation (25(25) $J[\rho ,\psi ]=\sum _{l=1}^L\frac{w_{\gamma }}{{\overline{P_{E_{\gamma ,l}}}}^2}\left(P_{E_{\gamma ,l}}\right[\rho ,\psi ]-\overline{P_{E_{\gamma ,l}}}{)}^2+\sum _{n=1}^N\frac{w_d}{{\overline{R_n}}^2}(R_n\left[\psi \right]-\overline{R_n}{)}^2.$(25) ), so that $J_{R,\min }=\left\{\rho |\eta \right(\rho )={\eta }_{\min }\}$. Now we can compute the intersection of the sets of minima of both terms in Equation (25(25) $J[\rho ,\psi ]=\sum _{l=1}^L\frac{w_{\gamma }}{{\overline{P_{E_{\gamma ,l}}}}^2}\left(P_{E_{\gamma ,l}}\right[\rho ,\psi ]-\overline{P_{E_{\gamma ,l}}}{)}^2+\sum _{n=1}^N\frac{w_d}{{\overline{R_n}}^2}(R_n\left[\psi \right]-\overline{R_n}{)}^2.$(25) ). Inserting ${\psi }_{\min }$ and using the homogeneity of the sample and the linearity of the integral, Equation (27(27) $P_{E_{\gamma ,l}}[\rho ,\psi ]={\int }_{V_{\mathrm{sam}}}{ϵ}_{E_{\gamma ,l}}\left(x\right){\rho }_l{\sigma }_{E_{\gamma ,l}}\frac{1}{\gamma }\psi \left(x\right)\mathrm{d}x.$(27) ) yields (35) $\overline{P_{E_{\gamma ,l}}}=\frac{{\rho }_l}{\gamma }{I}_{E_{\gamma ,l}}={\rho }_l\sqrt{3\sum _{l^{\prime }}{\sigma }_{t,l^{\prime }}{\rho }_{l^{\prime }}}{I}_{E_{\gamma ,l}}$(35) where (36) $I_{E_{\gamma ,l}}={\int }_{V_{\mathrm{sam}}}{ϵ}_{E_{\gamma ,l}}\left(x\right){\sigma }_{E_{\gamma ,l}}{\psi }_{\min }\left(x\right)\mathrm{d}x.$(36) From this system of L equations, where L is the number of elements, we can, using some algebra, derive linear relations for between the different elemental densities, which are valid for all pairs $l^{\prime },l\mathrm{with}1\le l^{\prime },l\le L$: (37) ${\rho }_{l^{\prime }}=\frac{\overline{P_{E_{\gamma ,l^{\prime }}}}{I}_{E_{\gamma ,l}}}{\overline{P_{E_{\gamma ,l}}}{I}_{E_{\gamma ,l^{\prime }}}}{\rho }_l.$(37) By inserting them in Equation (35(35) $\overline{P_{E_{\gamma ,l}}}=\frac{{\rho }_l}{\gamma }{I}_{E_{\gamma ,l}}={\rho }_l\sqrt{3\sum _{l^{\prime }}{\sigma }_{t,l^{\prime }}{\rho }_{l^{\prime }}}{I}_{E_{\gamma ,l}}$(35) ), we obtain for each element l the relation (38) $\overline{P_{E_{\gamma ,l}}}={\rho }_l^{3/2}\sqrt{3\sum _{l^{\prime }}{\sigma }_{t,l^{\prime }}\frac{\overline{P_{E_{\gamma ,l^{\prime }}}}{I}_{E_{\gamma ,l}}}{\overline{P_{E_{\gamma ,l}}}{I}_{E_{\gamma ,l^{\prime }}}}}{I}_{E_{\gamma ,l}}.$(38) Since all expressions apart from ${\rho }_l$ are constant, it is apparent that there is exactly one value for each ${\rho }_l,1\le l\le L$, that satisfies the equation. As a consequence, $J_{\min }$ contains just one point.

It should be noted that the argument would still be valid if a different model had been chosen for the gamma counting rates defined in Equation (6(6) ${ϵ}_{E_{\gamma ,l}}\left(x\right)={ϵ}_{\mathrm{geom},E_{\gamma ,l}}\left(x\right)\cdot {ϵ}_{\mathrm{att},E_{\gamma ,l}}\left(x\right)\cdot {ϵ}_{\det ,E_{\gamma ,l}},$(6) ). The generalization to heterogeneous media is discussed in the concluding Section 7.

5. Implementation

For the neutron transport computations, a neutron transport code named SPARC [6] is used. We apply the finite element method for the spatial discretization of Equation (3(3) $-\mathrm{\nabla }\cdot \left(\frac{1}{3{\mathrm{\Sigma }}_t\left(x\right)}\mathrm{\nabla \varphi }\right(x\left)\right)+{\mathrm{\Sigma }}_a\left(x\right)\varphi \left(x\right)=Q\left(x\right).$(3) ), to handle the geometric complexity of the underlying model. Furthermore, the finite element method is well suited for Poisson-like equations, as the SP₁-equation. The FEniCS library [15,16] with the automatic mesh generator Gmsh [17] is used for the implementation. We use PETSc [18] as linear algebra backend in conjunction with hypre [19], a library of high performance preconditioners. For this work we use a combination of a gmres solver [20] and an algebraic multigrid preconditioner [21]. Since the SP₁-equation is self-adjoint, a modification of the source term is sufficient to also compute the adjoint neutron flux. The gradient of the Lagrangian can then be computed as detailed out in Section 3.

The model used for the following simulation studies is shown in Figure 4. It depicts a zoom into the relevant part of the computational domain. The remaining parts are filled solely with more moderating material. It depicts the moderator (green), the sample (yellow), the gamma detector (purple), the neutron detectors (red) and on the left side the source (in blue). Figure 4(a) shows the computational domain, while Figure 4(b) shows the Finite Element Mesh used for the reconstruction including approx. 84.000 Finite Elements.

Figure 4. Computational Model used in SPARC including the different regions of the computational domain (left side) and the Finite Element mesh used for the reconstruction (right side). Depicted are the moderator (green), the sample (yellow), the gamma detector (purple), the neutron detectors (red) and on the left side the source (in blue). (a) Computational Domain and (b) Computational Mesh used for Reconstruction.

6. Numerical results

In this section, we present numerical results which show the applicability of the herein developed theory. In Subsection 6.1, we first present a comparison of the gradient computed with the Optimal Control formalism as introduced in Subsection 3.2 and a Finite Difference approach in order to verify our implementation and the adjoint calculus in Section 3. Thereafter, we show exemplary minimizations for which the computed gradient is used within an L-BFGS-B scheme. In Subsection 6.2 we study both the stability of the reconstruction under given measurement errors as well as the applicability of our method for different materials. Throughout our simulation studies we have used fixed numerical grids, but the investigation of the effect of measurement noise of varying intensities in Subsection 6.2 shows that we are not running into problems related to the so-called inverse crime in accordance with [22].

6.1. Verification of the implementation

In order to verify the implementation, the gradient of the objective function is computed with both Finite Differences and the Optimal Control formalism. Figure 5 depicts an exemplary case which compares the gradient components along the line segment defined by ${\rho }_2=1-{\rho }_1$, ${\rho }_1\in [0,1]$. In the upper subfigure the absolute value of the gradient components is taken so that a logarithmic scale can be applied for better visualization. The plot shows a very good agreement between the results for both ${\rho }_1$ and ${\rho }_2$. In the lower subfigure, the difference between the gradients computed with both methods is depicted. The relative error is smaller than $0.1\mathrm{\%}$ for both ${\rho }_1$ and ${\rho }_2$.

Figure 5. Comparison of the partial derivatives of the objective function computed with Finite Differences and the Optimal Control Formalism along the line segment defined by ${\rho }_2=1-{\rho }_1$, ${\rho }_1\in [0,1]$. (a) Comparison of the gradient components computed with Finite Differences and the Optimal Control Formalism and (b) Relative error of the gradient components computed with the Optimal Control formalism.

While Figure 5 depicts the verification of the gradient itself, Figure 6 depicts the gradient-based minimization process. The gradients computed with the Optimal Control formalism are used to carry out the minimization using the L-BFGS-B [23] algorithm implemented in SciPy [24]. Figure 6 depicts the objective function as well as three exemplary minimizations. The minimizations are represented by the blue, cyan and green lines, while the objective function is represented by a contour plot. It should be noted that each colour in the contour plot represents an order of magnitude, so a major challenge here is to minimize a function which values vary over many orders of magnitude. For each case, the minimum (purple) is found within only a few iterations.

Figure 6. The logarithmised objective function with blue, cyan and green lines showing three exemplary minimizations using the L-BFGS-B algorithm.

6.2. Simulation studies

In order to investigate the stability of the reconstruction under measurement errors, a simulation study has been carried out with normally distributed errors. The standard deviations were as large as 5, 10, and 20% of the simulated counting rates. For each error level 50 Simulations have been carried out. The reconstructed values are depicted in Figure 7. The error-free reconstruction is depicted as the black point in the center, while the reconstructed values for ρ are depicted in red, blue, and green, depending on the error level. It can be seen that the reconstructed values for ρ are scattered in a wider area for larger error levels. Table 1 shows the reconstruction errors for the different error levels. The errors are to be understood as mean absolute percentage errors. For all three cases the reconstruction errors are much smaller than the errors on the counting rates. At the same time the number of needed L-BFGS-B iterations is nearly independent from the error on the simulated counting rates.

In order to test the applicability of our method, four cases with different material properties have been investigated. The materials themselves are artificial, but using parameters that are realistic for the application. The first two cases are combinations of a highly scattering and a highly absorbing material in each case. The third and fourth cases investigate combinations of materials with similar neutronics properties, in one case two highly scattering materials and in the other case two highly absorbing materials. The upper half of Table 2 shows the errors of the reconstruction for the four different cases. It can be seen that the reconstruction of ρ works well for all investigated cases. The error on the counting rates is $10\mathrm{\%}$ , while the error on the reconstructed densities is typically 5–$6\mathrm{\%}$ . In one case the error on ${\rho }_2$ is as low as $2.6\mathrm{\%}$. This seems to be an exception, however.

Figure 7. Reconstructed values of ρ for normally distributed measurements errors. 50 Simulations were carried out for each error level. The average errors on the reconstruction and the number of needed L-BFGS-B iterations can be found in Table 1.

Table 1. Average absolute errors on the reconstruction of ${\rho }_1$/${\rho }_2$ for three different error levels as well as the average number of needed L-BFGS-B iterations.

Std. dev./total value measurement [%]	Reconstruction error ${\rho }_1$ [%]	Reconstruction error ${\rho }_2$ [%]	#iterations
5	3.1	1.4	6.3
10	6.1	3.4	6.3
20	11.7	8.0	6.6

For each case, 50 simulations have been carried out. The reconstructed values are depicted in Figure 7. The used material properties are the same as defined in Table A1 for the case Mixed Scattering/Absorption 1.

Table 2. Average absolute errors on the reconstruction of ${\rho }_1$/${\rho }_2$ for four different material combinations as well as the average number of needed L-BFGS-B iterations.

Reconstruction with regularization	Reconstruction error ${\rho }_1$ [%]	Reconstruction error ${\rho }_2$ [%]	#iterations
Mixed Scattering/Absorption 1	5.4	2.6	6.3
Mixed Scattering/Absorption 2	5.2	5.5	5.5
Scattering Dominated	6.1	5.1	5.8
Absorption Dominated	5.9	6.1	4.9
reconstruction without regularization
Mixed Scattering/Absorption 1	50.0	51.8	11.0
Mixed Scattering/Absorption 2	23.7	23.3	5.8
Scattering Dominated	63.1	66.1	9.0
Absorption Dominated	46.5	48.9	9.9

The measurement values have a normally distributed error of 10%. For each material combination, 50 simulations have been carried out. The material properties can be found in Table A1. The upper half of the table shows the reconstruction with the term punishing deviations from the total Mass $\overline{M}$ which acts as a regularization term. The lower half of the table shows the reconstruction without said term. For both cases, the material properties and the measurement values are identical.

The regularizing effect of the term punishing deviations from the total mass introduced in Equation (12(12) $J_m\left[\rho \right]=\frac{w_m}{{\overline{M}}^2}(c_1{\rho }_1+c_2{\rho }_2-\overline{M}{)}^2$(12) ) can be seen by comparing the upper and the lower half of Table 2. The lower half shows the reconstruction results with the weight of the term $w_m$ set to zero. The same initial and measurement values as before have been used. Without the regularization, the error of the reconstruction massively increases, from values between 2.6 and $6.1\mathrm{\%}$ to values between 23.3 and $66.1\mathrm{\%}$ . The average number of iterations also increased considerably, from approximately five to six to values between nine and eleven. The second test case is the sole exception. These results underline the importance of the regularization term in order to increase the accuracy of the proposed method.

7. Conclusion and discussion

In this work we have considered a new kind of parameter estimation problem related to PGNAA where two different types of detectors are used, one to record a gamma spectrum and additional detectors for the total neutron flux in direct proximity to the sample. This parameter estimation problem arises in a specific measurement facility, QUANTOM, where radioactive waste drums are analyzed for their elemental composition. We invert a combined neutron-photon transport model to obtain a vector of elemental densities. In a simplified setting considering a homogeneous sample and one-group neutron transport we proved uniqueness of a solution. To our knowledge, this is the first mathematically rigorous proof of uniqueness for a PGNAA problem.

To the current state of the art PGNAA spectra are evaluated in an iterative manner by solving a fixed point problem [7,8]. This naturally reflects the dependence of the elemental composition and the neutron flux distribution in the sample. The method proposed in this work offers an alternative. Formulating the parameter estimation problem as an optimization problem we make use of a gradient descent method to solve for the elemental composition as a parameter of the neutron transport model. As a side effect, our method yields the elemental composition as absolute quantities in grams or kilograms. Furthermore, this approach is especially useful when several measurements are combined, since all measurements can be combined in the objective function. For the measurements in QUANTOM the relative formulation in [7,8] does not apply. Hence, our method offers more flexibility to deal with a broader range of measurement tasks related to PGNAA. In situations where both methods may be applied (standard PGNAA with one measurement and one detector type) it is not obvious which methods yields the best performance in terms of convergence and computation times. This will be a topic of future research.

Furthermore, we see a close connection to the theory of Diffuse Optical Tomography (DOT). In both DOT and the PGNAA problems presented in this paper, the governing equation is a diffusion equation. The measurement of the flux alone is not sufficient to ensure a unique reconstruction. In DOT, it is not possible to reconstruct both the scattering and absorption coefficients from the photon flux measurements alone [14]. Similarly, we have seen that the neutron flux is insufficient to reconstruct the elemental densities, as discussed in Section 4. Only the combination with the energy resolved measurements of the gamma detector, which yields measurement results for each element individually, ensures uniqueness. While in DOT no general uniqueness result exists due to the range of different measurables and unknowns [25], it is still possible to obtain uniqueness results for certain scenarios using additional assumptions on the scattering and absorption coefficients that are to be reconstructed [26]. For PGNAA comparable studies do not exist yet since this work represents the first result in this direction. Therefore, a complete study of possible necessary and sufficient conditions for existence and uniqueness of solutions should be a focus of future work.

The measurement results of the additional 3He neutron counters are included in the objective function (Equation (11(11) $J_{\det }[\rho ,\varphi ]=\sum _{l=1}^L\frac{w_{\gamma }}{{\overline{P_{E_{\gamma ,l}}}}^2}\left(P_{E_{\gamma ,l}}\right[\rho ,\varphi ]-\overline{P_{E_{\gamma ,l}}}{)}^2+\sum _{n=1}^N\frac{w_d}{{\overline{R_n}}^2}(R_n\left[\varphi \right]-\overline{R_n}{)}^2.$(11) )). Monotonicity of the neutron flux with respect to the parameter η leads to monotonicity of the counting rates and this fact was used to prove the main uniqueness result Theorem 4.1. The proof works for any measurement device of the quantity ‘total neutron flux’, as long as its constituting functional, representing the measurement device, is linear and positive with respect to the neutron flux ϕ. This means, that in praxis devices like fission chambers are also adequate to support QUANTOM measurements. Furthermore, the positioning and the number of 3He detectors were not used in the proof. Nevertheless, they might influence both the accuracy and speed of evaluations and are worth to be further investigated in the future to optimize QUANTOM measurements.

The homogeneity of the sample is a key assumption for our proof of uniqueness in Section 4. This case is of practical importance for PGNAA, since most measurements are carried out with homogeneous samples, especially at research reactors [5]. Monitoring of the total neutron flux is applied by monitor materials in this case. Hence, our results may be adapted also for this situation.

In the case of an inhomogeneous sample our results do not guarantee a unique solution. In the QUANTOM facility, this situation is covered by performing numerous measurements at different measurement positions, corresponding to a complete surface scan of the drum by the HPGe detector. Measurements of the total neutron flux by the 3He detectors are also performed in every measurement position. A joint evaluation of all measurement positions shall account for inhomogeneities in the drum. Nevertheless, the rigorous examination of uniqueness of a solution in this situation needs to be performed in the future. It should be noted that the gradient based method in introduced in Section 3 holds for the inhomogeneous case. Since more measurements are necessary, the objective function will include more terms. The method relies then on the assumption of homogeneity in discretized parts of the drum.

Acknowledgments

The authors would like to thank Framatome for the permission to use the cut-view of the QUANTOM facility (Figure 1). We would also like to thank René van Geemert (Framatome) for helpful discussions on the topic of parameter estimation problems and neutron transport. We would also like to acknowledge the contribution of Jannick Wolters and Christopher Helmes in the development of the SPARC code, which was used to produce the numerical results for this paper.

Disclosure statement

The authors report there are no competing interests to declare.

Additional information

Funding

The work of A. J. was performed within the project RAPID, funded by the German Federal Ministry of Education and Research (BMBF) under the grant number 033RK094B. The work of A.J. was also funded by the European Regional Development Fund (ERDF) project ZEBRA, funded by the European Union and the federal state of North Rhine-Westphalia under the funding number EFRE-0800566. The work of K.K., when affiliated with AiNT, was performed in parts within the project QUANTOM, funded by the German Federal Ministry of Education and Research (BMBF) under the grant no. 15S9406B, and in project RAPID, funded by the German Federal Ministry of Education and Research (BMBF) under grant no. 033RK094A. The authors are solely responsible for the content of this publication.

Appendix A.

Proof of Lemma 4.2

Let D be a bounded domain with smooth boundary $\mathrm{\partial D}$. Let $u\in W^{1,2}\left(D\right)$ satisfy $-\mathrm{\Delta u}\left(x\right)+\eta \left(x\right)u\left(x\right)=Q\left(x\right)\mathrm{in}D$with bounded, non-negative η and Q that are strictly positive on a subset of D that is not a null set. Then u is strictly positive in D.

Proof.

We carry out the proof by constructing a contradiction. Let us first define $u=-\lambda$ so that Equation (29(29) $-\mathrm{\Delta u}\left(x\right)+\eta \left(x\right)u\left(x\right)=Q\left(x\right)\mathrm{in}D$(29) ) can be rewritten as (A1) $L\left(x\right)u\left(x\right)=Q\left(x\right)\mathrm{in}D.$(A1) with the differential operator defined by (A2) $L\left(x\right)\lambda \left(x\right)=(\mathrm{\Delta }-\eta (x\left)\right)\lambda \left(x\right).$(A2) The coefficients of L are bounded and the partial differential equation is strictly elliptic. We also have $\int -\mathrm{\eta v}\mathrm{d}x\le 0$ for $v\left(x\right)\le 0$ since $\eta >0$ and $\mathrm{Lu}\ge 0$ since Q is non-negative. Therefore we can apply the strong maximum principle (theorem 8.19 in [27]).

If $\underset{D}{sup}u\ge 0$, then according to the strong maximum principle u is constant in D and for $u\ne 0$ we have (A3) ${\int }_D\mathrm{\eta v}\mathrm{d}x=0\mathrm{\forall }v\ge 0,v\in C_0^1\left(D\right).$(A3) While u = 0 cannot satisfy Equation (A1(A1) $L\left(x\right)u\left(x\right)=Q\left(x\right)\mathrm{in}D.$(A1) ) since Q>0 on a non-null subset of D, $u\ne 0$ cannot satisfy Equation (A3(A3) ${\int }_D\mathrm{\eta v}\mathrm{d}x=0\mathrm{\forall }v\ge 0,v\in C_0^1\left(D\right).$(A3) ) since $\eta >0$ for non-null subset of D. Therefore, we have $\underset{D}{sup}u<0$ and thus $\underset{D}{inf}\lambda >0$.

Appendix B.

Proof of Lemma 4.3

Let D be a bounded domain with smooth boundary $\mathrm{\partial D}$. Then the problem $\begin{array}{rl}-\mathrm{\Delta u}\left(x\right)+\alpha \left(x\right)u\left(x\right)& =Q\left(x\right)\mathrm{in}D\\ u\left(x\right)+\text{β}\left(x\right)\frac{\mathrm{\partial }u}{\mathrm{\partial }n}\left(x\right)& =0\mathrm{on}\mathrm{\partial D}\end{array}$has at most one solution u if $\alpha ,\text{β}>0$.

Proof.

Let $u_1$ and $u_2$ be solutions of Equation  (30(30) $\begin{array}{rl}-\mathrm{\Delta u}\left(x\right)+\alpha \left(x\right)u\left(x\right)& =Q\left(x\right)\mathrm{in}D\\ u\left(x\right)+\text{β}\left(x\right)\frac{\mathrm{\partial }u}{\mathrm{\partial }n}\left(x\right)& =0\mathrm{on}\mathrm{\partial D}\end{array}$(30) ) and $v=u_1-u_2$. For $w=v\mathrm{\nabla v}$, the divergence theorem (A4) ${\int }_D\mathrm{\nabla }\cdot w\mathrm{d}x={\int }_{\mathrm{\partial }D}w\cdot n\mathrm{d}s$(A4) yields (A5) ${\int }_D(\mathrm{\nabla v}{)}^2+v\mathrm{\Delta v}\mathrm{d}x={\int }_{\mathrm{\partial }D}v\frac{\mathrm{\partial }v}{\mathrm{\partial }n}\mathrm{d}s.$(A5) Using $-\mathrm{\Delta v}+\mathrm{\alpha v}=0$ in D and inserting the boundary condition, we obtain (A6) ${\int }_D(\mathrm{\nabla v}{)}^2+\alpha v^2\mathrm{d}x={\int }_{\mathrm{\partial }D}-\text{β}{\left(\frac{\mathrm{\partial }v}{\mathrm{\partial }n}\right)}^2\mathrm{d}s.$(A6) Since the integral on the left side is non-negative and the integral on the right side non-positive, both integrands are identically zero. Therefore, $\frac{\mathrm{\partial }v}{\mathrm{\partial }n}=0$ and from the boundary condition also v is identically zero on $\mathrm{\partial D}$. From the weak maximum principle [27] follows that v is identically zero on D.

Appendix C.

Material properties used in Section 6

Table C1. Material properties.

	${\sigma }_{t1}$	${\sigma }_{a1}$	${\sigma }_{t2}$	${\sigma }_{a2}$	$c_1$	$c_2$
Mixed Scattering/Absorption 1	1.00	0.01	0.50	0.25	1	2
Mixed Scattering/Absorption 2	1.00	0.05	1.00	0.25	1	1
Scattering Dominated	1.00	0.01	1.00	0.05	1	1
Absorption Dominated	0.50	0.30	0.50	0.25	2	2