Path Planning for Localization of Radiation Sources Based on Principal Component Analysis

Abstract

In this paper, we propose a path planning method for the localization of radiation sources using a mobile robot equipped with an imaging gamma-ray detector, which has a field of view in all directions. The ability to detect and localize radiation sources is essential for ensuring nuclear safety, security, and surveillance. To enable the autonomous localization of radiation sources, the robot must have the ability to automatically determine the next location for gamma ray measurement instead of following a predefined path. The number of incident events is approximated to be the squared inverse proportional to the distance between the radiation source and the detector. Therefore, the closer the distance to the source, the shorter the time required to obtain the same radiation counts measured by the detector. Hence, the proposed method is designed to reduce this distance to a position where a sufficient number of gamma-ray events can be obtained; then, a path to surround the radiation sources is generated. The proposed method generates this path by performing principal component analysis based on the results obtained from previous measurements. Both simulations and actual experiments demonstrate that the proposed method can automatically generate a measurement path and accurately localize radiation sources.

1. Introduction

Although radiation sources are applied in various fields such as energy, medical, and manufacturing, gamma rays, which are emitted from high-dose radiation sources, such as those exceeding 5 Sv, can cause serious damage to the human body. Consequently, the application of radiation sources is strictly regulated by law. However, there are cases where radiation sources are stolen or leaked out, and this has become a critical social problem.

There have been some cases in which radiation sources were stolen. In Georgia, five suspects were arrested for the illegal trade of radioactive uranium [1]. Accordingly, illegally traded radiation sources possibly fall into the hands of terrorist organizations, and are harnessed to create dirty bombs comprising explosives and radiation sources, or silent source attacks that leave high-dose radiation sources to harm humans. In addition, radiation sources can cause great damage when used in urban areas or important facilities. However, the detection and removal of radiation sources currently relies on human workers, even though it is extremely dangerous to approach and investigate these sources [2]. Owing to the high risk of exposure to radiation, it is necessary to develop a radiation source exploration system with mobile robots in urban areas and various facilities [3,4].

Accordingly, a method for the localization of radiation sources in an environment using a radiation detector has been studied. Radiation caused by radioactive decay can be categorized into electromagnetic radiation, such as gamma rays, and particulate radiation, such as alpha and beta rays. Although, alpha and beta rays can be easily blocked by thin materials, gamma rays exhibit high penetration power. Hence, gamma-ray detection is the conventional method for detecting radiation sources. The conventional search approach, such as uniform search, is based on a predefined path that scans the area. Ziock and Goldstein proposed a uniform search method to detect a radiation source [5]. Klimenko et al. also proposed a strategy of searching along a predefined path while adaptively optimizing detection time [6]. In addition, information-driven approaches have been studied. Various information gains, such as entropy and Fisher information, have been proposed to determine the best next measurement position. Hutchinson et al. used the maximum entropy sampling concept to determine optimal measurement positions [7]. Ristic et al. designed a method for maximizing the information gain in the Renyi divergence sense to determine subsequent measurement positions and detection times [8]. Liu and Abbaszadeh proposed a double Q-learning algorithm [9], a data-driven approach. These studies were essentially designed to use a radiation detector more as a counter rather than as an imager. A counter-type detector counts the number of incident gamma rays or obtains the radiation dose rate. However, the incident gamma ray direction is not determined.

As an alternative to the aforementioned techniques, gamma camera-based searching, which allows for determining the incident gamma ray direction, has been studied. Kim et al. proposed a method for three-dimensional (3D) reconstruction of the radiation source distribution in an environment using a mobile robot equipped with a Compton camera [10]. Shikaze et al. proposed a method for estimating the radiation source distribution using an unmanned helicopter equipped with a Compton camera [11]. However, these techniques are based on a predefined path. Search path planning using a gamma camera has not been extensively studied. Because high-dose radiation sources exceeding 5 Sv can cause significant damage to the human body, it is necessary to detect and eliminate it at an early stage. Therefore, an investigation that sequentially scans the entire environment is ineffective when the target environment is large.

To enable autonomous radiation source localization, a mobile robot requires the ability to automatically determine the next location for gamma ray measurement instead of following a predefined path. Based on this requirement, we propose a path planning method for the localization of radiation sources using a mobile robot equipped with a gamma-ray detector, which has a field of view (FOV) in all directions. The number of incident events is approximated to be the squared inverse proportional to the distance between the radiation source and the detector. Therefore, the lower the distance to the source, the shorter the time required to obtain the same radiation counts measured by the detector. The proposed method is designed to reduce this distance to a position where a sufficient number of gamma-ray events can be obtained, which then generates a path that surrounds the radiation sources. Figure 1 presents a conceptual image of the measurement path, which is the objective of the proposed method. The proposed method generates the above measurement path by performing principal component analysis (PCA) based on the measurement results obtained, up to the current time. We demonstrate that the proposed method can automatically generate a measurement path and accurately localize radiation sources via simulations and experiments.

The remainder of this paper is organized as follows. Section 2 presents the schematic representation and details of the proposed method. Section 3 presents the experiments and results. Here, the results of the simulation under the three conditions and the experiments are provided. Finally, Section 4 presents the conclusions of this study.

2. Proposed Method

2.1. Overview

