573
Views
0
CrossRef citations to date
0
Altmetric
Technical Papers

Tsallis Divergence as Strategy for Radioactive Source Search and Location

, , , &
Pages 1114-1124 | Received 30 May 2021, Accepted 09 Mar 2022, Published online: 20 Apr 2022

Abstract

For the problem of searching radioactive sources in a certain area, a search method combining Tsallis divergence strategy with a particle filtering algorithm is proposed. The method this paper proposes for searching for radioactive sources is carried out by using a mobile platform equipped with a NaI(Tl) scintillator detector. The estimation model of the parameters of the radioactive source is constructed and is based on the inverse square law and the fact that the count values of the NaI(Tl) detector in nuclear decay processes obeys the Poisson distribution. The Tsallis divergence strategy is used as a reward function to control the movement of the platform. The posterior distribution of the parameters of the radioactive source is continuously and iteratively updated by using the particle filtering algorithm. The results of Monte Carlo simulations and practical experiments verify the effectiveness of the algorithm.

I. INTRODUCTION

Since the twentieth century, the nuclear industry and nuclear technology have grown rapidly over the world. Nuclear materials and radioactive sources have gradually been used in all aspects of production. A trend of cross-fertilization between the nuclear discipline and other disciplines has begun to emerge. While nuclear technology facilitates people’s lives, there is also a risk of losing control of radioactive sources, which could lead to environmental pollution and public health hazards. Therefore, it is necessary to conduct in-depth research on radioactive source search and location techniques.

The search of uncontrolled radioactive sources is an important prerequisite for emergency response. So, a great deal of research has been carried out on the problem of detection and location of radioactive sources. The corresponding methods are divided into three categories: (1) orientation-based methods for locating radioactive sources that mostly use collimatorsCitation1–3 or shieldingCitation4–6; (2) a method based on the variation of radiation field counts (rates) with detection position that can be used for estimation of the parameters of radioactive sources by using Bayesian estimation,Citation7–10 maximum likelihood estimation,Citation11,Citation12 and other mathematical statistical methods; (3) gamma cameras.Citation13–16 Orientation-based methods are generally suitable for scenes with a small area and clear radiation distribution. With gamma cameras, long imaging times are required for radioactive source search. The method based on the change of radiation field count (rate) with the detection position has the advantages of a large search region and low detector requirements and is the most practical method for searching and locating radioactive sources.

For the search problem, this paper proposes an estimation method based on Tsallis divergence. Through the recursive Bayesian estimation model, the particle filtering algorithm is used to estimate the parameters of the radioactive source. The Tsallis divergence is used as a search strategy to control the movement of the mobile platform to achieve the purpose of searching the radioactive source. Compared with other radioactive source search methods, the method proposed in this paper has the advantage of being able to search for an unknown number of sources in a certain area instead of searching for one radioactive source in the preset area.Citation10,Citation11

The rest of the paper is organized as follows. Section II introduces the mathematical formulation of the problem. The principle of radioactive source search, including the particle filtering algorithm and control strategy of the mobile platform based on Tsallis divergence, are presented in Sec. III. The Monte Carlo simulations and practical experiments are carried out in Secs. IV and V, followed by conclusions in Sec. VI.

II. BACKGROUND AND THEORY

II.A. Radioactive Source Locating Model

Radiation counting recorded by a detector is a random process governed by a Poisson distribution.Citation17 The probability that the detector measures in a certain time τ at the observation position and obtains a count z is

(1) P(z;λ)=λzz!eλ,(1)

where λ is a function of the source parameters.

Let A be a flat open-field rectangular-shaped search region. Assume that a radioactive source is present in A. Let θi=xi,yi,εiR2 × R+, and denote the parameter vector of the source, where (xi, yi) denotes the source position in Cartesian coordinates and εi denotes the source intensity, where i indicates the i’th source. The vector of the radioactive source parameters can be expressed as θ=θ1,,θs, where s indicates the number of radioactive sources. Then, the function of the radioactive source parameters λ(θ) can be given by EquationEq. (2)..

(2) λ(θ)=ub+i=1sεiSηeμidi4πri2,(2)

where

ub=

= average counting rate due to background radiation (from cosmic rays and environmental radiation background in earth); ub is assumed to be constant over the search region

ri=

= distance between the detector and the i’th source

S=

= effective area of detector

η=

