Multitarget Tracking by Particle Filtering Based on RSS Measurement in Wireless Sensor Networks

We propose an algorithm for multitarget tracking by particle filtering in wireless sensor networks based on received signal strength (RSS) measurement where we also localize a newly appearing target whose location and reference power are unknown. Therefore, the number, the reference power, and the initial locations of targets are unknown in this problem. At the initial localization step, we apply approximate least squares (LS) method to roughly estimate the target location. After the initial location is estimated, we estimate the reference power. This is possible because we can use multiple number of measurements for estimating multiparameters. The proposed approach is particularly emphasized on the initialization step that completes the whole multitarget tracking system by particle filtering in a challenging scenario. The proposed approach is validated by computer simulations for its effectiveness.


Introduction
In wireless sensor networks, we can take advantage of many sensors to estimate certain states of targets. In this paper, we propose an algorithm to detect and track targets based on the received signal strength (RSS) measurement where an unknown number of targets vary with time, and the "reference power" of the targets is also unknown and varies with time. In fact, RSS measurement model is not very popular in the literature regarding target tracking or target localization [1,2] because, whereas bearing and/or range measurement provides one-to-one mapping between a measurement and a parameter of interest, RSS measurement generally does not provide one-to-one mapping except for the cases when the measurement is obtained over the sensor fields by using "a grid spacing" and when "standard commercial radiofrequency (RF) prediction tools" are used. Therefore, we do not assume "one-to-one" mapping between a measurement and a target in this paper. A salient feature of the RSS measurement is that the parameter to be estimated is directly related with the measurement while bearing and/or range measurement is not. This means that the data fusion methods, such as time of arrival (TOA) [3], time difference of arrival (TDOA) [4][5][6], or angle of arrival (AOA) [7,8], require a preprocessor before acquiring the measurement (bearing or range) while RSS [9][10][11] is almost raw data from the target. Therefore, we might be able to apply maximum likelihood (ML) method directly to RSS measurement-based estimator.
The least squares (LS) method is successfully employed for localization methods using the data fusion methods, for example, TOA, TDOA, AOA, and RSS, specifically, in a mobile locating systems [12,13] or sensor node localization in a wireless ad hoc sensor networks [14]. Interestingly, a maximum likelihood (ML) method can be applied to the RSS model to improve the performance beyond the LS method [15]. Once a target is identified, particle filtering can be applied to track the dynamic states of the target. The problem of multitarget tracking is not an easy task while various approaches, for example, extended Kalman filter, data fusion, and joint probabilistic data association, were proposed in the literature [16][17][18][19][20][21].
We assume that the initializing step is applied to only a new single target identification. First, we directly estimate the distance between a possible new target and each sensor based on measurements and formulate the equation with respect to the squared-estimated-distance for each sensor. And then 2 International Journal of Distributed Sensor Networks we can cancel out the quadratic part of equations before applying linear LS method. After initializing the location of the new possible target, we estimate the reference power given the estimated location when the reference power is unknown. The contributions of the paper consist in proposing an algorithm for multitarget tracking in wireless sensor networks where initial positions of targets, the reference power of the targets, the number of targets are unknown and estimated. A number of statistical approaches are jointly employed in this algorithm such as least squares method, iterative maximum likelihood method, and particle filtering. The challenge consists in tracking multitargets based on a single signal strength transmitted from targets. This paper is organized as follows. In the following section, the problem is formulated in a mathematical model. Next, we show how to initialize targets. The initialization process is described for both cases of known reference power and unknown reference power. And then particle filtering is described for tracking varying number of targets. The computer simulations to verify the validity of the proposed algorithm follow. The paper is concluded at the final section.

