Simultaneous Reconstruction of Emission and Attenuation in Passive Gamma Emission Tomography of Spent Nuclear Fuel

TL;DR

This paper introduces a novel method for simultaneously reconstructing emission and attenuation maps in passive gamma emission tomography of spent nuclear fuel, improving image accuracy without relying heavily on prior information.

Contribution

The authors develop a constrained minimization approach that jointly reconstructs emission and attenuation maps, enhancing accuracy over traditional methods in nuclear safeguards inspections.

Findings

Significantly improved reconstruction accuracy over FBP.

Better classification of fuel rod status.

Effective in low prior information scenarios.

Abstract

The International Atomic Energy Agency (IAEA) has recently approved passive gamma emission tomography (PGET) as a method for inspecting spent nuclear fuel assemblies (SFAs), an important aspect of international nuclear safeguards which aim at preventing the proliferation of nuclear weapons. The PGET instrument is essentially a single photon emission computed tomography (SPECT) system that allows the reconstruction of axial cross-sections of the emission map of the SFA. The fuel material heavily self-attenuates its gamma-ray emissions, so that correctly accounting for the attenuation is a critical factor in producing accurate images. Due to the nature of the inspections, it is desirable to use as little a priori information as possible about the fuel, including the attenuation map, in the reconstruction process. Current reconstruction methods either do not correct for attenuation, assume a uniform attenuation throughout the fuel assembly, or assume an attenuation map based on an initial filtered back projection (FBP) reconstruction. Here, we propose a method to simultaneously reconstruct the emission and attenuation maps by formulating the reconstruction as a constrained minimization problem with a least squares data fidelity term and regularization terms. Using simulated data, we compare the proposed method to FBP, showing that our approach produces significantly better reconstructions by various numerical metrics and a much better classification of spent, missing, and fresh fuel rods.

Figures

Equations

$$H_0(\mu {)}_{i,p}=r_{i,p}\exp \left(-c_{i,p}{d}_{i,p}^T\mu \right).$$

$$r_{i,p}=\frac{1}{N_{p,\text{vox}}}\sum _{s=1}^{N_{p,\text{vox}}}{r}_s^{3D}.$$

$$c_{i,p}=\frac{N_{p,\text{vox}}}{r_{i,p}}\sum _{s=1}^{N_{p,\text{vox}}}\frac{r_s^{3D}}{\cos {\alpha }_s}.$$

$$H_{\varphi }(\mu {)}_{i,p}=r_{i,p}\exp \left(-c_{i,p}{d}_{i,p}^T{R}_{\varphi }\mu \right).$$

$$H(\mu )=\left[\begin{array}{c}{H}_{{\varphi }_1}(\mu )R_{{\varphi }_1}\\ ⋮\\ H_{{\varphi }_{N_{\text{ang}}}}(\mu )R_{{\varphi }_{N_{\text{ang}}}}\end{array}\right].$$

$${\tilde{H}}_{\varphi }(\tilde{\mu }{)}_{i,s}={\tilde{r}}_{i,s}\exp \left(-{\tilde{d}}_{i,s}^T{\tilde{R}}_{\varphi }\tilde{\mu }\right).$$

$$\begin{array}{rl}{\text{min}}_{(\lambda ,\mu )\in {\mathbb{R}}^{2N_{\text{pix}}}}& \left\{{\Vert F\left(\lambda ,\mu \right)-m\Vert }_2^2+\sum _i{\alpha }_i{P}_i\left(\lambda ,\mu \right)\right\}\\ & \text{subjectto}A\left[\begin{array}{c}\lambda \\ \mu \end{array}\right]\le b.\end{array}$$

$${\alpha }_{\lambda }{\Vert L\lambda \Vert }_2^2+{\alpha }_{\mu }{\Vert L\mu \Vert }_2^2,$$

$$\displaystyle\text{ker}_{L}=\left[\begin{array}[]{rrr}0&1&0\\ 1&-4&\>\>\>1\\ 0&1&0\end{array}\right].$$

$$B_{\lambda }=\left[\begin{array}{ccc}{r}_1& \cdots & r_{N_{\text{rod}}}\end{array}\right]\text{and}{B}_{\mu }=\left[\begin{array}{cccc}{r}_1& \cdots & r_{N_{\text{rod}}}& w\end{array}\right].$$

$${\alpha }_{\lambda }{\Vert P_{\lambda }\lambda \Vert }_2^2+{\alpha }_{\mu }{\Vert P_{\mu }\mu \Vert }_2^2.$$

$${\Vert r(\lambda ,\mu )\Vert }_2^2:={\Vert \begin{array}{c}F(\lambda ,\mu )-m\\ \sqrt{{\alpha }_{\lambda }}{M}_{\lambda }\lambda \\ \sqrt{{\alpha }_{\mu }}{M}_{\mu }\mu \end{array}\Vert }_2^2,$$