= detector (intrinsic) efficiency

di=

= thickness of the shield

μi=

= linear attenuation coefficient, which depends on the energy of gamma radiation between the detector and the source.

III. SEARCH PRINCIPLE

III.A. Recursive Bayesian Estimation

According to Bayes’ theorem, the recursive Bayesian estimation method can be used to update the posterior probability density function (PDF) of the system state using new data obtained from detection. The posterior PDF is then used to find the optimal estimation of the state parameters:

(3) pa|b=pb|apapb,(3)

where p(b) is the normalization constant. It can be given by EquationEq. (4):

(4) pb=pb|apadx.(4)

Suppose that the posterior PDF of the parameter of the radioactive source at time k  1, i.e., Pθ|z1:k1, is known, where z1:k1=[z1,,zk1] is the set of detector measurements. After receiving a measurement zk at time k, the posterior PDF of the parameter of the radioactive source can be updated as

(5) Pθ|z1:k=Pzk|θPθ|z1:k1Pzk|z1:k1,(5)

where Pzk|z1:k1 is the normalization constant. The parameter estimation of the radioactive source can be calculated after obtaining the updated PDF. The function Pzk|θ is the likelihood function defined by Eq. (1):

(6) pzk|θ=Pzk;λθ.(6)

III.B. Particle Filtering

For nonlinear recursive Bayesian estimation problems such as the estimation of the parameters of the radioactive source, particle filtering is a good method.Citation18 The core idea of particle filtering is to use a set of discrete particles to approximate the posterior distribution of the system state. It means that a set of particles with weights θk,i,wk,ii=1N are used to approximate the posterior distribution of the radioactive source parameters.

III.B.1. Particle Initialization

In the particle initialization phase, a particle set θ0,i,w0,ii=1N is created based on the priori information, where θ0,i=xs0,i,ys0,i,εs0,i. In this paper, it is assumed that the radioactive sources are uniformly distributed in the search region, i.e., xs0,i,ys0,i:U(A), where subscript s denotes the parameter of the s’th radioactive source and A denotes the area of the search region. The activity of the radioactive source obeys the gamma distribution, i.e., εs0,i:Γ(κ,φ), where κ and φ are artificially set parameters of the gamma distribution. The probability of radioactive sources in the search region is given by EquationEq. (7):

(7) P0s=smaxs+1s=0smaxs+1,(7)

where the maximum number of sources smax = 2. The initial particle weights are all equal, i.e., w0,i=1N, i=1,,N.

III.B.2. Update Step

After the zk is obtained, the posterior PDF of the radioactive source parameters can be updated according to Eq. (5). At time k  1, after completing resampling, the weights of each particle are set as 1/N. Therefore, the weight of each particle at time k can be updated as

(8) wk,i=Cpzk|θk,i,(8)

where C is the normalization factor.

If sampling operation was performed directly, particle degradation would inevitably occur.Citation10 The problem of particle degradation is manifested by the fact that few particles have large weights while most particles have almost zero weights. Computational resources would be wasted on useless particle updates.

To overcome the particle degradation problem, the technique of progressive correction is used in this paper.Citation19 Sampling operation is divided into different stages. In the early stages of the sampling operation, the weights of the particles are larger than the final particles. This operation allows more particles to be selected. Then, the weights gradually approach the real weights and complete the sampling operation. The detailed procedure can be found in CitationRef. 19. In the calculation phase t=1,,T, the weights of the particles can be expressed as

(9) wk,it=Cpzk|θk,it=1Tγt,(9)

where 0γt<1 and t=1Tγt=1.

III.B.3. Estimation of Radioactive Source Parameters

After the updated step, a new set of particles can be obtained. The number of radioactive sources in the region can be estimated by using the set of particles:

(10) sˆ=1Ni=1Nsk,i,(10)

where sk,i denotes the number of radioactive sources contained in the i’th particle at time k and sˆ denotes the estimation of the number of radioactive sources. Then, using EquationEq. (11), the radioactive source parameters can be estimated as the sample mean:

(11) θˆ=1Nsˆi=1Nθk,iδksk,i,sˆ,(11)

where Nsˆ indicates the number of particles and δk is the Kroneker delta.

III.C. Platform Control

