Radiological mapping and uncertainty quantification by a fast Microcanonical Langevin Monte Carlo sampler

Abstract

Radiological mapping plays a critical role in nuclear emergency response and environmental management activities. A radiation image, representing the spatial and intensity distribution of the radioactivity, is reconstructed from the radiation data and associated contextual information. Typical image reconstruction methods, such as Maximum Likelihood Expectation-Maximization (ML-EM), only provide point estimates of the pixel or voxel activities without associated uncertainties. Here, we apply a new Microcanonical Langevin Monte Carlo (MCLMC) sampler for radiation image reconstruction and uncertainty quantification. The MCLMC sampler properties are first tested with synthetic radiation images. Methods to obtain the radiation distribution estimate and the associated uncertainty from the samples drawn by MCLMC are discussed. Given sufficient measurement statistics, the radiation distribution estimate obtained from MCLMC results closely resembles the ground truth with less risk of over- or under-fitting compared to ML-EM. When MCLMC is run in parallel on a GPU, it can converge in about 10 seconds for an image with $10^3$--$10^4$ pixels, which is significantly faster than other comparable Markov Chain Monte Carlo (MCMC) samplers. We also tested MCLMC on a dataset from a real distributed source radiological mapping campaign, and the reconstructed results agree well with ground truth. The fast MCLMC sampler therefore enables improved imaging accuracy and prompt uncertainty quantification for reconstructed radiation images, which can better inform decision-making in response to radiological events.

Figures (15)

• Figure 1: Synthetic radiological mapping scenario in which a detector rasters over a field containing a distributed radiological source consisting of three plumes. Each dot represents a $t_i=5$ s duration measurement position and is colorized by the measured counts in that time interval. The pixel sizes are set to 1 m$^2$. The detector trajectory starts from the green dot near $(0, 0)$ and stops at the red dot near $(0, 24)$. The $\mathsf{X}$ and $+$ represent a typical pixel in the background region and in the source region, respectively.
• Figure 2: Reconstruction by ML-EM with different number of iterations. Total measurement time $T=10$ minutes. Top row: reconstructed intensity distributions. Bottom row: measured counts vs. mean counts calculated from Eq. \ref{['eq:lambda']} using the reconstructed intensities. PSNR/SSIM is 11.79/0.23, 16.61/0.59, 23.50/0.89, 19.21/0.64, respectively, for 10, 100, 1000, 10000 iterations.
• Figure 3: Performance metrics vs. number of ML-EM iterations for different total measurement times $T$. (a) Poisson loss calculated using the ML-EM reconstructed intensities. (b) SSIM and (c) PSNR from comparison of reconstructed image at different number of ML-EM iterations against the ground truth. The overall highest PSNR = 29.79 and SSIM = 0.97 are obtained at 3000 iterations for $T=1000$ minutes.
• Figure 4: Prior and posterior distributions of MCLMC for $T=100$ minutes. (a), (b) Maps showing the mean and standard deviation of the truncated Gaussian prior distribution for all pixels, from $10$ ML-EM iterations. The $\mathsf{X}$ and $+$ represent a typical pixel in the background region and in the source region, respectively (same location as in Fig. \ref{['fig:synthetic_data_ground_truth']}). (c), (d) Marginal posterior distributions for a pixel in background/radioactive source region shown in (a). HDI: Highest Density Interval. GT: ground truth.
• Figure 5: Intensity and uncertainty maps. Top 4 rows: intensity maps showing different properties of the posterior distribution. Row 5: $68\%$ Highest Density Interval of marginal distributions for individual pixels. Row 6: measured counts vs. mean counts calculated from Eq. \ref{['eq:lambda']} using the value of mean, marginal mode and joint mode in the top 3 rows. Image quality metrics are given in Table \ref{['tab:mean_vs_mode_at_different_measurement_time']}.