$$\displaystyle\left\|\begin{matrix}\left[\begin{array}[]{c}J_{r}(x^{(k)})\\ \sqrt{\beta^{(k)}}I\end{array}\right]x_{\text{step}}+\left[\begin{array}[]{c}r(x^{(k)})\\ 0\end{array}\right]\end{matrix}\right\|_{2}^{2},$$

$$A{x}_{\text{step}}\le b-A{x}^{(k)}.$$

$$2J_r(x^{(k)}{)}^T{J}_r(x^{(k)})+2{\beta }^{(k)}I,$$

$$\frac{{\Vert x_{\text{truth}}-x_{\text{recon}}\Vert }_2}{{\Vert x_{\text{truth}}\Vert }_2}\cdot 100\%,$$

Full text

1 Introduction

As part of an effort to deter the proliferation of nuclear weapons, various technical measures referred to as ”safeguards” are used to verify the declarations made by the signatories to the Treaty on the Non-Proliferation of Nuclear Weapons about their nuclear material and activities [1]. The monitoring of spent fuel assemblies (SFAs) from nuclear power plants (NPPs) is an important task within these safeguards, aiming at detecting any eventual diversion of spent nuclear fuel for non-declared purposes. Ideally, a single fuel pin missing from an SFA should be detected. For any safeguards investigation of SFAs, it is important to use a minimum amount of a priori information on the SFA under study, in order to avoid biasing and potentially misleading the investigation.

Each non-nuclear-weapon State Party to the Treaty is required to conclude a safeguards agreement with the International Atomic Energy Agency (IAEA). As such, safeguards activities are largely monitored and coordinated by the IAEA.

Since the 1980s, IAEA has developed, in collaboration with some of its Member States, gamma ray emission tomography (GET) for imaging SFAs [2]. GET was deemed attractive for detecting partial defects (part of the fuel of an SFA missing) and verifying the integrity of the SFA because it has the potential to directly image the spatial distribution of the active material and the relative locations of the pins in the SFA in a non-destructive way. This effort has culminated at the end of 2017 in IAEA approval to use the PGET instrument (Passive Gamma Emission Tomography) [3, 4, 5] in inspections.

The PGET instrument is able to reliably identify single missing or replaced pins in WWER-440, BWR and PWR SFAs with burnups in the range of 5.7–-57.8 GWd/tU and cooling times from 1.9-–26.6 years [3, 4, 5]. Images reconstructed using gamma ray energies higher than those of 137Cs (>700 keV) have better water-to-fuel contrast in fuel cooled for up to 20 years and thus have the most potential for missing fuel pin detection [6]. Following these first results, IAEA has expressed the need for new image reconstruction and processing methods for a more accurate assessment of the locations and count of missing pins and a more accurate calculation of the relative radioactivity levels of individual pins.

An SFA is a challenging object for tomographic imaging as it contains materials with very different emission and attenuation properties: strongly attenuating and emitting spent nuclear fuel (commonly uranium dioxide) and less attenuating material with zero emission (water or air). When considering diversion scenarios beyond missing pins, one could also consider, e.g., pins replaced with fresh nuclear fuel (strongly attenuating and zero emission) or with activated materials other than nuclear fuel (moderately attenuating and high emission).

An SFA consists of a regular lattice of about 100 to 300 fuel pins depending on fuel type, with most often one or several empty lattice positions (so-called water channels). Attenuation of the 662 keV gamma rays from 137Cs between the center and the edge of an SFA is of the order of a factor 100. This high gamma ray attenuation combined with the heterogeneous nature of SFAs makes detailed attenuation information essential for producing very realistic images.

Acquiring such attenuation information by means of a separate imaging procedure such as high energy CT is not practical from an operational point of view. Information on the geometry of the SFA, e.g., provided by the NPP, can in principle be used to obtain detailed information on attenuation. However, the requirement to make use of as little a priori information as possible precludes this.

On the other hand, partial defect detection does not necessarily require images with accurate (relative) intensities: it is more important to have good contrast between emitting and non-emitting regions. In practice, however, images with more accurate (relative) intensities typically have better contrast.

The development of GET for spent fuel verification has largely been conducted under the IAEA Support Program projects JNT 1510 and JNT 1955 (phase I). The JNT 1955 project [7, 8] and related work [9] used simulated data and investigated both analytic (filtered back-projection (FBP) wihout and with a posteriori attenuation correction) and algebraic image reconstruction techniques (which combined the Additive Simultaneous Iterative Reconstruction Technique with homogeneous attenuation information throughout the area covered by the SFA).

