Multi-robot cooperative search for radioactive sources based on particle swarm optimization particle filter

2.1. Problem description

Assume that the missing radioactive source $$ \theta=\left[x_{s}, y_{s}, I_{s}\right]^{T} $$ is located on the ground level in a mountainous area, where $$ \left(x_{s}, y_{s}\right) $$ is the position coordinate of the radioactive source, and $$ I_{s} $$ represents the activity of the source. In this study, a team of mobile robots was deployed to search for a lost radioactive source $$ \theta $$ in a complex mountainous environment. Each robot has the ability to move autonomously and perform local search operations while seamlessly cooperating with other robots to achieve a comprehensive and systematic global search.

Nuclear radiation decay is a random process that occurs in radionuclides. This randomness makes the decaying particles uncertain. According to previous studies, radiation counts from nuclear decay obey Poisson statistics ^[16]. Therefore, the probability that a radiation detector registers $$ z \in Z^{+} $$counts per second (CPS) from the source that emits on average $$ m $$ CPS is:

where $$ \lambda=m \varepsilon $$ denotes the actual average count of particles detected by the detector per second and the parameter of the Poisson distribution; $$ \varepsilon $$ is the intrinsic detection efficiency of the detector, which is the ratio of the number of particles detected by the detector to the total number of particles incident on the detector at the same time. $$ \lambda=m \varepsilon $$ represents the true average count measured by the detector per second; $$ \varepsilon $$ is the inherent detection efficiency of the detector, which is the ratio of the number of particles detected by the detector to the total number of particles incident on the detector at the same time. Given the source parameter vector $$ {\theta_k} $$, the likelihood function of the radiation count measurement $$ z_{k} $$ at the recorded measurement location $$ \left(x_{k}, y_{k}\right), k=1, 2, \ldots $$ is given by:

where $$ \mathcal{P} $$ is the Poisson distribution defined in Equation (1), and $$ \lambda( {\theta_k}) $$ is the mean radiation count recorded by the detector. When the distance between the radiation detector and the radiation source is $$ r $$, the calculation of $$ \lambda( {\theta_k}) $$ is performed as follows:

where $$ R_{i} $$ and $$ u_{i} $$ are the distance and radiation attenuation coefficient of the gamma rays through the ith medium, respectively; $$ I_{b} $$ is the background radiation count, which is the average count per second when $$ I_{s}=0 $$ or when $$ r $$ is infinite. The radiation detector measures the intensity of the radiation source in CPS and calculates the count value. Simultaneously, the energy response constant $$ \rho\left(C P S \cdot \mu S v^{-1} \cdot h\right) $$ is considered ^[17]. The count value $$ z $$ is expressed as the product of the mean radiation count $$ \lambda( {\theta_k}) $$ and the energy response constant $$ z=\lambda( {\theta_k}) \times \rho $$.

The activity of a radioactive source can be greatly affected by factors such as detector performance limitations, environmental absorption, multipath attenuation, and reflections. These factors often bring errors to detection data, thus affecting the accuracy of radioactive source positioning and parameter estimation. Figure 2 depicts the radiation profile attenuated by simulated radiation data according to Equation(3). It can be clearly observed from Figure 2 that due to the influence of the object's radiation attenuation coefficient, the radiation of the radioactive source model presents a complex distribution.

Figure 2. Heatmap of the radiation distribution in the simulated environment, where the color bars indicate the signal strength in CPS. In the map, the source is at the upper right with an activity of about $$ 2.8 \times 10^5 $$ CPS. CPS: Counts per second.

2.2. Particle filter

Within the framework of the PF algorithm, we use an ensemble of particles and their associated particle weights to estimate the posterior probability density of the system state. Each particle represents a potential system state, and the accuracy of state estimation increases with the number of particles. At the same time, each particle weight represents the likelihood probability associated with the system state of each particle. Compared with traditional Kalman filtering methods, particle filtering has proven to be more suitable for dealing with nonlinear and non-Gaussian systems ^[18].

