# Closed-Loop Estimation for Randomly Sampled Measurements in Target Tracking System

Authors: Jin Xue-bo; Lian Xiao-feng; Su Ting-li; Shi Yan; Miao Bei-bei
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

Many tracking applications need to deal with the randomly sampled measurements, for which the traditional recursive estimation method may fail. Moreover, getting the accurate dynamic model of the target becomes more difficult. Therefore, it is necessary to update the dynamic model with the real-time information of the tracking system. This paper provides a solution for the target tracking system with randomly sampling measurement. Here, the irregular sampling interval is transformed to a time-varying parameter by calculating the matrix exponential, and the dynamic parameter is estimated by the online estimated state with Yule-Walker method, which is called the closed-loop estimation. The convergence condition of the closed-loop estimation is proved. Simulations and experiments show that the closed-loop estimation method can obtain good estimation performance, even with very high irregular rate of sampling interval, and the developed model has a strong advantage for the long trajectory tracking comparing the other models.


---

## Body

## 1. Introduction

Target tracking is the most important preliminary step for many higher-level analysis applications. Nowadays, some new sensors have been used in the tracking systems, such as the radio frequency identification (RFID) readers. The RFID stores and retrieves data through the electromagnetic transmission to an RF compatible integrated circuit. Once the tag gets close to the readers, the distance between the readers and tags can be got and sent to the data processing center. The measurements of RFID are randomly sampled [1] because of the data-driven measurement mechanisms. Data-driven approach has been used in many applications [2], and the irregular sampling is one of the important issues in this approach.In general, the video tracking system has to extract the visual information at each frame [3], which costs much computing amount. In [4], the target is tracked by some selected frames to reduce the calculation cost and achieve the real-time tracking, which also results in the randomly sampled tracking problem. If the output measurements are obtained at a set of irregular sampling times, the traditional recursive estimation from K to K+1 may fail in general [5]. Both the model and the estimation method should be reconsidered.Reference [6] transformed the randomly sampled measurement tracking to some time-varying parameters and used the current model to describe the processing model [7–11], which assumes a priori probability density of the acceleration as Rayleigh density. Due to the randomly sampled measurement, this assumption is no longer satisfied.Except the current model, there were several other models used in the tracking, such as constant-velocity (CV) model, constant acceleration model (CA), and Singer model (zero mean first-order Markov model) [12, 13]. The CV models [1] emphasize that the accelerations are small. In maneuvering target tracking, the inclusion of acceleration in the state vector would degrade tracking performance. The main attractive feature of this model is its simplicity. It is sometimes used in the maneuvering target tracking techniques, such as the so-called noise-level adjustment, when the maneuver is quite small or random. It is also simply referred to as the CA model or more precisely the nearly CA model. The Singer model regards the target acceleration as a first-order semi-Markov process with zero mean, which is in essence a priori model since it does not use online information on the target maneuver. Again because of the irregular sampling time, the priori model does not meet the actual dynamic model of the target.The approach to update the system model online has attracted great interest of the researchers. For example, the interacting multiple model (IMM) [14, 15] method considers the change of the system dynamics as a Markovian parameter, whose transition probability is set based on the online estimation and then fusions several models for the tracking, while IMM suffers heavy computational burden on condition that the maneuvering target has complex motion. Moreover, the complex movement can also lead to frequent switch between different models, which can cause the tracking performance to decline. Another model [16, 17] estimated the state of a power system, where the bus voltages are transformed to a system parameter. But the works of this closed-loop estimation have not yet been involved in the randomly sampled tracking system.This paper will develop a joint state-and-parameter estimation method for the target tracking system with randomly sampled measurements, where the estimation problem is reformulated as two loosely coupled linear subproblems. This paper is organized as follows. Section2 derives the system dynamic model under the random sampling time and gives the estimation method based on Kalman filter. The convergence of the algorithm is proved in Section 3. The simulations and experiments are provided in Section 4. Finally, some concluding remarks are given in Section 5.


## 2. System Model and Closed-Loop Estimation Method

We begin with the continuous dynamic model of the moving target. Letx,x˙, and x¨ be the target location, velocity, and acceleration along a generic direction, and the state is expressed as x=[x,x˙,x¨]T. Assume the nonzero mean acceleration satisfies x¨(t)=a-(t)+a(t), where a-(t) is the mean of acceleration in the interval [0t] and a(t) is a zero mean first-order stationary Markov process with variance δa2. We have a˙(t)=-αa(t)+w(t); α is maneuver frequency and w(t) is zero mean processing white noise with variance δw2=2αδa2. The parameter α=1/τ is the reciprocal of the maneuver time constant τ and thus depends on how long the maneuver lasts. For example, for an aircraft τ≈60s for a lazy turn and τ≈10–20s for an evasive maneuver. The parameter δa2=E[a2(t)] is the “instantaneous variance” of the acceleration.Then, we can obtain the acceleration satisfyingx¨˙(t)=-αx¨(t)+αa-(t)+w(t) and the following state-space representation of the continuous time model can be obtained:
(1)x˙(t)=Ax(t)+Ua-(t)+Bw(t),
where
(2)A=[01000100-α],U=[00α],B=[001],
and w(t) is process noise with the covariance matrix given by w(t)~N(0,2αδa2). Assume the measurement data is obtained at the sampling time ti and the measurement equation is as follows:
(3)z(ti)=H(ti)x(ti)+v(ti),i=0,1,2,…,
where H(ti) is measurement matrix and v(ti) is measurement noise with the known variance R; that is, v(ti)~N(0,R).


