# IMM Iterated Extendedℋ∞ Particle Filter Algorithm

Authors: Yang Wan; Shouyong Wang; Xing Qin
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

In order to solve the tracking problem of radar maneuvering target in nonlinear system model
and non-Gaussian noise background, this paper puts forward one interacting multiple model (IMM) iterated extendedℋ∞ particle filter algorithm (IMM-IEHPF). The algorithm makes use of multiple modes to model the target motion form to track any maneuvering target and each mode uses iterated extended ℋ∞ particle filter (IEHPF) to deal with the state estimation problem of nonlinear non-Gaussian system. IEHPF is an improved particle filter algorithm, which utilizes iterated extended ℋ∞ filter (IEHF) to obtain the mean value and covariance of each particle and describes importance density function as a combination of Gaussian distribution. Then according to the function, draw particles to approximate the state posteriori density of each mode. Due to the high filter accuracy of IEHF and the adaptation of system noise with arbitrary distribution as well as strong robustness, the importance density function generated by this method is more approximate to the true sate posteriori density. Finally, a numerical example is included to illustrate the effectiveness of the proposed methods.


---

## Body

## 1. Introduction

In the field of manoeuvring target tracking, IMM [1, 2] algorithm is a popular method. It uses multiple modes to model the target motion form to track any maneuvering target and modes transition governed by a first-order Markov chain. However, in every mode a Kalman filter is running that the performance is deteriorated by nonlinearities and non-Gaussian. Although EKF [3, 4] and UKF [5] can deal with the nonlinear problem, they cannot avoid the limitation of Gaussian hypotheses.Recently, particle filter (PF) [6] has been paid attention to due to its capacity of dealing with the state estimation problem in any nonlinear non-Gaussian system and it has been widely applied [7–10]. It represents the required state posteriori density by a stochastic particles (or samples) with associated weights and to compute estimates based on these particles and weights. Particle filter methods for Markovian switching systems have also been proposed in [11–13]. However, a major drawback of these methods is that the number of particles in each mode is proportional to the mode probability. If a mode probability is close to 0, the particles in that mode will be lost in a large number. When the mode probability increases again, particles need to be reconfirmed, which certainly will result in a significant increase of estimation errors in certain period. This phenomenon is known to cause numerical problems and is also discussed in [14]. To overcome this problem, Boers and Driessen [15] combine IMM with particle filter (IMM-PF) and there is an interaction between modes. Each mode uses a fixed number of particles and runs a regularised particle filter. The state posteriori density of the system is a weighted sum of posteriori density of each mode. Subsequently, the approach has been extensively used in variety manoeuvring target tracking scenarios. For example, in [16, 17], IMM-PF was applied to the problem of tracking ground target, and in [18, 19], it was applied to solve the maneuvering target tracking problem in non-Gaussian noise background. The multirate IMM-PF developed in [20] placed emphasis on computation savings, which enabled the particles in each model to be updated at a different rate. An IMM Gaussian particle filter was proposed in [21], which made particle filter able to concurrently process and saved operation time.However, in IMM-PF, the particles in each mode are drawn according to the prior density without any consideration of current measurements, which have a large deviation with the particles generated from the true posteriori density. Particularly in the case of the measurement model with high accuracy or target maneuvering, this issue will be more serious. To import the current measurements into the particle sampling process, an EKPF and an UPF were applied in each mode of the IMM algorithm presented in [20] and in [22], respectively. They used EKF or UKF to construct an importance density function for approximating the posteriori density for the state in each mode, during which the measurements were utilized. And then the particles were drawn from the function. Nevertheless, due to the limitation of Gaussian hypothesis, their filter performance has no obvious improvement in non-Gaussian condition.On the other hand, over the past decades, considerable attention has been devoted toℋ∞ filtering, and many important results have been reported in [23–34], and references therein. One of its main advantages is the fact that it is insensitive to the exact knowledge of the statistics of the noise signals. The idea of substituting an extended ℋ∞ filter (EHF) for EKF or UKF in the PF framework (called PF-EHF) was first proposed in [24]. However, direct utilization of EHF can face the same problem of model linearisation error as EKF. To eliminate this error and further enhance the filter accuracy, this paper uses an iterative technique to improve the Jacobian linearisation of nonlinear system on basis of EHF, namely, linearize the dynamical model at smoothed values and linearize measurement model at filtered value, respectively (called IEHF). Meanwhile, the iterative technique can repeat many times at one step to improve the degree of linearisation approximation of nonlinear system. As the importance density function generated by IEHF, the resulting filter called IEHPF.In this paper, we propose to combine IMM with IEHPF (IMM-IEHPF) for the tracking problem of maneuvering target in nonlinear non-Gaussian system. The filter accuracy of IEHF is higher than that of EKF and UKF, and IEHF is insensitive to the exact knowledge of the noise processes, therefore the importance density function generated by it is more approximate to the true state posteriori density.The structure of this paper is as follows. Section2 presents the problem formulation of the dynamical system model for maneuvering target. Section 3 describes the IMM-IEHPF algorithm to be used for the modes. In Section 4, through simulation example, several algorithms are compared and analyzed. Finally, the conclusion is provided in Section 5.