Measurement and State Space Model.
We consider a 2D rectangular field of interest with uniformly distributed sensors. The RSS measurement at the sensor from source targets at the time step is described as follows [22]: where l , is the location of the target at the time step , is the sensor index, is the number of targets, Ψ is the reference power of the target , that is, the received power from the source at the reference distance 0 , r is the sensor location, is the attenuation factor ( ≥ 1), V is the background zero-mean Gaussian noise, and is the total number of sensors used in the field. The RSS measurement model does not provide much information about targets, for example, the number of targets, the distance from sensors, and the direction of the source. We use three or four best sensors, that is the three strongest received signals when we apply the LS method with known reference power while four best measurements are employed with unknown reference power for initialization of a newly appearing target. In the tracking system, we have to be able to estimate the location and the reference power of targets, in addition to the number of targets at every time step based on a single "signal strength" at each sensor, which is a challenging problem.
We model a single moving target system by a linear state space model as follows [1,18,23]: where x = [ 1, 2,1,2, ] ⊤ is the state vector which indicates the position and the velocity of a target, respectively, in a 2D Cartesian coordinate system. G and G are known matrices according to the classical dynamics defined by where u is the Gaussian noise-like acceleration perturbation and is sampling time-period (s). Targets maneuver with random acceleration based on classical dynamics (discrete-time version). Only a part of the state is related with the measurement in the RSS model because the state comprises location and velocity while only the location of a target contributes to the measurement. Nevertheless, once we estimate any one of them, we can relate them and find the rest of the states since all states are related by classical mechanics. The Doppler effect should be taken into account for a target with considerably high speed in a real scenario in order to obtain enhanced performance.

Varying
Pattern of the Number of Targets. We use particle filtering to track the variable and multiple number of targets in this paper. At any time step, particles are propagating according to the dynamic state space model. We assume a certain pattern by which the number of targets varies between two consecutive time steps. We assume that the number of targets varies according to three patterns as follows [24].
(1) The number of targets remains the same as at the previous time step with the same identities.
(2) The number of targets increases by a newly appearing target.
(3) The number of targets decreases by one, that is, in the tracked targets at the previous time step.
If the sampling time-period for the discretization is short enough, the above assumption will be satisfied.

Initializing a Single Target with Known Reference Power
The approximate LS method is not the best performing method for the localization [15]. However, if we use particle filtering for target tracking, a fast and simple initializing method can be more efficient (as LS method does) rather than a complicated and highly accurate method because particle filtering will retrieve the right track eventually even with a low-accuracy location-estimate. Therefore, we jointly employ the approximate LS method and particle filtering for initialization and tracking, respectively. Generally, the prior distribution of a target, that is, the probabilistic distribution of the target, is assumed to be known in the application of particle filtering.

Least Squares Method.
Let us suppose that we have m unknowns and n linear equations, and we can describe the problem by the equations as follows: where is a n-by-m, x is a m-by-1, and is a n-by-1 matrices, respectively. Depending on the sizes of m and n, we can consider the problem in three cases. If m > n, we have a smaller number of equations than that of unknowns. In this case the solution is underdetermined or incompletely specified. In this case, there are many solutions that satisfy the above equation. An approach to this case is known to be minimum norm solution that has the criterion as follows: min ‖x‖ such that x = .
If the rank of is n, then ⊤ is invertible, and a solution can be obtained as follows: where ⊤ ( ⊤ ) −1 is called pseudoinverse of for the underdetermined problem. If n = m, is a square matrix, and, unless is singular, we can obtain a unique solution by using the inverse of as follows: Nonetheless, if is singular, then we have no solution or many solutions.
Finally, if m < n, we have a bigger number of equations than that of unknowns. In this case the equations are inconsistent and the solution is overdetermined. In this case, no solution exists that satisfies all of n equations. Therefore, our concern is now to find x of the best approximation that is close to a solution in terms of least squared errors for all n equations. The criterion of this approach is to minimize the norm of errors, which is described as follows: The error vector satisfies the following orthogonal condition: Then, the following is satisfied: If has full rank, then ⊤ is invertible, and the least squares solution can be obtained as follows: where ( ⊤ ) −1 ⊤ is the pseudoinverse of the matrix for the overdetermined problem. Therefore, linear least squares solution can be obtained in the overdetermined problem. Please note that pseudoinverse has different forms depending on whether the problem is underdetermined or overdetermined. Furthermore, both formulations of pseudoinverse are identical with the exact inverse of a matrix when m = n (on condition that is not singular). Therefore, as long as we have two or more linear equations with two unknowns, least squares method can find either "exact solution" or "an approximate solution in terms of minimized norm of errors. " In this paper, we have two or more number of equations than the number of unknowns (i.e., two). We need three measurements of quadratic equations in order to obtain two linear equations. Besides, we need four quadratic equations of measurements to obtain two linear equations when unknown reference power is jointly estimated.
We can use sensor measurements as many as we want. However, the signal to noise ratio decreases drastically as the distance increases between a sensor and a target. On the other hand, there is an ambiguity about better information among a number of measurements, particularly, between the second, the third, and the fourth. This is why we may want to employ more linear equations than that of unknowns sometimes while we need two linear equations most of the time in this paper. In this case, we need to obtain pseudoinverse for overdetermined problem. As shown in [25], increased number of linear equations improves the performance significantly for the case of "time-of-arrival (TOA)" or "time-difference-of-arrival" measurement while it is less effective for the problem employing RSS measurement. As a matter of fact, we do not need more than two linear equations because the information beyond two is redundant and even degrades the performance due to low SNR most of the time. Therefore, we use three sensor measurements with known reference power and four sensor measurements with unknown reference power.
From (1), we estimate the squared distance between sensors and a target location before applying the LS method when we assume = 2 and 0 = 1 m (omitting the time index) as follows: If we denote the estimated distances by , , and and corresponding locations of three sensors by ( 1 , 2 ), ( 1 , 2 ), and ( 1 , 2 ), respectively, the LS method solves the equations as follows: and ( ⊤ ) −1 ⊤ is the pseudoinverse of , and x = [ ] ⊤ is the initialized target location. In order to obtain linear equations, we need to remove the quadrature part of equations by using the formulation of " 2 − 2 " and " 2 − 2 . " Therefore, we need three measurements to obtain two linear 4 International Journal of Distributed Sensor Networks equations for estimating 2D position. Similarly, we need four measurements in order to apply this approach for estimating 3D position. The performance assessment of the LS initial localization is presented in the following section.