In this section, the control problem of the mobile platform is considered as a partially observable Markov decision process (POMDP). The elements of the POMDP in the estimation problem contain the posterior PDF of the set of particles, a set of control vectors, and a reward function. The control vectors denote the set of actions that the mobile platform can take in the current state. The reward function represents the value of the reward that the mobile platform can obtain by reaching the next state. Let ukUk denote the control vector applied at time k, where Uk is the set of control vectors at time k and Uk=,,,,Z,[,],\^. It can be supposed that the posterior PDF of the set of particles at time k − 1 is pθ,s|z1:k,u0:k1. Thus, after obtaining zk, the control vector applied at time k can be expressed as

(12) uk=argmaxvUkERv,pθ,s|z1:k,u0:k1,(12)

where R(·) represents the reward function. The value of the reward function depends on the measurements obtained after applying the control vector.

Different strategies can be implemented by choosing different reward functions. For this kind of problem, researchers have considered methods such as entropyCitation20 and fisher information.Citation8 In this paper, Tsallis divergence is chosen as the reward function. This method would help us to control the platform with the maximum information gain through the entire search.

The Tsallis divergenceCitation21,Citation22 between two functions f0 and f1 can be defined as

(13) DTαf0f1=11α1f0α(x)f11α(x)dx,(13)

where the integral is of the dimension of the state space and α > 0 is the parameter of the function.

The reward function used in EquationEq. (12) is taken to be the Tsallis divergence between the posterior density at time k, f0=ps,θz1:k,u0:k1, and the future posterior density, f1=ps,θz1:k,zk,u0:k, where zk is the measurement detected by the detector after the application of the control vector. To reduce the complexity of the notation, in the following the past measurement z1:k and control u0:k1 will be omitted. Thus, ps,θ denotes the posterior density of the particle set at time k, and ps,θzk,uk denotes the posterior density of the particle set after the application of the control vector. The reward functionCitation23 in EquationEq. (12) can be expressed as

(14) Ruk,pθ,s=11α1s=0smaxps,θzk,ukαps,θ1αdx.(14)

Substituting the posterior distribution of particles (after resampling with equal weights), the Tsallis divergence between the posterior density at time k and the future posterior density can be approximated by EquationEq. (15) (CitationRef. 23):

(15) Ruk,pθ,s11α1γαzkukγ1zkukα,(15)

where for zN,

(16) γτzuk=1/Ni=1Npzθk,i,sτ.(16)

The likelihood function pzθk,i,s can be calculated by EquationEq. (6). Thus, the expectation value of the reward function in EquationEq. (12) can be approximated as

(17) ERv,pθ,s|z1:k,u0:k111αz=0z1γ1zuk1γαzukγ1zukα,(17)

where z1 is the settable value, which is set to 25 in this paper.

When using the Tsallis divergence method in practice, there may be scenarios where the detection position is too close to the radioactive source. The parameter estimation model is not valid for this scenario. In order to prevent the platform approaching too closely to the radioactive sources, this paper adds a multiplication factor to the formula. When the radiation field is strong, the multiplication factor is small and vice versa:

(18) ERv,pθ,s|z1:k,u0:k1e(ψλ(θ))1αz=0z1γ1zuk1γαzukγ1zukα,(18)

where λθ=1/Ni=1Nλθk,i is the sample mean of the mean radiation count defined in EquationEq. (2); ψ0 is a user-defined parameter, and in this paper ψ is set as 0.1.

IV. MONTE CARLO SIMULATIONS

IV.A. Environmental Settings

In this paper, a 20 × 10-m rectangular search region is constructed, and measurement points are set at 2-m intervals within the search region. The platform with a NaI(Tl) scintillator detector moves over the measurement points in the search region. The height of the detector is set at 0.5 m above the ground. The platform can only move on eight measurement points adjacent to itself at a time (excluding the case where the measured point position exceeds the search region boundary). The maximum number of radioactive sources present in the search region is two. In this paper, we have constructed a simulation scenario shown in using Geant4. In the Geant4, scenarios were simulated where one or two radioactive sources are present in the search region. The NaI(Tl) scintillator detector was selected to record the data in the simulation. And, then, the simulated measurement of each detection point could be obtained. After obtaining the detection data in the simulation, the proposed algorithm was used to estimate the parameters of radioactive sources.

Fig. 1. Search region and measurement point schematic.

Fig. 1. Search region and measurement point schematic.

