Journal of Systems Engineering and Electronics Vol. 19, No. 6, 2008, pp. 1127–1132
Posterior Cramer-Rao lower bounds for multitarget bearings-only tracking Guo Lei1,2 , Tang Bin1 , Liu Gang2 & Xiao Fei1 1. School of Physical Electronics, Univ. of Electronic Science and Technology of China, Chengdu 610054, P. R. China; 2. National Key Laboratory of Information Integrated Control, Chengdu 610036, P. R. China (Received September 27, 2007)
Abstract: Usually, only the Cramer-Rao lower bound (CRLB) of single target is taken into consideration in the state estimate of passive tracking systems. As for the case of multitarget, there are few works done due to its complexity. The recursion formula of the posterior Cramer-Rao lower bound (PCRLB) in multitarget bearings-only tracking with the three kinds of data association is presented. Meanwhile, computer simulation is carried out for data association. The final result shows that the accuracy probability of data association has an important impact on the PCRLB.
Keywords: multiple target tracking, bearings-only tracking, posterior Cramer-Rao lower bounder, data association.
1. Introduction As one of the modern electronic reconnaissance technologies, passive tracking uses the electromagnetic wave radiated from target to orient and track target. Multitarget tracking using bearings-only measurement is the most widely used method in practice. It usually includes two parts, i.e., state estimation and data correlation. The Kalman filtering algorithm or extent Kalman filtering algorithm are mainly used in state estimation. At present, there are some data association algorithms available, for example, the joint probabilistic data association (JPDA), multiple hypothesis tracking (MHT), and probabilistic multiple hypothesis tracking (PMHT). Target state estimate accuracy is an important parameter to evaluate the performance of tracking algorithms. Usually, the Cramer-Rao lower bound (CRLB) is used as the quantification index to analyze target state estimate accuracy, which provides state estimation accuracy lower bound for the target with unknown state. Up to now, the CRLB has been used to analyze estimate performance in some articles. For example, in Refs. [1–2], the CRLB under the condition of uncertain measurement is presented and those main parameters affecting the CRLB are also analyzed, but
their discussion is only limited in the case of single target tracking. As for the case of multitargets, the obtainment of the CRLB is much more difficult. This is because both the filtering algorithm and the data correlation probability should be considered together in order to obtain the lower bound of target state estimate. Unfortunately, both are correlated with each other, which make analysis more difficult. The posterior CRLB is defined as the inversion of posterior FIM, which is superior to the CRLB[3] . In this article, under reasonable condition, a recursive Riccatic-like formulation of PCRLB of the MTT using bearings-only measurement is derived. Furthermore, we can find that the data association algorithm has an important impact on target state estimate accuracy with other conditions kept same. These conclusions can help to optimize only the design of target track systems using bearings-only measurement.
2. General model Let M be the number of the targets being tracked. Suppose that it has already been known and fixed here. The index i designates one among the M targets and is always used as superscript. At the time t, the state vector is Xt = (Xt1 , . . . , XtM ), one element of which is written as
1128
Guo Lei, Tang Bin, Liu Gang & Xiao Fei i , Vti ), Xti = Fti (Xt−1
i = 1, . . . , M
(1)
where Fti is the state transition matrix. The noise Vti is a Gaussion zero-mean vector with covariance matrix ΣV . The observation vector collected at time t is denoted as yt = (yt , . . . , ytmt ). The index j is used as superscript to refer to one of the mt measurement. The vector yt is composed of detection measurement and clutter measurement. The false alarms are assumed to be uniformly distributed in the observation area[1]. Their numbers are assumed to be characterized by the Poisson distribution with the density μf , where μf = λV, V is the volume of the observation area, and λ the average number of false alarms per unit volume. Due to the fact that as the origin of each measurement is unknown, we have to introduce the vector Kt to describe the associations between the measurements and the targets. Each component Ktj is a random variable whose value is among {0,. . . ,M}. If i = 0, Ktj = i indicates that ytj is associated with the ith target. Otherwise, it is a false alarm if i = 0. In the first case, ytj is a realization of the stochastic process (2) Ytj = Hti (Xti , Wtj ) if Ktj = i Hti is the measurement transition matrix. The measurement noise Wtj is supposed to be white noise and independent. If only the bearing measurement is used, a measurement produced by the ith target is generated according to i yt − ytobs Ytj = arctan + Wtj (3) xit − xobs t obs are the Cartesian coordinates The variables xobs t , yt of the observer, which are known. The variable (Ktj )j=1,...,mt is supposed to be uniformly distributed. Their common probability (πti )i=1,...,M is defined as
πti = P (Ktj = i),
j = 1, . . . , mt
i=1
In order to solve the data association, some assumptions are made in the following[3] . A1 One measurement might originate from one target or the clutter. A2 One target can produce zero or one measurement at one time. A3 One target can produce zero or several measurement at one time.
3. PCRLB The CRLB or posterior CRLB provides a means of quantifying the achievable accuracy of target state esˆ timation. Let X(Z) be an unbiased estimator of a parameter vector X, which is based on the measurement vector Z. The CRLB for the error covariance matrix is defined to be the inverse of the finish information matrix (FIM), i.e. ˆ − X][X ˆ − X]T } J −1 C = E{[X
(5)
(6)
The inequality (6) means that the difference C − J is a positive semidefinite matrix. The matrix J is given by 2 ∂ ln p(Z, X) J(i, j) = E − (7) ∂X(i)∂X(j) −1
where p(Z, X) is the joint probability density function (PDF) of (Z, X), X(i) denotes the ith component of X. Note that the expectation E{.} is operated on both Z and X. In Ref. [4], such approach is studied and developed. It presents a recursive formulae to derive a sequence of posterior FIMs JXt for the unbiased estimation of the target state, i.e.
(4)
where πti is the prior probability of an arbitrary measurement, which is associated with target i. The πt vector is considered as a realization of the stochastic vector Π t = (Πt0 , Πt1 , . . . , ΠtM ). The prior distribution of the variable Πt is presented in the following p(Πt ) = p(Πt0 )p(Πt1 , . . . , ΠtM |Πt0 )
where p(Πt1 , . . . , ΠtM |Πt0 ) is uniform on the hyperM plane defined by Πti = 1 − Πt0 .
22 21 11 −1 12 JXt+1 = DX − DX (JXt + DX ) DXt t t
(8)
11 t DX = E[−ΔX Xt log p(Xt+1 |Xt )] t
(9)
where 12 = E[−ΔXt+1 log p(Xt+1 |Xt )] DX t t X
21 12 T t = E[−ΔX DX Xt+1 log p(Xt+1 |Xt )] = [DXt ] t 22 = E[−ΔXt+1 log p(Xt+1 |Xt )]+ DX t t X
(10) (11)
Posterior Cramer-Rao lower bounds for multitarget bearings-only tracking X
E[−ΔXt+1 log p(Yt+1 |Xt+1 )] t+1
(12)
where the operators and Δ denote the first partial derivative and the second partial derivative, respectively
X
∂ ∂ = ,..., ∂x1 ∂xn
T
,
ΔYX =
X
T Y
obtain log p(Yt = ytmt |Xt = xt , Kt = kt ) =
mt
4. PCRLB for multiple targets The PCRLB in Ref. [4] can be extended and used in the case of multiple targets filtering. In the following, the measurement uncertainty and the existence of multiple targets are simultaneously taken into consideration. Then, we can find that all the recursive formulas are known except the second term of the right-hand side of Eq. (11). This term is expressed as the second partial derivatives of the logarithm of the measurement’s posterior PDF. Therefore, how to obtain this term is a key problem to solve. In order to obtain this equation, we need to know the association probability between measurements and tracks. Here, three kinds of cases are considered in the following for which three lower bounds are computed, respectively: The notation B1 represents the case in which the PCRLB is computed under the assumption that the associations are known. B1 is an ideal condition in order to compare other case. The notation B2 represents the case in which the PCRLB is computed under the assumptions A1 and A2. B2 is used in JPDA and MHT algorithm. The notation B3 represents the case in which the PCRLB is computed under the assumptions A1 and A3. B3 is used in PHMT algorithm. In the next part of the article, those PCRLBs under B1, B2, or B3 will be computed, respectively. PCRLB B1
(14) The gradient of the log-likelihood function is not zero only if there exists j i . In this case i
∇Xti log p(ytmt |xt , kt )
=
∇Xti p(ytj |xit )
(15)
i
p(ytj |xit )
We can obtain for all i = 1, . . . , M Xi
JX it (p(ytmt |xt , kt )) = t
i
EXt EY ji |X And J
2
Xti
1 Xti
i
∇Xti p(ytj |xit )(∇Xti p(ytj |xit ))T i
p(ytj |xit )2
t
t
p(p(yt |xt , kt )) = 0 if i1 = i2
(16)
(17)
In the general case, an observation model with an additive Gaussian noise is defined as Ref. [5] 1 p(ytj |xit ) =(π ny det Σ )−1/2 exp − (ytj − H(xit ))T 2 Σ x−1 (ytj − H(xit )) (18) where ny represents the measurement dimension and Y ∈ Rny . Then, we have 1
1
∇X i1 p(ytj |xit ) =p(ytj |xit )∇X i1 · t
t
1
1
H T (xit )Σ −1 (ytj − H(xit )) (19) Finally, we obtain the PCRLB for the case of Gaussian noise Xi
JX it (p(Yt |Xt )) = EXti ∇Xti H T (xit )Σ−1 (∇Xti H T (xit ))T t (20) 4.2
PCRLB B2
B2 represents the case with the assumptions A1 and A2. We can write
log
log p(Yt = ytmt |Xt = xt ) p(yt =
A1−A2
=
(yt1 , . . . , ytmt )|xt , kt )p(kt )
=
kt
log Suppose that the association vector is known. We can
kj
log p(ytj |xt t )
j=1
(13)
−1 provides a lower bound of the meanThe matrix JX t+1 square error for estimating Xt+1 . From Ref. [4], we 11 analytically obtain the following equalities: DX = t −1 iT −1 iT −1 i 12 21 diag{F ΣV F }, DXt = DXt = diag{ΣV F ΣV }. Xt+1 In Eq. (12), JXt+1 (p(Xt+1 |Xt )) = diag{ΣV−1 }, where Jαβ (p) represents E[−Δβα log(p)].
4.1
1129
mt kt j=1
p(ytj |xt , kt )p(kt )
(21)
1130
Guo Lei, Tang Bin, Liu Gang & Xiao Fei We obtain for i1 , i2 = M
The gradient of the log-likelihood is ∇Xti log p(yt |xt ) =
∇Xti
mt
J
1 Xti
p(ytj |xt , kt )p(kt )
j=1
kt
(22) p(yt |xt ) The association of one measurement with the ith target is denoted as kt ⊃ i. Under the assumption A2, there exists at most one such measurement, which is denoted as j i , then
(23) p(yt |xt ) Finally, we can obtain the PCRLB for the case of the Gaussian noise J
2
Xti
1 Xti
(p(Yt |Xt )) = 1
EXt ∇X i1 H T (xit )
−1
t
EYt |Xt ·
i p(yt |xt , kt )p(kt )(ytj
−
·
p(yt |xt )2
i 2 p(yt |xt , kt )p(kt )(ytj − H(xit ))T · −1
2
(∇X i2 H T (xit ))T
(24)
t
4.3
PCRLB B3
To our knowledge, the algorithm using the assumption A3 need a joint estimation of Xt and πt . The estimation quality of joint strongly affects state. We consider the PCRLB for the estimation of the join vec
1:M−1 1:M tor ( t , Xt ). Let us define Φt = ( t , Xt ). To 22 evaluate the third term of DΦ , we can write t log p(Yt = yt |Φt = φt ) mt j=1
log
A1−A3
=
log
mt
p(ytj |φt ) =
j=1
m0t V
− πt0 p(ytj |xM t )+
M−1 j i j M j M i (p(yt |xt ) − p(yt |xt ))πt + p(yt |xt )
(25)
For i = M , the gradient of the log-likelihood is mt ∇X i p(ytj |xit ) t
j=1
1
(p(Yt |Φt )) =
2
1
t
1
2
p(ytj |φt )2 −
EY j |Φt · t
j=1
p(ytj |xit )p(ytj |xit ) 2 H(xit ))T ]
mt
·Σ
−1
1
(ytj − H(xjt ))·
(∇X i2 H t
T
2 (xit ))T
(28)
p(ytj |φt )
In the computer simulation, we deal with classical bearings-only experiments with three targets. For simplicity, a nearly constant-velocity model is supposed. Note that the state vector Xti represents the coordinate and the velocities in the x − y plane: Xti = (xit , yti , vxit , vyti ) with I= 1,2,3. For each target, the discretized state equation associated with time period is ⎛ ⎛ 1 ⎞ ⎞ I 0 I I 2 2 2 i ⎠ Xti + ⎝ 2 ⎠ Vti (29) =⎝ Xt+1 0 I2 0 I2 where I2 is the 2-D identity matrix and Vti is the Gaussian zero-mean vector with the covariance ma2 2 , σvy ]. A set of mt measuretrix ΣV = diag[σx2 , σy2 , σvx ments is available at discrete times. The initial coordinates of the targets and of the observer are (in meter and meter/second, respectively) X01 = (200, 1 500, 1, −0.5)T, X02 = (0, 0, 1, 0)T
i=1
∇Xti log p(yt |φt ) = πti
1 Xti
5. Computer simulation
1 H(xit )
kt ⊃i2
2
Xti
EΦt πti πti ∇X i1 H T (xit )Σ−1
(ytj
kt ⊃i
Then, finally, we can obtain the PCRLB for the case of the Gaussian noise
i
kt ⊃i j=j i
(27)
J
p(ytj |xt , kt )p(kt )∇Xti p(ytj |xit )
(p(Yt |Φt )) =
1 2 mt ∇X i1 p(ytj |xit )(∇X i2 p(ytj |xit )T i1 i2 t t EYti |Φt EΦt πt πt p(ytj |φt )2 j=1
∇Xti log p(Yt |Xt ) =
2
Xti
(26)
X03 = (−200, −1 500, 1, 0.5)T X0obs = (200, −3 200, 1.2, 0.5)T The observer is following a leg-by-leg trajectory. Its velocity vector is kept constant on each leg except at the following instants
Posterior Cramer-Rao lower bounds for multitarget bearings-only tracking 1131 ⎞ ⎛ ⎞ ⎛ ⎞ ⎛ ⎞ ⎛ and λ=1. The final results are plotted in the Fig. 2. −0.6 vxobs 2.0 vxobs ⎠ ⎝ 400,800 ⎠ = ⎝ ⎠ ⎝ 200,600,900 ⎠ = ⎝ From Fig.2, it can be found that B2>B3>B1 irobs obs vy200,600,900 vy400,800 0.3 0.3 respective of whatever trace and volume are. Meanwhile, the gap between B3 and B1 is larger than that The trajectories of the three objects and the observer between B2 and B3. First, it means that the optimal are plotted in Fig. 1(a). performance of the algorithm using the assumptions We have computed the trace and the volume of the A1 and A2 are worse than that of the algorithm using three bounds with the standard dynamic noise σx = −1 the assumptions A1 and A3. It should be pointed out σy =0.000 5 ms , the detection probability Pd =0.9
Fig. 1
Fig. 2
Trajectories and measurement of three targets
Trace and volume of the three PCRLB: B1 (dashed), B2 (solid), B3 (dashdotted)
target to the case of multitargets filtering problem is studied. The state accuracy lower bounds for the bearings-only are derived according to the three association assumptions between the measurements and the targets. Using theory and computer simulation means that the data association algorithm has an important impact on the target state estimate accuracy. This conclusion can help to analyze and improve of the target track accuracy.
that the JPDA algorithm and the MHT algorithm use the assumptions A1 and A2, whereas the PMHT algorithm uses the assumptions A1 and A3. Therefore, the performance of the PHTM algorithm is better than those of the JPDA algorithm and the MHT algorithm if other conditions are kept same. Second, the performance of the algorithm with the association known is far better than those of the two former cases. Therefore, the data association algorithm has an important impact on the target state estimate accuracy.
References
6. Conclusions
[1] Ristic B, Zollo S, Arulampalam S. Performance bounds
The extension of the PCRLB from the case of single-
for maneuvring target tracking using asynchronous multi-
1132
Guo Lei, Tang Bin, Liu Gang & Xiao Fei
platform angle-only measurement. Proceedings of 4th In-
[5] Bar-shalom Y, Li X. Multitarget-mulitsensor tracking:
ternational Conference on Information Fusion, Montreal,
principles and techniques storrs.
Quebec, 2001: 695–703.
1995, 62–85.
CT: YBS Publishing:
[2] Niu R, Willett P K, Bra-Shalom Y. Matrix CRLB scaling due to measurements of uncertain. Origin IEEE Trans. on SP, 2001, 49(7): 1325–1335. [3] Bar-shalom Y, Fortmann T E. Tracking and data association. New York: Academic Press, 1988: 126-129. [4] Tichavsky P, Muracchik C H, Nehorai A. Posterior CramerRao bounds for discrete-time nonlinear filtering. Trans. on SP, 1998, 46(5): 1386–1396.
IEEE
Guo Lei was born in 1971. He is a lecturer in University of Electronic Science and Technology of China (UESTC). He is a Ph. D. candidate. His research interests include multisensor data fusion, passive tracking, sensors management, etc. E-mail: leiguo@uestc. edu.cn