Uncertainty Estimation of the Dose Rate in Real-Time Applications Using Gaussian Process Regression

Major standard organizations have addressed the issue of reporting uncertainties in dose rate estimations. There are, however, challenges in estimating uncertainties when the radiation environment is considered, especially in real-time dosimetry. This study reports on the implementation of Gaussian process regression based on a spectrum-to-dose conversion operator (i.e., G(E) function), the aim of which is to deal with uncertainty in dose rate estimation based on various irradiation geometries. Results show that the proposed approach provides the dose rate estimation as a probability distribution in a single measurement, thereby increasing its real-time applications. In particular, under various irradiation geometries, the mean values of the dose rate were closer to the true values than the point estimates calculated by a G(E) function obtained from the anterior–posterior irradiation geometry that is intended to provide conservative estimates. In most cases, the 95% confidence intervals of uncertainties included those conservative estimates and the true values over the range of 50–3000 keV. The proposed method, therefore, not only conforms to the concept of operational quantities (i.e., conservative estimates) but also provides more reliable results.


Introduction
The concepts of equivalent dose and effective dose were first introduced by the International Commission on Radiological Protection (ICRP) in order to provide recommendations and guidelines for the protection of people and the environment in an integrated manner in all exposure situations [1][2][3]. However, given that these concepts are not measurable quantities, the International Commission on Radiological Units and Measurements (ICRU) defined a few measurable operational quantities to establish convenient and appropriate evaluations of an equivalent and effective dose [4,5]. In cases of strongly penetrating radiations, such as gamma rays and neutrons, an adequate operational quantity for monitoring a specific area is defined by the ambient dose equivalent H*(10) (hereafter referred to as "ambient dose rate" and used interchangeably with the term "dose rate"). The ICRU defined the ambient dose rate as, "The dose equivalent at a point in a radiation field that would be produced by a corresponding expanded and aligned field in the ICRU sphere at a depth of 10 mm on the radius vector opposing the direction of the aligned field." The response of the ambient dose rate is highly dependent on photon energy and the angle of radiation incidence. Therefore, the requirement of IEC 60846:2009 advises that the relative response (2) where ∅(E i ) is the fluence rate at the energy of E i . In real conditions, the integral of continuous energies should be changed to the sum of discrete energies (or equivalently, the number of channels N; in this case, N = 869).
Sensors 2020, 5, x FOR PEER REVIEW 3 of 12 Figure 1. Illustration of the proposed concept of dealing with uncertainty in dose rate estimation compared with the conventional method.

G(E) Function
A general description of the ( ) function in terms of the dose conversion coefficient ℎ( ) at mono-energy and the response function of a detector ( , ), which represents the photon of energy depositing energy into the detector, can be expressed as where and are the minimum and maximum detectable energies deposited in the detector, respectively. The total dose rate ( ) in multi-energy radiation conditions can, therefore, be represented as Consequently, the total dose rate can be directly estimated using the G(E) function and a measured spectrum. According to related studies [7][8][9]11], the G(E) function can be expressed as Sensors 2020, 20, 2884 4 of 11 where A(K) is a parameter, K max is the number of terms, and M is constant? The values for K max and M were set to 7 and 0, respectively. It should be noted that the optimization of these parameters is not the main concern of this study. Dose rate can be represented by combining Equations (3) and (4): To compute A(K), it is required to obtain spectra with known mono or multiple energies and the corresponding dose rates. The availability of energy sources limits actual experiments; however, Monte Carlo simulations allow for the use of any energy, so corresponding dose rates can be calculated, given that detector geometry has been properly defined. Finally, A(K) were obtained using the gradient descent method [11].

