In 2011, the Department of Homeland Security (DHS) National Technical Nuclear Forensic Center (NTNFC) hosted a panel of Plutonium (Pu) experts to develop a plan for advancing forensic analysis of Pu materials. During this conversation, the experts concluded that different production processes produce final products with different characteristics. They hypothesized that these different characteristics, or “signatures”, on the final product could potentially allow a forensic analyst to determine which processes were used to produce the Special Nuclear Material (SNM), and, in turn, make inferences on where the material originated.

As a result of these discussions, scientists at Pacific Northwest National Laboratory (PNNL) conducted an experiment designed by statisticians at Los Alamos National Laboratory (LANL) and Sandia National Laboratories (SNL) to replicate historical and modern Pu processing methodologies and conditions. This experiment consisted of 76 runs, where each run considered the same set of nine processing parameters, whose values intentionally varied from run to run. For each run, the resulting SNM was imaged with a scanning electron microscope (SEM) to generate images of the various particles of Pu. These images were then segmented using LANL’s Morphological Analysis for Material Attribution (MAMA) software (Porter et al., 2016). This post-processing segmentation extracts the different particles that are present in each SEM image, and generates measurements based on the physical and morphological characteristics of the particles. Such measurements include particle areas, aspect ratios, convexities, circularities, gradients, and shadings (Zhang et al., 2021).

This paper serves to statistically relate the processing conditions of the different runs to the morphological and physical characteristics of the resulting Pu particles. In particular, we consider Pu particles that were processed using a solid oxalate feed, and particles that were processed in a 0.9 M oxalate solution. We relate the observed MAMA characteristics (which include, but are not limited to, those listed above) to the processing conditions used to produce the samples of Pu (which include, but are not limited to, temperature, nitric acid concentration, and Pu concentration) using inverse prediction methods. We consider Bayesian methodologies to quantify our predictions on the processing conditions, and to naturally incorporate uncertainty into these predictions.

We discuss these statistical methodologies in Section 2; in Section 3, we outline the specifics of the problem from a statistical perspective; in Section 4, we apply the methodologies to a dataset obtained from studying actual Pu particles generated under intentionally varied processing conditions within a design of experiments framework (rather than a simulated dataset); and in Section 5, we discuss the implications of our results.

To study the relationship between the processing parameters and the particle characteristics, we consider an inverse prediction framework. The traditional regression problem involves making predictions about responses, given a set of explanatory or input predictors. As a simple example, we may be interested in determining the weight of an individual, given their height, sex, and nationality. Conversely, the inverse prediction framework reverses the problem and so considers making predictions about the explanatory predictors, given a set of responses. Following the simple example above, we would now be interested in determining the height, sex and nationality of an individual, based on their weight. Inverse prediction is used across a variety of disciplines, and, in particular, has been used in forensic science, and nuclear forensics [see, for example, Lewis et al. (2018), Ries et al. (2018), and Ries et al. (2022)].

Mathematically, consider a q-dimensional response vector, y = y 1 , y 2 , … , y q ′ , and a p-dimensional input vector x = x 1 , x 2 , … , x p ′ . We express the relationship between x and y as y = g x | θ + ε , (1) where g ⋅ is the true underlying physical phenomenon that is responsible for producing the results, which are typically represented by an emulator or surrogate that mimics the outputs in terms of the inputs, θ is a vector of model parameters, and ɛ is a random vector that captures noise present in the observed data.

As mentioned above, the goal of the inverse prediction framework is to predict the value of the input variables x ′ that produced a new observation y ′. There are two approaches by which we can learn about y ′: We can either 1) construct a model that directly predicts x ′ as a function of y ′, or we can 2) invert a traditional forward model that predicts y ′ as a function of x ′. We will refer to the first approach as the “direct” model, and the second approach as the “inverse” model. It should be noted that, while it is more convenient to construct the “direct” model, this approach may violate regression assumptions. For example, standard linear regression models assume that the input variables are measured with negligible error. By constructing a model that directly predicts x ′ as a function of y ′, we are treating x as the output variable and y as the input variable and so this assumption no longer holds. Additionally, the “direct” model may not be well-suited for optimal design of experiments, since the vast majority of literature in design of experiments does not consider the inverse problem in its formulation (Anderson-Cook et al., 2015; Anderson-Cook et al., 2016).