Suppose we have $$ N $$ mobile robots deployed in a radiation field, each maintaining $$ M $$ random particles. The estimated source parameter vector of the $$ m $$ th particle after obtaining the $$ k $$ th observation is expressed as $$ \theta_{k}^{m}=\left[x_{k}^{m}, y_{k}^{m}, I_{k}^{m}\right]^{T}, m \in M $$. At the $$ k $$ th moment, the $$ i $$ th robot is located at the location $$ \left(x_{k}^{i}, y_{k}^{i}\right) $$, and the obtained observation vector is expressed as $$ g_{k}^{i}= $$$$ \left\lceil x_{k}^{i}, y_{k}^{i}, z_{k}^{i}\right\rceil^{T}, i \in N $$. If the initial confidence is given as $$ p\left(\theta_{0} \mid g_{0}\right)=p\left(\theta_{0}\right) $$, then according to the Markov hypothesis, the prior probability density can be expressed as:

where $$ P\left(\theta_{k} \mid \theta_{k-1}\right) $$ represents the state transition probability, $$ p\left(\theta_{k-1} \mid g_{1: k-1}\right) $$ is the posterior probability density corresponding to the previous step. After obtaining the observation vector $$ g_{k} $$, the prior value is updated through Bayesian theory, and the posterior probability density is obtained as:

where $$ p\left(g_{k} \mid g_{1: k-1}\right)=\int P\left(g_{k} \mid \theta_{k}\right) P\left(\theta_{k} \mid g_{1: k-1}\right) d \theta_{k} $$ is a normalization constant, which depends on the likelihood function $$ p\left(g_{k} \mid \theta_{k}\right) $$ and the statistical properties of the measurement noise. In the PF algorithm, we use the Monte Carlo (MC) method to approximate the posterior probability density of the target system state by the weighted sum of a set of particles ^[9].

where $$ M $$ is the number of particles, $$ w_{k}^{m} $$ is the weight of the $$ m $$ th particle, and $$ \delta(\cdot) $$ is the Dirac function. Furthermore, the importance sampling (IS) technique ^[19] and sequential IS (SIS) ^[18] technique can be used to approximate the weights of these random particles.

where $$ p(\cdot) $$ is the target distribution, representing the posterior probability density we hope to obtain. In the context of radiation source localization, the fusion of measurements from multiple sensors is crucial for achieving accurate estimations of source parameters. Consider a scenario where an array of sensors is deployed to monitor radiation levels using a total number of sensors. To effectively use these measurements, we formulate the estimation of the radiation source parameters as a probabilistic exhibit minimal cross-correlation in their measurement errors. The unnormalized weight for the $$ i $$ th potential parameter value at time $$ k $$. is expressed as:

where $$ {w^{\prime}}_{(k)}^{(i)} $$ denotes the unnormalized weight associated with the approximated posterior probability density function. Consequently, the normalized weight for the $$ i $$ th potential parameter value is given by:

The PF is a non-parametric filtering algorithm extensively employed for state estimation tasks. Nevertheless, a major drawback of the PF is the issue of particle degeneracy (a few particles end up having non-zero weights). This degeneracy occurs when the likelihood function becomes highly localized, indicating precise observation information. Consequently, the overlap between the likelihood probability distribution and the prior probability distribution diminishes significantly, resulting in an inadequate number of particles to represent the posterior distribution effectively. Under such circumstances, the estimation outcome of the PF may exhibit considerable errors or even yield incorrect estimations. To address the issue of particle poverty, it becomes imperative to employ a larger number of particles to approximate the posterior distribution more accurately. However, this approach comes at the cost of increased computational complexity for the algorithm. Moreover, the effectiveness of the PF is contingent upon the accuracy and adaptability of the proposed distribution. When a substantial disparity exists between the proposed distribution and the target distribution, the performance of the PF is adversely affected. For instance, if the target distribution exhibits multiple peaks and the proposed distribution fails to accurately capture these peaks, the PF performance will be compromised. Furthermore, the PF also has the problem of high computational complexity. In practical applications, a large number of particles are often required to approximate the posterior distribution adequately. However, as the number of particles increases, so does the computational burden of the algorithm. This limitation restricts the real-time applicability of PFs in certain scenarios.

3.1. Improved particle swarm optimization algorithm based on multi-robot