### 2.1. Model Discretization

We can get the following by the differential equation (1):
(4)x(t)=eA(t-t0)x(t0)+∫t0teA(t-λ)Ua-(λ)dλ+∫t0teA(t-λ)Bw(λ)dλ.We can see that for any known integration interval [t0,t], x(t) can be gotten at any time t if the initial state x(t0), the parametersA, U, B, a-(t), and w(t) in [t0,t] are known.Consider the time interval fromti-1 to ti and assume
(5)a-(λ)=a-(ti-1),w(λ)=w(ti-1),λ∈[ti-1ti];
we can have
(6)∫ti-1tieA(ti-λ)Ua-(λ)dλ=∫ti-1tieA(ti-λ)Udλa-(ti-1)∫ti-1tieA(ti-λ)Bw(λ)dλ=∫ti-1tieA(ti-λ)Bdλw(ti-1).Setthi=ti-ti-1; we have the system matrix as Ad(ti-1)=eAthi, Ud(ti-1)=∫ti-1tieA(ti-λ)Udλ and the noise wd(ti-1)=∫ti-1tieA(ti-λ)Bdλw(ti-1) with the covariance Qd(ti-1)=E[wd(ti-1)wdT(ti-1)].Because the process matrixA in (2) is not a full-rank matrix, we cannot calculate the matrix exponential eAthi by the Lagrange-Hermite interpolation. Here, we use Laplace transform and have
(7)(sI-A)-1=[s-100s-100s+α]-1=[1s1s21s2(s+α)01s1s(s+α)001s+α].The matrix exponentialeAthi can be gotten by the inverse Laplace transform as
(8)Ad(ti-1)=[1thiαthi-1+e-αthiα2011-e-αthiα00e-αthi]
and by the similar approach, we can get the system parameter
(9)Ud(ti-1)=[1α(-thi+α·thi22+1-e-α·thiα)thi-1-e-α·thiα1-e-α·thi]
and the variance of wd(ti-1) as
(10)Qd(ti-1)=E[wd(ti-1)wdT(ti-1)]=2αδα2[q11q12q13q12q22q23q13q23q33]
with the parameters described as
(11)q11=12α5[2α3thi331-e-2α·thi+2α·thi+2α3thi33-2α2thi2-4α·thie-α·thi]q12=12α4[e-2α·thi+1-2e-α·thi+2α·thie-α·thi-2α·thi+α2thi2]q13=12α3[1-e-2α·thi-2α·thie-α·thi]q22=12α3[4e-α·thi-3-e-2α·thi+2α·thi]q23=12α2[e-2α·thi+1-2α·thi]q33=12α[1-e-2α·thi].Then, we get the discrete state-space model of the tracking system as(12)x(ti)=Ad(ti-1)x(ti-1)+Ud(ti-1)a-(ti-1)+wd(ti-1)z(ti)=H(ti)x(ti)+v(ti),
where x=[x,x˙,x¨]T is the state of the system to be estimated and whose initial mean and covariance are known as x0 and P0, wd(ti) and v(ti) are white noise with zero mean and independent of the initial state x0, z(ti) is the measurement vector, H(ti) is measurement matrices, and v(ti) is measurement noise with known variance R. Until now, the irregular sampling is turned to the varying-parameter system. We can see the same sampling interval is just a particular case of the random sampling problem. Therefore, the model of the randomly sampling tracking is a general one.


### 2.2. System Parameters Estimation

Here, we assume the maneuver frequencyα and the variance of the acceleration δa2 are not constant but variable and expressed as αi and δai2. From the processing model of (12), we have the discrete time equation of the acceleration as
(13)x¨(ti)=βix¨(ti-1)+(1-βi)a-(ti-1)+wa(ti-1),
where βi=e-αithi and wa(ti-1) is a zero mean white noise sequence with the variance
(14)δawi2=δai2(1-βi2).αi is the maneuver frequency at the sampling time ti. a-(ti-1) is the mean of one interval, so we have a-(ti)=a-(ti-1). Set a(ti)=x¨(ti)-a-(ti); then we can obtain
(15)a(ti)=βia(ti-1)+wa(ti-1).Consider the estimation of accelerationa^(ti) is a random process; we have
(16)a-(ti-1)=1i∑i=0i-1x¨^(ti),
where i is the number of data. For a first-order stationary Markov process (15), we have the statistics relation between the autocorrelation functions r(0), r(1) with the parameters βi and δawi2 by the Yule-Walker method [18]
(17)ri(0)=1i∑i=0i-1a^(ti)a^(ti)ri(1)=1i∑i=1i-1a^(ti)a^(ti-1)βi=ri(1)ri(0)δawi2=ri(0)-βiri(1).Next, we can getαi and δai2 by δai2=δawi2/(1-βi2),αi=ln⁡βi/-thi, and then get the system parameters Ad(ti-1), Ud(ti-1), and Qd(ti-1) in process function (12).


### 2.3. Algorithm Summary