To realize a mobile robot equipped with a gamma-ray detector that sequentially repeats gamma-ray measurements and localize radiation sources, the proposed method determines the position of the next measurement point based on the results of previous measurements. The radiation source distribution is estimated by considering large environments such as urban areas, public facilities, or shopping malls. Moreover, we assume a situation where strong radiation sources are located in the target environment because it would not be necessary to use an inspection robot if the radiation source intensity is low and human workers can be deployed. Our strategy is as follows. First, human workers investigate the areas with low intensity and limit the scope of the inspection. After the target area is determined, the robot is placed near the target area with high intensity. The robot conducts inspection and identifies the exact radiation source locations. Although the region with high count rates is avoided by human inspectors, rapid detection and elimination of radiation sources are still required in cases such as dirty bombs comprising explosives and radiation sources. In addition, it is undesirable to limit the area for long periods even if there is no threat to the human body. Therefore, an investigation that sequentially scans the entire environment is ineffective when the target environment is large. Based on these requirements, the proposed method was designed, as illustrated in Figure 2.

The proposed system adopts a Compton camera, which uses Compton scattering kinematics to localize a radiation source, as the gamma-ray detector. Compared with other types of gamma cameras, such as pinhole or coded-aperture cameras, the Compton camera has a wide FOV [12,13]. Particularly, our all-around-view Compton camera, which employs a 3D pixel array CdTe detector, enables measurement in all directions [14]. In [15], it was concluded that a gamma camera with a wide FOV would be preferable for wide-area exploration. The details of the all-around view Compton camera are described in Section 2.2.

After the gamma-ray measurement, the reconstruction of the radiation source distribution to a 3D space is conducted based on the measurement results. The proposed system applies a simple back-projection method, which is a general reconstruction method. Details on the reconstruction of the radiation source distribution are provided in Section 2.3.

Generally, the lower the distance to the source, the shorter the time required to obtain the same radiation counts measured by the detector. Therefore, if the number of gamma-ray measurement events is insufficient at the current position, a mobile robot is required to approach the radiation sources until the number of events reaches a sufficient value. Moreover, to accurately localize the radiation sources, it is necessary to follow a path that surrounds them. Hence, the proposed method is designed to generate a path that surrounds the radiation source after approaching a distance where a sufficient number of measurement events can be obtained. The proposed method applies PCA to generate this path based on the reconstruction results of the radiation source distribution, which is constructed by simple back-projection. Comprehensive details on path planning are presented in Section 2.4.

The above process is repeated until the measurement path is closed. Then, the localization of the radiation sources is conducted. The final reconstruction of the radiation source distribution is performed by solely adopting the data obtained from where the number of incident events is significant, and then the location of radiation sources is finally estimated as the weighted average of the values added to pixels in the projection space. Because of simple back projection, each projection space pixel has a value that indicates how many times it is overlapped. Details regarding the projection space pixel value are provided in Section 2.3 and Equation (2). When multiple radiation sources are positioned, the k-means algorithm is applied to separate the area of the reconstructed distribution corresponding to each radiation source. Then, to localize each source, the weighted average is calculated relative to each separated area. The details are described in Section 2.5.

2.2. Gamma-Ray Measurement Using All-Around View Compton Camera

Various gamma-ray detectors have been proposed; however, they can all be classified into three types: gas-filled, scintillation, and semiconductor detectors [16]. In [8], the Geiger-Muller counter, which is a gas-filled detector, was used to obtain the radiation dose rate. However, scintillation and semiconductor detectors have adequate sensitivity and energy resolution for radiation dose rate measurement and gamma spectrometry [12]. In [17], the authors used a CsI(Tl) scintillator and a Silicon Photomultiplier as a detector. By counting the number of incident gamma rays, measurements were conducted. When high-energy resolution measurements are required, semiconductor detectors should be considered [12]. Using CdTe or CZT ensures high-energy-resolution gamma spectroscopy, although such detectors are expensive.

A gamma camera consists of one or multiple detectors, and it can be classified into three types: pinhole, coded-aperture, and Compton [12]. Although pinhole cameras have high angular resolution, they are based on heavy collimators. Therefore, they are not appropriate for mobile robots. In [18,19], a coded-aperture camera was used for gamma-ray imaging. The authors used an active mask and achieved a wide FOV. In [10,20], Compton cameras were used to localize radiation sources because of their wide FOV. Thus, because of its advantages, the Compton camera is adopted in our system.

As illustrated in Figure 3, the Compton camera consists of two layers: a scatterer and an absorber. The gamma ray, which interacts with the elements in the first and second layers, is recorded as the measurement data. The data contain the positions of the elements interacting with gamma rays at the two layers, respectively, with the deposited energy of radiation. Let A and B be the positions where the gamma ray is deposited on the scatterer and absorber, respectively, while $\theta$ is the scattering angle of the gamma ray with respect to the straight line connecting A and B. According to the kinematics of Compton scattering, the scattering angle is calculated as follows:

(1)$cos\theta =1-m_e{c}^2(\frac{1}{E_B}-\frac{1}{E_A+E_B}),$

where $m_e$, $E_A$ and $E_B$, and c represent the mass of an electron, deposited energy at each layer, and speed of light, respectively. The positions at which the two elements interact with the incident radiation are considered the coordinates relative to the internal coordinate system of the Compton camera. A cone representing the positions of the radiation sources on the surface is estimated via two consecutive interaction positions of the elements and scattering angle. As the two interaction positions A and B and the scattering angle $\theta$ are given, the radiation source lies on the cone surface with an axis direction $\overrightarrow{BA}$ and a half opening angle $\theta$. This cone is called a Compton cone, and the radiation source position is estimated on the multiple overlapping Compton cone surfaces via back projection.