Performance of Initial Localization.
In this section, we assess the performance of the LS and the ML methods in initial localization. These two methods have different features. Whereas the LS method does not consider any probabilistic assumptions, the ML method solves the problem in a probabilistic manner. Therefore, the ML approach finds the solution that maximizes the probability density function (PDF) given the measurement as follows: ln (y; l) where Ψ 0 = , = 2, |r − l| = D , and y = [ 1 2 3 ] ⊤ because we use only 3 measurements. Consequently, the ML method finds the value by which the first-order derivative of the log likelihood function becomes zero as follows: Because a solution in the closed form does not exist for the above equation, we need to apply an iterative method (i.e., the Newton-Raphson method) to find l that satisfies (16) as follows: where d is defined in Appendix A, and the solution is derived therein. In the ML approach, we adopt the solution of the LS method as an initial guess (l 0 ). Therefore, the ML approach certainly requires more steps beyond the LS method. We also compare the performances of these methods with Cramer-Rao bound (CRB) which is derived in Appendix B. When we apply the iterative method, we have to be careful because it can diverge sometimes, especially when the second-order derivative of the log-likelihood function is small [26]. Therefore, in the simulation, we set the threshold for the second-order derivative such that when inverse of derivative is larger than the threshold, it stops the iteration and adopts the LS estimate as the solution in order to avoid the divergence. We select 10 as the iteration number. The number of iterations was selected after extensive simulations with various values of iteration numbers. The performance was improved clearly as we increased the iteration number up to 5 in the computer simulations while, sometimes, we have enhanced performance with values greater than 5. Therefore, we safely select 10 as the iteration number for assuring that we do not obtain improved performance any more. The 2D locating system field is depicted in Figure 1. ). We assume that the noise power is specified in dBW and that the unit of measure for the noise is volts. For power calculations, it is assumed that there is a load of 1 ohm. If we compare each pair of the result when the target is located at (20,20), Figures  2-4 show that the ML method compresses the distribution of LS estimates more densely around the true value. Some of ML iterative estimates have relatively larger errors. This is because the ML iterative estimates diverge sometimes. We stop the iteration and employ the LS estimate in this case.
The simulation results are shown in Figures 5-7 where the values are averaged for results with respect to various locations of the targets. Figure 5 shows the mean of distance error of two methods. As it shows, the iterative ML method outperforms the LS method. In Figures 6-7, the performances of two methods are compared to the CRB. The CRB is a lower-bound for the variance of estimates of an arbitrary unbiased estimator. When the noise power is low (0.01 dB), the ML iterative method obtains performance almost the same as the CRB. The exact values of the simulation results are summarized in Table 1      -coordinate for the LS method is 1.19 m 2 , and the variance of the realized estimates in direction for iterative ML method is 0.13 m 2 , respectively. It shows that the variance of the estimates for the iterative ML method is below the CRB when the noise power is 0.01 dB, which means that the iterative ML estimator is not unbiased. The bias is caused since we approximated the ML estimator by the numerical iterative method, and, furthermore, our ML estimator does not have a large data record for measurement to be asymptotically optimal. We can see the highly improved performance by the  ML approach compared to that of the linearized LS method as shown in the results.