In this paper, we use Bayesian Model Calibration (Kennedy and O’Hagan, 2001; Higdon et al., 2004, 2008; Walters et al., 2018; Lee et al., 2019; Nguyen et al., 2021) to approach the inverse prediction framework. Model calibration is a process to estimate, or calibrate, model parameters in the context of an input to output relationship, and falls under the second method (the “inverse” model) discussed above. More specifically, we are interested in using a calibration approach to make predictions about the processing conditions that were used to produce samples of Pu. Each of these samples consists of several particles whose morphological characteristics contain information pertaining to the conditions under which they were produced. By studying the information contained in these particles, we can infer the associated processing conditions.

In this paper, we use a fully Bayesian adaptive spline surfaces (BASS) framework to model or build an emulator to approximate the underlying relationship between the inputs and outputs (Francom et al., 2018; Francom et al., 2019; Francom and Sansó, 2020) ¹ . Suppose y is a vector of morphological characteristics associated with a sample of Pu, and x is a vector of processing conditions used to produce a sample of Pu. Without loss of generality, suppose that each x ₁, x ₂, …, x _p ∈ [0, 1]. The input to output relationship is modeled using. y = a 0 + ∑ m = 1 M a m B m x + ε , (2) where a ₀ is the intercept, B _m ( x ) is a basis function on the input variables x , a _m is the coefficient for the mth basis function, and ɛ ∼ MVN(0, Σ) (Francom and Sansó, 2020). We consider priors for each of the parameters, a , σ ², and M. Note that, because we define a prior over M, the number of basis functions is not fixed, and varies throughout the sampling process. These parameters are sampled using a Markov Chain Monte Carlo (MCMC) process. The samples obtained from this process constitute the posterior distribution of the input parameters, given the observed output parameters.

Upon defining the appropriate BASS model, we perform Bayesian Model Calibration, which allows us to estimate the input parameters that make the model best match the data provided by the output parameters. In this instance, we do not include a discrepancy term. Ideally, a discrepancy term would capture model bias error (i.e., how well a physical model approximates reality), which, in our case, is challenging to handle. Future work will include incorporating this discrepancy term, as well as exploring alternative emulators that are able to better capture the relationship between the input and the output parameters. For more details about Bayesian Model Calibration, see Kennedy and O’Hagan. (2001), Higdon et al. (2004), Higdon et al. (2008), Walters et al. (2018), and Nguyen et al. (2021).

Suppose that run r ∈ 1 , … , R results in n _r images and L _r particles. The number of images varies between particles, and the number of particles varies between images, so that the number of particles per run is not consistent. For our data, R = 76. Each of the runs in our experiment yielded between 5 and 176 images per run, and between 81 and 4,643 particles per run. By using information contained in these runs, we can make predictions about the processing conditions of a new run. In this scenario, our input parameters include processing conditions such as temperature, nitric acid, and Pu concentrations, and responses include morphological characteristics such as particle shape and size.

Figure 1 demonstrates two different particles from two different runs, each with different processing conditions. Note that, while size and shape are useful for determining processing conditions of a sample, far more information is required to accurately predict the processing conditions of these two particles. For example, we can see that, while these particles are produced under entirely different processing conditions, they do maintain physical similarities, alluding to the difficulty of this problem. While both particles are made up of flat sheets, we can see that the arrangement of and the angles between these sheets do differ from one another. This information is captured by segmenting the particles using LANL’s MAMA software, which allows for capturing information about the morphological characteristics, such as area, diameter, and aspect ratio. By including the information from the MAMA software in our model, we can better distinguish between particles from different runs.

From Section 2, we have that y is a p-vector of morphological characteristics associated with a sample of Pu, and x is the q-vector of processing conditions used to produce that sample of Pu. Given our sample of L _r particles for each of our 76 runs, Y _r is the L _r × p matrix of morphological characteristics, where each row of Y _r corresponds to the morphological characteristics of the L _r particles associated with run r ∈ {1, …, 76}. Our full matrix of morphological characteristics is then given by the ∑ r = 1 76 L r × p matrix Y , such that Y = Y 1 Y 2 ⋮ Y R where each  Y r = y r 1 y r 2 ⋮ y r L r = y r 11 y r 12 ⋯ y r 1 p y r 21 y r 22 ⋯ y r 2 p ⋮ ⋮ ⋱ ⋮ y r L r 1 y r L r 2 ⋯ y r L r p . As an example, consider row y _r2 in the matrix Y _r . This vector corresponds to the p observed morphological characteristics associated with the second particle in run r, so that y _r21 corresponds to the observed value for the first morphological characteristics for particle 2 in run r, y _r22 corresponds to the observed value for the second morphological characteristic for particle 2 in run r, and so on.