The PGET instrument recently approved by IAEA for inspections resulted from the JNT 1510 project [10]. In its implementation as approved for SFA inspections, the image from an initial FBP without attenuation correction (so without a priori knowledge) is used to determine (assumed) pin locations and fuel assembly location. This information is then used to construct a heterogeneous attenuation map that is included in a second image reconstruction using the Novikov inversion formula [11] and resulting in the final image [3].

In an effort to get closer to the goal of not using any a priori information on the materials and geometry in GET for SFA safeguards, we investigated an approach for the simultaneous reconstruction of emission and attenuation. In the context of medical single photon emission computed tomography (SPECT), the simultaneous reconstruction from emission data has been a topic of research since the late 1970s [12]. For a recent extensive review, we refer to [13]. Here, we mention only some iterative methods used for simultaneous reconstruction in the case of arbitrary attenuation maps [12, 14, 15, 16, 17, 18, 19, 20].

The approach we propose is close to the ones in [15] and [19]. Indeed, we also formulate the reconstruction as a minimization problem with a least squares (LSQ) data fidelity term and regularization terms. However, our choice for the regularizers and the minimization algorithm differ, and we also use linear bounds that are specific for the application.

In particular, we investigate two regularization terms. The smoothness prior has been extensively used in linear tomography problems yielding satisfactory results given that it is computationally efficient and simple to implement. On the other hand, the geometry aware prior, to the best of our knowledge, has not been proposed before and is tailored for this specific application. Compared to the smoothness prior, the geometry aware prior maintains the computational efficiency, but improves the reconstruction at the small cost of reasonable assumptions about the fuel assembly geometry, available, for instance, from an initial FBP reconstruction.

This paper is organized as follows. In Section 2 we describe the PGET instrument, introduce the discrete measurement model along with the simulation of data and the minimization problem. Results are presented in Section 3 and discussed in Section 4. We draw some conclusions and indicate future perspectives in Section 5.

2 Methods

2.1 Measurement with the PGET instrument

The PGET safeguards instrument measures the gamma ray emissions from the nuclear fuel. It is a 1D SPECT system, using a 1D linear collimator in front of a 1D array of gamma ray detectors. This geometry allows for the reconstruction of 2D cross-section images of the fuel. A simplified schematic representation of the instrument is shown in Figure 1. The PGET is made up of 2 detector banks with 87 CZT gamma ray detectors behind a tungsten collimator in each bank mounted on a plate inside a water-tight enclosure. During the course of each measurement, that plate rotates 360 degrees to measure data projections around the whole fuel assembly.

The CZT detectors have dimensions $2\text{\ mm}\times 4.8\text{\ mm}\times 4.8\text{\ mm}$. The detector and collimator pitch in each bank is 4 mm, and the position of the 2 banks on opposite sides of the rotating plate are offset by 2 mm so that the detectors of the banks can be interleaved to achieve an effective detector spacing of 2 mm.

The tungsten linear collimator slits are 10 cm deep, 1.5 mm wide, and taper from 70 mm tall at the front to 5 mm at the back. On all other sides except for the collimator opening, the detectors are shielded by at least 2 cm of tungsten. The water-tight enclosure is made out of 3-mm-thick stainless steel plating.

For each measured data projection, each CZT detector records the number of counts above 4 user-determined gamma-ray energy thresholds over a user-determined measurement time. These are used to calculate the number of gamma-ray counts in broad energy windows defined between consecutive energy thresholds. The thresholds are typically selected to enhance the contribution of a single gamma-ray emitting isotope (137Cs, 154Eu) in each window. More detailed descriptions of the PGET instrument can be found in [4, 6, 10].

The PGET instrument geometry was pixelized in order to create the measurement model used in our proposed image reconstruction algorithm.

2.2 Discrete measurement model

We assume that the volume contributing to the measurement is uniform in its emission and attenuation along the direction of the fuel rods. Further, we represent the volume by its 2D axial cross-section, which we divide into an $n$ by $n$ grid of pixels indexed from $1$ to $N_{\text{pix}}=n^{2}$. We denote by $\lambda=(\lambda_{1},\dotsc,\lambda_{N_{\text{pix}}})\in\mathbb{R}^{N_{\text{pix}}}_{+}$ the vector of the emission values of the cross-section, and likewise by $\mu=(\mu_{1},\dotsc,\mu_{N_{\text{pix}}})\in\mathbb{R}^{N_{\text{pix}}}_{+}$ the attenuation values. Here $\mathbb{R}_{+}$ is the set of non-negative real numbers. A scaled-down example of these discrete cross-section maps and the pixel indexing can be seen in Figure 2.

Only those elements of $\lambda$ and $\mu$ that correspond to pixels inside the maximal circular disk contained in the $n\times n$ images are within the region of interest in the measurement and thus are actually variables. This disk can be seen in the attenuation map in Figure 2. The elements outside the disk are always set to zero. This is easily implemented in practice and we ignore this consideration in the rest of the article to simplify the following descriptions.