In , the area within the black rectangle indicates the region where the radioactive source may be present, and the circles represent the measurement points. In this paper, it is assumed that the prior distribution of the position of the sources obeys a uniform distribution and that the prior distribution of the relative source intensity of the sources obeys the gamma distribution, e.g., x0,i:U20,20, y0,i:U10,10, and ε0,iΓ(1.5,140). And, the activity of the radioactive source is set to 2 × 107 Bq. The energy of the sources is set as 1.33 MeV, and the detector (intrinsic) efficiency is set as 0.3. The mean counting rate of the background radiation is μb = 0.18/s. The effective detection area of the detector S is a constant. In this study we restrict the search region to a relatively small search area, so the influence of air attenuation can be ignored. Therefore, EquationEq. (2) can be simplified as

(19) λ(θ)=0.18+i=1smaxνiri2,(19)

where νi is the relative intensity of the i’th radioactive source. In this paper, relative intensity is set to the detector measurement 1 m away from the radioactive source; ri is the distance from the i’th radioactive source to the detector.

IV.B. Experiments with One Source

First, we study the performance of the algorithm for different numbers of moves. One radioactive source is located at (665, −171) cm in the search region. The α value is set to 0.5, and the number of particles N is set to 5000; the rest of the settings are the same as in Sec. IV.A. The experiment is carried out by changing the number of moves (one measurement per move of the mobile platform). The results are shown in and . Not surprisingly, the error in the estimation decreases when the number of detections increases. However, when the number of detections exceeds 16, increasing the number of detections has little effect on the estimation results but increases the calculation time and reduces the calculation efficiency.

Fig. 2. Search path and radioactive source parameter estimation results: (a) moves are 5; (b) moves are 10; (c) moves are 15; (d) moves are 20.

Fig. 2. Search path and radioactive source parameter estimation results: (a) moves are 5; (b) moves are 10; (c) moves are 15; (d) moves are 20.

Fig. 3. Error of estimation of radioactive source parameters.

Fig. 3. Error of estimation of radioactive source parameters.

Next, using the same simulation setup as before, we study the performance of the search algorithm for different search starting points. The results are presented in .

Fig. 4. Experimental results: (a) the starting point is (2000, 1000) cm; (b) the starting point is (2000, −1000) cm.

Fig. 4. Experimental results: (a) the starting point is (2000, 1000) cm; (b) the starting point is (2000, −1000) cm.

It can be concluded from the experimental results that the choice of different search starting positions has little effect on the results. The control strategy only considers differences in the posterior density of present and future radioactive source parameters.

IV.C. Experiments with Two Sources

In this section, we study the performance of the search algorithm in the case of two radioactive sources. Two radioactive sources are located at (629, 345) cm, (−453, −310) cm in the search region. The α value is set to 0.5, the number of moves is set as 30, and the number of particles N is set to 5000; the rest of the settings are the same as in Sec. IV.A. The results are shown in .

Fig. 5. Results of parameter estimation for two radioactive sources: (a) search path; (b) estimated number of radioactive sources.

Fig. 5. Results of parameter estimation for two radioactive sources: (a) search path; (b) estimated number of radioactive sources.

As can be seen from , in contrast to the case of one source, the Tsallis divergence–based method is able to essentially complete the task of estimating the parameters of the two radioactive sources. The dimensionality of the algorithm’s search space has increased due to the addition of a radioactive source. As a result, the accuracy of the algorithm for source estimation is reduced.

Then, we study the performance of the search algorithm for different values of parameter α. The simulation setting is identical to those in Sec. IV.A, except that α takes the values 0.1, 0.5, and 0.999. The results are shown in .

TABLE I Search Performance for Different Values of α (30 Moves)

Using α = 0.1, the search is the quickest but is the least accurate. Using α = 0.5, the estimation result is better than if α = 0.1 while the search time is longer than if α = 0.1. Using α = 0.999, the estimation result is almost the same as α = 0.5, but the search time is the longest. According to EquationEq. (14), different values of α indicate different shares of f0 and f1 in the reward function. While searching for radioactive sources, the future and the current posterior PDF are similar, and α = 0.5 is chosen for the best source search.

IV.D. Search Path

Next, the set of Monte Carlo simulations compares the performance of the described Tsallis divergence–based strategy with the performance of a uniform search along a prespecified path. The number of movements is set to 30; the rest of the setup is the same as in Sec. IV.A. Experiments are taken with one and two radioactive sources in the search region. The results are shown in and and

