Physics-based forward model for near-real-time quantitative imaging of spent nuclear fuel assemblies

Abstract

The Passive Gamma Emission Tomography (PGET) instrument, authorized by the International Atomic Energy Agency (IAEA) for verification of spent nuclear fuel, aims to reconstruct 2D cross-sectional images of spent fuel assemblies (SFAs), identify missing or present fuel pins, and quantify fuel pin activities. Although the first two objectives are reliably achieved, accurate determination of fuel pin activities remains a challenge due to intense self-shielding and scattering effects. We have developed a linear inverse approach that addresses these effects and demonstrated superior image quality and identification accuracy in simulation studies. This approach frames the image reconstruction process as an inverse problem, relying on a physics-based forward model of the PGET system. We improved our forward model by incorporating collimator septal penetration and detector scattering effects. The enhanced forward model enables near-real-time sinogram simulation and system matrix calculation, which is> 100,000 times faster than 3D Monte Carlo simulations. The model was validated through simulations of VVER-1000 and VVER-440 SFA, and a relative difference of 3.7% in counts was achieved between MCNP and our forward model. Based on this enhanced model, we successfully reconstructed images from the simulated data, identified 100% of the fuel pins, and achieved an average uncertainty of 2.3% in activity quantification. We applied the reconstruction method to measured data of VVER-440 SFAs, successfully imaging all the pins, including the innermost ones, and identifying the water channel within the SFA. The high accuracy and low computational cost of our forward model demonstrate its potential for real-world inspection scenarios and enable future algorithm development.

Introduction

Spent nuclear fuel (SNF) assemblies discharged from nuclear power plants are temporarily stored in water pools at plant sites before being transferred to permanent repositories¹ or reprocessing facilities². Before final disposal, the International Atomic Energy Agency (IAEA) requires verification that no nuclear materials were diverted from the pool, in accordance with the Treaty on the Non-Proliferation of Nuclear Weapons (NPT). The passive gamma emission tomography (PGET) instrument³, shown in Fig. 1a, was developed to help fulfill this mission. Compared to other spent fuel inspection systems, PGET is currently the only instrument capable of verifying spent fuel up to the level of a single fuel pin⁴. PGET is equipped with two arrays of tungsten-collimated high-energy-resolution cadmium-zinc-telluride (CZT) detectors to detect gamma rays emitted by long-lived fission products, such as ¹³⁷Cs and ¹⁵⁴Eu, in the spent fuel. The detectors revolve around the SFA and record counts at different angles, forming a sinogram (Fig. 1b). By applying appropriate imaging techniques, we can reconstruct a two-dimensional (2D) image of the SFA (Figure 1e). Automated fuel pin identification is performed to detect present or missing fuel pins, and their activities are quantified (Fig. 1f–g).

Fig. 1

PGET inspection of a SFA. (a) PGET system with a mock-up fuel assembly inside⁵. (b) Detector counts at different angles are recorded to form a sinogram. (c) A linear forward model of PGET is developed. (d) The collection of sinograms in response to each individual fuel pin, known as the system response matrix, is computed. (e) A 2D image of the SFA is reconstructed from the sinogram. (f) Fuel pins identified from the 2D image. (g) Histogram of number of photons emitted per pin, which is proportional to the fuel pin activity.

The accuracy of fuel pin identification and activity quantification depends on the quality of the reconstructed image. Filtered back-projection (FBP) has been employed to perform a fast reconstruction^6,7, but the image quality is unsatisfactory, especially when rendering fuel rods or water channels around the center of the SFA. The inverse approach, which formulates image reconstruction as an inverse problem, has demonstrated a superior image quality and identification accuracy^8,9,10. This approach employs a forward model to calculate the detector response for any pixelated gamma-ray source (Fig. 1c), forming a system response matrix (Fig. 1d). Image reconstruction is performed by unfolding the two-dimensional gamma-ray emission pattern from the measured sinogram, with prior knowledge of the system response matrix.