As in [21], we first describe how the measurements are formed in the case that the detector array, with $N_{\text{det}}$ detectors, is located to the side of the cross-section images, as seen in Figure 2. This detector position is considered here to be the zero measurement angle. The measurements at other angles are then easily computed using this zero angle setup and rotating the contents of the emission and attenuation images.

The forward projection at the zero measurement angle can be expressed in the form $F_{0}(\lambda,\mu)=H_{0}(\mu)\lambda$, where $H_{0}(\mu)$ is a $N_{\text{det}}\times N_{\text{pix}}$ matrix depending on $\mu$. The element $H_{0}(\mu)_{i,p}$, which is the coefficient by which the emission value of pixel $p$, that is $\lambda_{p}$, contributes to the measurement at detector $i$, can be expressed as

[TABLE]

Here $r_{i,p}\in\mathbb{R}_{+}$ is the spatial response of pixel $p$ with regard to detector $i$, namely, it expresses the probability that photons emitted isotropically in the volume represented by pixel $p$ propagate towards the visible part of the detector $i$. How this is determined is described in detail below.

The $q$th component of the vector $d_{i,p}\in\mathbb{R}^{N_{\text{pix}}}_{+}$ is the distance that the line connecting the center of pixel $p$ and the center of detector $i$ travels inside pixel $q$. This is illustrated in Figure 3. The product $d_{i,p}^{T}\mu$ can be understood as a line integral of $\mu$ along the line from pixel $p$ to detector $i$.

The term $c_{i,p}>1$ is a correction factor for the distances $d_{i,p}$. This is to take into account the fact that the distance traveled by photons emitted from the volume represented by pixel $p$ (i.e. the vertical extent of the SFA seen through the collimator slits) is usually longer than the distance from the center of pixel $p$ to detector $i$. How this is formed is described in detail below.

The spatial response $r_{i,p}$ is computed similarly to [22]. First the volume that pixel $p$ represents is divided into voxels, indexed from $1$ to $N_{p,\text{vox}}$, as seen in Figure 4(a). Now the spatial response of each voxel $s$, denoted by $r_{s}^{3D}$ (we drop the dependence on $i$ and $p$ from the notation for simplicity), is just the probability that a photon emitted from the center of voxel $s$ starts off towards the visible part of detector $i$. This is equal to the solid angle spanned at the center of voxel $s$ by the visible part of detector $i$ divided by $4\pi$ (Figure 4(b)). The spatial response $r_{i,p}$ is then just the average of the spatial responses of individual voxels (Figure 4(c)):

[TABLE]

For the correction factor $c_{i,p}$ consider the angle $\alpha_{s}$ between the line from pixel $p$ to detector $i$ and the line from voxel $s$ to detector $i$ as illustrated in Figure 4(d). Multiplying the length of the former line by $1/\text{cos}\alpha_{s}$ gives the length of the latter line. The correction factor $c_{i,p}$ is now the weighted average of the factors $1/\text{cos}\alpha_{s}$ with the spatial responses $r_{s}^{3D}$ used as weights:

[TABLE]

If the image resolution $n\times n$ is low, then it is advantageous to compute the spatial responses $r$ and the correction factors $c$ using a higher resolution and then downsample them to $n\times n$ by averaging.

The forward projection at an arbitrary angle $\phi$ can now be expressed as $F_{\phi}(\lambda,\mu)=H_{\phi}(\mu)R_{\phi}\lambda$, where $R_{\phi}$ is a $N_{\text{pix}}\times N_{\text{pix}}$ matrix that rotates the contents of the cross-section images by angle $\phi$ using bilinear interpolation, and $H_{\phi}(\mu)$ is a $N_{\text{pix}}\times N_{\text{pix}}$ matrix defined similarly to $H_{0}(\mu)$ in (1):

[TABLE]

Finally, the whole forward projection with measurement angles $\phi_{1},\dotsc,\phi_{N_{\text{ang}}}$ can be expressed as $F(\lambda,\mu)=H(\mu)\lambda$, where $H(\mu)$ is the $N_{\text{det}}\cdot N_{\text{ang}}\times N_{\text{pix}}$ system matrix

[TABLE]

2.3 Simulation of data

Recovering both the attenuation and emission simultaneously is a nonlinear and ill-posed inverse problem. Therefore, it is important to avoid the so-called inverse crime [23, Ch. 2.3]. In other words, the simulated data should not be produced by exactly the same computational model that is used in the reconstruction algorithm, namely the model in Section 2.2.

For this purpose, we briefly introduce another model, which is essentially a fully 3D version of the previous one. Instead of considering only a cross-section and dividing that into pixels, we divide the whole volume that can be seen by the detectors through the collimator slits into voxels, indexed $1,\dotsc,N_{\text{vox}}$.

