Maneuver Detection Based on Information Geometry and Geodesic Shooting

Maneuver Detection Based on Information Geometry and Geodesic Shooting Marion Pilté Frédéric Barbaresco THALES AIR SYSTEMS, Limours, France THALES AIR SYSTEMS, Limours, France Email: [email protected] Email: [email protected] Abstract—This paper explores a new way of monitoring track- In the following, two models will be considered to describe ing quality based on information geometry methods. The method the target’s dynamics: proposed takes into account all parameters of the movement of • Rectilinear movement, with constant velocity the target with the use of an appropriate distance which reflects the mean and the covariance of the distribution we obtain with • Circular movement, with constant angular velocity the Kalman filter. To achieve that, we develop the distance in The problem of tracking targets can be addressed in two the manifold of multivariate gaussians, compute the Fisher-Rao ways: tracking based on one or on several models representing distance, and compare it with bounds. For more details, see [1]. Index Terms—Radar, Information geometry, Fisher metric, different behaviours of a target. The first solution is only Kalman filter. appropriate if we are able to determine in advance a best model followed by the target, whereas the second one uses a set of I. I NTRODUCTION possible models for the target at an instant k. We will use the second method, based on the first. The development of fixed antenna radars enables to adapt the update rate of the radar. It is then necessary to decide when A. Linear Kalman filter to generate more pulses, like during a maneuver. A performant The discrete Kalman filter (KF), as described in [3], gives an tracking quality monitor is thus needed. estimation and a prediction of the state of the target. A cycle The most used monitors are CUSUM (Cumulative Sum of the discrete Kalman Filter is given in [3]. The parameters based detector), GLR (Generalized Likelihood Ratio), de- for the KF are: the model m, the covariance matrices of the scribed in [2]. But these detectors measure only filtered noises of the system and the measure (Q and R respectively), normalized innovation. and the initial state x̂0|0 and the state covariance P0|0 . The problem that must be solved is the following: the target This is the elmentary filter, now we describe a filter which can follow different kinematic models which are predefined. copes with different models at the same time. We have to detect a degradation in tracking quality. The monitoring we propose uses the outputs of a Kalman filter: the B. Multiple model filter : The Interacting Multiple Model estimated and predicted state and the covariance matrices have The simple Kalman Filter supposes that the target maneuver to be computed before applying the detection. We proceed can be represented by a best model. Multiple model algorithms in two steps: first, the filtering step for which a multiple are able to track targets whose movement has an uncertainty, model filter is used, then the monitoring step for which we because they deal with more than one model within each consider the space of multivariate Gaussians, which describes iteration. The idea is first to define a set of models, then to run the space of estimates and predictions at the output of the a bank of elementary filters, each based on a unique model, Kalman filter whose outputs (state vector and covariance) can and at last generate global estimations of a process based on be seen as parameters of a distribution. We shall need a notion the results of these elementary filters. Let us call M the set of distance, which is the backbone of our new method. of the N models : M = {m1 , ..., mN }. Section II recalls the basics of the simple and multiple We use the Interacting Multiple Model (IMM) algorithm, model Kalman filter. Section III focuses on information ge- the most used in practice. The parameters of the IMM are: ometry tools, finally, an application is presented in section IV. the structure of the set of models M, the covariances and the noises of the system and the measure (Q and R respectively), II. R EMINDER ON THE K ALMAN FILTER the initial state x̂0|0 and the covariance of the state P0|0 , and The target is considered as a ponctual mass. We use carte- the jump structure and the transition probability πji between sian coordinates. The state vector equation is: the models of the selected set. The cycle of the IMM is presented in [3]. xk+1 = Fk xk + Gk uk + Ek ωk (1) After having performed filtering, we explain a new method The measurement equation is linear: to obtain a track quality estimator: we use manifolds in order to compute distances between predicted and estimated states zk = Hk xk + vk (2) coupled with their covariance matrices.  III. I NFORMATION GEOMETRY ( µ√1i2 , σ1i )−( µ√2i2 , −σ2i ), ρ2 = ( µ√1i2 , σ1i )−( µ√2i2 , σ2i ). For θ1 = To construct an appropriate distance, we present a metric (µ11 , σ11 , ..., µ1n , σ1n ) and θ1 = (µ21 , σ21 , ..., µ2n , σ2n ), the in the general case of Riemannian manifold, then we show distance between the two distributions is: v two methods to compute the distance in the manifold of u n u X |ρ1 | + |ρ2 | 2 multivariate Gaussians. dD (θ1 , θ2 ) = t2 (log ) (9) i=1 |ρ1 | − |ρ2 | A. Derivation of an appropriate metric Sub-manifold NΣ where Σ is constant: The distance between Information geometry allows to find a distance on the space two distributions θ1 = (µ1 , Σ0 ) and θ2 = (µ2 , Σ0 ) is: of probability distribution parameters. The length element on q the manifold of multivariate Gaussians, is obtained when we dΣ (θ1 , θ2 ) = (µ1 − µ2 )2 Σ−1 0 (µ1 − µ2 ) (10) develop the Kullback-Leibler divergence up to the order 2. It writes: Z We have to note that the sub-manifolds ND and NΣ are p(x) not entirely geodesic, i.e. dD and dΣ are only upper bounds K(p, q) = p(x). ln( ) (3) q(x) for the Fisher-Rao distance if we consider the global space of Let us call p1 = p(x|θ) and p2 = p(x|θ + dθ). The multivariate Gaussian distributions. development of K to the order 2 gives: Sub-manifold Nµ where µ is constant: The distance is: n n ∂ 2 ln(p(X|θ)) 1X d2µ (Σ1 , Σ2 ) = [log(λi )]2 = dF ((µ, Σ1 ), (µ, Σ2 )) (11) X 2 ds = −2K(p1 , p2 ) = E( )dθi dθj (4) 2 i=1 i,j=1 ∂θi ∂θj For the particular manifold of multivariate Gaussians, we where 0 < λ1 ≤ λ2 ≤ ... ≤ λn are the eigenvalues of Σ−1 1 Σ2 . can compute ds2 explicitly thanks to the probability density We can now use these expressions to deduce one lower −( n ) bound and three upper bounds for the distance in Np . (x−µ)t Σ−1 (xµ ) p(x; µ, Σ) = (2π) 2 √ exp(− 2 ), and we obtain: det(Σ) Lower bound: The authors of [4] have derived a lower bound 1 thanks to an isometric embedding of the space of parameters ds2 = dµt Σ−1 dµ + tr[(Σ−1 dΣ)2 ] (5) of Np into the manifold of positive definite matrices. Given 2 ! Let a curve θ(t) = (θ1 (t), ..., θn (t)), with t1 ≤ t ≤ t2 Σi + µi µti µti θ1 = (µ1 , Σ1 ) and θ2 = (µ2 , Σ2 ), let Si = verifying θ(t1 ) = θ1 , θ(t2 ) = θ2 and ∀t ∈ [t1 , t2 ], θ(t) ∈ Θ. µi 1 We also suppose that this curve is of class C 1 by piece in Rn . for i = 1, 2. A lower bound for the distance between θ1 and The length of the curve θ(t) is: θ2 is: v u n+1 Z t2 Z 1 X n u1 X ds LB(θ1 , θ1 ) = t [log(λi )]2 (12) | dt| = | { gij (θ(t))θ˙i θ˙j }| (6) 2 i=1 t1 dt 0 i,j=1 Among all the curves, we call geodesic a curve of minimal where λk , 1 ≤ k ≤ n + 1 are the eigenvalues of S1−1 S2 . length (which is not a priori unique). The distance is: First upper bound: In [5], the authors compute the Fisher- Z 1 X n Rao distance in Np between two points of the sub-manifold d(θ1 , θ2 ) = inf | { gij (θ(t))θ˙i θ˙j }| (7) Nα = {pθ ; θ = (µ, Σ); Σ = αΣ0 , Σ0 ∈ Pn (R), α ∈ R∗+ }. t→θ(t)∈∆ 0 i,j=1 Given two distributions θ1 = (µ1 , Σ) and θ2 = (µ2 , αΣ), √ We will now explicitely derive the equations of the 2 2 α 1 1 dF (θ1 , θ2 ) = dα (θ1 , θ2 ) = 2arccosh( + √ + √ δδ t ) geodesics for this metric. We cannot give a formal expression 2 2 α 4 α (13) in the general case. However, it is possible to determine a 1 where δ = Σ− 2 (µ2 − µ1 ). system of differential equations satisfied by these geodesics. Then, we just have to apply the triangular inequality and Let us note Γijk the Christoffel symbol of the first type. choose the value of α: For θ1 = (µ1 , Σ1 ), θ2 = (µ2 , Σ2 ) and n n θα = (µ2 , αΣ1 ), we have: gik θ¨i + Γijk θ˙i θ˙j = 0, k = 1, ..., n X X (8) i=1 i,j=1 dF (θ1 , θ2 ) ≤ dα (θ1 , θα ) + dµ (θα , θ2 ) (14) B. Bounds computation for the distance U Bα = dα (θ1 , θα ) + dµ (θα , θ2 ) is thus an upper bound A first method to evaluate the distance is to establish some for the Fisher-Rao distance in Np . So a first upper bound is −1 −1 1 bounds. We summarize the results obtained in [4]. The method obtained by choosing α = ||Σ1 2 Σ2 Σ1 2 || n . uses the fact that we know exact expressions of the distance Second upper bound: Another way to choose α enables for certain sub-manifolds of Np . We detail these particular us to define a second upper bound U B2: we choose α = sub-manifolds where we know the distance. minα {dα (θ1 , θα ) + dµ (θα , θ2 )}. Sub-manifold ND where the covariance Σ is diagonal: Third upper bound: This new upper bound is based on Let Σ = diag(σ12 , ..., σn2 ), σi > 0, ∀i. Let us call ρ1 = isometries in the manifold Np and on the distance in the manifold ND . First let us state the following proposition (the initial velocity, only with the starting and finishing points. proof is in [4]: We use the algorithm described in [8]. We first introduce the Proposition: For any two distributions in Np , we have: following notations: B = expA (V ) and V = logA (B). The parallel transport equations in the manifold enable to transport dF ((µ1 , Σ1 ), (µ2 , Σ2 )) = dF ((0, In ), (µ3 , D3 )) (15) the vector field V on a curve of Np . With 0 the null vector, In the identity matrix, D3 is the The Jacobi field is a vector field defined along a geodesic γ. −1 −1 diagonal matrix given by the eigenvalues of A = Σ1 2 Σ2 Σ1 2 , For example, we consider a geodesic γ between θ and θ1 with Q is the orthogonal matrix whose columns are the eigenvectors an initial tangent vector V , and we suppose that V undergoes −1 a perturbation to become V + W . The variation of the final of A (A = QD3 Q), and µ3 = Qt Σ1 2 (µ2 − µ1 ). Proposition: The Fisher-Rao distance between two multivari- point θ1 can be determined thanks to the Jacobi field, with J(0) = 0 and J(0)˙ = W . In terms of exponential map: ate Gaussian distributions θ1 = (µ1 , Σ1 ) and θ2 = (µ2 , Σ2 ) is upper bounded by (see [4]): d J(t) = expθ0 (t(V + τ W ))|τ =0 (20) v u n p p dτ u X |1 − D3i |2 + |µ3i |2 + |1 + D3i |2 + |µ3i |2 U B3 = t2 2 log ( p p ) We present the principle of the algorithm on figure 1. It is the |1 + D3i |2 + |µ3i |2 − |1 − D3i |2 + |µ3i |2 same algorithm as in [8], except that we have an exact formula i=1 (16) for the geodesic knowing the initial tangent vector. where µ3i are the coordinates of µ3 and D3i the diagonal terms of D3 , as established in the proposition. C. Geodesic shooting algorithm In order to use a shooting algorithm, we first need to derive the exact expression of the geodesic linking two distributions, knowing the initial speed. This is done in [6]. We summarize here the main steps. A geodesic on Np satisfies: ( Σ̈ + µ̇µ̇t − Σ̇Σ−1 Σ̇ = 0 (17) µ̈ − Σ̇Σ−1 µ̇ = 0 Fig. 1: Shooting algorithm principle Let ∆(t) = Σ(t)−1 and δ(t) = Σ(t)−1 µ(t). The system of differential equations partially integrated is then: The space we consider, Np is of negative curvature, so we have the existence and unicity of the minimal geodesic. ˙ = −B∆ + xδ t    ∆ The geodesic distance is then given by: δ̇ = −Bδ + (1 + δ t ∆−1 δ)x (18) r 1 d = µ̇(0)Σ−1 (0)µ̇(0) + tr((Σ−1 (0)Σ̇(0))2 ) (21)  ∆(0) = Ip , δ(0) = 0  2 We then have the initial velocities of the curves: We want to compare the behavior of the bounds and the one ˙ ( ∆(0) = −B(B ∈ Symp (R)) of the shooting algorithm. We thus perform two experiences (19) (as in [4]): we want to compute the distance between (0, I2 ) δ̇(0) = x(x ∈ Rp ) and a varying rotation matrix. In the first experience, we fix   −B x 0 its eigenvalues, and make the angle vary. In the second expe- Let A =  xt  0 −xt  ∈ M at2p+1 (R).  rience, we fix the angle and make the eigenvalues vary. The curves obtained are presented figure 2, only for comparison 0 −xt B between the lower bound and upper bound 3 (the upper bound Eriksen in [7] has proven the following result: that performs best, see [4]) and the geodesic distance. Proposition: For t ∈ R, let Λ(t) := exp(tA) = We observe that the bounds have the same behavior as the ∆ δ Φ P∞ (tA)n  t true geodesic distance, this will be important for the following =:  δ ǫ γ t . Then the geodesic that goes  n=0 n! since the computation of the bounds is very fast, whereas the Φt γ Γ algorithm reclaims iterations. through the origin (0, Ip ) with initial tangent (x, −B) is the curve given by (δ(t), ∆(t)). D. Gradient descent to find the medians In order to obtain the equation of the geodesic, we can When performing tracking quality monitoring, we compute compute the exponential as in [6]. the distance between the estimation and the prediction at time Algorihm to compute the distance between multivariate k, but the curves obtained are imperfect. A common method to Gaussians: solve this problem is to take the median over five estimations The preceeding formula needs to know the initial tangent and then compute the distance to the prediction. R The median vector. The algorithm we present here is able to compute this m is found by minimizing f : m ∈ Np → A d(m, a)da.  Fig. 4: Reference curve of the maneuver acceleration Fig. 2: Comparison between bounds and geodesic distance Fig. 5: Difference between the distances for 4 bounds (left), Fig. 3: Gradient descent on Np (left), Behavior of the algo- 2 bounds and some geodesic points (right) rithm, with geodesics between three points (right) V. C ONCLUSION AND FUTURE WORK Let us call expm the exponential map at point m, and γv the A new method to detect tracking quality degradation has geodesic with initial tangent v. With γ̇(0) = −∇f (mn ), ∇f = been presented. We have chosen a space that represents quite R exp−1 (a) − A ||expm da, the Karcher flow is defined by: well the data we have at the output of an IMM: the space of −1 m (a)|| multivariate Gaussians. On this space we use a distance that mn+1 = γn (tn ) = expmn (−tn .∇f (mn )) (22) represents both the state and the covariance of the error. We We perform a regular gradient descent with the Karcher flow. note also that the bounds are sufficient to compute the distance, The principle of the algorithm is presented on figure 3 left, since we detect the degradation in the figure were there are the behavior on figure 3 right. only bounds, but the shooting algorithm is essential to compute The drawback of this method is the need of the shooting medians. The method proposed detects the degradations quite algorithm, because it is computationally expensive. well. Future work will consist of using the median to suppress IV. A PPLICATION TO TRACKING QUALITY MONITORING part of the imperfections of the curves and develop a way to To perform tracking quality monitoring, we first take an compare the different methods: we want to be able to define IMM filter to make predictions. We then evaluate the distance a signal-to-noise ratio, or threshold values so that we are able (i) (i) between (xk , Pk|k ) and (x̂k|k−1 , Pk|k−1 ), ie the distance be- to perform relevant comparisons between methods. tween the estimation at time k and the predictions we make R EFERENCES for each model. We then evaluate the difference between [1] M. Pilté, F. Barbaresco, Tracking Quality Monitoring Based on Informa- the results obtained for each model. Degradation in tracking tion Geometry and Geodesic Shooting, Submitted at International Radar quality correspond to moments where the difference switches Symposium Conference, 2016 sign. In practice, the results are noisy, so we consider there is [2] J. Ru, A. Bashi, Z. Rong Li, Performance Comparison of Target Maneuver Onset Detection Algorithms, Conf. on Signal and Data Processing of a degradation when the difference changes abruptly. Small Targets, 2004 We give an example of a detection of the tracking quality [3] J. B. B. Gomes, An Overview on Target Tracking Using Multiple Model degradation when a target switches from a rectilinear move- Methods, Instituto Superior Tecnico, 2008 [4] J. Strapasson, J. Porto, S. Costa, On Bounds for the Fisher-Rao Distance ment to a circular movement. We thus have two models, and Between Maultivariate Normal Distributions, AIP Conference Proceed- compute and plot the difference of the distances. The reference ings, Volume 1641, Issue 1, p.313-320, 2015 curve is given on figure 4. The results for the four bounds are [5] M. Calvo, J. Oller, A distance between elliptical distributions based in an embedding into the Siegel group, Journal of Computational and Applied shown on figure 5, left. We also show the results only for LB Mathematics, 2001 and U B3 along with some points of the geodesic distance, this [6] T. Imai, A. Takaesu, M. Wakayama, Remarks on geodesics for multivari- is shown on figure 5, right. The reason why we plot only some ate normal models, Journal of Math-for-Industry, Vol. 3 (2011B-6), pp. 125130, 2011 points for the distance computed with the shooting algorithm [7] P. S. Eriksen, Geodesics connected with the Fisher metric on the multi- is that the computation time is quite long. The method is able variate normal manifold, Proceedings of the GST Workshop, 1987 to detect the quality degradation when the target changes mode [8] M. Han, F. C. Park, DTI Segmentation and Fiber Tracking Using Metrics on Multivariate Normal Distributions, J Math Imaging Vis, 2013 and the bounds are sufficient.