Although the all-around view Compton camera cannot clearly differentiate between the scatterer and absorber, the principle is the same as that of a general Compton camera. A Compton cone can be obtained from the two consecutive interaction positions of the elements, A and B, with the scattering angle. However, it is difficult to determine whether the two interacting elements are absorbers or scatterers because the time between the scatter and the subsequent absorption ranges from ps to ns. Therefore, for each incident event, two types of calculations are performed for each interacting element in the case of an absorber and a scatterer. The all-around view Compton camera projects Compton cones on both sides of a straight line connecting the two interacting elements, as illustrated in Figure 4. In addition, an energy window is set to identify the type of radionuclide based on the gamma-ray energy spectrum emitted from the radiation source. Through this approach, the effect of natural radioactive materials contained in the target area can be reduced. Our previous paper comprehensively describes a sensing mechanism that adopts the all-around view Compton camera [14].

2.3. Reconstruction of Radiation Source Distribution via Simple Back-Projection

This section describes how to use the all-around view Compton camera to reconstruct the measurement data into a 3D radiation source distribution via simple back-projection, which is a general reconstruction method [21]. Several studies have been conducted to propose a methodology for imaging the distribution of radioactivity [22,23,24]. Among them, maximum likelihood expectation maximization (MLEM) is widely used to reconstruct data obtained from a radiation detector [25,26]. The MLEM can estimate the position of the radiation source with high accuracy via iterative calculations. Hellfeld et al. proposed a system to enable real-time online reconstruction through MLEM [19]. This approach can be applied to our system and can significantly improve the performance. Future work may include employing the MLEM in our system. In this study, simple back-projection is adopted as a basic method for reconstructing radiation source distributions.

In the case of using a normal Compton camera with scatterer and absorber layers, let ${\theta }_i$ be the scattering angle at the measurement point index i when N denotes all the measurement points. Then, let $\mathit{k}$ be the vector connecting the interacting elements on the absorber and scatterer layers, as illustrated in Figure 5. In addition, let $\mathit{m}$ be the vector connecting the interacting element on the scatterer and voxels in the projection space. The pixel value $f\left(V\right)$ of the voxels can be defined by the delta function in the simple back-projection method as follows:

(2)$f\left(V\right)=\sum _{i=1}^N\delta (\frac{\mathit{k}·\mathit{m}}{\left|\mathit{k}\right|\left|\mathit{m}\right|}-cos{\theta }_i).$

As described in Section 2.2, it is difficult to distinguish between an absorber and a scatterer because the process, from the incident to their absorption, takes place at a very high speed. Therefore, for each incident event, two types of calculations are performed for each interacting element in the case of an absorber and a scatterer. The measurement data G can be represented by the following matrix:

(3)$G=\left(\begin{array}{cccc}{P}_{s1}& P_{a1}& E_{s1}& E_{a1}\\ P_{s2}& P_{a2}& E_{s2}& E_{a2}\\ ⋮& ⋮& ⋮& ⋮\\ P_{sN}& P_{aN}& E_{sN}& E_{aN}\end{array}\right),$

where each row represents the data measured in one incident event, while $P_s$ and $P_a$ are the positions of the interacting elements on the scatterer and absorber, respectively. $E_s$ and $E_a$ represent the deposited energies at $P_s$ and $P_a$, respectively. Because the intensity of the radiation source is constant, the measurement data that do not satisfy the following equation can be ignored:

(4)$|E_s+E_a-E_{rad}|<\Delta E,$

where $\Delta E$ is a margin that considers the energy resolution of the Compton camera. Based on the equations above, a simple back-projection method is implemented as described in Algorithm 1.

Algorithm 1: Simple Back Projection

After sensing the gamma rays at each measurement point, the measurement data are projected into the 3D voxel space via the simple back-projection method. The voxels in 3D space exhibit the number of times they overlap with the Compton cone, which occurs for each gamma-ray event. To compress the data, voxels with high values are extracted. However, the information with respect to the Z axis is not required for determining the movement direction of the robot, because the mobile robot can move in two-dimensional (2D) space. Therefore, the proposed method projects the voxel data onto a 2D plane. Finally, the radiation source distribution is obtained, as illustrated in Figure 6. The proposed method generates a measurement path based on the 2D information.

2.4. Path Planning via PCA

For the accurate localization of radiation sources, a mobile robot is required to closely approach radiation sources to acquire a sufficient number of incident events, which then moves along a path that surrounds the radiation sources. The proposed method generates the measurement paths via PCA. After the reconstruction of the radiation source distribution with the simple back-projection method, PCA of $f\left(V\right)$ being projected onto a 2D plane is conducted. The next measurement point is determined along the direction of principal components. The measurement path is defined as a set of the next measurement points, and the path is generated by repeating the above calculation.

PCA is a method for expressing variables distributed in space using specific indices. Specifically, it extracts the principal component axis that maximizes the variance of the data. The obtained directional component of this axis is called the first principal component. In this study, PCA is applied to 2D variables; thus, the axis orthogonal to the first principal component is uniquely determined. The directional component of this axis is called the second principal component.

The principal component axis can be obtained from the eigenvector of the variance-covariance matrix of the data, while the variance of the data in the direction along the principal component axis can be obtained as the eigenvalue. The procedure for performing PCA based on the pixels projected on the $xy$ plane is as presented follows:

The PCA is performed with the current measurement point as the coordinate origin. Therefore, it is necessary to perform a coordinate transformation because the reconstructed distribution is represented in world coordinates. Let $\mathit{X}=\left\{{\mathit{x}}_i,\cdots ,{\mathit{x}}_N\right\}$ be the coordinates of the voxels that have been translated. The coordinate origin is the current measurement point, and let $\mathit{X}$ be the coordinates of the pixels on the $xy$ plane. ${\mathit{x}}_k={\left[x_k{y}_k\right]}^T$ is defined by the index k, and the centroid $\overline{\mathit{x}}$ of these data can be represented as

(5)$\overline{\mathit{x}}=\frac{1}{N}\sum _{k=1}^N{\mathit{x}}_k.$

Then, the difference between $\overline{x}$ and each data is calculated:

(6)${\widehat{\mathit{x}}}_k={\mathit{x}}_k-\overline{\mathit{x}}\left(k=1,\cdots ,N\right).$

From here, the variance-covariance matrix $\mathit{C}$ is represented as:

(7)$\mathit{C}=\frac{1}{N}\sum _{k=1}^N{\widehat{\mathit{x}}}_k{\widehat{\mathit{x}}}_k^T.$

By performing an eigenvalue analysis of $\mathit{C}$, the eigenvectors, ${\mathit{e}}_1$ and ${\mathit{e}}_2$, corresponding to the eigenvalues of the variance-covariance matrices, ${\lambda }_1$ and ${\lambda }_2$, can be obtained. The eigenvalues are the variances of the data along the corresponding eigenvectors, and the eigenvectors are the principal component axes. ${\mathit{e}}_1$ represents the direction that provides maximum variance to the data and corresponds to the first principal component. In addition, ${\mathit{e}}_2$ is the second principal component orthogonal to ${\mathit{e}}_1$.

The proposed method determines the next measurement point along the principal component vectors. If the number of incident gamma-ray events at position i is less than the threshold, the next measurement point is determined along the first principal component. However, if the number of incident gamma-ray events at position i is greater than or equal to the threshold, the next measurement point is determined along the second principal component. Considering $\mathit{e}$ as the selected principal component vector, the next measurement point is determined as:

(8)${\mathit{x}}_{i+1}={\mathit{x}}_i+A\ast \mathit{e},$

where A represents the constant distance that the mobile robot moves in one step. The distance moved in one step is tuned according to the target area size and the inspection urgency. If the target area is large and the inspection must be completed in a short time, the step size is determined to be a large value. However, this can reduce the localization accuracy. Conversely, if the target area is small, the inspection is not urgent, and high localization accuracy is required, a smaller step size value is determined. The value is set once before the inspection begins. This process is repeated, and path planning is only completed when the loop is closed, as illustrated in Figure 7. Path planning is implemented as described in Algorithm 2. The next measurement position is determined based on the measurement at the current position, and the previous measurements do not affect the next measurement point. Notably, there are sign ambiguities in both eigenvectors, ${\mathit{e}}_1$ and ${\mathit{e}}_2$. The sign of ${\mathit{e}}_1$ is defined such that the robot moves toward the radiation sources. Meanwhile, the selected direction of ${\mathit{e}}_2$ is irrelevant as long as the same direction is kept during the path planning. Additionally, because the detector is fixed on the mobile robot, it can only rotate with respect to the Z axis. Consequently, the detector faces the movement direction of the robot, as shown in Figure 7. The red arrows represent the direction that the detector faces; the aforementioned direction at position ${\mathit{x}}_i$ is determined to be along $\overrightarrow{{\mathit{x}}_{i-1}{\mathit{x}}_i}$.

Although the proposed method primarily targets a single radiation source, it can support multiple radiation sources by tuning the threshold to switch from the first principal component direction to the second principal component direction. When specifying the position of each radiation source is required, a large threshold is set to generate a path that will approach each radiation source. However, when it is necessary to understand the overall distribution of the radiation sources, setting a small threshold will generate a path that greatly surrounds the radiation sources. Accordingly, the proposed method can generate an appropriate path by setting the threshold according to the objective. This is considered an advantage of the proposed method.

Algorithm 2: Order of measuring

2.5. Localization of Multiple Radiation Sources

After the final radiation source estimation, the k-means algorithm is applied to the 3D voxel values obtained through simple back-projection for localizing multiple radiation sources. In this study, the number of clusters in k-means clustering is manually determined based on the observation of the final radiation source estimation. However, the number of clusters can be automatically estimated using an information criterion, such as the Bayesian information criterion [27].

3. Experiments

Simulations and experiments were conducted to evaluate the proposed method. The performance of the proposed method significantly depends on the measurement results of radiation, and it is substantially influenced by the distributed conditions of the radiation sources. Therefore, the evaluation of the proposed method was conducted under multiple conditions, in which the distribution of radiation sources in each condition was different.

Simulation experiments were conducted under three conditions: (a) a single radiation source was positioned in the environment, (b) multiple radiation sources were concentrated in a specific area, and (c) multiple radiation sources were dispersed in the environment. In addition, experiments were conducted under the condition that a single radiation source was located in the environment.

3.1. Simulations

First, the all-around view Compton camera was implemented and the measurement data of gamma rays were obtained using the library called Geant4 [28,29]. Figure 8 illustrates the all-around view Compton camera implemented in Geant4. The camera consists of a 12 × 4 × 30 array of CdTe detectors whose width, height, and depth are 8, 12.5, and 2.21 mm, respectively. The width, height, and depth of the Compton camera with a cover are 138.6, 60, and 140 mm, respectively. Geant4 can simulate radiation tracks when interacting with matter based on the Monte Carlo method. The reliability of the library has been comprehensively validated and used in several studies [30]. Regarding the gamma-ray measurement using Geant4, we can obtain the data that present the position of the element interacting with the gamma ray and deposited energy at each element. Using the obtained data, a simple back-projection was performed to reconstruct the distribution of radiation sources using the previous measurements. The next measurement point was then determined via PCA. After all the measurements were completed, the final distribution of the radiation sources was reconstructed, and the positions of the radiation sources were localized. The simulation experiment is illustrated in Figure 9.