## 2. System Model and Problem Description

Consider the following discrete nonlinear time-varying systems for maneuvering target tracking:(1)x(k+1)=f(x(k),r(k))+g(x(k),r(k))v(k),k∈ℕ,z(k)=h(x(k))+w(k),k∈ℕ,
where f(·) and g(·) represent the system function and process noise function in mode r(k), respectively. The target state vector is x(k)=[x(k)x˙(k)y(k)y˙(k)]T, where x(k), y(k), and x˙(k), y˙(k) denote the target position, velocity in Cartesian coordinates. T represents transpose. z(k) is the measurements at time k. h(·) is the nonlinear measurement function in mode r(k). The process noise v(k) and the measurement noise w(k) are mutually independent. The mode transition of the system is modelled by a Markov chain with
(2)Pr{r(k)=j∣r(k-1)=i}=πij,i,j∈{1,2,…,U},
where πij is the transfer probability and U is the number of modes in the IMM.The state estimation problem of system described in (1)~(2) solved by IMM can be divided into four parts, which are interaction, filtering, mode probability update, and combination. Among them, the combination is a weighted sum of filtering estimation and mode probability of each mode. Therefore, filtering is the key element. Take model i as an example and in view of Bayesian estimation, the filtering estimation of this mode is to obtain the posteriori density p(xi(k)∣Z(k)) of state xi(k) in model i given all measurements Z(k)={z(1),z(2),…,z(k)}. Based on this a posteriori density, an estimate of the state is, easily obtained as
(3)x^i(k)=∫xi(k)p(xi(k)∣Z(k))dxi(k).
So, the final combination is given by
(4)x^(k)=∑i=1Uui(k)x^i(k),
where x^i(k) and ui(k) represent the state estimation and mode probability of mode i at time k, respectively.It is well known that Kalman filter can get the optimal state estimation in the linear Gaussian system. However, in practice, most cases belong to nonlinear non-Gaussian system, so the analytic expression of posteriori densityp(xi(k)∣Z(k)) is difficult to get. Even though we can get it, the integral operation involved in (3) is very complex. Therefore, it is urgent to find out an effective filtering method to solve the problem of state estimation of nonlinear non-Gaussian system.


## 3. Proposed Algorithm

IMM-IEHPF algorithm proposed in this paper, in fact, is the combination of IMM algorithm and IEHPF algorithm. Compared with IMM-PF, IMM-IEHPF has considered the latest measurements in the particle sampling process by utilizing the prediction and updating mechanism of IEHF. Therefore, the distribution of particles is more approximate to the true state posteriori density and the filter accuracy can be enhanced. We first briefly describe theℋ∞ filter [23] ahead of IEHF and then provide the mathematical model and the concrete implementation steps of IMM-IEHPF algorithm.


### 3.1.ℋ∞ Filter