GP Regression
This section briefly introduces the concept of GP regression employed for implementation purposes to deal with uncertainty in dose rate estimation. More details can be found in [19,20].
The GP is expressed as a distribution over functions for which any finite subset of variables has a joint multivariate Gaussian distribution. Since the GP is described by Gaussian distribution, it is parameterized by its mean function m(x) and positive definite covariance function k(x, x ), also known as a kernel function: Typically, m(x) is set to 0 to avoid expensive computations in posterior distribution and make inferences only via the kernel function. The kernel function takes two indices x and x and returns their corresponding modeled covariance. By choosing an appropriate kernel function, it is possible to incorporate assumptions such as smoothness and likely patterns that are expected in the data. A popular choice of the kernel is the radial basis function kernel, where two points are exponentially correlated, depending on the distance between them.
The main assumption in GP modeling is that output y is an observation of f (x) that has been corrupted by Gaussian noise : where noise term is assumed independent and identically distributed with zero means. Hence f (x) is a latent variable whose posterior distribution will be inferred after observing new samples at various locations in the domain. The resultant inference is called GP regression. Suppose training outputs y t have been observed and predictions for test outputs f * have been made. They then follow a joint normal distribution: where X t and X * are the design matrices for training and test data, respectively. K(X t , X t ) represents the covariance matrix between all points observed so far in the training data, which is similarly true for other covariance matrices of K(X t , X * ), K(X * , X t ), and K(X * , X * ). I is an identity matrix whose diagonal elements and off-diagonal elements are 1 and 0, respectively? Conditioning f * on the observation y t p f * |X t , y t , X * , the predictive distribution of test points with respect to the mean and covariance matrix can be written as: Sensors 2020, 20, 2884

of 11
Consequently, the estimation of the posterior mean and covariance are involved in calculating four different covariance matrices.

Monte Carlo Modeling and Simulation
A Monte Carlo N-Particle Transport Code (MCNP6) [21] was used to validate the proposed method. A schematic of the MCNP6 model used for simulations is illustrated in Figure 2. The 5.08 × 5.08 cm (diameter × height) NaI(Tl) crystal with a density of 3.6 g cm −3 was covered by a MgO reflector with a density of 2 g cm −3 , which was surrounded by aluminum with a density of 2.7 g cm −3 . A 20-mm-thick aluminum plate was placed behind the crystal to mimic a phenomenon where photons are scattered or backscattered in a photomultiplier tube [22]. This is a reasonable assumption because it considers the scattered effects of the photomultiplier tube on a spectrum as the angle of irradiation direction changes. A parallel beam of photons distributed over a circular source was irradiated on the NaI(Tl) detector. In order to obtain G(E) g functions under various directions of irradiation, the directions of the detector were rotated against the circular source by specific angle α (i.e., 0, 45, and 90), where subscript g is the irradiation geometry that determines the G(E) function. For isotropic geometry, the detector was centered inside the spherical source, emitting fully isotropic irradiation of photons. matrix whose diagonal elements and off-diagonal elements are 1 and 0, respectively? Conditioning f * on the observation y (f * |X , y , X * ), the predictive distribution of test points with respect to the mean and covariance matrix can be written as: m (x) = (x, X ) (X , X ) + σ Ι y (9) (x, x ) = ( , x ) ( , X ) ( , ) + ( , x ) Consequently, the estimation of the posterior mean and covariance are involved in calculating four different covariance matrices.

