Multi-sensor particle ﬁltering with multi-step randomly delayed measurements

This paper develops particle ﬁltering for multi-sensor systems with randomly delayed measurements, where the general case that random delay can be multi-step rather than one-step or two-step is considered. Moreover, different sensors can have different delay steps and delay probabilities. Random delays are assumed to be mutually independent for different sensors and modelled by a separate sequence of random variables obeying discrete dis-tributions. Since random delay leads to the actual measurements being dependent rather than independent given states, and this dependence becomes more complicated with the increase of random delay step, a new formula of the local likelihood density is proposed and then a new weighting scheme is adopted in particle ﬁltering to deal with these difﬁculties. The proposed method is applied to two examples to testify its effectiveness and superiority.

and diverge especially for high-dimensional state space. Moreover, due to the inherent Gaussian assumption, they will yield a poor filtering estimate in strongly non-linear or non-Gaussian cases.
Particle filtering (PF) is a promising alternative to Gaussian approximation filters since it is not subject to the Gaussian assumption and can provide a decent estimate with sufficient samples [13][14][15]. Thus it is attractive to utilise PF to deal with RMD. PFs are proposed to deal with one-step RMD and multistep RMD in [16,17] and [18], respectively. For the multi-sensor system with multi-step RMD, if randomly delays of all sensors occur in a synchronous manner, the PF proposed in [18] can be directly used by augmenting the measurement equation. However, different sensors generally have different delay steps and delay probabilities, that is to say, random delays are asynchronous for different sensors. To tackle this case, a PF is proposed in [19].
All aforementioned PFs assume that the current actual measurement is independent of the previous actual measurements conditioned on the state sequence, which is impractical in the presence of measurement delay and brings a theoretical problem to these PF algorithms (See Remark 1 in Section 4 for details). To solve this problem, an improved PF is developed for the onestep or two-step delay in [20] and [21], respectively, where a new local likelihood density is designed. We find that the calculation of the local likelihood density requires to take all possible combinations of a segment of actual measurements into consideration, and the number of all possible combinations is a power exponential function of delay step (See Remark 2 in Section 4 for details). When delay step is small, one can enumerate all possible combinations as has been adopted in [20] and [21]. But for multi-step RMD, which is more likely to occur in practical applications, combinations become extremely intricate, and enumeration method is obviously impractical. Hence, the motivation of this paper is to design a particle filtering which can deal with the measurement dependence for multi-step asynchronous RMD in multi-sensor systems.
The remainder of the paper is organised as follows. In Section 2, we present the multi-sensor measurement model with asynchronous random delays. In Section 3, we review the core of the generic PF and the standard multi-sensor PF. In Section 4, we study the dependence of multi-step randomly delayed measurements and then propose a new PF to deal with multi-step RMD. Numerical simulations are shown in Section 5. Section 6 concludes the paper.

PROBLEM FORMULATION
This paper considers the following non-linear dynamic system with multiple independent sensors: where x k is the latent state with known initial probability density p(x 0 ), z k,m the ideal (undelayed) measurement of the mth sensor at time step k ∈ ℕ, and M the number of sensors. The system noise {v k , k ≥ 0} and the measurement noise {e k,m , k ≥ 0} M m=1 are white noise processes with arbitrary probability densities p v k (⋅) and p e k,m (⋅), respectively. Now assume that the ideal measurements {z 0,m } M m=1 are always available for use but that for k ≥ 1, the ideal measurements {z k,m } M m=1 may be randomly unavailable. Then the multi-sensor measurement model with asynchronously multi-step random delays can be modelled as where y k,m is the actual measurement of the mth sensor at time step k, and k,m denotes the delay step of the mth sensor at time step k. Let d m be the user-defined maximum possible measurement delay for the mth sensor. Then k,m obeys a d m -valued discrete distribution where the random event k,m = s represents that an s-step measurement delay occurs for the mth sensor at time step k, and p s k,m is the s-step delay probability for the mth sensor at time step k. It is assumed that k 1 ,m is independent of k 2 ,m when k 1 ≠ k 2 , which means that delays occur at different time steps are mutually independent for the same sensor m. In addition, we assume that k 1 ,m 1 is independent of k 2 ,m 2 when m 1 ≠ m 2 . This indicates that delays between different sensors are mutually independent. We further assume the mutual independence of . For notational simplicity, we denote by z k = {z k,1 , … , z k,M } the ideal measurement set of M sensors, y k = {y k,1 , … , y k,M } the actual measurement set of M sensors, and a i: j the set {a i , a i+1 , … , a j } of any general sequence{a k }.
For the non-linear system (1) with randomly delayed measurements (2), our objective is to design particle filtering to approximate the filtering density p(x k |y 0:k ). To fulfil this goal, we quote two widely used properties in the Bayesian filtering framework as follows [22].