This 3D model can be described in much the same terms as the previous 2D one. We use notation with tilde for the concepts related to the 3D model. Let $\tilde{\lambda},\tilde{\mu}\in\mathbb{R}_{+}^{N_{\text{vox}}}$ denote the discrete emission and attenuation maps of the volume. Now the 3D forward projection at angle $\phi$ can be expressed as $\tilde{F}_{\phi}(\tilde{\lambda},\tilde{\mu})=\tilde{H}_{\phi}(\tilde{\mu})\tilde{R}_{\phi}\tilde{\lambda}$, where $\tilde{H}_{\phi}(\tilde{\mu})$ is a $N_{\text{det}}\times N_{\text{vox}}$ matrix depending on $\tilde{\mu}$ and $\tilde{R}_{\phi}$ is a $N_{\text{vox}}\times N_{\text{vox}}$ matrix that rotates the volume by $\phi$ degrees. The element $\tilde{H}_{\phi}(\tilde{\mu})_{i,s}$, which is the coefficient by which the emission value of voxel $s$, that is $\tilde{\lambda}_{s}$, contributes to the measurement at detector $i$, can be expressed as

[TABLE]

The term $\tilde{r}_{i,s}$ is the 3D spatial response. This is actually the same spatial response $r_{s}^{3D}$ that was used to compute the 2D spatial response in Section 2.2. The $u$th element of $\tilde{d}_{i,s}$ is the distance that the line connecting the center of voxel $s$ and the center of detector $i$ travels inside voxel $u$.

As in the 2D case, all the matrices $\tilde{H}_{\phi}(\tilde{\mu})$ can be composed so that the whole forward projection can be expressed as $\tilde{F}(\tilde{\lambda},\tilde{\mu})=\tilde{H}(\tilde{\mu})\tilde{\lambda}$, where the system matrix $\tilde{H}(\tilde{\mu})$ is similar to (5).

2.4 Minimization problem

We formulate the reconstruction from a measurement $m\in\mathbb{R}^{N_{\text{det}}\cdot N_{\text{ang}}}$ as a constrained minimization problem with a LSQ data fidelity term and regularization terms $P_{i}$:

[TABLE]

The purpose of the regularization terms is to compensate for the incomplete data by incorporating a priori knowledge about the unknowns in the reconstruction. Regularized inversion is robust against modelling errors and measurement noise. For more information on regularization of nonlinear ill-posed inverse problems, see [23, 24].

The regularization parameters $\alpha_{i}$ balance the effect of the data fidelity term and the regularization terms. Matrix $A$ and vector $b$ are such that the inequality in (7), understood to hold componentwise, defines a convex set.

The data fidelity term $\left\|F\left(\lambda,\mu\right)-m\right\|_{2}^{2}$, if seen only as function of emission $\lambda$, is convex, but as function of attenuation $\mu$, or both $\lambda$ and $\mu$, it is non-convex. It is a smooth function in all cases.

2.4.1 Regularization terms

We use two different choices for the regularization terms. In both cases the terms are convex quadratics.

The first choice, called here the smoothness prior, predisposes the algorithm toward reconstructions that are smooth in the sense that the changes in the emission and attenuation images are gradual when moving from one pixel to the next. It has the form

[TABLE]

where $L$ is the discrete Laplace operator, that is, a two dimensional convolution with the kernel

[TABLE]

The convolution is computed only at the points where the non-zero elements of the kernel stay inside the disk of variables.

The other choice for the penalty terms, called the geometry aware prior, assumes that the positions and the diameters of the possible rods, whether they are actually present or not, are known. In practice this information can be gained, for example, by identifying the assembly type and its position from an FBP reconstruction.

Let $r_{i}\in\mathbb{R}^{N_{\text{pix}}}$, $1\leq i\leq N_{\text{rod}}$, be vectorized images each displaying one of the $N_{\text{rod}}$ rods and nothing else. In addition let $w\in\mathbb{R}^{N_{\text{pix}}}$ be a vectorized image of water outside the rod positions. A scaled-down example of these images can be seen in Figure 5.

We wish to assert that our emission reconstruction is close to being a linear combination of the rod vectors $r_{i}$, i.e., that it is close to the subspace $S_{\lambda}=\text{span}\left(r_{1},\dotsc,r_{N_{\text{rod}}}\right)$. Similarly we wish that the attenuation reconstruction is close to being a linear combination of the rods $r_{i}$ and the water vector $w$, meaning it is close to $S_{\mu}=\text{span}\left(r_{1},\dotsc,r_{N_{\text{rod}}},w\right)$.

To achieve this, we define the matrices

[TABLE]