It is crucial for the forward model to take into account all relevant photon interactions, including both photoelectric absorption and Compton scattering, in the calculation. Absolute and accurate quantification of fuel pin activity is only achievable if the system response includes contributions from both unscattered and scattered photons. This is especially the case for fuel pins near the center of the assembly being inspected, where detector counts contributed by Compton-scattered photons are non-negligible due to the heavy attenuation by the uranium fuel^11,12,13. Real-time forward models have been developed based on ray-tracing techniques^8,14, accounting for attenuation but not scattering. It is of interest to examine how the reconstruction quality and ability to quantify single-pin activities improve by including the scattering effect¹⁵.

Another desirable feature of the forward model is a fast computation speed. Ideally, the model should generate the system response matrix within a minute to minimize the waiting time between sinogram acquisition and image reconstruction. Models developed with general-purpose Monte Carlo radiation transport codes are accurate but time-consuming and incompatible with this time frame. For example, Monte Carlo N-Particle (MCNP) and Serpent models^16,17 of PGET have been developed, taking days to weeks to simulate a single sinogram on a high-performance computing (HPC) cluster consisting of thousands of cores. Existing approaches also include performing simulations for a limited number of projection angles and augmenting the data with data interpolation¹⁴ or machine learning algorithms¹⁸. However, even a small subset of angular views can be time-prohibitive to simulate without access to HPCs, which is likely the case at a near-reactor SNF pool. These limitations highlight the need for a PGET model that can accurately account for the scattering effect and generate projections in real time without HPC processing. In addition, a limited number of measured datasets is available at the moment, and they do not cover all possible SFA types and fuel pin configurations. A fast forward model can augment the measured dataset by generating simulated sinograms of missing SFA configurations, facilitating the development and assessment of various reconstruction algorithms. This is particularly relevant for the development of machine learning algorithms^15,19, which require a variety of SFA configurations for training and testing. In our previous work^9,10, we have developed a fast forward model that accounted for the gamma-ray self-shielding and down-scattering effects. While our previous model allowed us to image accurately and render the fuel pin position, the uncertainty associated with the pins’ activity was up to 6%, and it was tested only on simulated data sets. In this work, we enhanced our model to improve the modeling of the collimator’s point spread function (PSF) and applied it to reconstruct both simulated and measured datasets. The paper is structured as follows. In the first section, we describe the details of our forward model and its validation against Monte Carlo simulations to demonstrate its accuracy. In the second section, we describe our image reconstruction approach based on the developed forward model and apply it to both simulated and measured datasets. The discussions and conclusions are presented in the last section.

Fast forward model based on convolutional forced detection

In this section, we describe the details of our forward model that calculates the detector response (sinogram) for any given SFA configuration. The forward model was validated by comparing against MCNP simulation, and a 3.7% relative difference in counts was achieved.

Overview

Figure 2 shows the schematic of the PGET system that we aim to model, which is divided into two regions by the 16.5-cm radius inner wall. The outer region (Region II) of the system encompasses two detector arrays, each of 91 CZT detectors with a 4-mm pitch, placed behind a tungsten collimator with a 1.5 mm slit width and 10 cm length. The two detector arrays are offset by half a pitch, resulting in an effective pitch of 2 mm. The CZT detector has a dimension of 3.5 mm × 3.5 mm × 1.75 mm, and the tungsten collimator has a thickness of 2.5 mm. The inner region (Region I) of the system houses the underwater SFA. During the measurement, the PGET is lowered into the water pool and a SFA is inserted into its central cavity by a fuel-handling system. The detectors are then rotated around the SFA, and the counts are registered at preset angular intervals. It should be noted that in our simulation, instead of rotating the detectors, we rotated the SFA in the opposite direction. Our goal is to calculate the number of counts recorded by each detector at each angle, given the the following input:

• 1. The material map, including the position of each fuel pin, the material composition and density. For simplicity, we assume all fuel pins are homogeneous along the axial direction.

• 2. The activity of each fuel pin. For simplicity, we assume no spatial variation of activity inside the pin.

• 3. The energy spectrum of the gamma rays emitted by the spent fuel. For simplicity, we assume the spectrum is the same for all fuel pins.

• 4. Number of photons to simulate (NPS).

Fig. 2

Imaging of a spent nuclear fuel assembly in water using PGET.