Generic particle filtering
In particle filtering context, the joint density p(x 0:k |y 0:k ) is usually approximated by a set of weighted random samples as where (⋅) is the Dirac delta function. The importance weight w (i ) k is computed recursively by where q(x k |x 0:k−1 , y 0:k ) is known as the importance density function.

Standard multi-sensor particle filtering
In the multi-sensor case, since sensors are assumed to be independent of each other, the joint likelihood density p(y k |x 0:k , y 0:k−1 ) can be computed as where p(y k,m |x 0:k , y 0:k−1,m ) is called the local likelihood density. In the undelayed case, the actual measurement sequence y 0:k is exactly the ideal measurement sequence z 0:k . Then by using Properties 1 and 2, the importance weight (5) is simplified into where the joint likelihood density (6) is simplified into Moreover, if we choose q(x k |x k−1 , y k ) = p(x k |x k−1 ) and conduct resampling at each time step, the standard multi-sensor particle filtering, also known as the central particle filter [23] can be obtained.

MULTI-SENSOR PARTICLE FILTERING WITH RMD
In the case of RMD, the actual measurements are disordered, repeated and missing in general. This fact leads to that Properties 1 and 2 cannot be directly applied to actual measurements, and resultantly the importance weight (5) cannot be simplified into (7). Hence, when the standard multi-sensor particle filtering is used by roughly ignoring the RMD, it will inevitably yield a poor state estimate.
To develop a new particle filtering which can capture the information of actual measurements correctly, the importance weight (5) is required to be reformulated. More specifically, we need to figure out expressions of p(x k |x 0:k−1 , y 0:k−1 ) and p(y k,m |x 0:k , y 0:k−1,m ). To this end, we shall give the following two lemmas.

Lemma 1.
Conditioned on x k−1 , x k is independent of x 0:k−2 and y 0:k−1 , that is, Proof. Since the actual measurement set {y 0 , y 1 , … , y k−1 } is a subset of the ideal measurement set {z 0 , z 1 , … , z k−1 }, (9) can then be obtained by using Property 1. □ For the simplified form (10), according to the assumption that elements of { s,m } k s=1 for the fixed m are mutually independent and also independent of the state sequence and measurement sequence, we have

Lemma 2. The local likelihood density is simplified into
where By dividing all possible combinations of j k−d m ,m , … , j k,m into the following d m + 2 cases, we can obtain explicit expressions of ( j k−d m ,m , j k−d m +1,m , … , j k,m ).
In this case, the number of combinations is (d m + 1) d m , and the probability of these combinations is p 0 k,m . Since y k,m takes the quantity z k,m , and none of y k−d m :k−1,m can take the quantity z k,m , by using Property 2 we have Case 2: j k,m = 1 and j k−1,m ≠ 0.
In this case, the number of combinations is d m (d m + 1) d m −1 , and the probability of these combinations is p 1 k,m (1 − p 0 k−1,m ). Since y k,m takes the quantity z k−1,m and none of y k−d m :k−1,m can take the quantity z k−1,m , y k,m does not depend on y k−d m :k−1,m conditioned on x k−d m :k,m . By using the Bayes formula and Property 2, we have Hence In this case, the number of combinations is d s−1 m (d m + 1) d m −s+1 , and the probability of these combinations is where we have where C j k,m denotes the indicator function of set C j k,m . In summary, we obtain the following theorem.
where a j k,m = In fact, the measurement y k,m can only occur a k-step delay at most when k < d m , that is, p j k,m = 0 for k < j ≤ d m . Thus Equation (22) in Theorem 1 can be simplified by discarding the combinations of j k−d m ,m , j k−d m +1,m , … , j k,m with zero probability.
In Theorem 1, the dependence of the actual measurements y k−d m :k,m is indeed considered. But the local likelihood density given by Theorem 1 has an unusual form because of the presence of Dirac delta functions and indicator functions, which brings a trouble for practical calculation. However, when y k and at least one of y k−d m :k−1,m take the same quantity, the second term in (22) is ∞, and resultantly the local likelihood density is ∞. This means that y k,m is uninformative about x k−d m :k and should be discarded. The above analysis can provide us with a solution for this trouble. Let k,m be the index satisfying whered m = min{k − 1, d m }. Then the local likelihood density can be assigned as Remark 3. k,m has the ability to distinguish the actual measurements which are repeated. Moreover, the probability that y k,m is repeated is ℙ( k,m = 1) = a k,m . Repeated measurements are regarded as valueless and will be discarded in (27).
Weighting: • Compute w (i ) k by (28), and normalisew In this paper, we use the transition density as the proposal density and conduct resampling at each time step. The recursive formula of importance weight (5) is thus simplified into (28) The resulting particle filtering for multiple sensors with RMD is presented in Table 1.