The expression $P_{\lambda}=I-B_{\lambda}(B_{\lambda}^{T}B_{\lambda})^{-1}B_{\lambda}^{T}$, where $I$ is the identity matrix, is the projection onto the orthogonal complement of $S_{\lambda}$. Define $P_{\mu}$ similarly. The geometry aware prior then has the form

[TABLE]

2.4.2 Bounds

The linear bounds that are used can be described as applying equally to all pairs of pixel values $(\lambda_{p},\mu_{p})$, that is, to all pairs of emission and attenuation values from the same pixel. Hence the bounds can be visualized in the emission-attenuation-plane, where they form a triangle, as seen in Figure 6. The values inside the triangle are feasible.

The bounds allow, in particular, the three materials relevant to us in this study: water (no emission, low attenuation), spent fuel rod (high emission, high attenuation) and fresh fuel rod (no emission, high attenuation), but they exclude the physically unlikely case of a material with high emission and low attenuation.

The triangle bounds can be described by giving an upper bound for emission values and both upper and lower bounds for attenuation values as these determine the three vertices of the triangle (the lower bound for emission is assumed to be zero). Some ways of estimating these upper and lower bounds are required in practice and this is discussed briefly in Section 4. Here we simply modify the true upper and lower bounds, as described in Section 3, to simulate error in estimating these values.

2.4.3 Minimization algorithm

The regularization terms that we use are such that the functional being minimized in (7) can be naturally written as non-linear LSQ term

[TABLE]

where the matrices $M_{\lambda}$ and $M_{\mu}$ depend on the choice of the penalty: they are either the discrete Laplace operator $L$ from (8) or the projection matrices $P_{\lambda}$ and $P_{\mu}$ from (14). We exploit this formulation of the problem and use a minimization method that is similar to the Levenberg-Marquardt algorithm (LMA) as described in [25].

Denote here by $x$ the combination of the emission and attenuation vectors, that is, $x=\left[\begin{matrix}\lambda^{T}&\mu^{T}\end{matrix}\right]^{T}$, and write $r(x)$ for the residual $r(\lambda,\mu)$ in (15). The LMA is a method of unconstrained optimization. At each iteration $k$, it minimizes, with regard to the next step $x_{\text{step}}$, a linear LSQ term that results from linearizing the residual $r(x)$ at the current iterate $x^{(k)}$ and from adding a regularization term:

[TABLE]

Here $J_{r}(x^{(k)})$ is the Jacobian matrix of the residual $r(x)$, $I$ is the $2N_{\text{pix}}\times 2N_{\text{pix}}$ identity matrix, and $\beta^{(k)}$ is the LM parameter modified at each step.

Differing from the usual LMA, we minimize (16) using linear constraints that keep the next iterate $x^{(k+1)}=x^{(k)}+x_{\text{step}}$ feasible:

[TABLE]

This minimization is done using the scaled gradient projection (SGP) method [26], where we use for scaling the inverse of the diagonal matrix that has the same diagonal as

[TABLE]

which is the Hessian of (16) with regards to $x_{\text{step}}$.

The Jacobian matrix $J_{r}(x^{(k)})$ is computed analytically, namely, based on an exact expression for it.

3 Results

The sinogram data used in this study was simulated with the 3D model described in Section 2.3. The fuel assembly phantoms used consist of $660\times 660\times 639$ voxels. Each axial column of voxels had uniform emission and attenuation properties. The measurements were simulated at 672 detector points per angle and then interpolated down to 167 detectors. Finally, Gaussian white noise, with a standard deviation of 2% of the maximum value of all the measurements, was added to form the data used in the reconstruction process.

Cross-sections of the phantoms used in simulating the data, downsampled by a factor of 4 to the PGET intrument’s reconstruction resolution of $165\times 165$ and cropped to include only the $69\times 69$ pixel area of the fuel assembly, can be seen in Figure 7. At this resolution, one pixel represent an area of $2\text{\ mm}\times 2\text{\ mm}$ when compared with the PGET instrument size.

The phantoms depict a GE12 assembly in water very near the center of the tomograph. An unmodified GE12 assembly consists of 92 UO2 rods with a diameter of 8.8 mm on a $10\times 10$ lattice with two $2\times 2$ regions without fuel. We modified this nominal GE12 assembly to include both missing rods and rods replaced by fresh UO2 rods at varying distances from the assembly center. The missing rods are in the top left half of the assembly and the replaced rods are in the lower right half. The attenuation phantom does not include the steel support structures of the assembly nor the steel interior wall of the tomograph. The attenuation values used are those corresponding to the 662 keV gamma-rays emitted by 137Cs, namely, 0.0085 mm-1 for water and 0.1356 mm-1 for UO2.