Now we summarize the closed-loop estimation algorithm for the randomly sampled measurements as follows.(
1
)   Initialization(i=0). Consider
(18)x^(t0∣t0)=x0,P(t0∣t0)=P0,α0,δa02,a-(t0),r0(t0)=x¨0·x¨0,r0(t1)=x¨0.(
2
)  Recursion(i:=i+1)(a)
System update: setthi=ti-ti-1 and the system parameter as
(19)A^d(ti-1)=[1thiαithi-1+e-αithiαi2011-e-αithiαi00e-αithi](20)U^d(ti-1)=[1αi(-thi+αi·thi22+1-e-αi·thiαi)thi-1-e-αi·thiαi1-e-αi·thi]and the variance of the wd(ti-1) as
(21)Q^d(ti-1)=E[wd(ti-1)wdT(ti-1)]=2αiδαi2[q11q12q13q12q22q23q13q23q33]
with parameters described as
(22)q11=12αi5[2αi3thi331-e-2αi·thi+2αi·thi+2αi3thi33-2αi2thi2-4αi·thie-α·thi]q12=12αi4[e-2αi·thi+1-2e-αi·thi+2αi·thie-αi·thi-2αi·thi+αi2thi2]q13=12αi3[1-e-2αi·thi-2αi·thie-αi·thi]q22=12αi3[4e-αi·thi-3-e-2αi·thi+2αi·thi]q23=12αi2[e-2αi·thi+1-2αi·thi]q33=12αi[1-e-2αi·thi].(b)
State prediction: consider(23)x^(ti∣ti-1)=A^d(ti-1)x^(ti-1∣ti-1)+U^d(ti-1)a-(ti-1)P(ti∣ti-1)=A^d(ti-1)P(ti-1∣ti-1)A^dT(ti-1)+Q^d(ti-1).(c)
State update: consider(24)x^(ti∣ti)=x^(ti∣ti-1)+K(ti)[z(ti)-H(ti)x^(ti∣ti-1)](25)K(ti)=P(ti∣ti-1)HT(ti)×[H(ti)P(ti∣ti-1)HT(ti)+R(ti)]-1(26)P(ti∣ti)=[I-K(ti)H(ti)]P(ti∣ti-1).(d)
Parameter adaptation: the mean of the acceleration(27)a-(ti-1)=1i∑i=0i-1x¨^(ti∣ti).Wheni≤K0, the maneuver frequency αi is set to α0 and the covariance of the noise δai2 is gotten by the following:
(28)δαi2={4-ππ[aM-x¨^(ti∣ti)]2whenx¨^(ti∣ti)>04-ππ[x¨^(ti∣ti)-a-M]2whenx¨^(ti∣ti)<0asmallpositiveconstantwhenx¨^(ti∣ti)=0.Wheni>K0, the parameter is updated by the following
(29)a^(ti)=x¨^(ti∣ti)-a-(ti)(30)ri(1)=ri-1(1)+1i[a^(ti)a^(ti-1)-ri-1(1)](31)ri(0)=ri-1(0)+1i[a^(ti)a^(ti)-ri-1(0)](32)βi=ri(1)ri(0)δaw2i=ri(0)-βiri(1)(33)δai2=δawi21-βi2αi=ln⁡βi-thi.The irregular sampling timeti-1, ti and the interval thi reflect in the time-varying parameters of the system, so we can conclude that the Kalman filter shown in (23)–(33) based on system (12) with system parameters (19)–(22) can obtain the same estimation performance as regular sampling Kalman filter.


## 2.1. Model Discretization

We can get the following by the differential equation (1):
(4)x(t)=eA(t-t0)x(t0)+∫t0teA(t-λ)Ua-(λ)dλ+∫t0teA(t-λ)Bw(λ)dλ.We can see that for any known integration interval [t0,t], x(t) can be gotten at any time t if the initial state x(t0), the parametersA, U, B, a-(t), and w(t) in [t0,t] are known.Consider the time interval fromti-1 to ti and assume
(5)a-(λ)=a-(ti-1),w(λ)=w(ti-1),λ∈[ti-1ti];
we can have
(6)∫ti-1tieA(ti-λ)Ua-(λ)dλ=∫ti-1tieA(ti-λ)Udλa-(ti-1)∫ti-1tieA(ti-λ)Bw(λ)dλ=∫ti-1tieA(ti-λ)Bdλw(ti-1).Setthi=ti-ti-1; we have the system matrix as Ad(ti-1)=eAthi, Ud(ti-1)=∫ti-1tieA(ti-λ)Udλ and the noise wd(ti-1)=∫ti-1tieA(ti-λ)Bdλw(ti-1) with the covariance Qd(ti-1)=E[wd(ti-1)wdT(ti-1)].Because the process matrixA in (2) is not a full-rank matrix, we cannot calculate the matrix exponential eAthi by the Lagrange-Hermite interpolation. Here, we use Laplace transform and have
(7)(sI-A)-1=[s-100s-100s+α]-1=[1s1s21s2(s+α)01s1s(s+α)001s+α].The matrix exponentialeAthi can be gotten by the inverse Laplace transform as
(8)Ad(ti-1)=[1thiαthi-1+e-αthiα2011-e-αthiα00e-αthi]
and by the similar approach, we can get the system parameter
(9)Ud(ti-1)=[1α(-thi+α·thi22+1-e-α·thiα)thi-1-e-α·thiα1-e-α·thi]
and the variance of wd(ti-1) as
(10)Qd(ti-1)=E[wd(ti-1)wdT(ti-1)]=2αδα2[q11q12q13q12q22q23q13q23q33]
with the parameters described as
(11)q11=12α5[2α3thi331-e-2α·thi+2α·thi+2α3thi33-2α2thi2-4α·thie-α·thi]q12=12α4[e-2α·thi+1-2e-α·thi+2α·thie-α·thi-2α·thi+α2thi2]q13=12α3[1-e-2α·thi-2α·thie-α·thi]q22=12α3[4e-α·thi-3-e-2α·thi+2α·thi]q23=12α2[e-2α·thi+1-2α·thi]q33=12α[1-e-2α·thi].Then, we get the discrete state-space model of the tracking system as(12)x(ti)=Ad(ti-1)x(ti-1)+Ud(ti-1)a-(ti-1)+wd(ti-1)z(ti)=H(ti)x(ti)+v(ti),
where x=[x,x˙,x¨]T is the state of the system to be estimated and whose initial mean and covariance are known as x0 and P0, wd(ti) and v(ti) are white noise with zero mean and independent of the initial state x0, z(ti) is the measurement vector, H(ti) is measurement matrices, and v(ti) is measurement noise with known variance R. Until now, the irregular sampling is turned to the varying-parameter system. We can see the same sampling interval is just a particular case of the random sampling problem. Therefore, the model of the randomly sampling tracking is a general one.


## 2.2. System Parameters Estimation

Here, we assume the maneuver frequencyα and the variance of the acceleration δa2 are not constant but variable and expressed as αi and δai2. From the processing model of (12), we have the discrete time equation of the acceleration as
(13)x¨(ti)=βix¨(ti-1)+(1-βi)a-(ti-1)+wa(ti-1),
where βi=e-αithi and wa(ti-1) is a zero mean white noise sequence with the variance
(14)δawi2=δai2(1-βi2).αi is the maneuver frequency at the sampling time ti. a-(ti-1) is the mean of one interval, so we have a-(ti)=a-(ti-1). Set a(ti)=x¨(ti)-a-(ti); then we can obtain
(15)a(ti)=βia(ti-1)+wa(ti-1).Consider the estimation of accelerationa^(ti) is a random process; we have
(16)a-(ti-1)=1i∑i=0i-1x¨^(ti),
where i is the number of data. For a first-order stationary Markov process (15), we have the statistics relation between the autocorrelation functions r(0), r(1) with the parameters βi and δawi2 by the Yule-Walker method [18]
(17)ri(0)=1i∑i=0i-1a^(ti)a^(ti)ri(1)=1i∑i=1i-1a^(ti)a^(ti-1)βi=ri(1)ri(0)δawi2=ri(0)-βiri(1).Next, we can getαi and δai2 by δai2=δawi2/(1-βi2),αi=ln⁡βi/-thi, and then get the system parameters Ad(ti-1), Ud(ti-1), and Qd(ti-1) in process function (12).


## 2.3. Algorithm Summary

Now we summarize the closed-loop estimation algorithm for the randomly sampled measurements as follows.(
1
)   Initialization(i=0). Consider
(18)x^(t0∣t0)=x0,P(t0∣t0)=P0,α0,δa02,a-(t0),r0(t0)=x¨0·x¨0,r0(t1)=x¨0.(
2
)  Recursion(i:=i+1)(a)
System update: setthi=ti-ti-1 and the system parameter as
(19)A^d(ti-1)=[1thiαithi-1+e-αithiαi2011-e-αithiαi00e-αithi](20)U^d(ti-1)=[1αi(-thi+αi·thi22+1-e-αi·thiαi)thi-1-e-αi·thiαi1-e-αi·thi]and the variance of the wd(ti-1) as
(21)Q^d(ti-1)=E[wd(ti-1)wdT(ti-1)]=2αiδαi2[q11q12q13q12q22q23q13q23q33]
with parameters described as
(22)q11=12αi5[2αi3thi331-e-2αi·thi+2αi·thi+2αi3thi33-2αi2thi2-4αi·thie-α·thi]q12=12αi4[e-2αi·thi+1-2e-αi·thi+2αi·thie-αi·thi-2αi·thi+αi2thi2]q13=12αi3[1-e-2αi·thi-2αi·thie-αi·thi]q22=12αi3[4e-αi·thi-3-e-2αi·thi+2αi·thi]q23=12αi2[e-2αi·thi+1-2αi·thi]q33=12αi[1-e-2αi·thi].(b)
State prediction: consider(23)x^(ti∣ti-1)=A^d(ti-1)x^(ti-1∣ti-1)+U^d(ti-1)a-(ti-1)P(ti∣ti-1)=A^d(ti-1)P(ti-1∣ti-1)A^dT(ti-1)+Q^d(ti-1).(c)
State update: consider(24)x^(ti∣ti)=x^(ti∣ti-1)+K(ti)[z(ti)-H(ti)x^(ti∣ti-1)](25)K(ti)=P(ti∣ti-1)HT(ti)×[H(ti)P(ti∣ti-1)HT(ti)+R(ti)]-1(26)P(ti∣ti)=[I-K(ti)H(ti)]P(ti∣ti-1).(d)
Parameter adaptation: the mean of the acceleration(27)a-(ti-1)=1i∑i=0i-1x¨^(ti∣ti).Wheni≤K0, the maneuver frequency αi is set to α0 and the covariance of the noise δai2 is gotten by the following:
(28)δαi2={4-ππ[aM-x¨^(ti∣ti)]2whenx¨^(ti∣ti)>04-ππ[x¨^(ti∣ti)-a-M]2whenx¨^(ti∣ti)<0asmallpositiveconstantwhenx¨^(ti∣ti)=0.Wheni>K0, the parameter is updated by the following
(29)a^(ti)=x¨^(ti∣ti)-a-(ti)(30)ri(1)=ri-1(1)+1i[a^(ti)a^(ti-1)-ri-1(1)](31)ri(0)=ri-1(0)+1i[a^(ti)a^(ti)-ri-1(0)](32)βi=ri(1)ri(0)δaw2i=ri(0)-βiri(1)(33)δai2=δawi21-βi2αi=ln⁡βi-thi.The irregular sampling timeti-1, ti and the interval thi reflect in the time-varying parameters of the system, so we can conclude that the Kalman filter shown in (23)–(33) based on system (12) with system parameters (19)–(22) can obtain the same estimation performance as regular sampling Kalman filter.


## 3. Proof of the Convergence

Based on the closed-loop estimation algorithm (18)–(33), we can see that the parameter used to estimate state is an estimated one and similarly the estimated states to calculate parameters αi and δai2 have estimation errors too. Therefore, it is important to guarantee the convergence of the estimation of the states and parameters.From (27), (29), and (33), we know if the estimation x¨^(ti∣ti) increased suddenly, a^(ti) will increase greatly because the mean changes less than x¨^(ti∣ti), and δawi2 becomes large too. Then, a very large positive δai2 will be obtained, and Q^d(ti) will also contain a large number of elements (here, we call it a big matrix). From the Riccati equation of Kalman filter
(34)P(ti+1∣ti)=A^d(ti){P(ti∣ti-1)-P(ti∣ti-1)HT(ti)SSSSS×[H(ti)P(ti∣ti-1)HT(ti)+R(ti)]-1SSS×H(ti)P(ti∣ti-1)HT}A^dT(ti)+Q^d(ti),
we find that P(ti+1∣ti) will be a big matrix if Q^d(ti) is a big one, and K(ti+1) will increase greatly. As a result, the estimation state x^(ti∣ti)=x^(ti∣ti-1)+K(ti)[z(ti)-H(ti)x^(ti∣ti-1)] will be a big matrix too. This trend results in positive feedback loops, which means x^(ti∣ti) will become larger and larger, and finally, divergence. We give the following theorem to guarantee the algorithm convergence.Theorem 1.
The estimationx^(ti+1∣ti+1) is bounded if the variance of the target acceleration δai2 has an upper bound; that is, there is a positive δ02 satisfying δai2≤δ02.Proof.
We firstly consider maneuvering frequencyαi. From (19), (21), and (22), we know if αi→0 and δai2≤δ02, the target has the constant acceleration maneuvering, and the system model is the constant acceleration model with the parameter as follows
(35)A^d(ti-1)⟶A-d(ti-1)=[1thithi2201thi001]Q^d(ti-1)⟶Q-d(ti-1)=δαi2[thi520thi48thi36thi48thi33thi22thi36thi22thi].
Ifαi→∞ and δai2≤δ02, we can get the system parameter matrix such as A^d(ti-1)→A--d(ti-1)=[1thi0010001] and Q^d(ti-1)→Q--d(ti-1)=δαi2[000000001]. Therefore, we can see that A^d(ti-1) and Q^d(ti-1) are the monotonic matrix with finite value elements.
Then, we consider the solution of Riccati equation (34) on the condition that the system parameter matrix has errors, such as A^d=A~d+ΔAd and Q^d=Q~d+ΔQd, where A~d and Q~d are the actual system parameters and ΔAd and ΔQdare the errors of the system parameter. Unlike the research about the uncertainty system, here we do not know the actual system matrices A~d and Q~d, but we can know the upper bound of the system parameters A^d(ti-1) and Q^d(ti-1), when δai2≤δ02, such as
(36)Aupper(ti-1)=A-d(ti-1)=[1thithi2201thi001],Qupper(ti-1)=δ02[thi520thi48thi36thi48thi33thi22thi36thi22thi].
The perturbed discrete algebraic Riccati equation is as follows:(37)P=A^d(ti)PA^dT(ti)-A^d(ti)PHT(ti)[H(ti)PHT(ti)+R(ti)]-1×H(ti)PA^dT(ti)+Q^d(ti).
We know that (37) is equal to
(38)P=A^d(ti)(P-1+HT(ti)R(ti)H(ti))-1×A^dT(ti)+Q^d(ti).
Then, for any vectors, we have
(39)sTPs=sT[A^d(ti)(P-1+HT(ti)R(ti)H(ti))-1×A^dT(ti)+Q^d(ti)(P-1+HT(ti)R(ti)H(ti))-1]s≤λ1(P-1+HT(ti)R(ti)H(ti))-1sTA^d(ti)A^dT(ti)s+sTQ^d(ti)s,
where λ1(X) is the maximum eigenvalue. By the relation of vector eigenvalue λm(X-1)=λM-m+1-1(X), where λ1(X)≥λ2(X)≥⋯≥λM(X), we have
(40)sTPs≤sTA^d(ti)A^dT(ti)sλM(P-1+HT(ti)R(ti)H(ti))+sTQ^d(ti)s.
That is,(41)P≤A^d(ti)A^dT(ti)λM(P-1+HT(ti)R(ti)H(ti))+Q^d(ti).
We have(42)λM(P-1+HT(ti)R(ti)H(ti))≥1λ1(P)+λM(HT(ti)R(ti)H(ti)).
Then, by (41) and (42), we have
(43)P≤A^d(ti)A^dT(ti)(1/λ1(P))+λM(HT(ti)R(ti)H(ti))+Q^d(ti)=λ1(P)A^d(ti)A^dT(ti)1+λ1(P)λM(HT(ti)R(ti)H(ti))+Q^d(ti)≤λ1(P)Aupper(ti)AupperT(ti)1+λ1(P)λM(HT(ti)R(ti)H(ti))+Qupper(ti).
Next, by the relation of Hermite matrix and its eigenvalue, we have(44)λ1(P)≤λ1(P)λ1(Aupper(ti)AupperT(ti))1+λ1(P)λM(HT(ti)R(ti)H(ti))+λ1(Qupper(ti)).
Then, we have(45)λ12(P)λM(HT(ti)R(ti)H(ti))+λ1(P)[1-λ1(Aupper(ti)AupperT(ti))sss-λM(HT(ti)R(ti)H(ti))λ1(Qupper(ti))]-λ1(Qupper(ti))≤0.
Assume thatλM(HT(ti)R(ti)H(ti))>0 and set
(46)1-λ1(Aupper(ti)AupperT(ti))-λM(HT(ti)R(ti)H(ti))λ1(Qupper(ti))=p1.
We have the following solution of (45):
(47)λ1(P)≤-p1+p12+4λM(HT(ti)R(ti)H(ti))λ1(Qupper(ti))2λM(HT(ti)R(ti)H(ti)).
Therefore, we can conclude that the maximum eigenvalue of estimation covarianceP has the upper bound shown as (47) if δai2≤δ02.
If one step predictive covariance is bounded, that is,|P(ti∣ti-1)|≤P0, then we know P(ti+1∣ti) must be bounded by (47) with the fact that |Q^d(ti)|≤Q0. And based on (25), we know K(ti+1) must be a bounded matrix and x^(ti+1∣ti+1) must be bounded too.


## 4. Simulations and Experiments

### 4.1. The Estimation by Different Extraction Rate and Irregular Rate

The method here is applied to a two-dimensional planar video tracking. Here, as a tracking problem, we just use the simple background and one target. The video gotten by the Image Capture Test Bed is shown in Figure1.Figure 1
The video with simple background and one target.We control the car maneuvering on the test bed and catch the images of target movement by a stationary camera. For every image of the video, the target is extracted based on the color and then we get the measurement data of maneuvering target on the Image Capture Test Bed like Figure2.Figure 2
The measurement of maneuvering target got from the video.We know that the camera catches the image under the same interval, and that will produce large amounts of image data. If we can use some of images in the video for tracking, the image storage and computation cost will greatly reduce. But “using some of images” means that the measurements no longer have the same sampling interval. Here, define the Extraction Rate as(48)EXrate=extractednumberofimagesfromthevideototalnumberofimagesinthevideo×100%
to describe the image compression rate. And define the Irregular Rate to measure the sampling interval as
(49)IRrate=∑i=1N|thi-∑i=1Nthi|N.The state for the target in the 2D space isx(k)=[x(k)x˙(k)x¨(k)y(k)y˙(k)y¨(k)]. The initial state estimate x0 and covariance P0 are assumed to be x0=[x(0)00y(0)00]T and P0=diag⁡(10,10,10,10,10,10).We extract 243 images from a video with 491 images where EXrate = 49.49% and IRrate = 0.1043 and by the algorithm developed with the initial parametersα0=1/20, δa02=10, a-0=0, αM=3, K0=3, we get the estimation of trajectory with estimation covariance 10.0881 along the horizontal axis and 8.1660 along the vertical axis, shown in Figure 3. The estimation trajectories of horizontal and longitudinal axis is shown in Figure 4 and the estimation error are shown in Figure 5.Figure 3
The real trajectory and the estimation trajectory.The estimations of horizontal and longitudinal axis.
(a)
(b)The location estimation errors.
(a)
(b)To illustrate how the irregular rate affects estimation performance, the algorithm is used to estimate the target trajectory under different Irregular Rate (shown in Table1) with the same Extraction Rate, EXrate=49.86%. The RMSE position is defined as RMSE2D=RMSEH2+RMSEL2, where RMSEH and RMSEL are the root-mean square errors (RMSE) of position for horizontal and longitudinal axis, respectively. The relation between RMSE2D and IRrate is shown in Figure 6. We can see that the Irregular Rate affects the estimation performance very little. The Irregular Rate changes 21 times almost from 0.0088 in Case 1 to 0.1928 in Case 10, but RMSE2D is about 14 for all IRrate. We can conclude that IRrate does not affect the tracking performance, when with the same EXrate.Table 1
The different irregular rate for 10 cases.
Case  1
Case  2
Case  3
Case  4
Case  5
Case  6
Case  7
Case  8
Case  9
Case  10
0.01
0.06
0.08
0.09
0.10
0.13
0.14
0.14
0.18
0.19Figure 6
The relation betweenRMSE2D and IRrate.


### 4.2. The Performance with Different Models

Next, we compare the model developed here with other dynamics model, such as CV model [12], CA model [12], Singer model [13], current model [8], and IMM [14]. We set the process noise covariance as Q=1 for the CV and CA model and σw2=1 and α=1/20 for Singer model. Because the current model is very sensitive to the priori parameters, we give several system parameters such as α=1/30, αmax⁡=3 (current model I), α=1/20, αmax⁡=30 (current model II), and α=1/20, αmax⁡=3 (current model III). After 100 Monte Carlo simulation runs, RMSE2D are calculated. For different trajectory with different EXrate and IRrate, the estimation results are shown in Figures 7(a)–7(f), where in order to show clearly, we use the black “O” to describe the actual trajectory at the sampling time in Figures 7(e) and 7(f).The tracking results for videos.
(a)
(b)
(c)
(d)
(e)
(f)Table2 and Figure 8 show RMSE2D under the different IRrate and EXrate. We can see that the model here can get the better estimation performance than CV, CA, Singer model, current model, and IMM for almost all EXrate and IRrate. We also note that the current model needs the right parameter, or else the performance will become worse.Table 2
The estimation covariance and irregular rate in Figure7.
The video (frames)
EXrate
IRrate
RMSE2D
Total
Extracted
The modelhere
CVmodel
CAmodel
Singermodel
Current model Iα=1/30 αmax⁡=3
Current model IIα=1/20 αmax⁡=30
Current model IIIα=1/20 αmax⁡=3
IMM
(a)
381
339
88.98%
0.0261
10.72
29.54
16.99
33.76
90.84
11.56
10.14
21.88
(b)
241
188
78.01%
0.0514
12.76
54.16
21.03
51.37
121.0
12.01
13.01
24.51
(c)
221
150
67.87%
0.0711
14.39
66.22
37.25
67.68
144.7
17.30
17.64
29.35
(d)
231
121
52.38%
0.1087
16.69
39.14
20.10
37.94
94.85
13.84
15.82
30.90
(e)
371
107
28.84%
0.2124
23.34
31.01
23.88
27.12
49.61
20.73
20.10
37.43
(f)
211
28
13.27%
0.4635
40.93
69.85
25.92
49.74
73.78
23.14
26.06
38.14Figure 8
RMS
E
2
D under different EXrate and IRrate.We note that in Figure7(f), the tracking error of the developed model is larger than current models II, III, and IMM, even CA. We find that there is a big estimation error at 5th second. The reason is that there are not enough data gotten to update the parameter at K0=4. Therefore, the estimation error is bigger. But we also note that the estimation error declined quickly, so the developed model has a strong advantage for the long trajectory tracking comparing the other models.Another fact we also noticed is that thoughIRrate almost does not affect the tracking performance, it is obvious that low EXrate can decline the tracking performance. This is because the lower EXrate means less measured data gotten and less useful information that can be provided; therefore, the estimate is more inaccurate.As to the sampling intervalthi, the lower EXrate means larger thi. If the sampling interval thi is large enough to break Shannon Sampling Theorem, the estimation performance will decline.


### 4.3. The Estimation of Video Target

At last, we use the developed method to track a target in real scene. In order to decrease the calculation cost, we select some frames from the video according to the characteristics of the movement. That is, if we find that the target is stationary or moves slowly, then we discard these frames. We use a threshold to test whether a target makes a big maneuver or not. Obviously, a large threshold can make the calculation cost lower, but lowerEXrate will make the performance decrease too.So the threshold should be carefully selected to balance the calculation cost and performance. Here, we select 95 frames from 245 frames;EXrate and IRrate are 38.77% and 0.1367, respectively. Figure 9 gives the tracking results of number 1, 27, 40, 65, 74, 97, 128, 129, 158, 181, 189, and 226 frames in the video. The estimation of target is marked by “black” dot. The estimation covariance of RMSE2D as 1.034 mm is obtained (the tracking area is 300 * 300 mm2).Figure 9
The tracking results in number 1, 27, 40, 65, 74, 97, 128, 129, 158, 181, 189, and 226 frames.


## 4.1. The Estimation by Different Extraction Rate and Irregular Rate

The method here is applied to a two-dimensional planar video tracking. Here, as a tracking problem, we just use the simple background and one target. The video gotten by the Image Capture Test Bed is shown in Figure1.Figure 1
The video with simple background and one target.We control the car maneuvering on the test bed and catch the images of target movement by a stationary camera. For every image of the video, the target is extracted based on the color and then we get the measurement data of maneuvering target on the Image Capture Test Bed like Figure2.Figure 2
The measurement of maneuvering target got from the video.We know that the camera catches the image under the same interval, and that will produce large amounts of image data. If we can use some of images in the video for tracking, the image storage and computation cost will greatly reduce. But “using some of images” means that the measurements no longer have the same sampling interval. Here, define the Extraction Rate as(48)EXrate=extractednumberofimagesfromthevideototalnumberofimagesinthevideo×100%
to describe the image compression rate. And define the Irregular Rate to measure the sampling interval as
(49)IRrate=∑i=1N|thi-∑i=1Nthi|N.The state for the target in the 2D space isx(k)=[x(k)x˙(k)x¨(k)y(k)y˙(k)y¨(k)]. The initial state estimate x0 and covariance P0 are assumed to be x0=[x(0)00y(0)00]T and P0=diag⁡(10,10,10,10,10,10).We extract 243 images from a video with 491 images where EXrate = 49.49% and IRrate = 0.1043 and by the algorithm developed with the initial parametersα0=1/20, δa02=10, a-0=0, αM=3, K0=3, we get the estimation of trajectory with estimation covariance 10.0881 along the horizontal axis and 8.1660 along the vertical axis, shown in Figure 3. The estimation trajectories of horizontal and longitudinal axis is shown in Figure 4 and the estimation error are shown in Figure 5.Figure 3
The real trajectory and the estimation trajectory.The estimations of horizontal and longitudinal axis.
(a)
(b)The location estimation errors.
(a)
(b)To illustrate how the irregular rate affects estimation performance, the algorithm is used to estimate the target trajectory under different Irregular Rate (shown in Table1) with the same Extraction Rate, EXrate=49.86%. The RMSE position is defined as RMSE2D=RMSEH2+RMSEL2, where RMSEH and RMSEL are the root-mean square errors (RMSE) of position for horizontal and longitudinal axis, respectively. The relation between RMSE2D and IRrate is shown in Figure 6. We can see that the Irregular Rate affects the estimation performance very little. The Irregular Rate changes 21 times almost from 0.0088 in Case 1 to 0.1928 in Case 10, but RMSE2D is about 14 for all IRrate. We can conclude that IRrate does not affect the tracking performance, when with the same EXrate.Table 1
The different irregular rate for 10 cases.
Case  1
Case  2
Case  3
Case  4
Case  5
Case  6
Case  7
Case  8
Case  9
Case  10
0.01
0.06
0.08
0.09
0.10
0.13
0.14
0.14
0.18
0.19Figure 6
The relation betweenRMSE2D and IRrate.


## 4.2. The Performance with Different Models

Next, we compare the model developed here with other dynamics model, such as CV model [12], CA model [12], Singer model [13], current model [8], and IMM [14]. We set the process noise covariance as Q=1 for the CV and CA model and σw2=1 and α=1/20 for Singer model. Because the current model is very sensitive to the priori parameters, we give several system parameters such as α=1/30, αmax⁡=3 (current model I), α=1/20, αmax⁡=30 (current model II), and α=1/20, αmax⁡=3 (current model III). After 100 Monte Carlo simulation runs, RMSE2D are calculated. For different trajectory with different EXrate and IRrate, the estimation results are shown in Figures 7(a)–7(f), where in order to show clearly, we use the black “O” to describe the actual trajectory at the sampling time in Figures 7(e) and 7(f).The tracking results for videos.
(a)
(b)
(c)
(d)
(e)
(f)Table2 and Figure 8 show RMSE2D under the different IRrate and EXrate. We can see that the model here can get the better estimation performance than CV, CA, Singer model, current model, and IMM for almost all EXrate and IRrate. We also note that the current model needs the right parameter, or else the performance will become worse.Table 2
The estimation covariance and irregular rate in Figure7.
The video (frames)
EXrate
IRrate
RMSE2D
Total
Extracted
The modelhere
CVmodel
CAmodel
Singermodel
Current model Iα=1/30 αmax⁡=3
Current model IIα=1/20 αmax⁡=30
Current model IIIα=1/20 αmax⁡=3
IMM
(a)
381
339
88.98%
0.0261
10.72
29.54
16.99
33.76
90.84
11.56
10.14
21.88
(b)
241
188
78.01%
0.0514
12.76
54.16
21.03
51.37
121.0
12.01
13.01
24.51
(c)
221
150
67.87%
0.0711
14.39
66.22
37.25
67.68
144.7
17.30
17.64
29.35
(d)
231
121
52.38%
0.1087
16.69
39.14
20.10
37.94
94.85
13.84
15.82
30.90
(e)
371
107
28.84%
0.2124
23.34
31.01
23.88
27.12
49.61
20.73
20.10
37.43
(f)
211
28
13.27%
0.4635
40.93
69.85
25.92
49.74
73.78
23.14
26.06
38.14Figure 8
RMS
E
2
D under different EXrate and IRrate.We note that in Figure7(f), the tracking error of the developed model is larger than current models II, III, and IMM, even CA. We find that there is a big estimation error at 5th second. The reason is that there are not enough data gotten to update the parameter at K0=4. Therefore, the estimation error is bigger. But we also note that the estimation error declined quickly, so the developed model has a strong advantage for the long trajectory tracking comparing the other models.Another fact we also noticed is that thoughIRrate almost does not affect the tracking performance, it is obvious that low EXrate can decline the tracking performance. This is because the lower EXrate means less measured data gotten and less useful information that can be provided; therefore, the estimate is more inaccurate.As to the sampling intervalthi, the lower EXrate means larger thi. If the sampling interval thi is large enough to break Shannon Sampling Theorem, the estimation performance will decline.


## 4.3. The Estimation of Video Target

At last, we use the developed method to track a target in real scene. In order to decrease the calculation cost, we select some frames from the video according to the characteristics of the movement. That is, if we find that the target is stationary or moves slowly, then we discard these frames. We use a threshold to test whether a target makes a big maneuver or not. Obviously, a large threshold can make the calculation cost lower, but lowerEXrate will make the performance decrease too.So the threshold should be carefully selected to balance the calculation cost and performance. Here, we select 95 frames from 245 frames;EXrate and IRrate are 38.77% and 0.1367, respectively. Figure 9 gives the tracking results of number 1, 27, 40, 65, 74, 97, 128, 129, 158, 181, 189, and 226 frames in the video. The estimation of target is marked by “black” dot. The estimation covariance of RMSE2D as 1.034 mm is obtained (the tracking area is 300 * 300 mm2).Figure 9
The tracking results in number 1, 27, 40, 65, 74, 97, 128, 129, 158, 181, 189, and 226 frames.


## 5. Conclusions

The main contribution of this paper is to model the real-time system dynamics at the random sampling points.(1) By calculating the matrix exponential with inverse Laplace transform, the irregular sampling interval is transformed to time-varying parameters matrix of the system. (2) Based on the statistics relation between the autocorrelation function and the covariance of Markov random processing, the system model with online parameter is developed. The proof and the experimental results show that the developed method can get good tracking performance.As an example, the developed method is used for the video tracking problem. According to the motion characteristics of the target, some frames are selected for the tracking purpose. The tracking results show that good tracking performance is obtained by a smaller amount of calculation.



---

*Source: 315908-2014-02-26.xml*