Monte Carlo Modeling and Simulation
A Monte Carlo N-Particle Transport Code (MCNP6) [21] was used to validate the proposed method. A schematic of the MCNP6 model used for simulations is illustrated in Figure 2. The 5.08 5.08 cm (diameter height) NaI(Tl) crystal with a density of 3.6 g cm -3 was covered by a MgO reflector with a density of 2 g cm -3 , which was surrounded by aluminum with a density of 2.7 g cm -3 . A 20-mm-thick aluminum plate was placed behind the crystal to mimic a phenomenon where photons are scattered or backscattered in a photomultiplier tube [22]. This is a reasonable assumption because it considers the scattered effects of the photomultiplier tube on a spectrum as the angle of irradiation direction changes. A parallel beam of photons distributed over a circular source was irradiated on the NaI(Tl) detector. In order to obtain ( ) functions under various directions of irradiation, the directions of the detector were rotated against the circular source by specific angle α (i.e., 0°, 45°, and 90°), where subscript is the irradiation geometry that determines the ( ) function. For isotropic geometry, the detector was centered inside the spherical source, emitting fully isotropic irradiation of photons. Since an actual spectrum is influenced by the broadening effect due to the statistical variation of the scintillation light signals and various electronic sources of noise, a simulated spectrum must be modified to accommodate such effects. This can be regarded as a convolution process of an ideal spectrum with the kernel of a broadening filter. A Gaussian-energy broadening filter is often applied, meaning that a delta function type of peak becomes a Gaussian function with full-width at halfmaximum value (FWHM = 2.36 sigma). In the MCNP, a non-linear function with three parameters regarding FWHM is specified to apply broadening effects on the ideal spectrum. The optimal values of parameters obtainable from measured spectra were found using a genetic algorithm [23]. Since an actual spectrum is influenced by the broadening effect due to the statistical variation of the scintillation light signals and various electronic sources of noise, a simulated spectrum must be modified to accommodate such effects. This can be regarded as a convolution process of an ideal spectrum with the kernel of a broadening filter. A Gaussian-energy broadening filter is often applied, meaning that a delta function type of peak becomes a Gaussian function with full-width at half-maximum value (FWHM = 2.36 × sigma). In the MCNP, a non-linear function with three parameters regarding FWHM is specified to apply broadening effects on the ideal spectrum. The optimal values of parameters obtainable from measured spectra were found using a genetic algorithm [23]. Figure 3 shows the determined G(E) functions of the NaI (Tl) detector as a function of energy deposited in the crystal under four different irradiation geometries. As expected, the G(E) 0 function, i.e., G(E) function for the angle of incidence of 0, tended to yield higher values over the entire energy range than those for other directions of photons with respect to the detector axis. Additionally, there were relatively small differences between the values of G(E) functions for energies above 200 keV, showing good agreement with previous results [9,12]. This could be ascribed to the fact that the energy deposition of relatively high energy depends primarily on the volume of the crystal. On the Sensors 2020, 20, 2884 6 of 11 other hand, the values of the G(E) functions obtained from different types of irradiation geometry tended to relatively disperse, especially for energies below 200 keV. This is because the photons in that energy range have high interaction probabilities with the crystal, so the energy deposition in the crystal becomes proportional to the projected area of the crystal incident surface. These results suggest that the estimated dose rate may drift from the true value, especially in the low energy range, depending on the G(E) functions calculated by different types of irradiation geometry. deposited in the crystal under four different irradiation geometries. As expected, the ( ) ° function, i.e., ( ) function for the angle of incidence of 0°, tended to yield higher values over the entire energy range than those for other directions of photons with respect to the detector axis. Additionally, there were relatively small differences between the values of ( ) functions for energies above 200 keV, showing good agreement with previous results [9,12]. This could be ascribed to the fact that the energy deposition of relatively high energy depends primarily on the volume of the crystal. On the other hand, the values of the ( ) functions obtained from different types of irradiation geometry tended to relatively disperse, especially for energies below 200 keV. This is because the photons in that energy range have high interaction probabilities with the crystal, so the energy deposition in the crystal becomes proportional to the projected area of the crystal incident surface. These results suggest that the estimated dose rate may drift from the true value, especially in the low energy range, depending on the ( ) functions calculated by different types of irradiation geometry.  Figure 4 shows the posterior mean of the GP model as well as its probabilistic nature in the form of a 95% confidence interval using the data from previously-determined ( ) functions; the ( ) functions are also illustrated in this figure for comparison. The result shows that the GP model no longer has a single value for energy but a distribution (i.e., Gaussian distribution) indexed by energy. In addition, the entire data points of ( ) functions determined under different types of irradiation geometry were found inside the 95% confidence region of the posterior. In particular, the relative vertical width of the confidence region with regard to the mean value tended to increase as it moved to the low energy range, especially for energies below 200 keV, to account for variations induced by radiation energy and direction of radiation incidence. It should be noted that the absolute values of the confidence region obtained from the GP model over the entire energy range are almost similar.  Figure 4 shows the posterior mean of the GP model as well as its probabilistic nature in the form of a 95% confidence interval using the data from previously-determined G(E) functions; the G(E) functions are also illustrated in this figure for comparison. The result shows that the GP model no longer has a single value for energy but a distribution (i.e., Gaussian distribution) indexed by energy. In addition, the entire data points of G(E) functions determined under different types of irradiation geometry were found inside the 95% confidence region of the posterior. In particular, the relative vertical width of the confidence region with regard to the mean value tended to increase as it moved to the low energy range, especially for energies below 200 keV, to account for variations induced by radiation energy and direction of radiation incidence. It should be noted that the absolute values of the confidence region obtained from the GP model over the entire energy range are almost similar.       Figure 5 shows an example of independent realization functions (i.e., G(E) GPR functions) randomly sampled from the GP model. As expected, each G(E) GPR function represented a different path because of the randomness of the stochastic process, fluctuating around the mean of the GP model. This implies that multiple dose rate values can be calculated using the G(E) GPR functions multiplied by the observed spectrum, resulting in not only the best dose rate estimate but also its uncertainty, which might contain the true value. The mean of the GP model was lower than that of the G(E) 0 function, which generally overestimates dose rates, so it nearly coincided with the G(E) ISO function (see Figure 4). That is, the mean dose rate values and the dose rates estimated by the G(E) ISO function might be in good agreement. This result is quite promising because isotropic geometry can be a reasonable assumption for irradiations often received from naturally occurring radioisotopes in homes or the surrounding environments [9,16]. . Gaussian process (GP) regression using the data from the previously-determined ( ) functions under the angles of incidence of 0° (black dashed line), 45° (green dotted line), and 90° (cyan short-dashed line), and isotropic geometry (gray dash-dotted line). The blue solid line represents the mean of the GP model. The blue shaded area denotes a 95% confidence interval. The previously determined ( ) functions are also illustrated for comparison. The inset shows the same graph on a linear scale. Figure 5 shows an example of independent realization functions (i.e., ( ) functions) randomly sampled from the GP model. As expected, each ( ) function represented a different path because of the randomness of the stochastic process, fluctuating around the mean of the GP model. This implies that multiple dose rate values can be calculated using the ( ) functions multiplied by the observed spectrum, resulting in not only the best dose rate estimate but also its uncertainty, which might contain the true value. The mean of the GP model was lower than that of the ( ) 0° function, which generally overestimates dose rates, so it nearly coincided with the ( ) function (see Figure 4). That is, the mean dose rate values and the dose rates estimated by the ( ) function might be in good agreement. This result is quite promising because isotropic geometry can be a reasonable assumption for irradiations often received from naturally occurring radioisotopes in homes or the surrounding environments [9,16].