Least Squares for Multitargets.
One of the most difficult scenarios in target tracking problem arises when the multiple number of targets varies. The LS method needs to be modified to be applied when the number of targets varies. We suppose that the number of targets varies following the varying pattern as specified at the beginning of Section 2.  we reformulate the measurement equations with respect to the distance between the new target and each sensor. Then, we apply the approximate LS method to the newly formulated equations. The whole procedure for initializing a new target is described in Algorithm 1.
Algorithm 1 (initializing a new target in multitargets via least squares method (with known reference power). (i) We suppose the number of sensors in the field is , and the number of continuing targets is . At any time step, we receive the measurements 1: as follows and we define = 10 log 10 ( where Γ = ∑ =1 (Ψ 0 /|r − l | ) which is an estimated or predicted part of measurement for the continuing targets, Υ = |r − l +1 | which is the estimated distance between the new target and each sensor, Ψ +1 is the reference power of the new target, and l 1: is the estimated or predicted locations of targets propagating from the previous time step.
(iii) Send the information of {r min , Υ min } to the LS algorithm, and obtain the initial location of the newly appeared target using additional two neighboring sensors of the best sensor (r min ). Find the best two neighboring sensors that have shorter distances than the rest of the neighbors (make sure that these three sensors form a "right triangle" but not straight line).

Initializing a Target with Unknown Reference Power
If the reference power of the target is unknown, the problem is more challenging. We initialize a new target location first and, then, estimate the power of the new target. The path loss exponent is a function of the environment. The value = 2 corresponds to a free-space environment, which is unrealistic while we can employ it without loss of generality in this section. From (12), we cancel out the common reference power (i.e., Ψ) using a couple of equations for the sensors and as follows: Then, we have new multiple equations to apply the linear LS method as follows: We use = 1 and = 2, 3, 4 (i.e., four strongest received signals) in order to take advantage of good measurement information (high signal to noise ratio) even though a different combination of sensor measurement can be employed. The approximate LS method is applied to the newly formulated linear equations instead of quadratic equations. Once the target location, l, is estimated, we can directly apply the ML method for estimating the power of the target from the multiple equations. If we rewrite the measurement equations, + 10 log 10 (Υ 2 ) = 10 log 10 Ψ + V .
If we set If a new target is initialized in multiple targets with unknown reference power, the solution needs to be slightly modified.
We have the measurement equation as After canceling out the power, The rest of the steps are the same as in the initialization with known reference power, and only is differently formulated.
In multiple targets case, we cannot apply the ML method directly for estimating the reference power of the new target because we cannot solve the ML function directly in this case due to the Γ factor, and the maximum log-likelihood function has the formula as follows: Therefore, we apply the LS method one more time for estimating the power of the new target in this case.

Particle Filtering as the Solution to Tracking System
The initialization algorithm is only a part of whole tracking system by particle filtering. When the reference power of the target is unknown, the generated particle will have one more element in the state vector, and the initialization algorithm will obtain the location and reference power of the new target. When we apply particle filtering, we need to select an appropriate model among all possible models. The model which has the maximum weight sum will be selected according to the predefined varying patterns of the number of targets. Any particle generates the offspring of a new target even if it can be removed after the "model selection" step or resampling step [27]. Therefore, the "initialization step" is applied to every single particle at every time step of the tracking system.

Sampling Importance Resampling Particle Filter.
There are many versions of particle filtering [28][29][30] and our choice for a tracking system is "sampling importance resampling (SIR)" particle filter. If we denote state function, measurement function, state, measurement, and the weight by f (⋅), h (⋅), , y , and (where is the time step index and is the particle index), respectively, the importance density, ( | −1 , y 1: ), will turn out to be the prior density, ( | −1 ), and because resampling is executed at every time step, the weight can be computed by ∝ (y | ), which is the likelihood function.

SIR Particle Filter Combined with Initialization Algorithm.
In a multiple and varying number of target-tracking system, the state space equation must include the state of the number of targets. We denote it by K at the time step . K has three patterns to propagate as the previous assumption; that is, We propagate equal number of particles for each model. The number of targets follows the model for the next time step. Therefore, every single particle will produce multiple descendants depending on the number of possible models. The number of models at the next time step depends on K −1 . For instance, if K −1 = 0, then the possible models will be either 0 or 1. If K −1 = > 0, then the possible K = − 1, , and + 1. However, when > 1, it will have −1 more multiple models because targets have equal possibility of disappearance. When K = 0, the state space will be empty except for the number of targets. The posterior function of interest will be (K 1: , S 1: , x 1: | y 1: ) (where S is the reference power of the target and x is the state vector), and the distribution is approximated by the "random measure"; that is, where is the total number of the particles. We can compute the weight when we use SIR particle filter; that is, If we use only the three best sensors and assuming that sensors are not correlated, then We need to apply the initialization algorithm when we generate a particle that has a newly appearing target (the algorithm generates the estimated location and reference power): this is the case when K = K −1 + 1. In Algorithm 2, we described how to apply initialization, combined with particle filtering. In Algorithm 2, we apply the initialization algorithm that is explained in Section 4 for estimating the distance between each sensor and the new target (Υ , ) and the reference power of the new target. If we cancel out the reference power, we do not use neighboring sensors anymore as shown in Algorithm 2. Note that when the reference power is known, the initialization algorithm has to use the sensors that are neighboring to each other (see Algorithm 1).
(c) Model 3: K +2+ = 2 (particle + 2 + is generated from particle ). We suppose , = 10 log 10 ( where Γ = ∑ which is predicted part of measurements for the continuing targets, Υ , = |r − l new, | which is the estimated distance between the new target and each sensor, Ψ new is the reference power of the new target, and l , is the predicted location of targets propagating from the previous time step, which is propagated from the previous particle.

Computer Simulations
In this section, we present computer simulations based on the proposed approach. The validity of the proposed algorithm is verified by using MATLAB version 7.12.0.635 (R2011a). All simulations were performed by -files coded by the authors. We apply the algorithm shown in Algorithm 2 at  every time step for the whole tracking system. We assume that the reference power of the target is unknown and varies as The initial reference power is given by 200 and 180 (J/s) to each target. The number of particles is 300 when we apply particle filtering. A target appears at the coordination of (0, 150) initially, and, two time steps later, another target appears at the coordination of (200, 0) and both targets remain together until the time step is 62. From the time step 63, the second target disappears and only first target lasts until the time step 80. Acceleration-noise follows zero mean Gaussian distribution with variance 1, background noise power is 0.01 (dB), the variance of the reference power is 0.001, and the number of sensors is 5-by-5 = 25. A single simulation result is depicted in Figure 8. We ran 1000 simulations of the scenario that generates the result of Figure 8. The results show that generally the algorithm performs well with accurate detection of a new target and disappearance of the target. When we initialize a new target, the target which has unduly low reference power and the target that is located out of the field of interest are removed even if it is initialized as a new target.  Table 2.

Summary and Conclusions
In this paper, we proposed a general multi-target-tracking algorithm by particle filtering based on the RSS measurement model where multiple targets are unidentified, the number of targets varies, and also unknown reference power of the target varies. We introduced the initialization step in the tracking algorithm where the LS and ML methods are applied as the substeps. The proposed approach can be applied to very crucial practical problems such as tracking mobile phone users or wearable devices in a wireless cellular system or tracking moving sensor nodes in ad hoc wireless sensor networks even though generally the reference power is known in these practical problems. This approach also can be applied for ubiquitous 3D positioning system in emerging visible-light communications systems based on lightemitting diodes. The proposed approach assumes generally a number of parameters are unknown, that is, the number of targets, the reference power of the target, and the initial location of the target, which provides possibilities of using the proposed approach in highly challenging environment.

A. Derivation of Newton-Raphson Method
From (1), if we set Ψ 0 = , = 2, and |r − l| = D , then the log likelihood function becomes ln (y; l) where y = [ 1 2 3 ] ⊤ because we use only 3 measurements. It can be simplified as where is a constant which does not depend on the parameter we want to estimate, and we can rewrite it as We need to apply iterative method to find l that satisfies (A.7) as follows: Before we proceed, let us make some abbreviations as follows: We use the initial guess of l 0 from the solution of the "linearized LS method" by using three measurements.

B. Cramer-Rao Bound (CRB)
We derive the CRB of the parameter, that is, the location of a "source target, " in this section. The parameter is denoted by