We compare reconstructions done using the two different regularization terms, the smoothness prior and the geometry aware prior, both using the same bounds, and also a reconstruction done using FBP. All the reconstructions used 90 measurement angles spaced evenly over the full 360-degree rotation.

The upper and lower bounds for emission and attenuation values that determine the triangle of the linear bounds used in the iterative reconstructions were slightly extended from the minimum and maximum values used in simulating the data. Namely, the emission upper bound was 10% larger than the value used for UO2 in simulating the data and the lower bound was 0; the attenuation upper bound was 5% larger than the value used for the UO2; and the attenuation lower bound was 5% smaller that the value used for water.

The initial guess for the iterative reconstructions consisted of only water everywhere. The regularization parameters $\alpha_{\lambda}$ and $\alpha_{\mu}$ were chosen heuristically by sampling several values and choosing the ones yielding the best reconstructions according to the quality metrics in Table 1.

The iterative reconstruction algorithm was stopped when the decrease in the objective function being minimized (the function in (7)) dropped below 0.1% between iterations. This resulted in 9 and 10 iterations for the smoothness and the geometry aware priors, respectively. The reconstruction process took 6 and 7 minutes, respectively, using a Matlab R2018a implementation of the algorithm on a laptop with i5-5300U CPU at 2.3 GHz and 16 GB of RAM.

Reconstruction images, cropped to include only the fuel assembly, can be seen in Figure 7, and in Table 1 are values of different metrics comparing these cropped reconstructions to the cropped ground truth. The metrics used are the relative error (RE), computed as

[TABLE]

the structural similarity index (SSIM) [27] and the Haar wavelet-based perceptual similarity index (HaarPSI) [28]. For the first metric, smaller is better. The values of the SSIM and HaarPSI metrics range from 0 to 1 and larger is better.

The tout court FBP reconstructed image contains pixels with negative values and its scale is entirely different from the ground truth. To better compare the methods, the FBP image is modified before applying the metrics and Figure 7 shows the modified version of the image. First, the negative values in the FBP image are set to zero, then the image is scaled so that the average of the pixel values in the cropped section matches the average of the cropped ground truth emission image. Finally, the pixel values that are larger than the emission upper bound that was used in the iterative reconstruction are set to the value of the upper bound.

Figure 8 shows plots constructed like those the IAEA uses for the classification of spent fuel rods and missing rods [4]. For every fuel array position in the reconstructed emission image, the average value over the central $2\times 2$ pixels is computed to represent the emission value of that position. In the plots, the difference between a position’s emission value and the average value of its neighbors is plotted against the distance of the position from the assembly center. The location of the water channels is not assumed to be known here, that is, they are not excluded when computing the average of the neighboring positions.

In Figure 9, emission and attenuation values of rod positions, computed again as averages of the $2\times 2$ pixel centers, are plotted in the emission-attenuation-plane. These plots present an alternate classification tool for fuel array positions in an SFA. However, they require that both emission and attenuation are reconstructed and hence are not applicable to the FBP reconstruction.

4 Discussion

By all three metrics in Table 1 and by visual comparison of the reconstruction images in Figure 7, the iterative reconstruction methods proposed here produced a significant improvement over the FBP reconstruction. Again by all three metrics, of the two regularization choices for the iterative method, the geometry aware prior produced better results than the smoothness prior, although the difference is not as significant as between the FBP and the iterative methods.

From the plots in Figure 8 it is clear that compared to FBP the iterative reconstruction methods produce a better separation between the rod positions with and without spent fuel. Therefore, the proposed methods should allow for easier classification of the rods, if the classification is based on similar images. In this comparison, both choices for the prior term produce equally good results.

Reconstructing the attenuation and emission simultaneously with an iterative method offers two advantages. First, the attenuation correction makes the relative emission values more accurate. Second, the rod positions can be classified using the combination of their emission and attenuation values, that is, based on the plots in Figure 9. Here, both iterative methods produce a clear separation between positions with spent fuel rods, fresh fuel rods and water. The geometry aware prior produces a somewhat tighter grouping of the positions with water in them than the smoothness prior.

The smoothness prior is a trade-off between model accuracy, computational performance and the amount of prior information used. Both the emission and attenuation values make a jump at the boundaries of rod and water, especially for higher reconstruction resolution. Therefore, assuming smoothness of those coefficients is not entirely accurate. However, the smoothness prior is well-understood, computationally efficient, and simple to implement and explain. Moreover, it does not assume any information about the geometry of the assembly. These are strong benefits in a real-world safeguards imaging task.

The geometry aware prior maintains the computational performance of the smoothness prior and improves on the model accuracy at the cost of making assumptions about the fuel assembly geometry. However, these assumptions are not unreasonable as they amount to identifying the fuel assembly type and its position, which is something that can be done from an FBP reconstruction. Such an identification is part of the current method used by IAEA [3]. This identification is also very relevant for the second objective of GET, which is the quantitative assessment of individual fuel rod properties (e.g. the activity of key isotopes, cooling time, relative burnup), as knowledge of basic fuel parameters (e.g., assembly type and nominal fuel composition) is considered necessary for it [7, 8].

