Data collected at more than one sensor can be combined or considered separately. The former is often referred to as the “data fusion” step. Algorithms without data fusion make decisions based on readings from a single detector only and thus require no sensor-to-sensor or sensor-to-server communication. On the other hand, algorithms that fuse sensor data require the data to be transferred to one or a few centralized servers (base stations) for processing, creating more complexity in system architecture and introducing higher costs. Therefore, it is essential to understand and quantify the benefits of sensor data fusion when designing a radiation sensor network. For example, if the independent analysis is almost as good as algorithms with data fusion, communication between sensors can be simplified.

Data fusion is undoubtedly advantageous for discovering and recognizing a source, but it is unclear how much if at all, it improves detection. While simply integrating data (e.g., summing all data without filtering) improves signal strength, it also introduces noise, lowering SNR. The combined SNR may be lower than considering SNR from each sensor individually. For example, view a grid of widely-spaced sensors and a weak static source near one detector. The detector closest to the source will have a high SNR, whereas the other sensors receive mostly noise. Combining the data blindly from all four sensors will always decrease the overall SNR. More quantifying the benefits from data fusion can be found in (Chin et al., 2010; Liu, 2010).

In a noise- and error-free situation, determining the location of a radioactive source problem can be solved in closed form employing four ideal sensors and the Apollonius circle (Chin et al., 2008). When uncertainties and noise like background radiation are considered, more sensors are required to give correct results, especially for extremely low source intensities. In real-world scenarios, noise and measurement errors are unavoidable. Hence, the fusion of measurements from n sensors, n > 3, is necessary to increase the solution’s reliability.

The localization of a point radiation source of unknown strength A _u expressed in the unit of micro-Curie (μCi) called the source rate will be considered. The source is located at an unknown location (x _u, y _u ). The source gives a radiation intensity of I ( p ) = I ( x , y ) = A u ×   E × 2.22 × 10 6 ( x u − x ) 2 + ( y u − y ) 2 (1) (expressed in counts per minute or CPM) when measured by a sensor at point p = (x, y), where E is an efficiency constant unique to the sensor. The radiation count induced by the source and observed at the sensor i per unit time is a Poisson random variable with parameter λ = I(x _i , y _i ), not accounting for the background radiation B(x, y) (Rao et al., 2008) denoting the background radiation strength at (x, y) expressed in CPM, called the background rate. The radiation count measurement (due to the background radiation) at a sensor i located at (x _i , y _i ) is given by the Poisson random variable with parameter B(x _i , y _i ). For comparison purposes, the radiation model, actual measurements, and calibration process (Rao et al., 2008) are followed.

Apollonius Circle describes the locus of point p, which has a fixed distance from point A equal multiple of its distance from another fixed point B. If the m:n equals one, the locus is the perpendicular bisector of segment AB. If m ≠ n, then it is a circle (Kim, 2021).

In Figure 1, AP:BP = m:n. When the point P moves, keeping this ratio, the locus of p is a Circle, which is called the circle of Apollonius. This circle connects A and B’s interior and exterior division points (the interior bisector PP1 and the exterior bisector PP2 in Figure 2).

Particle filter is a sequential Monte-Carlo approach used to estimate the dynamic state parameters of nonlinear and/or non-Gaussian systems (Fox et al., 1999; Marimon et al., 2007). The essential idea is to approximate the probability density functions (PDFs) of the state of a dynamic model by random samples (particles) with associated weights, then propagating them across iterations based on a probabilistic model of the state update and the measurements (Arulampalam et al., 2002).

J. Cook et al. (Cook et al., 2020) proved the theoretical convergence of the estimate to the state posterior probability density function of Particle Filter. So, weak convergence in radioactive source localization applications is not a problem.

A. Kyriakis and K. Karafasoulis (Kyriakis and Karafasoulis, 2020) pointed that Particle Filter algorithm complexity increases exponentially with the number of sources. Therefore, the next section of this paper presents a H/W implementation of the proposed Particle Filter to speed up the processing that can help deal with this situation and real-time scenarios.

The Particle Filter Range Search algorithm works as follows. Initially, the 1,000 particles are evenly distributed, i.e., the initial state of the source is entirely unknown. Then, it converges to a source position estimation after using the radius of the Apollonius circle in the particle filter likelihood function. Figure 3 illustrates the flow of the proposed Algorithm. First, two consecutive sensors (e.g., S1 and S2) readings are used to obtain a candidate position for the radioactive source. Next, two points (p₁, p₂) on an Apollonius circle are calculated. Then the radius of the circle r calculated according to Eq. 2 is used in the likelihood function of the particle filter. Finally, a probability is assigned to each particle according to the standard Gaussian formula, and their weights are calculated as in Eq. 3. r = ( x s − x ) 2 + ( y s − y ) 2 (2) Where xs and ys are positions of the sensors in the environment. w ( q r m ) = 1 σ 2 π e − 1 2   ( r ∧ − r m σ ) 2 (3) Where r_m is the range measurement (Apollonius)

q = [x, y] is the state vector

r^ is the range estimate from the particle to the sensor