NUMERICAL SIMULATIONS
In this section, the efficiency of the proposed particle filtering is shown by two typical non-linear examples. One is the nonstationary growth model, and the other is the two-dimensional radar tracking model. The following filtering algorithms will be used.
Here, newMSPFRMD is our proposed method, while others are used for comparison. It is worth noting that

Example 1: Non-stationary growth model
The non-stationary growth model has been widely used as a benchmark problem to testify the performance of non-linear filtering. It is described as where the initial state follows x 0 ∼  (0, 25), and noise processes {v k , k ≥ 0} and {e k,m , k ≥ 0} 3 m=1 are mutually independent, white and Gaussian distributed according to v k ∼  (0, 10) and e k,m ∼  (0, 1) for m = 1, 2, 3, respectively.
Suppose that the actual measurement y k,1 suffers from a twostep random delay, and y k,2 and y k,3 suffer from a three-step random delay. The concrete delay probabilities are given in Table 2.
To make a fair comparison, we conduct K = 1000 independent Monte Carlo runs with T = 500 time steps. In each Monte Carlo run, N = 200 particles are employed for all particle filtering methods, and the initial state are randomly chosen from its initial probability distribution.
We firstly use the single sensor 2 to investigate the performance of all filtering methods. For clarity, we only present the RMSE curves of the first 100 and last 50 time steps for all filtering methods in Figures 1 and 2, respectively. Table 3 displays the TARMSEs and the total running time. Clearly, SMSPFideal performs best in terms of RMSE as it uses ideal measurements. By contrast, SMSPFignored has a very low estimation accuracy. The remaining three methods, CQKFRMD, MSPFRMD and newMSPFRMD, have mechanisms to handle the multi-step RMD. But they show different performances. Both MSPFRMD       and newSMSPFRMD can deal with RMD well, but the latter performs slightly better and consumes less time than the former. This is because the dependence of measurements is considered, and repeated measurements are discarded. CQKFMRD gives the worst results due to its inherent Gaussian assumption, and it requires the most running time due to state augmentation.
Since the multi-sensor CQKFRMD has not been developed yet, we use four other methods to estimate the state in the multisensor case. The curves of RMSEs for the first 100 and last 50 time steps are plotted in Figures 3 and 4, respectively, and simulation results are summarised in Table 3. Similarly to the single sensor case, SMSPFideal gives the best result, followed by newMSPFRMD. By contrast, SMSPFignored has the worst performance. Even though MSPFRMD costs more running time, its estimation error is bigger than newMSPFRMD.

Example 2: Two-dimensional radar tracking
Assume that the target is moving in the x − y plane according to the following dynamics The target is observed by two radar stations, and the corresponding measurement equations are given by where z k,m is the ideal measurement at time step k for the mth radar, which contains the radial distance r k,m and the bearing k,m for m = 1, 2. The measurement noises e k,1 and e k,2 are mutually independent white noise processes and obey e k,1 ∼  (0 2 , R k,1 ) and e k,2 ∼  (0 2 , R k,2 ), respectively, where R k,1 = R k,2 = diag(100 m 2 , 10 −6 rad 2 ∕s 2 ). It is assumed that the measurements from these two radars suffer from asynchronously random delays, and the actual measurement y k = [y k,1 , y k,2 ] satisfy (2). The maximum possible measurement delays for two sensors are set to d 1 = d 2 = 3, and the corresponding delay probabilities are given in Table 4.
The following results are for K = 500 Monte Carlo simulations with T = 500 time steps and N = 2000 particles in each    Table 5 lists TARMSEs and the total running time. It is clear that SMSPFideal has the smallest estimation error, while SMSPFignored has the largest estimation error. Our proposed method, newMSPFRMD, outperforms MSPFRMD in terms of RMSEs for all four state components, position and velocity. Moreover, it can save almost 15% of the total running time compared with MSPFRMD.

CONCLUSION
A new multi-sensor particle filtering is proposed to deal with the dependence of multi-step randomly delayed measurements. The simulation results validate the superiority of this method in terms of estimation accuracy and running time. It is very general and can be applied to various cases including: (a) The single sensor with multi-step random delay. (b) Multiple sensors with different delay steps and delay probabilities. (c) Delay probabilities can be time varying for the single sensor or multiple sensors.