Our model is based on the Convolutional Forced Detection (CFD) algorithm²⁰, which has been used for response simulation in Single Photon Emission Computed Tomography applications. The calculation consists of two parts: stochastic photon transport in Region I accounting for the gamma-ray down-scattering and self-shielding effects, and deterministic detector response calculation in Region II using the collimator’s point spread function (PSF). While the photon transport in Region I requires on-site calculation as it depends on the material map of the SFA, the PSF can be pre-calculated and retrieved when needed, thereby significantly reducing the computation time. The detector counts are obtained by convolving the photon distribution in Region I with the PSF.

Step 1: photon transport in region I

Monte Carlo photon transport is performed in Region I, with the scheme in Fig. 3. At the beginning, a new photon with an initial weight of 1 is created inside the fuel region. The photon’s energy is sampled from the emission spectrum. A copy of this photon that originates from initial emission (Type I) is appended to a photon list that saves all photons that can possibly reach a detector and deposit its energy. The photon direction is sampled isotropically. The photon is then transported through Region I using the delta-tracking algorithm. At each interaction site, we force the interaction to be Compton scattering, as photoabsorption in Region I makes no contribution to detector counts, and the cross-sections for coherent scattering and pair production are small in our energy range of interest (400–1500 keV). The photon’s weight is multiplied by the ratio of Compton scattering cross section and total cross section $\frac{\sigma _{CS}}{\sigma _{total}}$, i.e., the probability that the interaction is Compton scattering, to avoid biasing the simulation (see Supplementary Fig. S1 online). A copy of this photon that originates from Compton scattering (Type II) is saved to the photon list. The scattering angle and photon’s energy are sampled from the Klein-Nishina formula. The photon is then transported to a new position and the process repeats until the photon exits Region I or its energy falls below the detection threshold. The process is repeated for NPS photons. After the simulation, we obtain a photon list that saves the position, energy, and direction of all photons that can possibly contribute to the detector counts.

Fig. 3

Photon transport in region I. $\vec {r}$: position. $\vec {v}$: traveling direction. E: energy. w: photon’s weight.

Step 2: deterministic calculation in region II

The next step is to calculate the detector counts averaged over all possible angles for each photon in the list that is about to undergo emission (Type I) or Compton scattering (Type II). Each photon can be treated as a photon source, either isotropic (Type I) or anisotropic (Type II). Due to the high collimation of the tungsten collimator, only photons traveling almost perpendicularly to the detection plane have a high chance of reaching a detector (Fig. 4). The probability of a photon reaching Region II without attenuation and the probability density of the photon being emitted/scattered in any direction can be both approximated as angular-independent within this small acceptance angle. Accounting for the self-shielding effect, a Type I source in Region I is equivalent to an isotropic source in vacuum of intensity:

$$\begin{aligned} S = w e^{-\int \mu (E) dx}, \end{aligned}$$

(1)

and a Type II source is equivalent to an isotropic source of intensity:

$$\begin{aligned} S = w e^{-\int \mu (E) dx} \times \frac{d\sigma /d\Omega }{\int _{4\pi }d\sigma } \end{aligned}$$

(2)

where w is the photon’s weight, $e^{-\int \mu (E) dx}$ represents the penetration probability along the most contributing path, and $d\sigma /d\Omega$ is the differential cross section of Compton scattering evaluated at the scattering angle that enables the photon to be scattered along the most contributing path. Therefore, the detector counts are the sum of counts from a list of isotropic sources, each with an intensity S, position, and energy determined by a photon saved in Step 1.

Fig. 4

Schematic of the CFD algorithm. An emission or scattering event in Region I is treated as an isotropic source. The most contributing path (colored in purple) is the one that enables the photon to travel perpendicularly to the detection plane.

The detector counts in response to isotropic sources, referred to as the point spread function (PSF), can be pre-computed and stored in a table because they are independent of SFA properties. With the PSF known, we can calculate the detector counts contributed by each photon in the saved list by convolving the PSF at the photon’s position. The calculation scheme is shown in Fig. 5. For a projection angle $\alpha$, we first rotate all the photons by $\alpha$ in the opposite direction of the detector rotation. We then sum over all photons to obtain the total detector counts.

Fig. 5

Detector counts calculation in Region II.

PSF calculation and validation

For parallel-hole collimators, PSF is usually approximated by a Gaussian function²⁰:

$$\begin{aligned} \begin{gathered} \textrm{PSF}(x,y) = \frac{1}{\sqrt{2\pi }\sigma (y)}e^{\frac{-x^2}{2\sigma (y)^2}}, \sigma (y) = \frac{1+\frac{y}{l-2/\mu (E)}}{2\sqrt{2\ln {2}}}b \end{gathered} \end{aligned}$$

(3)

where b is the collimator aperture, l is the collimator length, $\mu (E)$ is the total attenuation coefficient in the collimator, x and y are the horizontal and vertical distance between the source and detector, respectively. Equation (3) does not sufficiently account for septal penetration and scattering from neighboring detectors for high-energy gamma-rays from SFA. In this work, we developed a custom MC code in C++ to simulate the collimator’s PSF by performing photon transport in the collimator. We simulated the PSF as a function of source position and energy while balancing accuracy and computational efficiency. A single collimator cell was discretized into 8 bins of 0.5 mm width along the X-axis, 21 bins along the Z-axis with an angular step size of $2^{\circ }$ covering the full opening angle, and 69 bins along the Y-axis with non-uniform step sizes (Fig. 6), totaling 11,592 source positions. The Y-axis step size was smaller near the collimator to account for more rapid variations. The choice of a 50 keV energy step from 650 to 1500 keV ensures a reasonable trade-off between resolution and computational cost, capturing energy-dependent variations in the PSF without excessive sampling. This step size is also compatible with the full energy peak width (20–30 keV) of the CZT detector within the energy range of interest⁶.We simulated the PSF for 1518 sources distributed in the Y-Z plane with energy ranging from 650 to 1500 keV in steps of 50 keV. 10⁸ photons were simulated for each source configuration. Two different detector responses were modeled:

• F8 tally: the photon energy deposited in the detector is recorded.

• F4 tally: the photon track length in the detector is recorded. This tally is used for validation purposes only and provides better statistics than F8 tally for the same number of photons.

Fig. 6

Positions of sources on the Y-Z plane where the PSF was simulated.

We used the analytical formula Eq. (3) and our MC code to calculate the PSF when a point source isotropically emitting 1.4 MeV gamma rays is placed 12.75cm away from the detector. This energy was selected because within our energy window of interest, photons of this energy can undergo several scattering reactions, which is more challenging to model. The detectors recorded counts in the 700-1500 keV energy window. The results are shown in Fig. 7. Compared to the ground truth (MCNP), the full width at half maximum (FWHM)FWHM of the PSF based on the analytical formula is larger. The long tail due to septal penetration and in-detector scattering is missing in the analytical PSF. In contrast, the PSF calculated by our custom MC code is in good agreement with the ground truth. It took approximately 30s for our MC code to simulate a PSF without parallelization, and 1.5 hours for MCNP with 16 cores.

Fig. 7

Comparison of PSFs obtained using the analytical formula, our custom MC code (NPS = $10^8$), and MCNP (NPS = $5\times 10^8$). The F4 tally is the average track length in the detector divided by its volume, which is an estimate of the average photon flux in the detector, while the F8 tally records the energy deposited by the photons in the detector volume.

Validation against MCNP simulation

We applied our forward model to simulate the sinogram for four different SFA configurations and compared the results to MCNP simulations. SFAs in Case 1 and Case 2 are VVER-1000 assemblies with all 331 positions filled (Fig. 8a). SFAs in Case 3 and Case 4 are VVER-440 assemblies with all 127 positions filled (Fig. 8b). All fuel pins have the same activity. The density, source emission spectrum, and detector energy window of each case are listed in Table 1. For Case 4, an emission spectrum (Fig. 9) displaying characteristic gamma lines of fission products such as ¹³⁷Cs, ¹³⁴Cs, and ¹⁵⁴Eu was adopted to simulate the gamma-ray emission of spent fuel. This spectrum was generated using SCALE 6.1 by modeling a Boiling Water Reactor $8\times 8$ fuel assembly with a burnup of 45 GWd/t and a cooling period of 8 years²¹. The energy window was set to 650-700 keV to capture the dominant ¹³⁷Cs gamma rays in the emission spectrum. MCNP simulations of the exact same configurations were performed to obtain the ground truth, and the NPS was set so that the error associated with the sinogram counts is below 5%.

Fig. 8