PSO is an advanced global optimization algorithm that draws inspiration from the collective behavior of organisms such as birds or fish. It was originally introduced by Kennedy and Eberhart in 1995 ^[20]. This algorithm efficiently explores multi-dimensional search spaces by simulating the movement and interaction of a swarm of particles. In PSO, every potential solution within the search space assumes the form of a particle, possessing both a position and a velocity. As the particles diffuse, each particle adjusts its position and velocity based on its current state while retaining knowledge of its own historical best position and the overall best position found so far. Through continuous iterative refinement, the ultimate global optimal solution gradually reveals itself. Where the velocity and position of each particle are updated using the following formula:

where $$ x_{t}^{i} $$ signifies the position of the ith particle at time t; $$ v_{t}^{i} $$ symbolizes the velocity of the ith particle at time t; $$ \text{pbest}^{i} $$ reflects the best historical position reached by the ith particle, and gbest represents the best global position discovered so far. $$ w $$ is the inertia weight, $$ c_{1} $$ and $$ c_{2} $$ are learning factors, and $$ r_{1} $$ and $$ r_{2} $$ are random numbers between 0 and 1.

The algorithm presented in this paper integrates the PSO algorithm into the PF algorithm, leveraging it as the generator of the proposed distribution. The algorithm uses an adaptive weighted PSO algorithm, which dynamically adjusts inertia weights and acceleration coefficients to strike a balance between global and local search capabilities. The adaptive weight of the algorithm can be calculated as follows:

where $$ w_{\max } $$ is the maximum value of the inertia weight, $$ w_{\min } $$ is the minimum value, and $$ T $$ denotes the maximum number of iterations. Furthermore, an adaptive fitness function is introduced to fuse the latest observations and recommendation distributions from multiple robots, defined as follows:

where $$ y_{\text {new }}^{i} $$ and $$ y_{\text {pre }}^{i} $$ denote the most recent observed and predicted values of the j-th robot, respectively; $$ \varGamma $$ represents the absolute deviation value between the observed value of the robot and the predicted value of the particle set at time $$ k $$. This characteristic is evident in the likelihood function $$ p\left(z_{k} \mid {\theta_k}\right) $$ for the radiation measurement value $$ y_{\text {new }}^{i} $$, wherein the weight increases as the measurement value approaches the peak and decreases as it deviates away from the peak. This paper presents an improved PSO algorithm based on multi-robot systems as a comparison to the conventional PSO method ^[21–23]. The proposed algorithm incorporates the most recent observation and prediction outcomes from every robot, allowing the particles of all robots to relocate towards regions with higher posterior probability density distributions. This adjustment aims to overcome the issue of particle degradation, which occurs when the weight distribution of particles becomes highly imbalanced. Algorithm 1 outlines the process of particle optimization using the multi-robot improved PSO algorithm.

3.2. Dynamic Window Approach dynamic window avoidance strategy

The Dynamic Window Approach (DWA) algorithm is a widely used method for obstacle avoidance ^[24]. Assuming the current state of the robot is denoted as $$ x_{t}=[x, y, \theta, v, \omega]^{T} $$, where $$ \mathrm{x} $$ and $$ \mathrm{y} $$ represent the coordinates of the robot, $$ \theta $$ represents the orientation, and $$ \mathrm{v} $$ and $$ \omega $$ represent the linear and angular velocities, respectively. The DWA algorithm generates a set of speed commands, which can be expressed as follows:

where $$ v_{\min } $$ and $$ v_{\max } $$ represent the lower limit and upper limit of the robot's linear speed, respectively; $$ w_{\min } $$ and $$ w_{\max } $$ present the lower limit and upper limit of the robot's angular speed, respectively; $$ v_{a} $$ and $$ w_{a} $$ represent the current linear speed and angular velocity of the robot, respectively; $$ \operatorname{dist}(v, w) $$ represents the shortest distance between the simulated trajectory and the obstacle at the current speed; $$ \Delta t $$ represents the time interval required for the robot to execute the speed command.