3.1.1. Condition (a) Single Radiation Source

A simulation experiment was conducted by arranging a single radiation source, which can be regarded as a point source, as illustrated in Figure 10. We positioned 1 MBq of cesium-137 15 m from the all-around view Compton camera. The threshold for switching the principal component vector was set to 400. Triplets and higher multiplicities of Compton events were not considered; only Compton event pairs were considered. The moving distance of the mobile robot at each time step between the measurement points was set to 5 m. The measurement at each point was performed for 60 s, and the robot stopped during the measurement. The background radiation was not considered.

Figure 11a illustrates the simulation result. SP represents the first measurement point where the robot started to move, while EP is the position at which the gamma-ray measurements were completed. The color bar on the right represents the value added to each pixel. The black and red arrows represent the directions of the first and second principal components, respectively. The ground truth of the position of the radiation source is represented as a red dot. The detector was fixed on the mobile robot; therefore, the detector orientation at position ${\mathit{x}}_i$ was determined along the robot’s movement direction, i.e., $\overrightarrow{{\mathit{x}}_{i-1}{\mathit{x}}_i}$. Figure 11b shows the detector orientation during the inspection. The detector orientation was zero when the X axis was facing upward, and clockwise rotation was defined as positive.

Because the number of incident events was insufficient at the first and second measurement points, the robot moved toward the first principal component, thus approaching the radiation source. Subsequently, the robot altered the moving direction toward the second principal component as a sufficient number of incident events was later obtained. Accordingly, a path was generated to surround the radiation source. Because the number of incident events fell below the threshold at the 7th and 11th measurement points, the robot moved toward the first principal component. After that, the measurement was completed at the 14th measurement point.

The final reconstruction of the radiation source distribution was performed by solely adopting the data obtained from the cases in which the number of incident events exceeded the threshold. It can be verified that the reconstructed distribution of the radiation source was estimated around the ground truth. In addition, the path was successfully generated, as intended. As the weighted average, the location of the radiation source was estimated as (−0.26 m, 0.23 m), while that of the ground truth was (0 m, 0 m). The estimation error was approximately 0.35 m. Based on this result, it was demonstrated that the proposed method can successfully generate the measurement path and accurately localize a single radiation source.

3.1.2. Condition (b) Multiple Radiation Sources Concentrated in a Specific Area

A simulation experiment was conducted in an environment where multiple radiation sources were locally concentrated, as illustrated in Figure 12a. The distance between the radiation sources was close; therefore, the identification of one location enables the elimination of all the radiation sources simultaneously. Three sources of cesium-137 (1 MBq) were positioned on a straight line, and the distance between the radiation sources was set to 1.5 m. The mobile robot started the measurement at a position 7 m from the central radiation source. The threshold for switching the principal component vector was set to 5000.

Figure 12b presents the simulation results obtained when the radiation sources are locally concentrated. Because the number of incident events fell below the threshold up to the 6th measurement point, a path approaching the radiation sources was generated along the direction of the first principal component. Subsequently, a path that moved toward the second principal component and surrounded the radiation sources was generated. At the 11th, 17th, and 24th measurement points, the number of incident events fell below the threshold, such that the robot moved toward the first principal component again to acquire sufficient events. The measurement was completed at the 26th measurement point.

Similar to the case of a single radiation source, the reconstruction of the radiation source distribution was performed by the simple back-projection method using only the data obtained when the number of incident events per interval exceeded the threshold. Here, it can be confirmed that the reconstructed distribution of radiation sources was estimated around the ground truth. In this case, it is difficult to separate the area corresponding to each radiation source using the k-means algorithm because the radiation sources were too close to each other. However, all the radiation sources can be eliminated simultaneously, even if the location of each radiation source is not specified. The estimated location of the sources was (0.10 m, 0.26 m). The localization error relative to the central radiation source was approximately 0.28 m. Based on this result, we demonstrated that the proposed method can successfully generate a measurement path in an environment where multiple radiation sources are concentrated in a specific area.

3.1.3. Condition (c) Multiple Radiation Sources Dispersed

A simulation experiment was conducted in an environment in which multiple radiation sources were dispersed. In this case, the location of each radiation source should be specified. Figure 13 presents the simulation conditions. Two radiation sources were placed in two locations: (−0.5 m, −2 m) and (0.5 m, 2 m). Two 1 MBq cesium-137 radiation sources were used. Then, the robot was positioned 7 m from the origin of the coordinates.

Under this condition, two cases were evaluated according to the threshold required to switch the principal component vector. Figure 14a,b present the results obtained when the threshold was set to 5000 and 2500, respectively.

When the threshold was set to 5000, a path was generated that moved in the direction of the first principal component up to the third measurement point, and in the direction of the second principal component, except for the 5th and 12th measurement points. In this case, a path was generated to solely surround radiation source 1. However, no enclosing path was generated for radiation source 2, and a path similar to the single radiation source condition was obtained.