Simulated SFAs with all positions filled with fuel pins. (a): Case 1 and 2. (b) Case 3 and 4.

Table 1 Simulation parameters for the four SFA configurations.

Fig. 9

Gamma emission spectrum for a fuel assembly with 8-year cooling and 45 GWd/t burnup²¹. The ¹³⁷Cs peak appears broader than the others due to the uneven energy bin width.

The simulated sinograms and the absolute difference between the MCNP and our forward model are shown in Fig. 10. Overall, the simulated sinograms are in good agreement both visually and quantitatively for all cases. The difference is the most significant at projection angles where the diagonal fuel pins are aligned and the attenuation is the greatest. The relative Root-Mean-Sqaured-Difference (RMSE) between two sinograms is calculated:

$$\begin{aligned} \text {RMSE} = \frac{\sqrt{\frac{1}{MN}\sum _{i=1}^{N}\sum _{j=1}^{M}(S_{ij}^{MCNP} - S_{ij}^{CFD})^2}}{\frac{1}{MN}\sum _{i=1}^{N}\sum _{j=1}^{M}S_{ij}^{MCNP}}, \end{aligned}$$

(4)

where M and N represent the number of detectors and angles, respectively. Table 2 shows the RMSE and time performance for each case. The RMSE is below 5% for Case 1, 3, and 4. In Case 2, the self-shielding was the heaviest due to the high fuel density and large number of rods, and the probability for photons emitted by the inner fuel pins to reach the detectors was exceptionally low, leading to low statistics for these inner pins in MCNP simulation. In terms of computational time, our forward model is over 100,000 times faster than MCNP simulation. It took approximately 15 seconds for our model to simulate a sinogram of VVER-440 assembly and 45 seconds for a VVER-1000 assembly, on an i9-7920X CPU @ 2.90GHz with 20 cores in parallel.

Fig. 10

Comparison of sinograms simulated by MCNP and our CFD-based forward model. Colorbar: counts normalized to NPS = 10⁹.

Table 2 Comparison of the CFD-based forward model and MCNP simulation.

Image reconstruction and activity quantification

In this section, we applied our forward model to reconstruct images from both simulated and measured datasets. Our reconstruction approach is based on a linear-inverse method, which has been shown to provide superior image quality compared to FBP. The reconstructed fuel activity was compared to the ground truth, and a 2.3% relative uncertainty was achieved.

Sinogram generation

Two sets of data were prepared: simulated sinograms generated by MCNP for algorithm validation and measured sinograms for demonstrating applicability to real data. The details of the datasets are as follows:

Simulated sinograms

Figure 11 shows the MCNP-simulated sinograms for the four SFA configurations described in the previous section. These sinograms were input to our image reconstruction algorithm. The reconstructed images and activity distributions were compared with the ground truth to assess the accuracy of the reconstruction algorithm.

Fig. 11

MCNP-simulated sinograms of four SFA configurations. Colorbar: counts normalized to NPS = 10⁹.

Measured sinograms

A measurement campaign was conducted in Loviisa Nuclear Power Plant, Finland in 2021. Sinograms of five VVER-440 SFAs with different burnups and cooling periods were acquired, as shown in Fig. 12. Counts in the 650–700 keV energy window were recorded. At the center of each SFA, a water channel was present instead of a fuel pin, which should exhibit the lowest activity. The separation, defined as the difference between the activity of two pin types, is a measure of the algorithm’s ability to distinguish between the water channel and fuel pins¹⁵, and a higher separation indicates a better reconstruction quality.

Fig. 12

Measured sinograms of five VVER-440 SFAs. Colorbar: counts.

Image reconstruction method

Our reconstruction algorithm consists of the following steps:

• 1. Coarse FBP reconstruction was performed, from which the SFA geometry was extracted to build a material map.

• 2. The material map was fed to our forward model to calculate the system response matrix.

• 3. Enhanced reconstruction was performed by unfolding the measured singoram from the response matrix by solving an inverse problem.

Initial FBP reconstruction

First, we used the ASTRA Toolbox²² to perform the FBP reconstruction. Malfunctioning detectors shown as dark stripes were removed from the sinogram data before reconstruction. The FBP-reconstructed image allows us to determine the geometry of the SFA and the centers of fuel cells where fuel pins may be present⁹.

