Progress on Associate-Particle Imaging Algorithms, 2020

ORNL/SPR-2020/1766 Progress on Associate-Particle Imaging Algorithms, 2020 P. A. Hausladen M. A. Blackston A. J. Gilbert J. Gregor J. K. Mattingly November 2020 Approved for public release. Distribution is unlimited. DOCUMENT AVAILABILITY Reports produced after January 1, 1996, are generally available free via US Department of Energy (DOE) SciTech Connect. Website www.osti.gov Reports produced before January 1, 1996, may be purchased by members of the public from the following source: National Technical Information Service 5285 Port Royal Road Springfield, VA 22161 Telephone 703-605-6000 (1-800-553-6847) TDD 703-487-4639 Fax 703-605-6900 E-mail info@ntis.gov Website http://classic.ntis.gov/ Reports are available to DOE employees, DOE contractors, Energy Technology Data Exchange representatives, and International Nuclear Information System representatives from the following source: Office of Scientific and Technical Information PO Box 62 Oak Ridge, TN 37831 Telephone 865-576-8401 Fax 865-576-5728 E-mail reports@osti.gov Website http://www.osti.gov/contact.html This report was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise, does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof. ORNL/SPR-2020/1766 Physics Division PROGRESS ON ASSOCIATE-PARTICLE IMAGING ALGORITHMS, 2020 P. A. Hausladen M. A. Blackston A. J. Gilbert J. Gregor J. K. Mattingly November 2020 Prepared by OAK RIDGE NATIONAL LABORATORY Oak Ridge, TN 37831-6283 managed by UT-BATTELLE, LLC for the US DEPARTMENT OF ENERGY under contract DE-AC05-00OR22725 CONTENTS LIST OF FIGURES .......................................................................................................................................v LIST OF TABLES ....................................................................................................................................... vi ACRONYMS ............................................................................................................................................. viii ABSTRACT ...................................................................................................................................................1 1. INTRODUCTION .................................................................................................................................1 2. PROGRESS ON ITERATIVE GAMMA-RAY TOF RECONSTRUCTION ......................................2 2.1 RECONSTRUCTION FRAMEWORK .......................................................................................3 2.2 GEOMETRY AND PHYSICS MODELING ..............................................................................4 2.3 PRELIMINARY RESULTS ........................................................................................................5 2.4 SUMMARY .................................................................................................................................7 3. PROGRESS ON PROJECTION OF API OBSERVABLES INTO A MATERIALS BASIS ..............7 3.1 METHODS ..................................................................................................................................8 3.1.1 Physics model .................................................................................................................8 3.1.2 The inverse problem .....................................................................................................10 3.1.3 Assumptions and limitations .........................................................................................11 3.2 MATERIAL RECONSTRUCTION ALGORITHM APPLICATION STUDY .......................11 3.2.1 Gap determination .........................................................................................................12 3.2.2 Quantifying low-Z material behind high-Z material.....................................................13 3.2.3 Discrimination between materials of similar Z .............................................................13 3.3 ADDITIONAL API PHYSICS MODELING............................................................................14 3.3.1 Modeling gamma emission ...........................................................................................15 3.3.2 Time spectrum as a material signature..........................................................................16 3.4 OUTLOOK FOR MATERIAL IDENTIFICATION .................................................................17 4. PROGRESS ON MATERIAL IDENTIFICATION USING FULL SPECTRAL ANALYSIS ..........18 4.1 INTRODUCTION .....................................................................................................................18 4.2 FY 2020 PROGRESS ................................................................................................................19 4.3 FY 2021 PLANS ........................................................................................................................30 5. ACKNOWLEDGEMENTS .................................................................................................................30 6. REFERENCES ....................................................................................................................................30 iii LIST OF FIGURES Figure 1. Example 3D reconstruction of gamma rays (in debugging phase). ................................................6 Figure 2. Projection data for reconstructed images........................................................................................7 Figure 3. Representative cross sections for (a) photons and (b) neutrons. ....................................................8 Figure 4. Example neutron and x-ray incident flux. ......................................................................................9 Figure 5. Example input (a) radiography data and (b) algorithm output. ....................................................10 Figure 6. Example simulated radiograph and reconstructed material areal densities for the two annuli...............................................................................................................................................12 Figure 7. Simulated (a) x-ray and (b) neutron radiographs of the AT400R nuclear storage container in which the top sphere is Pu and the bottom sphere is W..............................................14 Figure 8. Reconstruction outputs using radiography data from (a) x-rays only and (b) both x-rays and neutrons. ...................................................................................................................................14 Figure 9. (a) A schematic diagram of neutron interrogation producing gamma radiation, and (b) a comparison of the output of the analytical model given in Equation (3.3) and MCNP simulation for a small cube of Fe. ...................................................................................................15 Figure 10. (a) Geometry of detector, objects, and source for proof-of-concept model, and (b) time spectra for (black) full detector average, (blue) an 8 × 8 region of interest centered on the object on the left, and (dashed) an 8 × 8 ROI centered on the object on the right. ........................17 Figure 11. Illustration of a 1D voxel transfer function H relating emergent current 𝑱𝒐𝒖𝒕to incident current 𝑱𝒊𝒏.......................................................................................................................................18 Figure 12. Coupling of multiple voxel transfer functions. ...........................................................................19 Figure 13. Neutron and photon trajectories emerging from a 1D, 1 cm thick slab of polyethylene (left), Fe (center), and HEU metal (right) for a monodirectional source of 14 MeV neutrons incident on the left face of each slab. ...............................................................................20 Figure 14. Histograms of position, direction, energy, and time-delay for neutrons emergent from a 1D slab of polyethylene in response to a monodirectional source of 14 MeV incident neutrons. ..........................................................................................................................................21 Figure 15. Histograms of position, direction, energy, and time-delay for neutrons emergent from a 1D slab of Fe in response to a monodirectional source of 14 MeV incident neutrons. ..................22 Figure 16. Histograms of position, direction, energy, and time-delay for neutrons emergent from a 1D slab of HEU metal in response to a monodirectional source of 14 MeV incident neutrons. “Transmitted” neutrons emerge from the side of the voxel opposite the source, and “reflected” neutrons emerge from the same side as the source; uncollided neutrons (those that passed through the slab without interaction) are omitted. The “ridges” labeled in the histograms of emergent energy versus direction result from (a) elastic scatter by uranium. ..........................................................................................................................................23 Figure 17. Histograms of energy, direction, and time-delay for neutrons emergent from a 1D slab of (a) polyethylene, (b) Fe, and (c) HEU metal in response to a monodirectional source of 14 MeV incident neutrons. ..............................................................................................................24 Figure 18. Distribution of source neutron energy and direction used to generate 1D voxel transfer functions. .........................................................................................................................................25 Figure 19. Transfer functions relating emergent and incident energy for 1D slabs of (a) polyethylene, (b) Fe, and (c) HEU metal. .......................................................................................27 Figure 20. Distribution of energy and direction for the test case using an isotropic fission spectrum neutron source. ................................................................................................................28 Figure 21. Comparison of transmitted and reflected neutron current calculated directly using (black) MCNP to (red) the current predicted using transfer functions for 1D voxels of (a) polyethylene, (b) iron, and (c) HEU metal......................................................................................29 v LIST OF TABLES Table 1. Summary of quantitative gap determination using the reconstruction algorithm output...............13 vi ACRONYMS 1D 2D 3D API D-T DU FY HEU MCNP MVp ORNL PNNL SIRT TOF TV WLS One-dimensional Two-dimensional Three-dimensional associated-particle imaging deuterium-tritium depleted uranium fiscal year highly enriched uranium monte-carlo n particle megavolt peak Oak Ridge National Laboratory Pacific Northwest National Laboratory simultaneous iterative reconstruction technique time of flight total variation weighted least squares viii ABSTRACT The present work describes the progress on developing imaging algorithms that use fast neutron signatures acquired using the associated-particle imaging (API) method. The present work complements ongoing work to develop neutron source and detector hardware to enable field inspection by investigating algorithms that are capable of discriminating among critical materials or extracting three-dimensional (3D) geometrical information from single-sided or transmission measurements. The present work is divided into three approaches: 1. Iterative reconstruction of inelastic gamma-ray emissions to perform 3D time-of-flight (TOF) imaging in a single view in either transmission or backscatter configurations. Iterative reconstruction enables image resolution better than the inherent TOF resolution. 2. Decomposition of registered neutron and x-ray radiographs into an assumed material list for each pixel in the image. 3. Material identification using full spectral analysis that includes the emergent neutron and gamma ray energies, times, and angles. Progress for each approach is summarized for fiscal year (FY) 2020. 1. INTRODUCTION For the past decade, Oak Ridge National Laboratory (ORNL) has developed 3D tomographic imaging techniques using fast neutrons from a deuterium–tritium (D-T) neutron generator produced via the 𝑑 + 𝑡 → 𝛼 + 𝑛 reaction. This development has focused on the API technique, in which the time and location of detected alpha particles (ascertained by a detector embedded in the D-T neutron generator) determine the time and direction of the associated 14.1 MeV neutrons. Using API techniques, transmission imaging with excellent contrast using a wide-cone beam can be achieved. Induced-reaction imaging techniques that can reconstruct images of the probability of induced fission and hydrogen elastic scattering are also possible. In these induced-reaction imaging modalities, the lines of response are determined by the initial directions of the interrogating neutrons as determined by their associated alpha pixels. Progress has also been made on developing imaging techniques whose contrast originates from small-angle scattering and inelastic neutron scattering (reconstructing based on either detected gammas or on detected neutrons). To date, effort for developing imaging techniques has concentrated on tomographic measurements using laboratory imaging systems with gantries that position the source and detector with respect to each other. For these imaging techniques, image reconstruction does not correspond to uniquely identifying materials and densities for each voxel in the image. Rather, image reconstruction corresponds to reconstructing a single parameter for each voxel of the image. In transmission images, that parameter is the attenuation coefficient, which can be interpreted to be the total probability of a neutron interacting per centimeter. For induced-fission imaging, the parameter is the number of neutron pairs produced per centimeter per incident neutron. A further complication of reconstructing induced-reaction images is that inferring physical values for voxels in the image requires knowledge of the reconstructed 3D transmission through the volume to correct for efficiency. The present work is part of an effort to transition imaging using D-T neutrons from laboratory use to field use to enable a portable inspection capability. In this use case, the source and detector will be handpositioned, and imaging will be limited to a single or few views. In some instances, transmission imaging may not be possible. Other efforts address the source and detector hardware to reduce their size, power, 1 and complexity to the operator. The present work addresses algorithm development to extract maximum information from the available API signatures in the context of single or few-view measurements. This algorithm work can be categorized into the following two general themes. 1. Algorithms that directly calculate images from data with a minimum of assumptions to give operators the best estimate of a 3D object geometry from a single measurement. In this case, images are voxel maps in which each voxel value corresponds to the probability of a given interaction and may be calculated from data in either transmission or backscatter measurement geometries. 2. Algorithms that directly reconstruct for materials and densities (rather than a parameter that encodes the probability of a particular reaction) but may have constraints on the materials or geometry to make the optimization problem tractable. The present work investigates the following three distinct approaches that address these themes: • Iterative TOF reconstruction of inelastic gamma-ray emissions. This approach reconstructs the probability of inducing inelastic gamma rays per incident neutron in each voxel of an image. Because the speeds of 14 MeV neutrons and gamma rays are known, performing 3D TOF reconstruction in a single view in either transmission or backscatter configurations is possible. Iterative reconstruction enables image resolution better than the inherent TOF resolution. • Decomposition of x-ray and neutron transmission into a materials basis. This approach uses transmission imaging of D-T neutrons and high-energy x-rays to reconstruct a best-fit decomposition to an assumed material list for each pixel in registered neutron and x-ray radiographs. Including additional API observables is planned for future efforts. • Material identification using full spectral analysis. The goal of this approach is to make a best estimate of the materials in an inspected object of known geometry by simultaneously fitting all possible API signatures. The initial thrust of this effort is to identify a means of implementing a forward model that is sufficiently fast and calculates all of the desired signatures. This document reports progress on the development of these three approaches to analysis of API imaging data for FY 2020. This report also satisfies the FY 2020 deliverables, “Progress Report on Inelastic Gamma TOF Reconstruction, FY20,” “Progress Report on Material Identification Using Full Spectral Analysis, FY20,” and “Progress Report on Projection of API Observables into Materials Basis, FY20” of the project “Associated-Particle Imaging Algorithms for Material Identification” supported by the US Department of Energy’s Office of Defense Nuclear Nonproliferation Research and Development. 2. PROGRESS ON ITERATIVE GAMMA-RAY TOF RECONSTRUCTION Imaging the distribution of inelastic-scattering gamma rays produced in an object by the passage of interrogating D-T neutrons is desirable because the known speeds of 14 MeV neutrons and gamma rays make 3D TOF reconstruction possible in a single view when the source and detectors are positioned in either a transmission or a backscatter configuration. This form of imaging has been investigated previously by a number of authors using spectroscopic gamma ray detectors primarily to quantify the C, N, and O in cargo to screen for the presence of explosives [1]. ORNL began investigating iterative reconstruction of TOF images of inelastic gamma rays induced by 14 MeV D-T neutrons in the “Fast TOF Correlation” project performed in collaboration with Lawrence Livermore National Laboratory. In this work, gamma ray detection was performed by fast organic scintillator neutron detectors used primarily for multiplicity counting, but which are sensitive to gamma 2 rays and allow fast timing. Significantly, this work demonstrated that iterative reconstruction using the knowledge of the system response could achieve better spatial resolution than the intrinsic resolution given by TOF [2]. However, after these initial investigations, the ORNL effort on the Fast TOF Correlation project shifted to concentrate on developing a fast-timing neutron imaging panel. Investigation of iterative reconstruction shifted to the present work. In other previous work, ORNL collaborated with the University of Tennessee, Knoxville, on the “3D Tomography and Image Processing Using Fast Neutrons” project to combine modern, parallel iterative reconstruction techniques with novel associated-particle, fast-neutron–induced reaction (induced fission and elastic scatter) and transmission imaging methods [3]. The resulting reconstruction code for novel imaging techniques achieved sufficient speed for practical 3D image reconstruction, approached favorable solutions using regularization and constraints, and weighted data by appropriate errors. The present work uses the building blocks of this existing 3D tomographic reconstruction code to implement 3D reconstruction of inelastic gamma ray emissions that is sufficiently fast and selects “good” images that have smooth sections with sharp edges. 2.1 RECONSTRUCTION FRAMEWORK To reconstruct smooth images while preserving sharp edges consistent with machined parts, image reconstruction is modeled by the following weighted least squares (WLS) minimization problem that is subjected to a total variation (TV) constraint: 𝒙∗ = arg min !" ‖𝑨𝒙 − 𝒚‖&% + 𝛽‖𝒙‖&& s. t. TV(𝒙) ≤ 𝜖 "#$ (2.1) Here, matrix 𝑨 = D𝑎'( F incorporates knowledge about the system geometry and the physics associated with the induced gamma reactions; matrix 𝑾 = diag[𝑤' ] implements statistical weighting; and vectors 𝒙 = D𝑥( F and 𝒚 = [𝑦' ], respectively, represent the image being reconstructed and the acquired projection data. Each equation in the linear system models what happens to neutrons as they travel in the direction of an alpha ray from the D-T source toward the detector array. TV refers to the L1 norm of the L2 norms of the image gradients given by 𝑇𝑉(𝒙) = P‖∇' 𝒙‖ . (2.2) ' The effect of constraining the TV value of the image is sparsification of the gradient magnitude image. This sparsification leads to solutions that exhibit smooth regions with well-defined edges. The WLS-TV problem is solved using a relaxed, incremental, proximal gradient scheme that consists of a two-step iteration [4]. First, the WLS term is minimized using a proximal mapping that keeps the image somewhat close to the one produced in the previous iteration. Next, the TV constraint is satisfied by mapping the image to the closest point on the surface of an L1 ball. This mapping is implemented using the Chambolle–Pock algorithm for solving convex optimization problems [5]. Using a less mathematical description, the WLS computation initially converges to an image that has relatively sharp edges, even if it appears grainy or otherwise noisy. After that, the closest image that satisfies the TV constraint is selected. The process is then repeated. In each iteration, the WLS computation is prevented from deviating too far from the image produced by the TV constraint and vice versa. After several of these two-step iterations, a solution emerges for which the WLS residual error is relatively small, and the TV value is close to or smaller than the given constraint. 3 The WLS minimization is implemented using a modified version of a preconditioned gradient decent method known as the simultaneous iterative reconstruction technique (SIRT). During the 𝑘)* iteration, the image is forward projected to produce model-based projection data. The difference between the measured and the forward-projected data is then backprojected to update the image using appropriate scaling to both ensure and accelerate convergence. Mathematically [6], V W 𝒙(,) − 𝛼𝑪 V 𝑨0 𝑾S𝑨𝒙(,) − 𝒚W, 𝒙(,-.) = S𝑰 − 𝛼𝛽𝑪 (2.3) V has the following values: where preconditioning matrix 𝑪 V ≡ diag Y 𝑪 1 \. ∑' 𝑤' (∑* 𝑎'* )𝑎'( (2.4) Parameters 𝛼 and 𝛽 control the behavior of the algorithm. Both are set automatically using mathematical bounds, but they can be overwritten by the user. Parameter 𝛼 controls the step size for the update. Parameter 𝛽 controls the degree to which minimization of the data fidelity term drives the computation versus the desire to regularize the solution using the minimum-norm Tikhonov term. Parameter 𝜖 represents an upper bound on the TV value. This is set automatically based on a user-defined full width at half maximum value for the image resolution. The user also has the ability to set this bound to an explicit value, should that be desired. A more detailed descriptions of this specific WLS-TV framework can be found elsewhere [7, 8]. Notably, those implementations exploit the concept of “ordered subsets” to accelerate convergence, which is not applicable here because data was only acquired for only a single projection. Future work will look into alternative acceleration methods as well as use of the ray-based WLS algorithm known as the algebraic reconstruction technique [4], which updates the image using the following iteration scheme: ∀𝑖 ∶ 𝒙(,-.) = 𝒙, − 𝑨1' 𝒙, − 𝑦' 𝑨, |𝑨' |&& + (𝑤' 𝛽)2. ' (2.5) where 𝑨' denotes the corresponding row of the matrix. The algebraic reconstruction technique update equation is not regularized. 2.2 GEOMETRY AND PHYSICS MODELING The neutron source is modeled as a set of alpha normal vectors 𝒏3 emanating from a single-point location 𝒑4 . The detectors are modeled as a set of point locations 𝒑5 and surface normal vectors 𝒏5 . The image voxels are modeled by their center-point locations 𝒑( . An arbitrary number of neutron source locations is supported. All alpha and detector variables are indexed relative to a specific neutron source location to allow for multiple projections (views), possibly based on different detector configurations. The authors omitted the source indexing here for the sake of simpler notation. Each acquired projection data entry 𝑦35) represents the number of gamma rays detected in a specified time window based on neutrons that were emitted in the direction of a particular alpha normal vector. Matrix 𝑨 models the number of neutrons emitted in the direction of an alpha normal vector times the probability that they interact with material at any one of the voxels in a manner that contributes to gamma rays being detected by each of the detectors within the mentioned windows. Because the alpha normal vectors do not line up with the voxel centers, image-space interpolation is used to obtain effective voxels 4 located at incremental distances from the source along the direction of each alpha normal vector. Thus, the indices of matrix elements 𝑎'( refer to 𝑖 ≡ (alpha direction, detector, time window) and 𝑗 ≡ (interpolation − based voxel location). At each incremental distance from the source, the probability of a gamma ray being detected is calculated as the product of the solid-angle Ω(5 , detector efficiency 𝜖5 , and the probability of the gamma ray reaching the detector within a given time window 𝑝(5) where the voxel subscript indicates that the variable depends on the location from which the gamma ray was emitted. Although most of the acquired projection data will represent true coincidences, some will represent accidental coincidences. The accidental coincidences are modeled as the expected rates calculated from the number of neutrons emitted times the number of gamma rays detected and the fraction of total scan time each time window corresponds to. Combined, forward projection can therefore be mathematically expressed as 𝑦35) = 𝑁3 P Ω(5 𝜖5 𝑝(5) 𝑥( + 𝑟35 , ( where 𝑟35 = 𝑁3 𝑁65 Δ𝑡 . 𝑇789: The direction of the emitted neutrons is associated with uncertainty, which is modeled by an alpha cone that has a Gaussian cross section. That is, the interpolated voxel center located at an incremental distance along an alpha normal vector mentioned above is actually the weighted sum of many interpolated voxel centers obtained as follows. Represented using spherical coordinates, 𝒑( = (𝑟, 𝜃, 𝜙) is tilted to produce 𝒑0( = (𝑟, 𝜃 + 𝑑𝜃, 𝜙), which is then spun around the former using a unit quaternion-based representation of rotation. The image is sampled at regular intervals during this process, each time using interpolation. The tilt determines the probability associated with the sampled voxels. Backprojection is the reverse of the forward projection. The above interpolations to compute the forward and backprojection are calculated on-the-fly at relatively high cost because of the need for upward of hundreds of interpolation operations for each incremental step along an alpha ray (i.e., fewer interpolations closer to the source, more interpolations further away as the alpha uncertainty cone opens up). 2.3 PRELIMINARY RESULTS Two example reconstructed images using the preliminary code are shown in the left and right columns of Figure 1. The brightest object that appears in both images is a depleted uranium (DU) cylinder. For one data set, the interior contains a DU ball. For the other, the interior is empty (i.e., no ball). Figure 1 shows 𝑥𝑦 and 𝑧𝑦 slices of the reconstructed image volumes. Although the DU ball is visible in the 𝑧𝑦 slice and the line profile through the middle thereof, it appears more like an isthmus in the 𝑥𝑦 slice. As desired, the DU ball is absent in the no-ball images and the line profile. The DU cylinder is cylindrical, but some 5 geometric distortion is apparent near the middle of both 𝑥𝑦 slices. The edges of the support cone take on artificially large values in several places. 0 0 60 -60 30 0 60 -30 60 0.01 0.008 0.006 0.004 0.002 0 20 40 60 80 100 120 140 No ball: WLS-TV image, Alpha pixels 1-512 0 60 -60 30 0 z coordinate (cm) 0.012 0 0 x coordinate (cm) Ball: Line plot for x=3cm, z=1.5cm 0.014 -60 0 Induced gammas per alpha per cm Induced gammas per alpha per cm z coordinate (cm) No ball: WLS-TV image, Alpha pixels 1-512 -60 y coordinate (cm) y coordinate (cm) y coordinate (cm) 0 60 -30 Ball: WLS-TV image, Alpha pixels 1-512 -60 y coordinate (cm) Ball: WLS-TV image, Alpha pixels 1-512 -60 160 180 y coordinate (cm) 60 x coordinate (cm) No ball: Line plot for x=3cm, z=1.5cm 0.014 0.012 0.01 0.008 0.006 0.004 0.002 0 0 20 40 60 80 100 120 140 160 180 y coordinate (cm) Figure 1. Example 3D reconstruction of gamma rays (in debugging phase). Each scan acquired data for 512 alpha rays and 72 neutron and gamma ray detectors. The image was set to be 90 × 180 ×180 voxels with dimensions of 0.67 cm on a side. SIRT was run for 100 iterations, which appeared to more than needed to reach convergence. Using a small four-core MacBook Pro, initialization took 32 s, establishing the TV constraint took 1,000 s, and running SIRT took 2,400 s. The total run time was thus about 60 min per reconstruction. Figure 2 shows the measured and forward projections for alpha pixel 239, which is oriented toward the DU ball. Good agreement exists between the two projections for both data sets. However, a nonlinear time shift can be seen for gamma rays that were induced near the source. This shift indicates the presence of either a geometry or a timing bug in the system model. 6 No ball: Measured projection, Alpha pixel 239 10 20 20 Detector number Detector number Ball: Measured projection, Alpha pixel 239 10 30 40 50 60 50 70 20 40 60 80 100 120 140 20 40 60 80 100 Time of flight (ns) Time of flight (ns) Ball: Forward projection, Alpha pixel 239 No ball: Forward projection, Alpha pixel 239 10 10 20 20 Detector number Detector number 40 60 70 30 40 50 60 120 140 120 140 30 40 50 60 70 70 40 60 80 100 120 140 20 80 100 Time of flight (ns) No ball: Measured vs forward projection comparison Measured Forward 500 400 300 200 100 0 20 60 Ball: Measured vs forward projection comparison 600 0 40 Time of flight (ns) Induced gammas per alpha per ns 20 Induced gammas per alpha per ns 30 40 60 80 100 120 140 600 Measured Forward 500 400 300 200 100 0 0 Time of flight (ns) 20 40 60 80 100 120 140 Time of flight (ns) Figure 2. Projection data for reconstructed images. Although the ball appears more like an isthmus in the upper image (source-to-detector plane), the ball is clearly visible in the transaxial plane (intersected by the alpha rays). Changing the voxels sizes and the sampling rates for the image-space interpolation had little effect. 2.4 SUMMARY A preexisting code for 3D neutron tomography was repurposed for the inelastic gamma scatter imaging. To implement these changes, the system geometry and system matrix code were rewritten entirely because no overlap existed. Several data structures were reimplemented to support indexing based on alpha direction, detector, and time window. The SIRT-based WLS solver was adapted to match the new indexing scheme. File I/O was also extended to include handling data in the hierarchical data format version 5 (HDF5) file format. The code was not written with computational efficiency in mind, but several methods can make it run faster. Some of these methods will be explored in the upcoming year, including a more efficient way of handling the alpha uncertainty cone. The possibility of accelerating SIRT in a manner similar to ordered subsets will also be investigated. However, extensive validation and debugging should be completed first to ensure overall correctness of the approach. Statistical weighting will be studied to see if the reconstruction quality improves when higher-count alpha rays carry more weight than lower-count rays (i.e., for the data shown above, a 50× difference between the minimum and maximum counts existed). 3. PROGRESS ON PROJECTION OF API OBSERVABLES INTO A MATERIALS BASIS Pacific Northwest National Laboratory (PNNL) previously developed a method to perform quantitative material reconstructions using radiography data for x-rays [9] and for x-rays and neutrons [10]. In this methodology, each pixel of the transmission image is decomposed into the optimal quantities of “basis materials,” in which the number of basis materials equals the number of distinct interrogations. This method requires an accurate understanding of the response of the radiography system and combines a physics-informed calculation of the expected response with inverse methods to quantitatively reconstruct the areal densities of a short list of known materials for the inspected object. The method also requires appropriate basis materials to be known. 7 Recent interest in using neutron radiography to noninvasively interrogate objects provided an opportunity to adapt this algorithm to combined neutron and x-ray radiography data. Data acquired using API systems are not only of particular interest because the API method enables reduced neutron scattering on the imaging plane for transmission measurements, but also because observation of additional, inducedreaction signatures such as secondary gamma emission and neutron scattering allow new forms of contrast [11]. Therefore, PNNL seeks to further develop this method to include additional observables from API systems to enable high-accuracy material reconstructions. 3.1 METHODS 3.1.1 Physics model For an inverse algorithm to be able to estimate a material composition from observed radiography data, a physical model must be provided that describes how the radiography systems respond to a given material composition. The present physics model includes exponential attenuation through the object, the energy distribution of the source, and the energy-dependent efficiency of the detector. This baseline model describes how a set of material densities will attenuate a beam of x-rays or neutrons between the particle source and a radiographic detector pixel. The baseline model for a given pixel ij and interrogation k, rijk, is defined as 𝑟'(, Sρ w⃗'( W = y 𝜀, (𝐸) 𝜙",, (𝐸) expD−µ w⃗(𝐸) ∙ wρ⃗'( F dE. (3.1) < In Equation (3.1), the response 𝑟 is calculated as a function of the incident spectral flux 𝜙" , the energydependent detector efficiency 𝜀, the attenuation coefficient µ, and the areal density 𝜌'( of each material along the path from the source to the detector pixel located at (𝑖, 𝑗). The attenuation and areal density are written as vectors 𝜇⃗ and wρ⃗'( because each has a length equaling the number of basis materials L. The index k denotes a given interrogation, such as x-ray or neutron. Notably, the integration over energy for each interrogation indicates that the radiographic detectors cannot determine particle energy, but instead, the integrated response increases linearly with both increasing particle energy and particle flux, which is typical of digital radiographic detectors. Example material attenuation coefficients for x-rays and neutrons are shown as a function of energy in Figure 3. (a) 1 10 2 (g/cm2 ) µ (g/cm2 ) 10 (b) steel poly Pb Pu 100 2 4 6 8 steel poly Pb Pu 10 1 10 2 2 10 4 6 8 10 12 14 energy (MeV) energy (MeV) Figure 3. Representative cross sections for (a) photons and (b) neutrons. 8 16 18 The attenuation coefficients currently used are based on tabulated cross sections from available compilations, though they could be modified to reflect experimentally observed attenuation that includes in-scattering. In Figure 4, example source energy spectra are shown for 6 and 9 megavolt peak (MVp) energy x-rays and 14.1 MeV D-T neutrons. Note that the megavolt peak value specifies the accelerated electron energy and corresponds to the maximum (or endpoint) energy of x-rays generated by the Bremsstrahlung x-ray source. In contrast to the monoenergetic D-T neutron source, the Bremsstrahlung x-ray spectra are very broad in energy. The unique way in which materials are attenuated by an interrogated object and detected by the imaging system allows discrimination and reconstruction of the material composition from radiography data. The x-ray spectra have additionally been filtered by 2 cm of Pb to remove low-energy flux that deleteriously affects material identification. Neutron radiography data can be particularly useful as an addition to x-ray radiography data because neutrons and x-rays interact with materials in unique ways, putting different signatures in the radiography data that can be used in the material reconstructions. For example, hydrogenous materials are particularly opaque to neutrons compared to x-rays. 6 MVp X-ray 9 MVp X-ray DT neutron 103 102 Scaled 0 (particles/cm2 -s) 104 101 100 0 2 4 6 8 10 12 14 energy (MeV) Figure 4. Example neutron and x-ray incident flux. A matrix for each of the interrogations k is organized according to the relative pixel locations: 𝑟.. 𝑟.& 𝑟&. 𝑟&& 𝐑, = ‚ ⋮ 𝑟?. 𝑟?& 𝑟.> ⋯ 𝑟 &> ⋱ ⋮ † ⋯ 𝑟?> To account for a finite spatial resolution of the imaging systems, convolution matrices are used in the physics model to represent this effect. Given the two-dimensional (2D) matrix Rk, a convolution matrix C can be applied to model the spatial blurring. The present work assumes a 2D Gaussian kernel, which is separable, so it can be applied in a single dimension at a time. The image data model for interrogation k is then given by 1 𝐃, = 𝐂., 𝐑 , 𝐂&, , where the 1 and 2 indices on C denote the two dimensions of application of the convolution. For the current image model, C1k = C2k. As a final step before use in the inverse problem, the image data matrices are vectorized into an IJK long vector according to 9 𝑑⃗(ρ w⃗) = 𝑣𝑒𝑐(𝐃. , 𝐃& , ⋯ , 𝐃, ) = D𝑟... ⋯ 𝑟?.. , 𝑟.&. ⋯ 𝑟?&. , ⋯ , 𝑟?>. ⋯ 𝑟?>@ F, where the vector ρ w⃗ is of length IJL, comprising all materials L at all image pixels IJ. 3.1.2 The inverse problem Given a model of how the radiography systems respond to a set of materials, the inverse problem seeks to take an observation and the physics model to estimate the best-fit material composition that produced the observed data. The baseline inverse problem is stated as a least-squares fit of the data model vector 𝑑⃗(ρ w⃗) to the equivalent observed data vector 𝑑⃗AB4 with TV regularization. & 𝜌⃗∗ = arg min !" Œ𝑑⃗(ρ w⃗) − 𝑑⃗AB4 Œ% + 𝛼 TV(𝜌⃗) (3.2) DC⃗ The optimization problem is to find the best set of material densities 𝜌⃗∗ that fit the observed data. The TV regularization term is added to the inverse problem to stabilize solutions that can be very sensitive to noise. Using the parameter α, the TV regularization term holds the expectation that the solution is piecewise smooth, or that variations of image composition between neighboring pixels are expected to be small. Balancing the regularization term with the data misfit is important for accurate material reconstructions. Previous work has shown that α can be picked adaptively [10]. Equation (3.2) is an optimization problem over the parameter set of material areal densities 𝜌⃗. In this work, this equation is minimized using the Gauss–Newton method. The Newton method is an iterative optimization method that approximates the minimization function, Equation (3.2) as quadratic at each iteration and minimizes that approximation by varying the parameter vector 𝜌⃗. This repeats until some convergence criteria is met, such as very small changes in 𝜌⃗. Since calculating the second derivatives required for the quadratic approximation can be computationally expensive, the Gauss–Newton method approximates the second derivative as a product of first derivatives of the minimization function. This has been shown previously to be an effective method that balances computational complexity with quick convergence [10]. An example of the input radiographic data and the output material reconstructions given by the algorithm is shown in Figure 5. (a) (b) 6 MVp X-ray Est. Steel (g/cm2 ) 0 80 60 1 40 2 20 0 3 Est. Poly (g/cm2 ) DT neutron 0 10.0 1 7.5 5.0 2 2.5 3 Figure 5. Example input (a) radiography data and (b) algorithm output. 10 0.0 3.1.3 Assumptions and limitations The accuracy of the output of an inverse algorithm depends on knowledge of the system response, including the input source spectrum, transport through the object, and the detector response to a given energy of particle flux. Assumptions: • Radiography data are well registered, such that pixels in one radiograph can be matched to any other radiography data used in the analysis. • A good conjecture of the material composition is available to make a set of materials (or representative materials) to reconstruct. • The response of the system to a given set of materials is thoroughly understood. Spatial resolution, scattering, and more complex interactions in the object or its surroundings are important contributors to the system response. These effects can be modeled or approximated. The present algorithm includes modeling of the spatial resolution and energy response of the detector. Scatter approximations have been implemented previously [10], although flexibility to the different scatter profiles from a variety of objects has not been completed. • The interrogation systems are generally stable. This aligns with the well-understood system responses. An x-ray or neutron generator that varies significantly in spectrum or intensity over time would cause inaccuracies in the material reconstructions. The authors of this report have experience building flexibility to detector instabilities into the algorithm, similar to a flat-field correction with radiography data. This work accounts for instabilities in flux magnitude and endpoint in experimental hyperspectral x-ray radiography data using the observed spectrum. For energy-integrating detectors, an understanding of the energy response of the detectors could be used for similar corrections. Note that others have shown that material reconstructions may not be sensitive to the x-ray endpoint [12]. Limitations: • The number of materials in the material list cannot exceed the number of interrogations, and the interrogations must each carry unique information on the material composition. For example, scanning an object with a 6 and 6.5 MV Bremsstrahlung x-ray source would make material reconstruction difficult because the two spectra are largely overlapping. However, neutron and x-ray radiography carry largely unique material information because the particles interact differently with many materials. • The system models must be able to be computed quickly. The reconstruction is iterative, so full Monte Carlo models cannot be run at each iteration. The aim for reconstructions is in the order of minutes on a typical laptop computer. Reconstruction time scales with the number of pixels, materials, and iterations. Current reconstructions of three 100 × 100 radiographs for three-material reconstructions can currently be completed in about 5 mins for 30–50 iterations. 3.2 MATERIAL RECONSTRUCTION ALGORITHM APPLICATION STUDY Assessing the performance of radiography systems for particular applications is complicated. However, abstracting a desired capability, such as high-resolution or material penetration, and then devising a test to compare different systems by a common metric or process is often possible. In the present study, many desirable capabilities were abstracted and evaluated to determine the extent to which quantitative material 11 reconstructions could aid these tasks. In particular, the application study focused on the following three capabilities: 1. Quantifying the gap size between concentric annuli 2. Quantifying the amount of low materials of similar atomic number (-Z) material shielded by high-Z material 3. Discriminating among materials of similar Z For each capability, a material geometry was chosen for evaluation, radiation-transport simulations performed to produce radiographs, and analysis performed using the material reconstruction algorithm. The present report summarizes this work, but more information can be found in the presentation Material quantification with dual-mode radiography, PNNL-SA-154031. 3.2.1 Gap determination To evaluate quantification of gaps, x-ray and neutron radiographs were calculated via the Monte-Carlo NParticle (MCNP6) radiation transport simulations for an object composed of two annuli. In the absence of a gap, the object consisted of a stainless-steel inner annulus having an inner diameter of 13.4 cm and an outer diameter of 18.4 cm and an outer high-density polyethylene annulus with an inner diameter of 18.4 cm and an outer diameter of 23.4 cm. The presence of a gap was created by decreasing the thickness of each annulus equally but maintaining the overall outer and inner diameter of the object. Transmission images were calculated for 6 and 9 MVp energy Bremsstrahlung x-rays and 14.1 MeV D-T neutrons to obtain three radiographs built from 1 cm2 face F4 tallies arranged on a 100 × 100-pixel cylindrical detector. This number of pixels matched the expected neutron detector pixel pitch and eased the computational burden of investigating the accuracy of material reconstruction. Materials were defined from natural isotopic fractions to reference the latest cross section libraries. Flux convergence of <1% was obtained for 108 histories. Statistical noise was added to the simulated transmission neutron radiographs by assuming either 100 s (high noise) or 1,000 s (low noise) acquisitions of nominally 100 neutron counts per second per pixel in an unattenuated image pixel and sampling a Poisson distribution after scaling the simulated image data. Neutron radiographs were blurred with a 2D Gaussian kernel with a standard deviation of two pixels. The x-ray radiographs had lower noise by a factor of 10 and better spatial resolution by a factor of two. Gaps between 0 and 10 mm were simulated by modifying material thicknesses. Example simulated data and material reconstruction outputs are shown in Figure 6. (a) (b) (c) Figure 6. Example simulated radiograph and reconstructed material areal densities for the two annuli. The (a) radiograph is shown with simulated noise of a 100 s measurement without blurring along with reconstructed (b) steel and (c) high-density polyethylene areal densities. Using the reconstructed material compositions, gaps between the material layers were determined by fitting the reconstructed material outputs to a ray-traced model of an annulus to determine the inner and outer radii of each material. Results from this analysis are summarized in Table 1, with varying levels of image noise and blurring. The following are key conclusions from these results: 12 1. Material reconstructions allow quantitative material measurements, even in the presence of image noise and blurring, provided that a good system response model is provided. 2. Although neutron radiography data have higher noise and greater image blurring, higher-accuracy material reconstructions can be made compared to those with x-rays alone. Table 1. Summary of quantitative gap determination using the reconstruction algorithm output. Case High noise, no blur Low noise, blur Low noise, blur w/o blur model High noise, no neutron Est. gap at 90% confidence (mm) 4.0 4.3 8.1 >10 The material quantification algorithm does not presently appear to support a more accurate determination of edges than the direct analysis of the images. However, on an encouraging note, the determination of material areal densities appears to be minimally affected by the presence or absence of gaps. In the present work, accuracy in determining the gap size was limited by the coarse pixel pitch (corresponding to neutrons) implemented in the material quantification algorithm. Future work will address implementing the algorithm on the high-resolution pixel grid of the x-ray detector. Implementation at high resolution may require variable grid spacing, with increased sampling near rapidly varying regions of the image close to object edges, and reduced sampling in slowly varying areas of the image to reduce computational burden. Quantification of gaps in measured data can also be affected by small-angle scattering that, when scattered down open paths, is preferentially transmitted to the detector. 3.2.2 Quantifying low-Z material behind high-Z material The second desirable capability that has been identified is the ability to quantify low-Z, low-density material shielded by high-Z, high-density material. This quantification is challenging when the high-Z material attenuates most of the x-ray flux. Using MCNP, radiographs of overlapping lead and polyethylene slabs were simulated with a range of material thicknesses that ranged from 0.7 to 7 mean free paths of Pb and 0.1 to 1 mean free paths of polyethylene (to 6 MeV x-rays). Again, 6 and 9 MVp xray radiographs were simulated along with 14.1 MeV D-T neutron radiographs. For this study, the key benefit of neutron radiography over x-ray radiography for accurate reconstruction of the polyethylene thickness is the greater penetration through the thick Pb layer. This benefit is shown by lower errors in the polyethylene thickness estimated by the reconstruction algorithm when including the neutron information. For example, the estimated thickness extracted from a number of radiographs of 1.7 cm polyethylene thickness behind 7 cm of lead had an error of 22 ± 51% when using x-rays alone, but 12 ± 7% when using x-rays and neutrons. 3.2.3 Discrimination among materials of similar Z The third capability that has been identified is the ability to discriminate among materials of similar atomic number (Z) and density. Materials of similar Z exhibit similar x-ray attenuation, making it difficult to tell them apart even with a dual- or multi-energy x-ray interrogation. Neutron attenuation coefficients are distinct from x-ray and offer a unique material signature that, when combined with x-rays, can be used for accurate material reconstructions. 13 To investigate this capability, interrogation of high-Z spheres in an AT400R nuclear materials storage container was simulated using MCNP. The container has layers of polyethylene and steel and contained either two Pu spheres, or one Pu sphere and one W sphere. Three x-ray radiographs at 5, 6, and 9 MVp and a 14.1 MeV D-T neutron radiograph were simulated, and the resulting radiographs are shown in Figure 7. This object is highly attenuating, especially for x-rays, which exhibit a thickness of about 7 mean free paths through the thickest portion of the object. The simulated image data provided to the reconstruction algorithm here assumed a 1,000 s neutron acquisition and no image blurring. Again, the xray radiographs are assumed to be ten times lower noise. (a) (b) Figure 7. Simulated (a) x-ray and (b) neutron radiographs of the AT400R nuclear storage container in which the top sphere is Pu and the bottom sphere is W. Outputs from the reconstruction algorithm are shown in Figure 8, using a three-material list of Pu, steel, and HDPE for the reconstruction. Given the high x-ray attenuation of the object, the material reconstructions using x-rays alone are poor and have difficulty quantifying the material in the thickest portion of the object. Adding the neutron data with its greater penetration and unique material signature allows material reconstruction behind those thickest portions. Notably, the estimations with neutron radiography data also show a variation in material estimated in the region of the bottom sphere as an underestimate of Pu and steel along with an overestimate of poly. This is an anomaly that indicates that the top sphere is not Pu. These results provide an example of how neutron radiography in addition to xray radiography can be useful for both increased penetration and quantitative analysis. (a) (b) Figure 8. Reconstruction outputs using radiography data from (a) x-rays only and (b) both x-rays and neutrons. 3.3 ADDITIONAL API PHYSICS MODELING Beyond the analysis of transmission radiography data, PNNL has begun work exploring how to use additional information that is available from the API method. In particular, this effort has concentrated on the following two API observables: 1. Secondary gamma emissions from API neutron interrogation. For material signatures, the outgoing gamma rays have an energy signature, such as the unique energy dependence of inelastic scattering 14 emissions, and an intensity signature, specifically with interactions such as fission [13] that result in considerable gamma generation. 2. The observed time spectrum from the API neutron interrogation. Because a multitude of interactions occur with a beam of neutrons passing through an object, the time spectrum can be seen as an aggregate measure that carries information about these interactions, including gamma generation and neutron scattering. This section introduces how these observables are being modeled and used, though this work is still continuing. 3.3.1 Modeling gamma emission Work on modeling the gamma-ray emissions from an interrogated object are summarized here. A more complete description can be found in the December 2019 report to NA-22, Prompt gamma emission for material identification, PNNL-29512. To quantify the material composition of an inspected object using its prompt gamma-ray emissions, a simple and accurate physics model is required that describes the generation of gamma rays from the interacting fast-neutron beam and the attenuation of these gamma rays as they exit the object. An initial model presently under development will use a pencil-beam attenuation model to study fast-neutron reaction rates for interactions that create secondary gammas in simple object models, as seen in Figure 9. (a) (b) MCNP sim analytical model 104 γ emissions counts/keV 103 102 101 100 Fast neutron interrogation 10 1 0 1000 2000 3000 4000 5000 energy (keV) Figure 9. (a) A schematic diagram of neutron interrogation producing gamma radiation, and (b) a comparison of the output of the analytical model given in Equation (3.3) and MCNP simulation for a small cube of Fe. An analytical model can be posed to determine fast neutron reaction rates 𝐹 at a location 𝑥 in an object that would create gammas at a variety of energies: 𝑑, = ∫< #$% (𝑟, 𝐹(𝐸F )𝐷(𝐸6 ) 𝑒 2G(<& )C + 𝐵𝐺(𝐸6 ))𝑑𝐸6 , (3.3) where 𝐹(𝐸F , ∆𝑧) = ∫∆I 𝜙(𝐸F , 𝑧)Σ(𝐸F ) 𝑑𝑧. Here, an observable 𝑑 of index 𝑘 is a function of the gamma production rate 𝑟 and the neutron reaction rate 𝐹 from which secondary gamma rays are generated from the interrogating neutron flux 𝜙 of energy 15 𝐸F and cross section Σ. The detector response function 𝐷S𝐸J W captures the gamma detector efficiency, finite energy resolution, and partial energy deposition. The exiting gammas are exponentially attenuated according to the object density 𝜌 and attenuation coefficient µ depending on the energy 𝐸J of the secondary gamma ray. The integration over z to get F may be over the whole object, or a small δ𝑧, depending on whether the interrogating neutrons are continuous or tagged/pulsed. The background detection rate from, for example, scattered gamma and neutron radiation, is given by BG. The OpenMC software [14] was used to retrieve evaluated nuclear data file (ENDF) data. This software puts the data into HDF5 files, which eases accessibility to the wide range of data available in ENDF. Although OpenMC also includes Fortran code to execute Monte Carlo criticality calculations, the present work uses the Python API to provide interpretable nuclear data. Notably, the OpenMC code was developed for neutron transport, though the most recent versions also include photon interaction data. With this tool, nuclear data can be quickly queried to find incident neutron data of interest, including secondary gammas from inelastic and capture interactions. Using the reaction type (MT) numbers as defined in the ENDF Manual, neutron interaction cross sections and gamma ray emission cross sections can be determined for a variety of interaction types. A simple MCNP simulation was made to examine the photon production physics and to determine how well Equation (3.3) can be used to anticipate the gamma emissions that could be observed with 14.1 MeV neutrons. In this simulation, a pencil beam of 14.1 MeV neutrons impinge on a small 5 × 5 × 5 mm3 material volume. A 4π-flux tally and 4 in. right-circular cylinder high purity germanium detector tally are placed 30 cm away to capture the spectrum of gammas that are emitted from the sample. Whereas a real object interrogated by a real D-T neutron source would require a more complex model, the objective was to use a “minimal” model to extract the physics processes generating photons, such as inelastic scattering, fission, (n, 2n), and so on. From the model, it was determined that most of the gammas produced on 208Pb are from (n, 2n) reactions generating the continuum in the spectrum in Figure 9, although how they could be used for material identification is unclear. As in previous work [11], this study focuses on the discrete emission lines of the inelastic (n,ng) excitations, although any gamma-rayemitting interaction could potentially be accommodated for in algorithm development. Figure 9 (b) shows the comparison of the output from the small-cube MCNP model and the analytical model, Equation (3.3). Agreement on the intensity of the discrete gamma lines is good, and the gamma rays emitted from reactions that result in a gamma-ray continuum are shifted to lower energies. In threebody interactions like (n, 2n), the energy of the outgoing gamma will be affected by the energy of the outgoing neutrons, which is not yet captured in the analytical model. Nonetheless, this result shows that an analytical model can capture many of the physical attributes of gamma-producing neutron interactions in a simple geometry. Work is continuing on developing realistic physics models and implementing the data into the material reconstruction algorithm. 3.3.2 Time spectrum as a material signature PNNL is currently evaluating the possibility of performing material quantification from pixel-based measurements of time spectra. Spatially dependent time spectra can be a signature of both the material and geometry of the object and could supplement the neutron/gamma-based image data beyond simply the reduction of scatter. Figure 10(a) shows the geometry of the simple proof-of-concept case of a plane detector with 1 cm resolution at 85 cm from a D-T neutron source. The object geometry includes four 3 cm-long cylinders 16 with radii of 4 cm. The objects are defined with a small, 1 cm, gap between the pairs of cylinders in the beam direction. Further, each pair is offset by 8 cm from the center line of the beam line. For this proofof-concept demonstration, we suppose a nonphysical interaction of the neutron that can only result in capture and emission of a prompt gamma ray, and that the emission angle of the gamma rays is uniform within 12° of the original neutron’s direction. The material specific (n,𝛾) probability for the cylinders on the right is twice the probability of the cylinders on the left. Although the geometry and physics are nonphysical at this point, this model is presented to begin to generate efficient physics models that capture the timing of the various neutron interactions in an object. (a) (b) Figure 10. (a) Geometry of detector, objects, and source for proof-of-concept model, and (b) time spectra for (black) full detector average, (blue) an 8 × 8 region of interest centered on the object on the left, and (dashed) an 8 × 8 ROI centered on the object on the right. For this model, a Monte Carlo framework for particle tracing with (n,𝛾) kinematics has been set up to evaluate the expected magnitude of differences in the time spectrum that would be anticipated for a given geometry and object. This framework samples both interaction probabilities in the cylinders and the gamma emission angle to calculate the time spectra that would be observed on the image plane for this given object. Time bins are 0.5 ns wide, and no further timing resolution (blurring) was applied. Figure 10(b) shows the time spectra for two 8 × 8 pixel patches centered in the object regions (blue) along with the full detector average (black). The left side (solid blue) shows more events at shorter times, from the faster moving gamma rays, and fewer slower neutrons at the longer times than the does the right side (dashed blue). The shorter times (~10 ns) indicate captures in the objects nearer to the source where the distance traveled by the neutrons is the shortest. Also, the central dip at about 11 ns is due to the 1 cm gap between each pair of cylinders. Isotropic photon emission processes can make reconstruction algorithms more complicated. Work is in progress to add realism to this model to determine the capability of using time spectra in reconstructing the material composition of an object. 3.4 OUTLOOK FOR MATERIAL IDENTIFICATION In the upcoming year, PNNL expects to split the effort between two thrusts that are both intended to increase the accuracy of material reconstruction. The first thrust will be to revise the forward model for transmission to account for additional physical effects such as x-ray and neutron scattering. The second thrust will be to include additional API observables beyond transmission. In particular, secondary gamma emission and the time spectrum on the API detector will both be explored in the context of high-accuracy material reconstruction. 17 4. 4.1 PROGRESS ON MATERIAL IDENTIFICATION USING FULL SPECTRAL ANALYSIS INTRODUCTION Most approaches to analyzing API neutron data isolate single observables such as neutron transmission, elastic scattering, inelastic collisions and emitted neutrons or gamma rays, (n,xn) reaction neutrons, or induced-fission neutrons or gamma rays. When attempting to identify materials and densities in an inspected object, simultaneously fitting all possible signatures and evaluating all combinations of materials, densities, and geometries of the inspected object is desirable. The present effort assumes a known geometry but unknown materials and densities. Evaluation of many thousands of combinations of materials will require approximate forward transport that is fast enough to perform the ensemble of calculations in minutes. The initial thrust of this effort is to identify means of implementing this fast forward model. During FY 2020, North Carolina State University began developing a method for accelerating transport calculations using precomputed voxel “transfer functions.” Figure 11 illustrates this approach; the voxel, represented in the figure as a one-dimensional (1D) slab of material, is modeled using a transfer function (also known as an impulse response or Green’s function) that predicts the particle current emergent from the voxel given the incident particle current: 𝐽AK) (𝑦) = y 𝐻(𝑦; 𝑥)𝐽'F (𝑥) 𝑑𝑥, (4.1) $ where the incident current 𝐽'F is a function of energy 𝐸'F , and direction 𝜇'F , the cosine of the particle’s direction with respect to the surface normal of the 1D voxel, 𝑥 = {𝐸'F , 𝜇'F } ∈ 𝑋 = (0, ∞) × [−1,1], (4.2) and the emergent current 𝐽AK) is a function of energy 𝐸AK) , direction 𝜇AK) , and time delay 𝜏 = 𝑡AK) − 𝑡'F , 𝑦 = {𝐸AK) , 𝜇AK) , 𝜏} ∈ 𝑌 = (0, ∞) × [−1,1] × (0, ∞). (4.3) In operator notation, the distribution of emergent particle energy, direction, and time-delay is related to the distribution of incident energy and direction by 𝐽AK) = 𝐇𝐽'F . (4.4) Figure 11. Illustration of a 1D voxel transfer function H relating emergent current 𝑱𝒐𝒖𝒕 to incident current 𝑱𝒊𝒏 . 18 If energy, direction, and time are histogrammed into discrete “bins,” then the operation is equivalent to matrix multiplication. For a series of coupled voxels, the current emerging from the last voxel can be predicted given the current incident on the first voxel, as Figure 12 illustrates. Note that coupling multiple voxels will require accounting for the current emergent from each face of the 1D voxel, because some particles will emerge from the same face through which the source particle entered. Additionally, this effect will complicate the calculation of current emergent from coupled voxels; the simple multiplication illustrated in Figure 12 will have to be modified to iteratively sweep forward and backward to account for particles that were “reflected” by the voxel. Figure 12. Coupling of multiple voxel transfer functions. 4.2 FY 2020 PROGRESS During FY 2020 Q3 and Q4, MCNP6.2 was used to model 1D, 1 cm thick slabs of polyethylene, Fe, and highly enriched U (HEU) metal and accumulate the corresponding voxel transfer function 𝐻 using the ptrac output-recording incident and emergent particle energy, direction, and time. First, the voxel transfer function was accumulated for a monodirectional beam of 14 MeV incident neutrons to visualize the relationship between emergent and incident particles. Figure 13 illustrates the particle trajectories emergent from each 1D slab. Neutron trajectories are shown in red, and photon trajectories are shown in blue; the length of each arrow is proportional to the particle’s energy. The authors also accumulated the histogram of emergent neutron position, energy, direction, and time-delay for each material; these histograms, shown in Figure 14–Figure 16, are 2D projections of the transfer function. Figure 17 shows the transfer functions’ dependence on energy and direction versus time-delay. 19 (a) (b) (c) Figure 13. Neutron and photon trajectories emerging from a 1D, 1 cm thick slab of polyethylene (left), Fe (center), and HEU metal (right) for a monodirectional source of 14 MeV neutrons incident on the left face of each slab. Neutron trajectories are shown in red, and photon trajectories are shown in blue. The length of each arrow is proportional to the emergent particle energy. 20 (b) (b) (c) (c) (c) (c) (a) Figure 14. Histograms of position, direction, energy, and time-delay for neutrons emergent from a 1D slab of polyethylene in response to a monodirectional source of 14 MeV incident neutrons. “Transmitted” neutrons emerge from the side of the voxel opposite the source, and “reflected” neutrons emerge from the same side as the source; uncollided neutrons (those that passed through the slab without interaction) are omitted. The “ridges” labeled in the histograms of emergent energy versus direction result from (a) elastic neutron scatter by 1H nuclei; (b) elastic scatter by 12C; (c) inelastic scatter by excited states of 12C. 21 (a) (a) (b) (b) Figure 15. Histograms of position, direction, energy, and time-delay for neutrons emergent from a 1D slab of Fe in response to a monodirectional source of 14 MeV incident neutrons. “Transmitted” neutrons emerge from the side of the voxel opposite the source, and “reflected” neutrons emerge from the same side as the source; uncollided neutrons (those that passed through the slab without interaction) are omitted. The “ridges” labeled in the histograms of emergent energy versus direction result from (a) elastic neutron scatter by Fe nuclei; (b) inelastic scatter by excited states of Fe. 22 (a) (a) Figure 16. Histograms of position, direction, energy, and time-delay for neutrons emergent from a 1D slab of HEU metal in response to a monodirectional source of 14 MeV incident neutrons. “Transmitted” neutrons emerge from the side of the voxel opposite the source, and “reflected” neutrons emerge from the same side as the source; uncollided neutrons (those that passed through the slab without interaction) are omitted. The “ridges” labeled in the histograms of emergent energy versus direction result from (a) elastic scatter by uranium. Figure 14 through Figure 17 omit uncollided neutrons. Additionally, Figure 14–Figure 16 partition the emergent neutron current into collided neutrons that emerged from the side of the slab opposite the source (“transmitted” neutrons) and those that emerged from the same side of the slab as the source (“reflected” neutrons). This terminology has been adopted throughout the remainder of this report. Notably, for the 14 MeV neutron source, the number of transmitted neutrons is greater than the number of reflected neutrons by approximately 10:1. 23 (a) (b) (c) Figure 17. Histograms of energy, direction, and time-delay for neutrons emergent from a 1D slab of (a) polyethylene, (b) Fe, and (c) HEU metal in response to a monodirectional source of 14 MeV incident neutrons. The discrete “ridges” that appear in the histograms of neutron energy versus direction in Figure 14–Figure 17 are caused by elastic neutron scattering and inelastic scatter by excited states, as labeled in the figures. 24 In Figure 17, these features show that neutrons scattering at a greater angle (and therefore traveling a greater distance) emerge later. As Figure 14–Figure 17 illustrate, accumulating the voxel transfer function for a single-incident energy and direction is relatively straightforward. In this case, 𝐽'F (𝑥) = 𝑄𝛿(𝐸 − 𝐸" )𝛿(𝜇 − 𝜇" ), (4.5) where 𝑄 denotes the total number of source particles, so the emergent current 𝐽AK) (𝑦) = 𝑄 ∙ 𝐻(𝑦; 𝐸 = 𝐸" , 𝜇 = 𝜇" ), (4.6) such that the transfer function for discrete incident energy 𝐸" and direction 𝜇" is simply1 𝐻(𝑦; 𝐸 = 𝐸" , 𝜇 = 𝜇" ) = 𝐽AK) (𝑦)/𝑄. (4.7) Consequently, the transfer function for a given voxel composition and geometry can be accumulated from a Monte Carlo simulation of the transport of source neutrons with randomly selected energy and direction. Each individual source neutron is started with an initial energy and direction randomly chosen from a distribution that spans the range of neutron energies and directions relevant to the problem, with probability selected to adequately sample all relevant neutron energies and directions. On a neutron-byneutron basis, 𝑄 from Equation (4.7) is equal to one, and each neutron history predicts one realization of the transfer function. For the transfer functions shown later in this report, the authors used a neutron source with directions isotropically distributed over 0° to 45° with respect to the 1D voxel surface normal and energies distributed uniformly on a logarithmic scale over 10-11 to 20 MeV (the minimum and maximum neutron energies tabulated for most cross sections). Figure 18 shows the histogrammed source neutron energy and direction distributions of the MCNP source term. Figure 18. Distribution of source neutron energy and direction used to generate 1D voxel transfer functions. The expected value of the transfer function can be estimated from a large number of source neutron histories by histogramming occurrences of neutron histories with given emergent energy, direction, and time delay for given incident energy and direction: 1 Note, the histograms shown in Figure 14–Figure 16 were not normalized by the total number of 14 MeV source neutrons; approximately 500,000 source neutrons were simulated to accumulate these histograms. 25 𝐻(𝑦 = {𝐸AK) , 𝜇AK) , 𝜏}; 𝑥 = {𝐸'F , 𝜇'F }) = 𝐽AK) (𝐸AK) , 𝜇AK) , 𝜏; 𝐸'F , 𝜇'F ) , 𝐽'F (𝐸'F , 𝜇'F ) (4.8) where the division of 𝐽AK) by 𝐽'F is performed elementwise on the rows of 𝐽AK) . Figure 19 shows the energy-to-energy transfer function for neutrons transmitted through and reflected from 1D voxels of polyethylene, Fe, and HEU metal. Observe that upscattering is evident at thermal energies for neutrons reflected from the polyethylene and Fe voxel, and fission spectrum neutrons are emergent from the HEU voxel for all incident neutron energies. Again, these transfer functions do not include uncollided neutrons; they were filtered from the Monte Carlo histories by removing those histories where the incident energy and direction were unchanged as the neutron passed through the voxel. The authors think the diagonal ridges evident in each transfer function are caused by low-angle scattering and plan to review the logic for excluding uncollided particles. 26 (a) (b) (c) Figure 19. Transfer functions relating emergent and incident energy for 1D slabs of (a) polyethylene, (b) Fe, and (c) HEU metal. The accuracy of the voxel transfer functions was tested by simulating an isotropic, fission spectrum neutron source incident on each voxel. Figure 20 shows the energy and direction distribution of the MCNP source term. 27 Figure 20. Distribution of energy and direction for the test case using an isotropic fission spectrum neutron source. Figure 21 compares the energy-dependent current of transmitted and reflected neutrons predicted by direct MCNP simulation (black) to the current predicted using the voxel transfer functions (red) shown in Figure 20. The error in the transfer function prediction is shown by the residuals (blue). In each case, the residuals are a small fraction of the current, but for reflected particles, the errors are higher because of the small sampling of reflected particles when generating the transfer functions. The figures also include uncertainties estimated from the count of neutrons contributing to each element of the transfer function; these are shown as high/low shaded regions, but they are too small to see clearly without magnifying the plot. 28 (a) (b) (c) Figure 21. Comparison of transmitted and reflected neutron current calculated directly using (black) MCNP to (red) the current predicted using transfer functions for 1D voxels of (a) polyethylene, (b) iron, and (c) HEU metal. 29 Each transfer function shown in Figure 19 was accumulated from approximately 500,000 neutron histories. The MCNP computational time ranged from a few seconds for polyethylene to roughly 3 h for HEU. All the MCNP calculations used analog transport, which is likely to remain necessary; ultimately, the simulations must include correlated neutrons and gammas from inelastic scatter and induced fission. Parsing the MCNP ptrac records takes approximately 1 min to extract the incident and emergent neutron and photon state variables (position, energy, direction, and time) and collate them into records relating incident and emergent particles. Loading and histogramming these records to accumulate the transfer function takes a few seconds (the majority of this time is spent loading the records; accumulating the transfer function takes about 20 ms). All of these operations can be performed in advance. The actual matrix multiplication to calculate the emergent current from the incident current takes less than a microsecond. This initial demonstration indicates that the voxel transfer function approach should enable sufficiently accelerated transport calculations necessary to analyze neutron-coded aperture images to identify materials present in the target. 4.3 FY 2021 PLANS In FY 2021, NCSU will extend the demonstration presented in this report to include time delay and direction in voxel transfer functions. This should not present a significant technical challenge, aside from requiring longer Monte Carlo simulations to adequately populate the elements of the transfer function with a sufficient number of particle histories. More importantly, NCSU will begin coupling voxels composed of different materials to (a) develop the procedure to iteratively sweep forward and backward to account for particles that were transmitted and reflected by the coupled voxels and (b) evaluate the accuracy of the transfer function for multiple coupled voxels. Based on the observations illustrated in Figure 14 through Figure 16 and Figure 19, the number of transmitted neutrons is greater than the number of reflected neutrons by approximately 10:1; consequently, the iterative backward sweep to account for reflected neutrons is likely to be a correction on the order of 10%. Implementing reduced-order modeling of voxel transfer functions is a stretch goal for FY 2021. In FY 2020 Q3, the potential to implement reduced-order modeling of the transfer function using truncated Karhunen–Loève expansions was reported. In FY 2021, implementing truncated Karhunen–Loève expansions and other reduced-order modeling techniques and their influence on the accuracy of the voxel transfer functions will be investigated. 5. ACKNOWLEDGMENTS This work is supported by the US Department of Energy, Office of Defense Nuclear Nonproliferation Research and Development in the National Nuclear Security Administration (NA-22). 6. REFERENCES [1] B. Perot et al., “Development of the EURITRACK Tagged Neutron Inspection System,” Nuclear Instruments and Methods in Physics Research Section B, vol. 261, pp. 295–298, 2007. [2] L. Nakae et al., Fast Time-of-Flight Correlation System—End-of-Year Report FY 2018, LLNLAR-763682, Lawrence Livermore National Laboratory, November 2018. 30 [3] P. Hausladen, M. Blackston, and J. Gregor, Progress Update on Iterative Reconstruction of Neutron Tomographic Images, ORNL/SPR-2019/1353, Oak Ridge National Laboratory, 2019. [4] S. Rose, M. Andersen, E. Sidky, and X. Pan, “Noise Properties of CT Images Reconstructed by Use of Constrained Total-Variation Data-Discrepancy Minimization,” Medical Physics, vol. 42, pp. 2690–2698, 2015. [5] A. Chambolle and T. Pock, “A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging,” Journal of Mathematical Imaging and Vision, vol. 41, no. 1, pp. 120– 145, 2011. [6] J. Gregor, P. Bingham, and L. F. Arrowood, “Total Variation Constrained Weighted Least Squares Using SIRT and Proximal Mappings,” 4th International Conference on Image Formation in X-ray Computed Tomography, Bamberg, Germany, July 2016. [7] J. Gregor and J. A. Fessler, “Comparison of SIRT and SQS for Regularized Weighted Least Squares Image Reconstruction,” IEEE Transactions on Computational Imaging, vol. 1, no. 1, pp. 44–55, 2015. [8] M. A. Blackston and P. A. Hausladen. “Fast-Neutron Elastic-Scatter Imaging for Material Characterization,” 2015 IEEE Nuclear Science Symposium and Medical Imaging Conference (NSS/MIC), San Diego, CA, 2015. [9] A. J. Gilbert, B. S. McDonald, S. M. Robinson, K. D. Jarman, T. A. White, and M. R. Deinert, “Non-Invasive Material Discrimination Using Spectral X-Ray Radiography,” Journal of Applied Physics, vol. 115, no. 15, 2014. [10] A. J. Gilbert, B. S. McDonald, and M. R. Deinert. “Advanced Algorithms for Radiographic Material Discrimination and Inspection System Design,” Nuclear Instruments and Methods in Physics Research Section B, vol. 385, pp. 51–58, 2016. [11] B. Canion, S. McConchie, and S. Landsberger, “Source Correlated Prompt Neutron Activation Analysis for Material Identification and Localization,” IEEE Transactions on Nuclear Science, vol. 64, no. 7, pp. 1725–1732, 2017. [12] M. Klasky, “Advanced Radiography Software: Algorithm for Radiographic Reconstructions,” NA84 Mid-Year Review: LANL, 2020. [13] P. A. Hausladen et al., “Induced-Fission Imaging of Nuclear Material,” INMM 51st Annual Meeting, pp. 1–10, 2010. [14] P. K. Romano, N. E. Horelik, B. R. Herman, A. G. Nelson, B. Forget, and K. Smith, “OpenMC: A State-of-the-Art Monte Carlo Code for Research and Development,” Annals of Nuclear Energy, vol. 82, pp. 90–97, 2015. 31