However, the path to surround both radiation sources was generated when the threshold was set to 2500. The robot moved in the direction of the first principal component from the SP to the third measurement point. Then, it moved in the direction of the second principal component. Subsequently, at the 11th, 17th, and 22nd, the movement was performed in the direction of the first principal component to acquire a sufficient number of incident events. At the other measurement points, movement was performed in the direction of the second principal component. It was confirmed that an appropriate path was generated as intended.

In both cases, to reconstruct the radiation source distribution, the simple back-projection method was adopted using only the data obtained when the number of incident events exceeded the threshold. Consequently, when the threshold was set as 5000, a distribution was obtained only for a radiation source 1. The location of radiation source 1 was estimated as (−0.45 m, −1.8 m), which is compared to that of the ground truth, (−0.5 m, −2 m). Hence, the estimation error was approximately 0.21 m. The location of radiation source 2 was not estimated. In contrast, when 2500 was set as the threshold, a distribution was successfully obtained around the ground truth. After the k-means algorithm was performed to separate the reconstructed distribution into two areas, the location of radiation source 1 was estimated as (−0.48 m, −1.5 m), which is compared to that of the ground truth, (−0.5 m, −2 m). Hence, the estimation error was approximately 0.50 m. In addition, the location of radiation source 2 was estimated as (0.63 m, 2.3 m), which is compared to that of the ground truth, (0.5 m, 2 m) with an estimation error of approximately 0.33 m.

Gamma rays emitted from radiation source 1 significantly influence the reconstruction result more than the gamma rays emitted from the radiation source 2. This is a limitation of simple back-projection. By setting a large value as the threshold, the robot approached very closely to the radiation source 1. Consequently, gamma rays from radiation source 2 were ignored when the distribution reconstruction was performed via simple back-projection. For the first time, both radiation sources were measured at the 10th measurement point. However, the number of incident events exceeded the threshold at the 10th measurement point, and the second principal component was selected as the direction of the next measurement point. At the 11th measurement point, the impact of radiation source 1 was significantly greater than that of radiation source 2. Accordingly, a path was generated to solely surround radiation source 1.

In an environment where multiple radiation sources were dispersed, it was confirmed that path planning significantly depends on the setting of the threshold required to switch the principal component vector. This challenge can be addressed by devising a procedure for adopting the proposed method. In the initial stage, the location of the radiation sources is approximately estimated by setting the threshold to a small value. If a specific location of radiation sources is required, an additional investigation can be performed by increasing the threshold value.

To address this challenge, additional simulations were conducted under the same conditions. In the initial stage, the robot moved and measured gamma rays with the threshold set to 1000. The movement distance of the mobile robot at each time step between the measurement points was set to 1 m. Figure 15a presents the determined path and the reconstructed distribution map with the proposed method in the initial stage. The red dots denote the two radiation sources, and the black line represents the determined path with the proposed method. It was confirmed that a path primarily surrounding the two radiation sources was generated. By setting a small threshold value, measurement points away from the sources were determined. It was confirmed that two location around the sources showed high values; however, there was no clear distinction between these locations.

To perform a more specific investigation, the threshold value was changed from 1000 to 3000. The robot started from the initial stage end point and continued the investigation. The movement distance of the mobile robot at each time step between the measurement points remained unchanged. Figure 15b depicts the determined path and the reconstructed distribution map with the proposed method. The red line denotes the determined path in the second stage. It was confirmed that the robot passed between the two sources and circumnavigated one source at a short distance. Through this additional investigation, the distribution map was updated, and the area around the two sources was clearly divided. These results confirm that an additional investigation can be performed by increasing the threshold value if a more specific investigation is required.

3.1.4. Performance Comparison

To evaluate the effectiveness of the proposed method, a performance comparison was conducted with an information-driven search and the uniform search methods described in [8]. These methods are based on the use of the GM counter as a gamma detector. The simulation condition was determined as being the same as that in [8]. A single source was placed in a 100 m × 100 m rectangular search area. The source strength was 5.7 MBq. The distance between measurement points was assumed to be 15 m, and the threshold for switching the principal component vector was set to 30.

Figure 16 presents the simulation results with the proposed method. The red dot is the radiation source, and the black line represents the determined path. It was confirmed that the robot accurately estimated the direction of the existing source despite the large area. The robot approached the source and measured gamma rays while circumnavigating the source. Table 1 shows a performance comparison with previous methods.

Note that the localization errors of the previous methods are the average errors over 50 Monte Carlo runs. Conversely, the error of the proposed method is the result of a single trial. It was confirmed that the proposed method achieves better localization accuracy with fewer measurement points than the previous methods do. However, the proposed method counted at each location for a longer duration overall. Therefore, the need for fewer measurement points is advantageous if it is time-consuming for the robot to navigate between points.

3.1.5. Discussion

Herein, a strategy to approach and circumnavigate radiation sources was proposed. To evaluate its validity, the performance on two paths was compared: (a) the robot did not approach the source and followed a certain path, as in a uniform search; (b) the robot approached the source and measured gamma rays while following a path that surrounds the sources.

Figure 17 presents the results. The red dot is a radiation source, and the black line denotes the movement trajectory of the robot. The paths were manually determined. The 5.7 MBq source was placed, and the robot continued the measurement until the number of incident events exceeded 30 counts. The robot’s movement speed was assumed to be 1 m/s.

Table 2 lists the comparison results. It was confirmed that path (b) led to better accuracy in a shorter time than path (a) did. This result demonstrates that approaching the source and having multiple perspectives near the source are effective for radiation source inspection. The number of incident events is approximated to be the squared inverse proportional to the distance between the radiation source and the detector. Therefore, the path approaching the radiation source can reduce the measurement time until sufficient incident events are obtained. The measurement time at P1 along path (b) was 32 s while those at P2, P3, and P4 were 1 s.