Response matrix simulation

The CFD-based forward model described in the previous section was employed to calculate the system response matrix. We first created a material map based on the FBP-reconstructed fuel pin positions. The diameter of the fuel pin, if not available, can be estimated based on the FWHM of the pin on the reconstruction⁹. We assigned a uniform source of unit intensity to each pin and applied the forward model to calculate the sinogram. The response matrix was constructed by stacking the sinograms for all fuel pins.

Gain correction for measured data

The bright/dark stripes in the measured sinogram (Fig. 12) were caused by the nonuniform gains of the detectors, which was not taken into account in our forward model. A gain correction procedure is therefore needed before performing the reconstruction. We adopted the gain correction algorithm in⁷, where it’s assumed that the mean projection averaged over all projection angles is a slowly varying function of detector position, regardless of the SFA configuration. We summed all the sinograms in the response matrix and calculated the mean projection. We then fitted the simulated mean projection to the measured one, and the ratio of the two is the gain correction factor.

FBP and inverse reconstruction

In this step, two different approaches were used to reconstruct the fuel pin activity: FBP and the linear-inverse approach. It should be noted that for comparison between simulation and reconstruction, we expressed the activity in terms of NPS (number of photons simulated/emitted). The activity in becquerel can be derived from the branching factor and acquisition time.

The FBP reconstruction was performed on the gain-corrected sinogram using the ASTRA Toolbox. We also applied the Chang’s method to compensate for self-shielding²³. After correction, we obtained a pixelated image of the SFA, with the pixel value proportional to the activity. By summing the pixels in each fuel pin, we obtained the relative activity of the fuel pins. Absolute activity quantification with this approach was not possible due to the approximations introduced in FBP.

The inverse approach relies on the linear relationship between the sinogram $\textbf{y}$ and fuel pin activity map $\textbf{x}$:

$$\begin{aligned} \begin{aligned} \textbf{y}= \textbf{A}\textbf{x}+ \textbf{n}, \end{aligned} \end{aligned}$$

(5)

where $\textbf{A}$ is the system response matrix, with each column being the sinogram in response to one fuel pin, and $\textbf{n}$ represents the noise. Knowing the response matrix $\textbf{A}$ and the sinogram $\textbf{y}$, we solved for $\textbf{x}$ using the fast iterative shrinkage thresholding algorithm (FISTA)²⁴. The algorithm converged within 200 iterations in 20 second, due to the small problem size. With this approach, we obtain the absolute activity of each fuel pin directly, without the intermediate step of summing pixel values as in FBP.

Reconstruction results

In this section, we present the reconstruction results of both simulated datasets and measured datasets.

Reconstruction on simulated dataset

Figure 13 shows the FBP-reconstructed images without and with Chang’s attenuation correction. The FBP images without attenuation correction show a significant amount of artifacts. The Chang’s attenuation correction improved the image contrast, especially at the center of the assembly. The corrected images were used to calculate the fuel pin activity.

Fig. 13

Comparison of the FBP images without and with Chang’s attenuation correction. Colorbar: relative activity in each pixel.

Figure 14 shows the fuel pin activity maps based on FBP and the linear-inverse reconstructions, as well as the ground truth. Since absolute activity was not obtainable from FBP, we applied a scaling factor to the FBP reconstruction such that the mean of the activities equals the ground truth. The activity reconstructed by the linear-inverse approach is plotted as is with no corrections. From Fig. 14, the FBP-based activity map displays a strong variation, especially at the center of the assembly, even with attenuation correction applied. In contrast, the inverse approach resulted in a more uniform activity map, in close agreement with the ground truth. This can also be observed in Fig. 15, which shows the histogram of pin activities. The ground truth is shown as the red dashed line. We observed that the linear-inverse approach produced a much narrower distribution than FBP in all cases. The relative uncertainty $\sigma$ of the fuel pin activity distribution, defined as the standard deviation of the distribution divided by the mean, was calculated for each case. The linear-inverse approach achieved a 2.3% uncertainty on average, as opposed to the 10.1% uncertainty with FBP.

Fig. 14