Consider the following state space model of linear time-varying system:(5)x(k)=F(k)x(k-1)+G(k)v(k),z(k)=H(k)x(k)+w(k),s(k)=L(k)x(k)k=0,1,2,…,
where F(k)∈𝒞n×n, G(k)∈𝒞n×m, H(k)∈𝒞a×n, and L(k)∈𝒞b×n are known and the definition of other parameters is the same as the above. The statistical properties of v(k) and w(k) are unknown and mutually independent. The system initial state x(0) is unknown. s(k) is the signal we intend to estimate. If we want to directly estimate x(k) then we set L(k)=In(In represents the n×n identity matrix). Let x^(k)=ℱ(z(0),z(1),…,z(k)) denote the estimate of x(k) for a given set of measurements Z(k). The estimation error e(k) can be described as
(6)e(k)=s(k)-s^(k)=x(k)-x^(k),
and the ℋ∞ norm [23] can be defined as
(7)J=supx(0),w,v∈h2((∑p=0keT(p)e(p))×(∑p=0kwT(p)w(p)(x(0)-x^(0))Π0-1(x(0)-x^(0))+∑p=0kvT(p)v(p)+∑p=0kwT(p)w(p))-1),
where Π0 is a positive-definite matrix which reflects a priori knowledge as to how close x(0) is to the initial guess x^(0). The ℋ∞ norm can thus be interpreted as the maximum energy gain from the input {Π0-1/2(x(0)-x^(0)),{v(p)}p=0k,{w(p)}p=0k} to the output {e(p)}p=0k. However, the closed-form solution of the ℋ∞ optimal estimation problem is available only in specific cases, so the realization of ℋ∞ optimal filter is a suboptimal filtering.Suboptimalℋ∞ filtering problem: given a scalar γ>0, find an ℋ∞ suboptimal estimation strategy ℱ(z(0),z(1),…,z(k)) that achieves J<γ2. This clearly requires checking whether
(8)γ≥infJ.
The solution of the optimal ℋ∞ filtering problem can be obtained by iterating on γ of the suboptimal ℋ∞ problem.For a givenγ>0, if the [F(k),G(k)] is full rank, then an estimator that achieves J<γ2 exists if, and only if [24],
(9)P-1(k+1k)+HT(k+1)H(k+1)-γ-2In>0,
where P(0)=Π0 is a given positive definite matrices and P(k+1/k) satisfies the Riccati recursion [24]
(10)P(k+1k)=F(k)P(kk-1)FT(k)+G(k)GT(k)-F(k)P(kk-1)[HT(k)InT]×Re-1(k)[H(k)In]P(kk-1)FT(k),
where
(11)Re(k)=[Ia00-γ2In]+[H(k)In]P(kk-1)[HT(k)InT],Ia represents the a×a identity matrix.The covariance updating is(12)P(k)=P(kk-1)-P(kk-1)[HT(k)InT]×Re-1(k)[H(k)In]P(kk-1).ℋ
∞ filter gain is given by
(13)K=P(kk-1)HT(k)(Ia+H(k)P(kk-1)HT(k))-1.So, the state updating equation ofℋ∞ filter is given by
(14)x^(k)=F(k)x^(k-1)+K(z(k)-H(k)F(k)x^(k-1)).ℋ
∞ filter can be extended to nonlinear systems by using Jacobian linearisation technique. The resulting filter is called EHF. Compared to EKF and UKF, the main advantage of EHF is that there is no limitation about system noise distribution and it treats the system noise as an energy-bounded signal. Through adjusting parameter γ, the state estimation error can be minimized.


### 3.2. IEHF

Although EHF can adapt to the case of system noise with non-Gaussian distribution, it also will face the problem of model linearisation error. Therefore, to further enhance the degree of linearisation approximation of nonlinear system, this paper makes improvements on model linearisation technique and proposes IEHF algorithm. Its fundamental idea is as below.Given the filtered estimationx^(k-1) at time k-1. Take the general nonlinear system described in (1) as an example, linearize the dynamical model at x^(k-1), and retain first-order small quantity (for simplicity, omit model variable r(k))
(15)x(k)≈f(x^(k-1))+∂f(x^(k-1))∂x^(k-1)×(x(k-1)-x^(k-1))+v(k).
Similarly, linearize measurement model at state predictive value x^(k/k-1)(16)z(k)≈h(x^(kk-1))+∂h(x^(k/k-1))∂x^(k/k-1)×(x(k)-x^(kk-1))+w(k).It is not difficult to see thatx^(k-1) and x^(k/k-1) have nothing to do with the current measurement z(k). However, the smoothed value x^(k-1/k) and filtered value x^(k) obtained from measurement z(k) are obviously superior to x^(k-1) and x^(k/k-1). If we linearize the dynamical model and linearize measurement model at the smoothed value and filtered value, respectively, the approximation degree of the model linearisation will be further improved. Therefore, this paper adopts the linearisation form as follows.Linearize the dynamical model at the smoothed valuex^(k-1/k)(17)x(k)≈f(x^(k-1k))+∂f(x^(k-1/k))∂x^(k-1/k)×(x(k-1)-x^(k-1k))+v(k).Linearize the measurement model at the filtered valuex^(k)(18)z(k)≈h(x^(k))+∂h(x^(k))∂x^(k)(x(k)-x^(k))+w(k).In order to further improve the filtering accuracy, the linearisation technique can be repeated many times at one step. All these contribute to IEHF algorithm. However, for many nonlinear filtering problems, the performance is not significantly improved after repeating iterations. Usually it ends after two or three times of iterations. According to the above thoughts, the implementation steps of IEHF algorithm are as follows.Assume that the iterative filtering valuex^d(k-1) and covariance Pd(k-1) at time k-1 are known, the iterative filtering value x^c(k) in c(1≤c≤d) iteration process at time k can be solved according to the following expression:
(19)x^c(k)=x^c(kk-1)+Kc{(x^c(kk-1)-x^c-1(k))z(k)-h(x^c-1(k))KcKcKc-hc-1(x^c(kk-1)-x^c-1(k))},
where
(20)x^c(kk-1)=f(x^c-1(k-1k))+fc-1(x^d(k-1)-x^c-1(k-1k)),Pc(kk-1)=fc-1Pd(k-1)f(c-1)T+Q(k-1),Kc(kk-1)=Pc(kk-1)h(c-1)T×{h(c-1)TPc(kk-1)h(c-1)T+Ia}-1,Pc(k)=Pc(kk-1)-Pc(kk-1)×[h(c-1)TInT]Re-1(k)[hc-1In]Pc(kk-1),fc-1=∂f(x^c-1(k-1/k))∂x^c-1(k-1/k),hc-1=∂h(x^c-1(k))∂x^c-1(k).The smoothing value required byc+1 iteration is
(21)x^c(k-1k)=x^d(k-1)+Pd(k-1)f(c-1)T×(Pc(kk-1))-1(x^c(k)-x^c(kk-1)).