Furthermore, the robot requires multiple perspectives to obtain the depth information because images acquired using a Compton camera are devoid of source-camera distance information. Having multiple perspectives near the radiation source improves the localization accuracy. Figure 18a presents a series of simulation results in Section 3.1.2. It was confirmed that measurement from one side was not sufficient to localize the sources accurately. However, a circumnavigational path is not absolutely necessary for radiation source localization, especially when a single source is placed in the environment. Figure 18b presents a series of simulation results in Section 3.1.1. It was confirmed that the source location was almost specified from one side. This aspect will be part of our future work defining the termination condition for the inspection.

3.2. Experiments

An actual experiment was conducted to evaluate the proposed method. Figure 19 illustrates the mobile robot and the all-around view Compton camera used in the experiment. Pioneer-3DX was adopted as the mobile robot [31], and an all-around view Compton camera was mounted above the mobile robot. The Compton camera has an energy resolution of approximately 60 keV; therefore, the margin to consider the energy resolution $\Delta E$ was set as 30 keV. The energy window to identify the radionuclide type was set as 630–690 keV for the source.

The position of the mobile robot was measured using motion capture systems. Reflection markers were installed on top of the all-around view Compton camera. In addition, a reflection marker was installed on the radiation source and used to obtain the ground truth of the radiation source position. Cesium-137 (2 MBq) was adopted as the radiation source.

The mobile robot was stationary for 10 min at each measurement point, and measurements were performed using an all-around view Compton camera. After measurement was completed at each point, the reconstruction of the radiation source distribution was performed via simple back-projection. Then, PCA was performed based on the reconstruction results to determine the next measurement point. When the number of incident events did not exceed the threshold, the next measurement point was determined according to the direction of the first principal component. However, the second principal component was selected to determine the next measurement point when the number of incident events exceeded the threshold. The moving distance of the mobile robot at each time step between the measurement points was set to 0.8 m. Figure 19 illustrates the experimental setup. In the right figure, the black arrow represents the direction that the mobile robot faces, and the red arrow denotes the movement direction to the next measurement point, as determined by the proposed method.

Figure 20 presents the results obtained from the environment in which a single radiation source was located. The red dot in the figure indicates the ground truth of the location of the radiation source. The threshold required to switch the principal component vector was set to 4000. The mobile robot moved from SP to the sixth measurement point in the direction approaching the radiation source, and then followed the path surrounding the radiation source. At the 9th measurement point, the number of incident events fell below the threshold, such that the path approaching the radiation source was selected again. After the 10th measurement point, a path was generated to surround the radiation source, and the measurement was completed at the 15th measurement point.

The final reconstruction of the radiation source distribution was performed via simple back-projection using only the data obtained from the case where the number of incident events exceeded the threshold. In Figure 20, the radiation source distribution is depicted in yellow. It can be confirmed that the radiation source distribution was successfully estimated around the ground truth. The location of the radiation source was estimated as ( 0.02 m, 0.07 m), which was compared to that of the ground truth, (0.23 m, 0.21 m). Hence, the estimation error was approximately 0.29 m. Based on this result, it was verified that the proposed method can be implemented to real systems and can successfully generate the measurement path for the localization of radiation sources. Moreover, the localization of the radiation source was accurately determined.

4. Conclusions

In this study, we proposed a path planning method for the localization of radiation sources using a mobile robot equipped with an all-around view Compton camera. The proposed method was designed to bring the radiation sources to a distance at which a sufficient number of incident events can be obtained, which then generates a path to surround the radiation sources. Via simulations and experiments, it was demonstrated that the proposed method can successfully generate a measurement path for the localization of radiation sources. In addition, the locations of the radiation sources can be estimated with high accuracy.

It was verified that path planning significantly depends on setting the threshold required to switch the principal component vector. To address this challenge, it is necessary to consider the gradient of the radiation intensity. Moreover, if a strong attenuator exists in the environment and the radiation sources are located over the object, the inspection difficulty increases. Additional path planning algorithms or inspection strategies should be utilized in such situations. We will focus on these points in our future work, as it can contribute to improving the reliability of the system.

Author Contributions

Conceptualization, T.K., H.W., R.K., Y.T., H.T., K.S., A.Y. and H.A; methodology, T.K., H.W. and R.K.; software, T.K., H.W., and R.K.; validation, T.K. and H.W.; formal analysis, T.K., H.W. and R.K.; investigation, T.K. and H.W.; resources, Y.T., H.T., K.S. and A.Y.; data curation, T.K. and H.W.; writing—original draft preparation, T.K. and H.W.; writing—review and editing, H.W. and R.K.; visualization, T.K. and H.W.; supervision, Y.T., H.T., K.S. and A.Y.; project administration, K.S. and H.A.; funding acquisition, K.S. and H.A. All authors have read and agreed to the published version of the manuscript.

Funding

This work was in part supported by JSPS KAKENHI Grant Number 19H00881 and the Nuclear Energy Science & Technology and Human Resource Development Project (through concentrating wisdom) from the Japan Atomic Energy Agency/Collaborative Laboratories for Advanced Decommissioning Science.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Conceptual image of measurement path: the proposed method is designed to approach the radiation sources to a distance where a sufficient number of gamma-ray events can be obtained, then, it generates a path to surround the radiation sources.