and σ: is the standard deviation from the PDF for that range measurement.

A pool of candidate particle positions is calculated, and an estimated source position is plotted as shown in Figure 3A. For the 2^nd round, the next consecutive sensor (S3) reading with the previous sensor reading (S2) is used to form another Apollonius circle, as shown in Figure 3B. Finally, the particles are resampled from the previous round’s Apollonius circle perimeter particle pool. Figures 3C,D show the following rounds, which will continue till utilizing all 18 sensors.

The pseudo-code of the proposed algorithm can be stated as follows.

The proposed architecture, shown in Figure 8, uses 2MxP dual-port FIFO to store the particles, where M is the number of particles and p is a two-component state vector, given by (x, y) to represent the X-position and Y-position. Two registers files are used: The first Mxw register file stores the particle weights where w represents the length of weight bits. The second Mxn register file keeps the replication factor for each particle (M = 2 ⁿ ), where n is the replication number. Figure 8 illustrates the data flow through processing stages of the proposed radioactive source localization H/W architecture. The operations of this PF H/W architecture can be summarized as follows.

1) In sample engine, the particles are generated from samples drawn using the random number generator (Mahmoud et al., 2016). Then, in each moving step, all particles are drifted by the same amount. The generated particles are stored in a FIFO for routing when it has a replication value. As all particles are generated, they can be used for weight calculation.

The straightforward implementation of the output calculation is accomplished by employing two accumulators to sum each of x and y values for all particles, which is then multiplied by its weight and multiplied by the inverse value of the total weight.

4) After sampling the total number of particles is done, the resample process starts and gives a replication factor for each particle according to its weight using one of the resampling methods. The basic idea of resampling is to eliminate small-sized particles and concentrate on particles with large weights. The resampling process calculates sequentially the replication factor for each particle for the next instant. The resample calculates one replication factor for every clock cycle.

In Residual Resampling (RR), resampling is performed where the number of replications of a specific particle is determined in one loop by truncating the product of the number of particles, and the simple truncation may result in a total number of replicated particles less than M. In general, R = ∑r ^(m) may not be equal to M. For solving this problem, for example, if the sum of all replicated particles is less than, M the hardware must decide which particles to replicate additionally so that the total number of particles is constant.

In the RR algorithm, for M input and N output (resampled) particles, the summation of all particles replication factors are less than the input particles M (N < M). The difference between input and output is obtained using the mechanism (Abd El-Halym et al., 2012), in which we added this remainder number to the replication factor of the highest weight particle.

The replication of each particle is calculated as follows: R i = ( w i ∗ M ∑ 1 M w i )   (4)

The simple structure of the resampling is executed by employing the lookup Table, which is used to calculate the inverse value of the weight summation, which lies between 1 and 64, then multiply this value by M. Accordingly, the number D = ( M / ∑ 1 M w i ) is calculated. The replication factor R i = ( w i ∗ D ) where i = 1, 2, 3 ….M values are calculated and stored in the replication factor memory.

The architecture is captured using VHDL and synthesized on a Xilinx Virtex 5 FPGA chip. Finally, the design is verified using ModelSim simulator. Table 4 summarizes the total resources utilization of the implemented architectures.

The proposed architecture uses two memory arrays: one MxW to store the weights of the particles where W equals to 8 bits used in our fixed-point representation and the other Mxn memory to store the replication factor where M = 2 ⁿ FIFO is implemented as a circular queue, where data can be stored and retrieved (in the order of its entry only) and thus has a read pointer and a write pointer. A synchronous FIFO uses the same clock for reading and writing. Moreover, total memory requirements for this architecture are: one dual-port FIFO of 2M words to store the particles vectors of width equals 24 bits.

The timing diagram of our proposed architecture operations is shown in Figure 14. Ti is the initial latency in particle generation, which equals two clock cycles. Therefore, the total cycle time required is TSIRF = (TW + Tres + αM) Tclk. The particle generation and computation weight operations can be overlapped in sequence such that TW is the startup latency of the weight calculation unit, which is equal to one clock. Therefore, the total latency cycles are three. Tres is the resampling time equal to M + 1 clock periods (one clock to add the M-N to the highest weight replication factor). In the routing of the particles, the value of α lies between 1 and 2 and depends on the replication factors of the particles. The best value α is 1, where every clock routs one particle such that the total routing time is M. The worst-case takes place when the first M−1 particles have zero replication factors. Only the last stored particle in the FIFO has a replication factor that equals M in this case; the FIFO needs to read M−1 particles first, till FIFO reaches the particle number M. Then the sampling, weighting, and the output calculations are started and repeated for those particle M times. For all other cases (when a considerable number of the particles have non-zero replication factors), the value of α is slightly more than one.

From the timing summary, the total time consumed for SIRF with RSR is 28.958 ns (maximum frequency: 34.533 MHz) using Xilinx FPGA chip xc5vlx50t with speed grade: -3.

Several architectures are implemented and compared to explain the tradeoff between speed up and resource utilization. The design occupies a small part of the targeted Xilinx chip. So, smaller and more affordable chips (e.g., Xilinx Zynq-7010) will be sufficient.