Comparison of fuel pin activity maps. First row: ground truth. Second row: reconstructed using FBP with attenuation correction applied. Third row: reconstructed using the linear-inverse approach. In the ground truth, all fuel pins have the same activity. Colorbar: number of photons emitted (absolute activity).

Fig. 15

Comparison of histograms of reconstructed fuel pin activities. First row: FBP with attenuation correction. Second row: Inverse approach. The red line represents the ground truth. $\sigma$ is the relative standard deviation of the activity distribution with respect to the mean.

Reconstruction on measured dataset

Figure 16 shows the FBP reconstructions without and with attenuation correction. The images with correction applied display a higher contrast, and were used to estimate the fuel pin activity.

Fig. 16

Comparison of the FBP images without and with Chang’s attenuation correction. Colorbar: relative activity in each pixel.

We compared the pin activity maps in Fig. 17. It is known that at the center of the VVER-440 assembly is a water channel where no fuel is present, hence the activity should be the lowest. This is accurately reproduced in the inverse reconstructions but not in FBP. Figure 18 shows the histogram of reconstructed pin activities, where the green line represents the activity at the water channel and the red line represents the mean activity of present fuel pins. The separation between the water channel and fuel pins is more pronounced in the inverse approach than in FBP, and exceeds 3-$\sigma$ for all measured cases (Fig. 19). Our separation is also consistent across different assemblies, indicating a robust performance against variations of detector gains and SFA properties (burnup and cooling time).

Fig. 17

Comparison of fuel pin activity maps. First row: reconstructed using FBP with attenuation correction applied. Second row: reconstructed using the linear-inverse approach. Colorbar: number of photons emitted (absolute activity).

Fig. 18

Comparison of histograms of reconstructed fuel pin activities. First row: FBP with attenuation correction. Second row: Inverse approach. Green line represents the activity at the water channel with the ground truth value of 0. The red line represents the mean activity of other fuel pins. $\sigma$ is the relative standard deviation of the activity distribution with respect to the mean.

Fig. 19

Separation between reconstructed activities at the central water channel and fuel pins. Different markers represent different assemblies, and different colors represent different reconstruction methods.

Discussion

It’s known that the strong self-shielding by high-density fuel materials must be corrected to ensure accurate activity quantification. However, the impact of Compton scattering on the reconstruction quality has not been well characterized. With our forward model, we calculated the counts without including Compton scattering in the response matrix, and compared them to the ones with scattering, shown in Fig. 20. The results show that the Compton scattering effect is significant for Cases 1-3, especially for the inner fuel pins. Excluding Compton scattering resulted in up to 50% underestimation of detector counts in the response matrix, which in turn led to an overestimation of pin activity up to 70% (Fig. 21). This indicates that the Compton scattering effect should be taken into account in the forward model to achieve accurate activity quantification. In the context of fuel diversion detection, this overestimation is particularly concerning because it can potentially result in missing fuel pins being falsely identified as present. Since the reconstructed activity levels would appear higher than they actually are, a diversion scenario-where fuel pins have been removed-could go undetected. This highlights the necessity of accurately modeling scattering effects to ensure reliable and sensitive fuel verification. The impact of Compton scattering on the reconstruction quality depends on the gamma emission spectrum of the fuel and the energy window used. For Case 4, with a narrow energy window around the 661 keV peak and a negligible amount of higher energy photons emitted by the fuel, the Compton scattering effect is of the order of a few percent.

Fig. 20

Decrease in detector response when no scattering was included in the response matrix. Colorbar: Counts(no scattering)/Counts(with scattering)-1.

Fig. 21

Increase in reconstructed fuel pin activity when no scattering was included in the response matrix. Colorbar: Activity (no scattering)/Activity (with scattering)-1.

An important input to our forward model is the gamma-ray emission spectrum of the fuel assembly, from which the photon’s energy is sampled. This spectrum can be inferred based on the burnup and cooling-time of the assembly using existing software such as ORIGEN²⁵. In this work, we assumed an emission spectrum dominated by the ¹³⁷Cs gamma line for all measured fuel assemblies, due to the lack of burnup information. We also assumed uniform burnup across all fuel rods, leading to a consistent gamma-ray emission spectrum for different rods and simplifying the simulation. However, in fuel assemblies with significant burnup variation, this assumption may introduce biases in the reconstructed activity distribution and potentially affect the identification of individual fuel rods. To account for this, a set of rod-dependent gamma-ray emission spectra can be incorporated into the simulation model to refine the response matrix and improve accuracy.