The matrix X is analogously defined for the processing conditions, where X _r is the L _r × q matrix of processing conditions, where each row of X _r is the q-vector of processing conditions associated with the corresponding run r. For dimensional consistency, we consider X _r as the L _r × q matrix, where x _r is merely repeated for each of the L _r rows in X _r , since all particles in run r are produced under the same set of processing conditions.

Before proceeding with calibration, we must first define our emulator. We choose to fit a BASS model using the associated BASS package in R (Francom and Sansó, 2020). We then perform model calibration to obtain the best set of processing conditions that are associated with a new set of observed morphological characteristics, captured by the matrix Y _R+1, that results from observing a set of particles from a new run.

Since run R+1 results in L _R+1 particles, we calibrate on each of the l _R+1 vectors of morphological characteristics. Each of these calibrations results in N _MCMC samples for each of the p processing conditions. That is, for a set of L _R+1 particles, we obtain L _R+1 N _MCMC × p matrices of predicted processing conditions, X R + 11 , … , X R + 1 L R + 1 . We can use the results from this set of L _R+1 matrices to determine a point estimate, and highest posterior density intervals for the p processing conditions. This process is outlined in Algorithm 1.

Algorithm 1.

Methodology for predicting processing conditions of newly observed particles.

Data: The matrix X of processing conditions for the set of observed particles; The matrix Y of morphological characteristics associated with each particle in the set of observed particles; A set of L _R+1 newly observed particles; The matrix Y _R+1 of morphological characteristics associated with each of the newly observed particles.

In this section, we apply the above methodology to our dataset. The data considered in this experiment consists of runs that were processed using two types of oxalate feed. Out of the 76 total runs, 24 of these runs were processed using a solid oxalate feed, and 52 were processed in 0.9 M oxalate solution. We analyze these two sets separately (Sections 4.1, 4.2), and jointly (Section 4.3). Considering the different types of runs separately and jointly allows us to study the effects of training models on data that are either a) more representative of an interdicted sample (i.e., when we consider the different types of oxalate feeds separately and train separate models for the different types of oxalate feeds), or b) trained on more data, and thus able to learn more (i.e., when we consider the different types of oxalate feeds together and train the model on the joint data). For example, classification techniques may allow us to distinguish between solid and solution runs. In this case, given a test run that we determine to be either a solid or solution run, we can train our model on runs that are more representative of the run we are studying.

In this paper, we looked at the ability of a BASS model to predict the processing conditions of particles of Pu, given a set of morphological characteristics by Bayesian Calibration via a Bayesian Adaptive Spline Surface model. Not only does this model allow for providing point estimates of a sample of SNM, but it also incorporates uncertainty into these predictions. This model proved to be sufficient at predicting processing conditions within a standard deviation of the observed processing conditions on average, when applied to a dataset of particles of Pu. In particular, we found that this model is best able to predict processing conditions 1, 2, and 4, and struggles most with predicting processing conditions 3, 7, and 8. By comparing Figures 3, 5, 7, we see that using both solid and solution runs to train the model (versus analyzing solid and solution runs separately) does not affect the ability of the model to predict the processing conditions of a new run, indicating that the BASS model is flexible enough to incorporate this information into its predictions.

By predicting the processing conditions of samples of SNM, we can begin to understand where material may have been produced, or by whom. While it remains that this the results of this method should not be used as the sole evidence for reaching a forensic conclusion, this analysis demonstrates that statistical models can aid the nuclear forensic community. In particular, these models help provide insight into the various sources of uncertainty and bias in the chemical processes. In addition, they allow forensic decision makes to numerically bound their observations, and quantitatively support their inferences and predictions. Future studies will aim to improve these results by considering the effects of particle sizes on the predictions, and will incorporate functional shape and texture data. Additionally, we would like to determine whether incorporating the discrepancy term in our calibration step would allow us to further decrease the RMSE of our predictions. Nevertheless, this methodology provides a solid foundation for future work in predicting processing conditions of particles of SNM using only their morphological characteristics.