Figure 2. Schematic of proposed method: the proposed method consists of three parts: gamma-ray measurement, reconstruction of the radiation source distribution, and path planning. The robot repeats the gamma-ray measurement, reconstruction, path planning, and movement until the path is closed.

Figure 3. Principle of Compton camera: the Compton camera comprises two layers: scatterer and absorber. Let A and B be the positions at which the gamma ray is deposited on the scatterer and absorber, respectively, while θ is the scattering angle of the gamma ray relative to the straight line connecting A and B.

Figure 4. Obtaining Compton cones using an all-around view Compton camera: the all-around view Compton camera projects Compton cones on both sides of a straight line connecting the two interacting elements.

Figure 5. Simple back-projection method: let k be the vector connecting the interacting elements on the absorber and scatterer layers, and m the vector connecting the interacting element on the scatterer and the voxels in projection space. The pixel value of the voxels, depicted in orange, can be calculated using the delta function in Equation (2).

Figure 6. Reconstruction of radiation source distribution: the measurement data at each point is projected into the 3D voxel space via the simple back-projection method. Because the mobile robot can move in 2D space, the voxel data are projected onto a 2D plane.

Figure 7. Loop closure of measurement path: the measurement path generation is completed when the loop is closed.

Figure 8. All-around view Compton camera implemented in Geant4: the camera consists of a 12 × 4 × 30 array of CdTe detectors. The red, green, and blue arrows represent the X-, Y-, and Z- axes, respectively.

Figure 9. Process of simulation experiments: the measurement data of gamma rays were produced using Geant4. Then, the radiation source distribution was reconstructed via simple back-projection. Finally, the next measurement point was determined using the PCA based on previous measurements.

Figure 10. Simulation condition in which one radiation source exists: a single radiation source was positioned 15 m from the mobile robot. Cesium-137 (1 MBq) was adopted as the radiation source.

Figure 11. Simulation result of the condition that single radiation source exists: (a) SP represents the first measurement point at which the robot started to move, while EP is the position at which gamma-ray measurements were completed. The color bar on the right represents the value added to each pixel via the proposed method. The blue and red arrows represent the directions of the first and second principal components, respectively. The ground truth of the position of the radiation source is represented by a red dot. It was confirmed that the distribution of the radiation source was estimated around the ground truth. In addition, the measurement path was successfully generated as intended. (b) The orientation of the detector was zero when X axis was facing up, and clockwise rotation was defined as positive.

Figure 12. Simulation condition and result of multiple radiation sources concentrated in a specific area: (a) Three 1 MBq cesium-137 were positioned on a straight line, and the distance between the radiation sources was set to 1.5 m. (b) The blue and red arrows represent the directions of the first and second principal components, respectively. The ground truth of the position of the radiation source is represented by a red dot. The path approaching the radiation source was generated along the direction of the first principal component. Subsequently, a path was generated, which moved toward the second principal component and surrounded the radiation sources.

Figure 13. Simulation condition in which multiple radiation sources were dispersed: two radiation sources were positioned in two locations: (−0.5 m, −2 m) and (0.5 m, 2 m). Cesium-137 (1 MBq) was adopted as the radiation source. The robot was positioned 7m from the origin of the coordinates.

Figure 14. Simulation results of dispersed multiple radiation sources: (a) The results obtained when the threshold was set to 5000 (b) The results obtained when the threshold was set as 2500. It was confirmed that the measurement path significantly depends on the setting of the threshold required to switch the principal component vector.

Figure 15. Simulation results of additional investigation: (a) The result obtained when the threshold was set to 1000 (b) The result obtained when the threshold was set as 3000. The black line depicts the determined path when the threshold was set to 1000, and the red line is the path with the threshold set to 3000. It was confirmed that the area around the sources was clearly divided through the additional investigation.

Figure 16. Simulation results in large area: the red dot is the radiation source, and the black line depicts the determined path. It is confirmed that the robot accurately estimated the direction of existing the source despite of the large area. The robot approached to the source and measured gamma rays while circumnavigating the source.

Figure 17. Comparisons of two paths: (a) the robot did not approach the source and followed a certain path. (b) the robot approached the source and followed a path that surrounds the source. The red dot is a radiation source, and the black line depicts the moving trajectory of the robot. The paths were determined manually. The source with 5.7 MBq was placed, and the robot continued the measurement until the number of incident events exceeded 30 counts.

Figure 18. Series of simulation results: (a) the results in Section 3.1.2. (b) the results in Section 3.1.1. The red dot depicts a radiation source, and the black line is the determined path with the proposed method.

Figure 19. Experimental setup: pioneer-3DX was adopted as the mobile robot, and the all-around view Compton camera was mounted above the mobile robot. A single radiation source is expected in the environment. In the right figure, the red line represents the measurement points determined by the proposed method.

Figure 20. Experimental results: the blue and red arrows represent the directions of the first and second principal components, respectively. The ground truth of the location of the radiation source is represented by a red dot. It can be confirmed that the radiation source distribution was successfully estimated around the ground truth.

Table 1

Performance comparison.

	Information-Driven [8]	Uniform Deterministic [8]	Proposed
Number of measurement points	24	20	14
Localization error in X [m]	0.75	3.55	0.10
Localization error in Y [m]	0.72	2.18	0.69
Total measurement time [s]	89.9	90.1	242

Table 2

Performance comparison of two paths.

	Path (a)	Path (b)
Localization error in X [m]	3.10	0.33
Localization error in Y [m]	0.27	0.09
Total measurement time [s]	80	35
Total moving distance [m]	80	83.4
Total search time [s]	160	118.4