The DWA algorithm generates multiple sets of feasible speed command combinations, represented as $$ u=[v, w]^{T} $$. These combinations are determined by the intersection of velocity boundary limits $$ v_{m} $$, acceleration limits $$ v_{d} $$, and environmental obstacle limits $$ v_{a} $$. Figure 3 visually illustrates the predicted trajectories derived from these combinations of velocity commands using the kinematic model of the mobile robot. In addition, weights are assigned to these trajectories through the evaluation function $$ g(u) $$ to select the driving speed corresponding to the optimal trajectory.

Figure 3. Mountain multi-robot trajectory prediction model. The mobile robot's historical search trajectories are depicted in different colors: blue, orange, yellow, and purple dots. The solid green line represents the predicted trajectory, while the green triangles along the trajectory represent the expected robot position.

where $$ \operatorname{dist}(u) $$, heading $$ (u) $$, and speed $$ (u) $$ represent the distance, angle, and speed, respectively, between the robot's next position $$ x_{k+1} $$ and the target position under decision $$ u $$. The weight coefficients, $$ w_{1}, w_{2} $$, and $$ w_{3} $$, determine each factor's relative importance in the speed command evaluation function. When a trajectory intersects with an obstacle at any point in the future, that trajectory is disregarded. Finally, the optimal decision $$ u^{*} $$ is decided to minimize $$ g(u) $$ as:

where the feasible decision set $$ U $$ refers to the set of permissible speed instructions, as defined by Equation (13). Algorithm 2 outlines the process of computing the reward function for the nth robot at time $$ k $$.

4.1. An illustrative run

The simulation results of multi-robot collaborative radioactive source search using PSO-PF algorithms are shown in Figure 4. The illustrative run uses the following parameters for the simulation:

Figure 4. An illustrative run of the source searching process using a multi-robot system. (A-C) Trajectories of the four robots and the estimated source locations at steps $$ k=1 $$, $$ k=5 $$, and $$ k=18 $$, respectively; (D) the PSO iteration and global fitness; (E) the observables of agents during the search. PSO: Particle swarm optimization.

● Multiple robots search for lost radioactive sources in a $$ 150 \mathrm{\; m} \times 150 \mathrm{\; m} $$ mountainous environment.

● The radioactive source is located at coordinates $$ (x_{s}=135 \mathrm{\; m}, y_{s}=90 \mathrm{\; m}) $$, and its activity is set to $$ I_{s}=2.94 \times 10^{6} \mathrm{\; CPS} $$.

● The intrinsic detection efficiency of the radiation detector $$ \varepsilon=60\% $$, the mountain radiation attenuation coefficient $$ u=0.05 $$, the energy response constant $$ \rho = 200 \mathrm{\; CPS/\mu Sv/h} $$, and the background radiation count is $$ I_{b}=20 \mathrm{\; s}^{-1} $$.

● The robot's movement is constrained by various physical parameters, including maximum linear velocity $$ v_{\max }=3 \mathrm{\; m/s} $$, maximum angular velocity $$ w_{\max }=90^{\circ}/\mathrm{s} $$, linear acceleration $$ \dot{v}=1 \mathrm{\; m/s}^{2} $$, angular acceleration $$ \dot{w}=30^{\circ}/\mathrm{s}^{2} $$, linear velocity resolution $$ 0.5 \mathrm{\; m/s} $$, and angular velocity resolution $$ 15^{\circ}/\mathrm{s} $$.

Figure 4A-C visually shows the source search area, the search trajectories of the robot at times $$ k $$ = 1, 5 and 18 using this algorithm, and the contour plots of the radiation distribution from the source. It is worth noting that the robot's trajectory, position, observation values, and particles are represented by different colors. Furthermore, the observations are described in the form of a staircase plot, and the global fitness and number of iterations are represented as a stacked plot of two variables with a common $$ \mathrm{x} $$-axis. During the search process, the fitness function of the PSO algorithm is calculated based on the latest observation values obtained by each robot, which makes the particle set move to a high-likelihood area and overcomes particle degradation. At step $$ k $$ = 1, the measurements obtained by the detector are small, and the positions of the particles are evenly distributed on the map. As the particle swarm algorithm iterates to find the optimal particle set, the particles gather around the estimated source location at step $$ k $$ = 5, indicating that the distance between the estimated source location and the true source location is decreasing. Finally, all robots complete the source search at step $$ k $$ = 18.