The bounds deliver a significant part of the reconstruction quality of the proposed method. They require some way of estimating the upper and lower bounds for the attenuation values and the upper bound for the emission values. The reconstruction quality is quite sensitive to getting these estimates somewhat correct and this is likely to be a challenge when moving to real data.

The attenuation bounds could be estimated based on the knowledge of the measurement energy window and assumptions about the materials being imaged, i.e., that water is the least and UO2 the most attenuating material present. One way of estimating the emission upper bound could be to again identify the fuel assembly type and its location from an FBP reconstruction, and then quickly simulate a sinogram using the assumed attenuation values for water and UO2 and some constant emission value for all the rods. The ratio between, e.g., the average value of the simulated sinogram and the constant emission value used in the simulation should be somewhat close to the ratio of the average value of the real sinogram and the upper bound for emission values, even if the real data comes from an assembly that is missing a few rods.

The way the regularization parameters $\alpha_{\lambda}$ and $\alpha_{\mu}$ were chosen here, by simply sampling several values and picking the ones that produced the best reconstruction, is time consuming. Yet these choices, once found, should work relatively well at least for some other reconstruction tasks. Since the forward model is linear in the emission $\lambda$, the input sinogram can be normalized so that neither the measurement time nor the intensity of the radiation should affect these choices much. Change of the assembly type or the measurement energy window(s) might have a larger effect. One may consider calibrating the regularization parameters using an imaging setup in a controlled environment before deploying the system to operative use. Another option is to use an automatic parameter choice method such as cross validation or discrepancy principle. The problem of parameter choice is outside the scope of this initial feasibility study.

The details of the current IAEA method are not publicly available. Although FBP is not the state-of-the-art reconstruction method for this application, it provides a well-known method for comparison. One could likely find a different way to scale the FBP reconstruction that would somewhat improve its performance by the metrics in Table 1, but the shape of the plot in Figure 8 is independent of this scaling. Also, using more measurement angles would enhance the FBP reconstruction to some extent, but the fact that the iterative method does not need more angles for a good quality reconstruction is a benefit as this allows for shorter measurement times.

With more optimized software and hardware, the reconstruction times would likely drop to a level that is acceptable from an operational point of view, i.e., below a minute or two.

This study is based on simulated data only. While we put effort into modelling the geometry of the measurement setup and the radiation physics, some real world phenomena are still left out. These effects include, for example, the fact that the detected radiation is not actually monoenergetic and that some part of the detected photons are scattered either in the fuel, in the instrument or in the detectors themselves. This will affect the accuracy of the forward model, but also the method by which the bounds are estimated, and will likely lead to further challenges when taking the method into practice. However, the computational model can be extended to include or approximate at least some of these additional features, so we believe that our regularized reconstruction approach offers a flexible framework for developing a real-world method that considerably improves image reconstruction in GET of SFAs. In particular, the reconstruction of an attenuation image of the SFA opens new possibilities for classification and inspection criteria for fuel.

5 Conclusion

In this paper, we propose an iterative reconstruction method for simultaneous reconstruction of emission and attenuation maps of axial cross-sections of SFAs from PGET measurements. The performance of the method, with two different regularization choices, was compared to FBP using simulated data with 90 measurement angles.

The proposed method shows significant improvement over FBP across multiple metrics that compared the reconstructions to the ground truth. It produces a better separation than FBP between spent fuel rods and rods that are missing or replaced with fresh fuel rods when classifying the rod positions in the emission reconstructions. Furthermore, the proposed method allows for a different, enhanced, approach to classification by using also the reconstructed attenuation information.

Of the two regularization choices for the proposed method, the geometry aware prior performs somewhat better than the smoothness prior, but the difference is not large. The geometry aware prior assumes some information about the geometry of the SFA being imaged, but this information can be estimated from an initial FBP reconstruction.

We expect further challenges in taking the method into practice with real data, but believe that this framework of regularized iterative reconstruction is a good starting point for a real world method that improves image reconstruction in GET of SFAs.

6 Acknowledgment

RB, CBC and PD were supported by Business Finland under Grant 1845/31/2014. RB also acknowledges partial support by the Academy of Finland through the Finnish Centre of Excellence in Inverse Modelling and Imaging 2018-2025. TAB and SS acknowledge support by the Academy of Finland through the Finnish Centre of Excellence in Inverse Modelling and Imaging 2018-2025, decision number 312339, and Academy Project 310822. TH was supported by the Academy of Finland via the decision number 314879.