In practice, it’s desirable to minimize the reliance on pre-existing information to avoid introducing unwanted bias to the verification results. Instead of relying on burnup data, a more robust approach for reconstructing the gamma-ray emission spectrum involves using the high-resolution spectrum measured by the CZT detector. Similar to the presented linear-inverse approach, we can simulate an energy-dependent detector response matrix with our forward model and reconstruct the unknown gamma-ray source emission spectrum by unfolding the measured CZT spectrum from the response matrix. By reconstructing the gamma-emission spectrum directly from measured data using the response matrix, we eliminate the need to assume or simulate the spectrum based on reactor-specific burnup histories. This approach enhances the objectivity of the verification process, as it minimizes dependence on operator-declared burnup data, which may be subject to uncertainties or biases. Instead, the reconstructed spectrum will be derived from fundamental physical interactions and measurements, offering a more independent means of fuel characterization. Further investigation into this approach will be a key focus of our future work.

In the gain correction step, we assumed that detector non-uniformity is angle-independent, an assumption supported by several experimental factors. First, the measurement system operates underwater with well-controlled environmental conditions, such as temperature, minimizing temporal variations in detector response. Second, before the measurement, we calibrated the detectors using the spent-fuel assembly as a strong ¹³⁷Cs source, ensuring consistent energy calibration⁷. Finally, while the gamma-ray energy spectra recorded at individual measurement angles were not saved for intercomparison, the summed spectrum over all angles was recorded and exhibits well-resolved characteristic peaks with good resolution¹². This indicates stable detector performance across all measurement angles, reinforcing the validity of the angle-independent assumption.

The assemblies studied in this work are VVER-440 types with hexagonal geometries, but the methodology is applicable to assemblies of arbitrary shapes and burnups. The assembly geometry can be initially approximated through the FBP reconstruction and then used as input for generating the response matrix. This enables the method to accommodate different fuel types, including PWR and BWR assemblies, without requiring prior knowledge of their geometry. Additionally, the model can handle fuels with varying burnups, provided that a gamma-ray emission spectrum is supplied as input. This spectrum can be determined either through burnup calculations or the proposed energy-dependent response matrix.

The distribution of reconstructed fuel pin activities in Fig. 18 resembles a Gaussian distribution. We compared the distribution to a Gaussian distribution in the quantile-quantile (Q–Q) plot in Fig. 22, and the two distributions are in good agreement, except the outlier that resulted from the low-activity water channel. The Q–Q plot thus provides another way of localizing the water channel by identifying the outlier on the plot. We also performed D’Agostino-Pearson normality test²⁶ for all measured cases, and the p-values are below 10^–4, suggesting that the null hypothesis be rejected. With the $\sigma$ deduced from the data, the probability that the activity of a fuel pin falls below 3$\sigma$ is 0.135%. For a VVER-440 assembly that encompasses 126 fuel pins, there is a chance that 0.17 fuel pins are below 3$\sigma$. The chance of non-observation of the central water channel is thus small.

Fig. 22

Quantile-quantile plot comparing the distribution of reconstructed activities to the Gaussian distribution. The points represent the activity distribution, while the red line represents the Gaussian distribution.

Conclusions

In conclusion, we have developed a linear inverse approach for reconstructing images of SFAs from both simulated and measured datasets. This approach relies on a fast and accurate forward model of the PGET system to enable near-real-time computation of the system response matrix. A key component of the forward model was the PSF of the collimator, which was calculated using a custom MC code to account for septal penetration and in-detector scattering. Our refined PSF closely aligns with the results from MCNP simulations, outperforming analytical calculations. The linear forward model was validated by comparing it to MCNP simulation, and a 3.7% relative difference in detector counts was achieved, which is within the uncertainty of the MCNP simulation. Reconstruction on simulated datasets demonstrated that our reconstruction algorithm can reliably reconstruct the fuel pin activity with a 2.3% relative uncertainty. Furthermore, when applied to measured datasets, our algorithm was able to identify both the fuel pins and water channels within the spent fuel assembly with high confidence.