Dose Rate Uncertainty Estimation
To validate the proposed method, various spectra were obtained for mono-energy over the range of 50-3000 keV at certain intervals with various geometries, e.g., the angles of incidence of 0 • , 45 • , and 90 • , and isotropic geometry. To calculate the uncertainty of the dose rate (i.e., 95% confidence interval), 100 G(E) GPR functions were randomly sampled from the GP model each time. Figure 6 shows a comparison of the energy response normalized to the energy of 622 keV emitted by Cs-137, estimated with G(E) GPR functions and the G(E) 0 function for the spectra obtained at the angle of incidence of 0 • . Here, the energy response was calculated by having the ratio of the estimated dose rate to the true dose rate at certain energy divided by the same ratio at the energy of 662 keV. That is, an increase in the value of the energy response suggests that the dose rate is overestimated, or vice versa. As we can see from the figure, the energy responses for the G(E) 0 function were reasonably close to one, which means that the estimated values of the dose rate and true values were in good agreement. This is because the test spectra were acquired under the same condition used for the G(E) 0 function calculation. These results are not as good as those that were reported by previous studies, especially for the low energy range. Nonetheless, it is worth emphasizing that the purpose of this study was to propose a concept that would make it possible to deal with uncertainty existing in the dose rate by taking into account the relative response of radiation energy and the direction of radiation incidence. For the dose rates estimated with G(E) GPR functions, the mean values deviated slightly more from the reference value 0 for energies below 200 keV. The reason is that there might be a discrepancy between the mean values of G(E) GPR functions sampled from the GP model and those from the G(E) 0 function, which is better suited with respect to the test spectra. In contrast, the proposed method was able to provide the uncertainty and the mean of the dose rate. As expected, the relative uncertainty tended to increase with a decrease in energy. In particular, the relative uncertainty increased sharply for energies under 200 keV because the GP model had relatively wide intervals for that energy range. In addition, Sensors 2020, 20, 2884 8 of 11 the 95% confidence interval of relative uncertainty mostly included the true value and the conservative values obtained with the G(E) 0 function. radiation incidence. For the dose rates estimated with ( ) functions, the mean values deviated slightly more from the reference value 0 for energies below 200 keV. The reason is that there might be a discrepancy between the mean values of ( ) functions sampled from the GP model and those from the G(E) ° function, which is better suited with respect to the test spectra. In contrast, the proposed method was able to provide the uncertainty and the mean of the dose rate. As expected, the relative uncertainty tended to increase with a decrease in energy. In particular, the relative uncertainty increased sharply for energies under 200 keV because the GP model had relatively wide intervals for that energy range. In addition, the 95% confidence interval of relative uncertainty mostly included the true value and the conservative values obtained with the G(E) ° function.  Figures 7-9 show the same comparison of energy response for the spectra obtained at the angles of incidence of 45° and 90°, and isotropic geometry. Although similar trends were observed in the non-zero angle of incidences, the estimated dose rate values obtained with the G(E) ° function were overestimated, especially for energies below 600 keV for all geometries in which the test spectra were acquired, which was expected. This shows the reason that the G(E) ° function is generally used to provide conservative dose estimates. In particular, the dose rate overestimation for the test spectra, assuming that photons were irradiated under the angle of incidence of 45°, is as high as 50% at 70 keV (see Figure 7). This is because the largest projected area of the crystal incident surface is generated at that specific angle, which increases interaction probabilities, especially for photons at low energies. show the same comparison of energy response for the spectra obtained at the angles of incidence of 45 and 90, and isotropic geometry. Although similar trends were observed in the non-zero angle of incidences, the estimated dose rate values obtained with the G(E) 0 function were overestimated, especially for energies below 600 keV for all geometries in which the test spectra were acquired, which was expected. This shows the reason that the G(E) 0 function is generally used to provide conservative dose estimates. In particular, the dose rate overestimation for the test spectra, assuming that photons were irradiated under the angle of incidence of 45 • , is as high as 50% at 70 keV (see Figure 7). This is because the largest projected area of the crystal incident surface is generated at that specific angle, which increases interaction probabilities, especially for photons at low energies. Likewise, the mean values of the dose rate estimate, with G(E) GPR functions, showed similar trends but were less overestimated. Furthermore, the uncertainties of those estimates included not only the conservative values estimated by the G(E) 0 function but also the true values in most situations. The uncertainty calculated by the proposed method is much more reliable than those that provide only a point estimate, because in the real world, it is not possible to know how far an estimate is from the true value.
Sensors 2020, 5, x FOR PEER REVIEW 9 of 12 but were less overestimated. Furthermore, the uncertainties of those estimates included not only the conservative values estimated by the G(E) ° function but also the true values in most situations. The uncertainty calculated by the proposed method is much more reliable than those that provide only a point estimate, because in the real world, it is not possible to know how far an estimate is from the true value.