Figure 4D and E shows that as the number of search steps increases, the robot's observations become increasingly larger, and the number of iterations of the PSO algorithm gradually decreases. It should be emphasized that when the robots are far away from the source location, the observations obtained by the robots from each other tend to be relatively small. Therefore, as the difference between the observations obtained by the robot and its predicted observations becomes smaller, the global optimal particle fitness value also becomes smaller. On the contrary, when the robot approaches the source position and obtains relatively large observed values, the global optimal particle obtains higher fitness values that increase as the gap between the observed and predicted values increases.

4.2. Parametric analysis

In this experiment, we aim to evaluate the performance comparison of the proposed PSO-PF algorithm and the Extended KalmanParticle Filter algorithm (EPF) in the multi-robot radioactive source search problem. The termination conditions for the source search process are defined as follows: The process concludes when the Euclidean distance between all robots and the estimated source position is less than 5 m and the variance of the samples in the PF is less than 1. Once the search process ends, it is deemed successful if the estimated source locations have converged to the true source location and the distances between all robots and the true source location are less than 5 m. Otherwise, it is considered a failure. Additionally, if the search steps reach 30 and the source has not been found, the source search process is terminated as well, indicating an unsuccessful search task. Figure 5 shows the results of single-step search time (SSST) and relative estimation errors (REE) for both algorithms performing a search task using the same simulation parameters as in the example run. It is worth noting that SSST represents the time required to estimate the source position and perform a move at each step, while a REE represents the ratio of the absolute error to actual value.

Figure 5. PSO-PF vs. PF for multi-robot radiation source search: (A) Relative Estimation Error of PSO-PF; (B) Relative Estimation Error of PF; (C) Single-step search time comparison. EPF: Extended Kalman Particle Filter; PSO-PF: Particle Swarm Optimization-Particle Filter.

Figure 5 presents a clear demonstration of the superior performance of the PSO-PF algorithm in comparison to the EPF algorithm, specifically in terms of the REE and SSST performance indicators. It is important to highlight that the PSO-PF algorithm exhibits a REE below 0.1 for the x-axis coordinate, y-axis coordinate, and activity degree (Ⅰ) of the source term. Conversely, while the EPF algorithm achieves a REE for the source location coordinates below 0.1 at $$ k $$ = 18, the REE for the source activity remains consistently around 0.2. While the PSO-PF algorithm excels in source parameter estimation, it initially incurs additional overhead in terms of single-step search time.

In order to understand the impact of particle numbers on the PSO-PF algorithm, 100 search tasks were also performed using the same parameters as in the example run. Figure 6 shows the search success rate and average search time when varying the number of particles. As the number of particles increases, the search success rate of the free energy strategy increases significantly, climbing from 56% to 98%. At the same time, the average search time increased from 2.75 to 7.04 s. In contrast, the DWA dynamic window obstacle avoidance strategy always maintains a high search success rate as the number of particles increases without significantly affecting the search time. The comparison between the DWA dynamic window obstacle avoidance strategy and the free energy strategy in Figure 6 shows that using the DWA dynamic window obstacle avoidance strategy in the proposed algorithm improves the efficiency and accuracy of multi-robot collaborative search work.

Figure 6. Performance comparison of the PSO-PF algorithm as the number of particles changes (A) The mean search time; (B) The search success rate. PSO-PF: Particle Swarm Optimization-Particle Filter.

Acknowledgments

The authors would like to thank the editor-in-chief, the associate editor, and the anonymous reviewers for their valuable comments.

Authors' contributions

Made substantial contributions to the research methodology, conceptual, and simulation and wrote and edited the original draft: Luo M, Huo J, Liu M

Performed critical review, comments, and revisions and provided administrative, technical, and material support: Zhou Z, Liu M

Availability of data and materials

Not applicable.

Financial support and sponsorship

This work was supported by the National Natural Science Foundation of China (No. 12205245) and the Natural Science Foundation of Sichuan Province (No. 2023NSFSC1437).

Conflicts of interest

All authors declared that there are no conflicts of interest.

Ethical approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Copyright

© The Author(s) 2023.