### 3.3. IMM-IEHPF

On the basis of IMM-PF [15] algorithm, this section provides the implementation process of IMM-IEHPF algorithm.Step 1 (interaction stage).
Calculate the mixing probability(22)uij(k-1)=1c-jπijui(k-1),c-j=∑i=1Uπijui(k-1),
where uij(k-1) represents the mixed probability of the mode at time k-1. The interactive input state x^oj(k-1) and the corresponding covariance P^oj(k-1) are given by
(23)x^oj(k-1)=∑i=1Ux^i(k-1)uij(k-1),i,j=1,2,…,U,P^oj(k-1)=∑i=1U‍uij(k-1)×{(x^i(k-1)-x^oj(k-1))tP^i(k-1)+(x^i(k-1)-x^oj(k-1))×(x^i(k-1)-x^oj(k-1))T},
where x^i(k-1) and P^i(k-1) are the state estimate and covariance estimation in mode-matched filter i at time k-1, respectively.Step 2 (filtering stage).
(
1
) According to the Gaussian distribution constructed by x^oj(k-1) and P^oj(k-1), each mode draws particles randomly,
(24)xil(k-1)~N(x^oj(k-1),P^oj(k-1)),l=1,2,…,Ns,i=1,2,…,U,
where Ns is the number of particles.
(
2
) Importance sampling. Use IEHF algorithm in Section 3.2 to update every particle xil(k-1) and obtain a combination of a bank of Gaussian distribution N(x^il,d(k),P^il,d(k)),l=1,2,…,Ns,i=1,2,…,U. Then draw particles according to this distribution
(25)x~il(k)~N(x^il,d(k),P^il,d(k)),i=1,2,…,Ns.
Calculate the particle weights [6]
(26)wil(k)=p(z(k)∣x~il(k))p(x~il∣xil(k-1))N(x^il,d(k),P^il,d(k)),
where p(z(k)∣x~il(k)) is likelihoods function, p(x~il∣xil(k-1)) is the priori density, and N(x^il,d(k),P^il,d(k)) represents importance density function generated by IEHF.
Normalizing(27)w~il(k)=wil(k)∑l=1Nswil(k).
(
3
) Resampling. According to the normalized weight, resample to get the new particles {x-il(k),(1/Ns)}i=1,2,…,Ul=1,2,…,Ns.
(4) State estimation of each mode(28)x^i(k)=1Ns∑l=1Ns‍x-il(k).Covariance estimation(29)P^i(k)=1Ns∑l=1Ns(x^i(k)-x-il(k))(x^i(k)-x-il(k))T.Step 3 (update the mode probability).
Particle residual of each mode(30)δil(k)=z(k)-h(x^il(k)),
where h(x^il(k)) is measurement predicted output.Covariance of residual(31)Si(k)=1Ns∑l=1Ns‍(z-i(k)-h(x^il(k)))(z-i(k)-h(x^il(k))T,
where z-i(k)=(1/Ns)∑l=1Ns‍h(x^il(k)) is mean of predicted output over the sample set.Mode likelihood function(32)Li(k)=1Ns∑l=1Ns‍N(δil(k);0,Si(k)).Mode probability(33)ui(k)=1CLi(k)c-i,C=∑i=1U‍Li(k)c-i.Step 4.
Combination(34)x^(k)=∑i=1Uui(k)x^i(k),P^(k)=∑i=1U‍ui(k)(P^i(k)+(x^i(k)-x^(k))(x^i(k)-x^(k))T).


## 3.1.ℋ∞ Filter

Consider the following state space model of linear time-varying system:(5)x(k)=F(k)x(k-1)+G(k)v(k),z(k)=H(k)x(k)+w(k),s(k)=L(k)x(k)k=0,1,2,…,
where F(k)∈𝒞n×n, G(k)∈𝒞n×m, H(k)∈𝒞a×n, and L(k)∈𝒞b×n are known and the definition of other parameters is the same as the above. The statistical properties of v(k) and w(k) are unknown and mutually independent. The system initial state x(0) is unknown. s(k) is the signal we intend to estimate. If we want to directly estimate x(k) then we set L(k)=In(In represents the n×n identity matrix). Let x^(k)=ℱ(z(0),z(1),…,z(k)) denote the estimate of x(k) for a given set of measurements Z(k). The estimation error e(k) can be described as
(6)e(k)=s(k)-s^(k)=x(k)-x^(k),
and the ℋ∞ norm [23] can be defined as
(7)J=supx(0),w,v∈h2((∑p=0keT(p)e(p))×(∑p=0kwT(p)w(p)(x(0)-x^(0))Π0-1(x(0)-x^(0))+∑p=0kvT(p)v(p)+∑p=0kwT(p)w(p))-1),
where Π0 is a positive-definite matrix which reflects a priori knowledge as to how close x(0) is to the initial guess x^(0). The ℋ∞ norm can thus be interpreted as the maximum energy gain from the input {Π0-1/2(x(0)-x^(0)),{v(p)}p=0k,{w(p)}p=0k} to the output {e(p)}p=0k. However, the closed-form solution of the ℋ∞ optimal estimation problem is available only in specific cases, so the realization of ℋ∞ optimal filter is a suboptimal filtering.Suboptimalℋ∞ filtering problem: given a scalar γ>0, find an ℋ∞ suboptimal estimation strategy ℱ(z(0),z(1),…,z(k)) that achieves J<γ2. This clearly requires checking whether
(8)γ≥infJ.
The solution of the optimal ℋ∞ filtering problem can be obtained by iterating on γ of the suboptimal ℋ∞ problem.For a givenγ>0, if the [F(k),G(k)] is full rank, then an estimator that achieves J<γ2 exists if, and only if [24],
(9)P-1(k+1k)+HT(k+1)H(k+1)-γ-2In>0,
where P(0)=Π0 is a given positive definite matrices and P(k+1/k) satisfies the Riccati recursion [24]
(10)P(k+1k)=F(k)P(kk-1)FT(k)+G(k)GT(k)-F(k)P(kk-1)[HT(k)InT]×Re-1(k)[H(k)In]P(kk-1)FT(k),
where
(11)Re(k)=[Ia00-γ2In]+[H(k)In]P(kk-1)[HT(k)InT],Ia represents the a×a identity matrix.The covariance updating is(12)P(k)=P(kk-1)-P(kk-1)[HT(k)InT]×Re-1(k)[H(k)In]P(kk-1).ℋ
∞ filter gain is given by
(13)K=P(kk-1)HT(k)(Ia+H(k)P(kk-1)HT(k))-1.So, the state updating equation ofℋ∞ filter is given by
(14)x^(k)=F(k)x^(k-1)+K(z(k)-H(k)F(k)x^(k-1)).ℋ
∞ filter can be extended to nonlinear systems by using Jacobian linearisation technique. The resulting filter is called EHF. Compared to EKF and UKF, the main advantage of EHF is that there is no limitation about system noise distribution and it treats the system noise as an energy-bounded signal. Through adjusting parameter γ, the state estimation error can be minimized.


## 3.2. IEHF

Although EHF can adapt to the case of system noise with non-Gaussian distribution, it also will face the problem of model linearisation error. Therefore, to further enhance the degree of linearisation approximation of nonlinear system, this paper makes improvements on model linearisation technique and proposes IEHF algorithm. Its fundamental idea is as below.Given the filtered estimationx^(k-1) at time k-1. Take the general nonlinear system described in (1) as an example, linearize the dynamical model at x^(k-1), and retain first-order small quantity (for simplicity, omit model variable r(k))
(15)x(k)≈f(x^(k-1))+∂f(x^(k-1))∂x^(k-1)×(x(k-1)-x^(k-1))+v(k).
Similarly, linearize measurement model at state predictive value x^(k/k-1)(16)z(k)≈h(x^(kk-1))+∂h(x^(k/k-1))∂x^(k/k-1)×(x(k)-x^(kk-1))+w(k).It is not difficult to see thatx^(k-1) and x^(k/k-1) have nothing to do with the current measurement z(k). However, the smoothed value x^(k-1/k) and filtered value x^(k) obtained from measurement z(k) are obviously superior to x^(k-1) and x^(k/k-1). If we linearize the dynamical model and linearize measurement model at the smoothed value and filtered value, respectively, the approximation degree of the model linearisation will be further improved. Therefore, this paper adopts the linearisation form as follows.Linearize the dynamical model at the smoothed valuex^(k-1/k)(17)x(k)≈f(x^(k-1k))+∂f(x^(k-1/k))∂x^(k-1/k)×(x(k-1)-x^(k-1k))+v(k).Linearize the measurement model at the filtered valuex^(k)(18)z(k)≈h(x^(k))+∂h(x^(k))∂x^(k)(x(k)-x^(k))+w(k).In order to further improve the filtering accuracy, the linearisation technique can be repeated many times at one step. All these contribute to IEHF algorithm. However, for many nonlinear filtering problems, the performance is not significantly improved after repeating iterations. Usually it ends after two or three times of iterations. According to the above thoughts, the implementation steps of IEHF algorithm are as follows.Assume that the iterative filtering valuex^d(k-1) and covariance Pd(k-1) at time k-1 are known, the iterative filtering value x^c(k) in c(1≤c≤d) iteration process at time k can be solved according to the following expression:
(19)x^c(k)=x^c(kk-1)+Kc{(x^c(kk-1)-x^c-1(k))z(k)-h(x^c-1(k))KcKcKc-hc-1(x^c(kk-1)-x^c-1(k))},
where
(20)x^c(kk-1)=f(x^c-1(k-1k))+fc-1(x^d(k-1)-x^c-1(k-1k)),Pc(kk-1)=fc-1Pd(k-1)f(c-1)T+Q(k-1),Kc(kk-1)=Pc(kk-1)h(c-1)T×{h(c-1)TPc(kk-1)h(c-1)T+Ia}-1,Pc(k)=Pc(kk-1)-Pc(kk-1)×[h(c-1)TInT]Re-1(k)[hc-1In]Pc(kk-1),fc-1=∂f(x^c-1(k-1/k))∂x^c-1(k-1/k),hc-1=∂h(x^c-1(k))∂x^c-1(k).The smoothing value required byc+1 iteration is
(21)x^c(k-1k)=x^d(k-1)+Pd(k-1)f(c-1)T×(Pc(kk-1))-1(x^c(k)-x^c(kk-1)).


## 3.3. IMM-IEHPF

On the basis of IMM-PF [15] algorithm, this section provides the implementation process of IMM-IEHPF algorithm.Step 1 (interaction stage).
Calculate the mixing probability(22)uij(k-1)=1c-jπijui(k-1),c-j=∑i=1Uπijui(k-1),
where uij(k-1) represents the mixed probability of the mode at time k-1. The interactive input state x^oj(k-1) and the corresponding covariance P^oj(k-1) are given by
(23)x^oj(k-1)=∑i=1Ux^i(k-1)uij(k-1),i,j=1,2,…,U,P^oj(k-1)=∑i=1U‍uij(k-1)×{(x^i(k-1)-x^oj(k-1))tP^i(k-1)+(x^i(k-1)-x^oj(k-1))×(x^i(k-1)-x^oj(k-1))T},
where x^i(k-1) and P^i(k-1) are the state estimate and covariance estimation in mode-matched filter i at time k-1, respectively.Step 2 (filtering stage).
(
1
) According to the Gaussian distribution constructed by x^oj(k-1) and P^oj(k-1), each mode draws particles randomly,
(24)xil(k-1)~N(x^oj(k-1),P^oj(k-1)),l=1,2,…,Ns,i=1,2,…,U,
where Ns is the number of particles.
(
2
) Importance sampling. Use IEHF algorithm in Section 3.2 to update every particle xil(k-1) and obtain a combination of a bank of Gaussian distribution N(x^il,d(k),P^il,d(k)),l=1,2,…,Ns,i=1,2,…,U. Then draw particles according to this distribution
(25)x~il(k)~N(x^il,d(k),P^il,d(k)),i=1,2,…,Ns.
Calculate the particle weights [6]
(26)wil(k)=p(z(k)∣x~il(k))p(x~il∣xil(k-1))N(x^il,d(k),P^il,d(k)),
where p(z(k)∣x~il(k)) is likelihoods function, p(x~il∣xil(k-1)) is the priori density, and N(x^il,d(k),P^il,d(k)) represents importance density function generated by IEHF.
Normalizing(27)w~il(k)=wil(k)∑l=1Nswil(k).
(
3
) Resampling. According to the normalized weight, resample to get the new particles {x-il(k),(1/Ns)}i=1,2,…,Ul=1,2,…,Ns.
(4) State estimation of each mode(28)x^i(k)=1Ns∑l=1Ns‍x-il(k).Covariance estimation(29)P^i(k)=1Ns∑l=1Ns(x^i(k)-x-il(k))(x^i(k)-x-il(k))T.Step 3 (update the mode probability).
Particle residual of each mode(30)δil(k)=z(k)-h(x^il(k)),
where h(x^il(k)) is measurement predicted output.Covariance of residual(31)Si(k)=1Ns∑l=1Ns‍(z-i(k)-h(x^il(k)))(z-i(k)-h(x^il(k))T,
where z-i(k)=(1/Ns)∑l=1Ns‍h(x^il(k)) is mean of predicted output over the sample set.Mode likelihood function(32)Li(k)=1Ns∑l=1Ns‍N(δil(k);0,Si(k)).Mode probability(33)ui(k)=1CLi(k)c-i,C=∑i=1U‍Li(k)c-i.Step 4.
Combination(34)x^(k)=∑i=1Uui(k)x^i(k),P^(k)=∑i=1U‍ui(k)(P^i(k)+(x^i(k)-x^(k))(x^i(k)-x^(k))T).


## 4. Simulation Results and Analysis

In order to validate the filtering performance of IMM-IEHPF algorithm, this paper compares standard IMM, IMM-PF, and IMM-IEHF algorithm. The experiment scene is designed as follows. Radar scanning interval isT=1s and repeat for 100 times. The target moves in straight line in the former 34s, then it turns anticlockwise with turning speed of ω=0.1rad/s until 67 s and moves in straight line in the last 33 s. Adopt the system model in (1)~(2) and the target state x(k)=[x(k)x˙(k)y(k)y˙(k)]T, which includes the position (m) and velocity component (m/s) in x-axis and y-axis. Mode r(k)∈{1,2,3}, r(k)=1 corresponds to the uniform rectilinear model. The r(k)=2 and r(k)=3 correspond to the anticlockwise turning mode and clockwise turning mode, respectively. With the known initial state x(0)=[200m15m/s150m6m/s]T, the target state was evoluted according to different system functions as follows.Mode1 (CV mode, r(k)=1):
(35)f(x(k),1)=[1T000100001T0001],g(x(k),1)=[T220T00T220T],Q1(k)=diag([52225222]).Mode2 (CT mode, r(k)=2):
(36)f(x(k),2)=[1sinωTω0(cosωT-1)ω0cosωT0-sinωT0(1-cosωT)ω1sinωTω0sinωT0cosωT],g(x(k),2)=[T220T00T220T],Q2(k)=diag([1025210252]),
where ω is the turning coefficient.Mode3 is the same as Mode 2; however, in mode 2, ω>0 corresponds to the anticlockwise manoeuvre, and in mode 3, ω<0 corresponds to the clockwise manoeuvre. The initial mode probabilities are set as u1(0)=u2(0)=u3(0)=1/3, and the mode transition matrix [πij] is given by [15]
(37)[πij]=[0.960.020.020.020.960.020.020.020.96].We measure the range and bearing, and the measurement functionh(·) is mode- and time-independent and is given by
(38)h(x(k))=[x2(k)+y2(k)arctany(k)x(k)]T.Measurement noise standard deviations areσr=30m for the range measurement, σθ=2mrad for the bearing, and the number of particles is Ns=1000 [15]. In IMM-IEHPF algorithm, parameter γ=1.6. Usually, in order to ensure the robustness of the algorithm, parameter γ should be as small as possible when meeting the existence conditions of (9). At present, it can only be got through trial and error. Iteration times d=3. We carried out Monte Carlo simulation for M=100 times and the position root mean square error (RMSE) at time k is defined as
(39)RMSE(k)=(1M∑mot=1M[(y^mot(k)-y(k))2(x^mot(k)-x(k))2+(y^mot(k)-y(k))2]1M∑mot=1M[(y^mot(k)-y(k))2(x^mot(k)-x(k))2)1/2.Figure1 compares the state estimation generated from a single run of IMM, IMM-IEHF, IMM-PF and IMM-IEHPF. For observing the state estimation diversity of the four methods, we have magnified Figure 1.Figure 1
Position estimate.As seen from Figure1, these four methods do well in the target state estimation in the two straight-line movements. However, in the anticlockwise turning period, the deviation from the true values obtained from IMM is larger than those of the other three methods. Among them, IMM-IEHPF has the best state estimation in the maneuvering stage. We can see that the trajectory of IMM-IEHPF is in accordance with the true trajectory.To compare the filtering accuracy of the previous four methods, Figure2 shows the position root mean square error of the four methods.  From Figure 2, we observe that the IMM filter diverges after t=34s until t=67s. The reasons for this is that in the IMM filter the speed is wrongly estimated. This is due to the nonlinearity in the measurements in the maneuvering mode. Compared with IMM, the position RMSE of IMM-IEHF, and IMM-PF also increase to certain extent, but both are less than that of IMM. IEHF improves the degree of approximation of the nonlinear system through iteration and PF can deal with nonlinear filtering problem by describing posteriori density through a series of particles, which achieve better filtering accuracy than EKF. However, IEHF is just a linear approximation for nonlinear system and improvement of filtering accuracy is limited. PF adopts the prior density as the importance density function without considering the latest measurement information. When the target maneuvering and the measurement information change greatly, the particles generated from the function may deviate from the true posteriori density, resulting in a bad filtering performance.Figure 2
Position RMSE.The IMM-IEHPF algorithm, however, can deal with target maneuvering and we see that the filter also performs well aftert=34s until t=67s. IEHPF takes advantages of IEHF and PF together in itself. In the filtering process, IEHPF not only considers the latest measurements in the particle sampling process but also reduces the influence of error caused by linearizing measurement model. When these algorithms are directly applied to the IMM framework, posteriori density of each mode is described much more well by IEHPF than EKF, IEHF, and PF. Therefore, the filtering accuracy of IMM-IEHPF is the highest.To quantitatively analyse the filtering performances of the four methods, Table1 lists the mean value and variance of position RMSE of each method in different observed bearing standard deviation, which stands for different disturbances of receiver noise and atmospheric. From Table 1, we can see that the mean value and variance of position RMSE of all four methods grow with the increase of the observed bearing standard deviation σθ, nevertheless, which of IMM-IEHPF are less than those of the others (which is in accordance with the analytical results in Figure 2), especially when the bearing standard deviation is considerable high, our proposed method still demonstrates higher tracking precision (i.e., σθ=10mrad, taking the mean of RMSE as the metric, it reduces by 57.86% compared with IMM, and 7.53% compared with IMM-PF).Table 1
Mean and variance of RMSE under different measurement bearing deviations.
Algorithm
Bearing standard deviationsσθ (mrad)
Mean
Variance
IMM
2
9.7293
98.2215
IMM-IEHF
2
4.0394
5.6625
IMM-PF
2
4.1397
2.4750
IMM-IEHPF
2
2.7593
0.3633
IMM
5
12.794
159.49
IMM-IEHF
5
6.5694
7.5583
IMM-PF
5
5.2413
2.9841
IMM-IEHPF
5
5.3559
1.1192
IMM
10
18.421
333.87
IMM-IEHF
10
8.6552
10.3592
IMM-PF
10
8.3949
8.2179
IMM-IEHPF
10
7.7631
3.4886


## 5. Conclusions

This paper has investigated the problem of maneuvering target tracking in the nonlinearities of the dynamic state-space and the non-Gaussian measurement noise. We have presented a new IEHPF method in the IMM framework. The performance of the proposed method is evaluated via simulations and compared with IMM, IMM-IEHF, and IMM-PF, which illustrates that the tracking performance of our proposed method is superior to the others. We have also provided insight into the influence of various bearing standard deviations. Three bearing standard deviations are given for comparing the performance of our proposed method with the other methods. We have found that IMM-IEHPF is weakly sensitive to noise levels. The simulation results have well verified the validity of the proposed method.This paper only considers a maneuvering target appearing in the scanned region of radar at whole time step. However, in actual radar target tracking, target may disappear at any time step. Therefore in future work, we will study how to construct a suitable model to describe target appearance or disappearance in maneuvering target tracking under a low SNR environment.



---

*Source: 970158-2013-04-22.xml*