Fig. 6. Uniform search of a radioactive source.

Fig. 6. Uniform search of a radioactive source.
.

TABLE II Search Performance for Different Search Strategies

Fig. 7. Uniform search of two radioactive sources: (a) 60 moves; (b) 120 moves.

Fig. 7. Uniform search of two radioactive sources: (a) 60 moves; (b) 120 moves.

It can be seen from and and that when the number of moves is almost equal, the accuracy of the radioactive source parameter estimation based on the Tsallis divergence method is better than the uniform search method. For a given number of moves, the search method based on Tsallis divergence can collect more information and get more accurate estimation results. To achieve the same estimation, the uniform search method requires more moves.

V. PRACTICAL EXPERIMENTS

A radiological experiment trial was conducted on a large, flat, and open area without any obstacles. The following experiment was designed. A mobile robot with a 2 × 2-in. NaI scintillator detector was used to detect a 137Cs source with an activity of about 17.95 MBq in the search region. The mobile robot is about 20 cm high. The NaI scintillator detector is mounted on the top of the robot. The Masero software was used to read the detector data to obtain the count value. The experimental scenario and apparatus are shown in , and the schematic diagram of the search experiment process is shown in .

Fig. 8. The experimental scenario and NaI detector.

Fig. 8. The experimental scenario and NaI detector.

Fig. 9. Search for radioactive sources.

Fig. 9. Search for radioactive sources.

The floor tiles in are 60 × 60 cm. A rectangular area of 10 cells × 10 cells between the two columns is chosen for the experiment. The measurement points are set on the four vertices of the floor tiles. The detector is used to detect at each measurement point for 60s, and the average value is taken as the count rate of that measurement point. First, four points in the search region are randomly selected for detection, and the average radiation background count rate of the environment can be obtained as 2.18/s. Then, the mobile robot is placed on a vertex of the search region and starts running the algorithm. It is assumed that the prior distribution of the position of the sources obeys the uniform distribution and that the prior distribution of the relative source intensity of the sources obeys the gamma distribution, e.g., x0,i:U3,3, y0,i:U3,3, and ε0,i\~Γ(1.5,100). The number of particles is set to 5000, the number of detections is set to 15, the value of α is set to 0.5, the parameter σ is set to 2, and the radioactive source is set at (48 cm, −175 cm). The experimental results are shown in .

Fig. 10. Search results.

Fig. 10. Search results.

As can be seen from , in the initial stage of the search, the robot is far away from the radioactive source, and the move path planned by the control algorithm module has a high randomness. Then, as the robot approaches the source, the count value of the detector rises. The control algorithm could plan the search path. When the robot is close to the source, the multiplication factor in the control algorithm works to control the measurement points from getting too close to the source. Finally, the counting information of the measurement points is used to complete the task of searching and locating the radioactive source through the particle filtering algorithm.

Next, experiments are conducted to compare the performance of the described Tsallis divergence–based strategy with the performance of a uniform search along a prespecified path. The number of movements is set to 60; the rest of the setup is the same as in Sec. IV.A. Experiments are taken with one radioactive source in the search region. The results are shown in and . From and and , it can be seen that the accuracy of the estimation of the radiation source parameters based on the Tsallis divergence is better than that of the uniform search.

TABLE III Search Performance for Different Search Strategies

Fig. 11. Performance of uniform search: (a) 30 moves; (b) 60 moves.

Fig. 11. Performance of uniform search: (a) 30 moves; (b) 60 moves.

VI. CONCLUSION

In this paper, a search algorithm based on the Tsallis divergence for searching the radioactive source is proposed. The number of radioactive sources and their location could be estimated by the algorithm. The search region is a flat open-field rectangular-shaped area, where measurements are collected by a robot equipped with a NaI(Tl) scintillator detector. The control strategy is done sequentially by computing the Tsallis divergence. Through Monte Carlo simulations and practical experiments, it is verified that the radioactive source search algorithm based on Tsallis divergence can complete the radioactive source parameter estimation.

In future work, the existence of obstacles in the search region should be considered to make the hypothetical conditions more in line with the real scene. In the search conditions, increasing the number of mobile platforms carrying detectors to cope with larger search regions by linking multiple detectors should be considered; this can save calculation time.

Disclosure Statement

No potential conflict of interest was reported by the author(s).

Data Availability Statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.

References