Discussion
This paper presented how GP regression can be applied to deal with dose rate uncertainty in real-time applications. The results demonstrated that the proposed approach is much more reliable and robust in comparison with existing methods. For conventional methods, a way to determine uncertainties is the statistical analysis of a series of observations (i.e., Type A uncertainties). In this case, however, they ignore other components of uncertainty determined by scientific judgment based on published data (e.g., G(E) function). Furthermore, in cases of real-time applications, the estimation of Type A uncertainties is practically impossible, so the uncertainty associated with dose rate can, therefore, not be reported. Although previous studies attempted to estimate dose rate uncertainty, they simply averaged G(E) functions for the angles of incidence of 0 • and 90 • over the entire energy range and neglected the combined effects of the direction of radiation incidence and photon energy on the detector response [24]. In contrast, this study constructed a GP model that considers the relative responses to various irradiation geometries and energy to provide the mean and its uncertainty for the estimated dose rate based on a single spectrum. In addition, the mean dose rate values were not heavily overestimated under various irradiation geometries, and the 95% confidence interval of uncertainties included the conservative estimates obtained with the G(E) 0 function and true values. An estimated value without a statement about its associated uncertainty is less informative because it does not make it possible to quantify the potential risk arising from radiation exposure and does not indicate the precision of the estimate. Lastly, the calculated uncertainty allows for a quantitative comparison with results reported by other investigators, which enables its assessment.

Conclusions
This work presented a new method for dose rate estimation that uses GP regression. The presented results confirmed that numerous G(E) GPR functions that account for the relative responses to radiation energy and irradiation directions could be randomly sampled from a GP model, making it possible to deal with uncertainty in dose rate estimation for real-time applications. While the conventional method overestimates the dose rate by as much as 50% under different irradiation geometries, the mean values of the dose rate estimated with G(E) GPR functions were closer to the true value. Furthermore, the overestimated values obtained with the G(E) 0 function and the true values were mostly found within the 95% confidence interval of uncertainty. Therefore, the proposed method conforms to the concept of operational quantities present in conservative estimates and provides a more reliable dose rate estimation.