Boosting Target Tracking Using Particle Filter with Flow Control Nima Moshtagh, Moses W. Chan Advanced Technology Center Lockheed Martin Space Systems Company 1111 Lockheed Martin Way, Sunnyvale CA 94089 ABSTRACT Target detection and tracking with passive infrared (IR) sensors can be challenging due to significant degradation and corruption of target signature by atmospheric transmission and clutter effects. This paper summarizes our efforts in phenomenology modeling of boosting targets with IR sensors, and developing algorithms for tracking targets in the presence of background clutter. On the phenomenology modeling side, the clutter images are generated using a high fidelity end-to-end simulation testbed. It models atmospheric transmission, structured clutter and solar reflections to create realistic background images. The dynamics and intensity of a boosting target are modeled and injected onto the background scene. Pixel level images are then generated with respect to the sensor characteristics. On the tracking analysis side, a particle filter for tracking targets in a sequence of clutter images is developed. The particle filter is augmented with a mechanism to control particle flow. Specifically, velocity feedback is used to constrain and control the particles. The performance of the developed “adaptive” particle filter is verified with tracking of a boosting target in the presence of clutter and occlusion. Keywords: P article Filtering, Target Tracking

1. INTRODUCTION Passive infrared (IR) sensors could be used for detection and tracking of ballistic missiles. Missile detection and tracking with passive IR sensors can be challenging due to significant degradation and corruption of target signature by the atmospheric transmission and clutter effect. We are interested in detecting and tracking a target in the boost stage, where its IR signature can be detected using IR sensors. Our main objective is to track a boosting target in the cluttered background, and accurately estimate its state (position, velocity, and acceleration). The estimates of the boosting target are used to predict its position during the ballistic phase, where the target’s motion is only governed by gravity. Traditional approaches to target detection and tracking in image sequences involve signal processing and thresholding of the data to detect candidate targets in the image. This approach may result in a large number of detections, including many false detections. Data association methodologies are developed to associate detections to the estimated tracks. When the background clutter is high (low signal-to-noise ratio), the probability of dropping track increases. Hence a number of tracking methodologies are developed that do not throw away older images, and work with information accumulated over the sequence of input images. Particularly, in the Bayesian tracking approach, the entire image is treated as the measurement, and probability of target existence in every pixel is computed. In target tracking literature, the approaches that avoid the explicit detection stage are known as track-before-detect methods.1–3 Particle filtering3–7 provides an approximate solution to a Bayesian filter by representing the posterior function∗ as a set of random samples and their associated weights. The particles and weights are updated iteratively as new observations are obtained, and they are used to obtain an estimate of target state. By constraining the motion of particles, one can control the particle flow in particle filtering. Particle flow moves the particles to a region of the state space with a higher probability of target existence. Recently, there has been a number of work on particle flow control.8, 9 Papi et al.9 show that prior knowledge on target’s state can be incorporated in either the prediction or the update steps of the Bayesian particle filter. However, the preceding works only consider fixed (hard) constraints given by prior knowledge of the target. Hard constraints limit the applicability of this ∗

Probability density of the state vector conditioned on all the previous measurements

200 180 160 140 120 100 80 60 40 20

Figure 1. A synthetic scene image, simulating the view from an infrared sensor. Pixels with larger values represent warmer areas.

approach. Tracking a maneuvering target with varying velocity and acceleration requires an adaptive constraint that can be updated as needed. In this work, we develop an adaptive particle filter algorithm that uses a feedback mechanism in the state update of the filter to control the flow of the particles. In particular, a velocity feedback is used to constrain the motion of the particles. The difference of our approach with previous methods is that using feedback we dynamically update the constraints. The advantage of such adaptive approach is that it can be applied to tracking maneuvering objects such as a boosting target, and it results in better tracking performance when target is occluded due to background clutter.

2. END-TO-END SIMULATION TESTBED A high fidelity end-to-end simulation testbed was used to generate a set of realistic clutter images by modeling atmospheric clutter and solar reflections. The dynamics and intensity of a boosting target are modeled and injected onto the background scene. Pixel level images are then generated with respect to the sensor characteristics. Figure (1) shows a sample clutter image generated by the toolbox. A sequence of such clutter images is used as the input to tracking algorithms. Figure (2) shows the architecture of the simulation testbed. The main components of the simulation toolbox are: • Boosting Target Modeling: Section 3 presents the motion model of a typical boosting target. The motion of a two-stage boosting target is simulated, and its motion in the image is obtained by projection of the three dimensional trajectory onto the image plane. Target dynamics in image plane are also derived in Section 3. • Sensor Modeling: The brightness of the target in the image plane is a mapping of the target radiance onto the image plane, which is affected by the sensor model and the background clutter. As the target moves across the background clutter, its appeared brightness could change due to the additive background noise. In this work, the target’s radiance is assumed to be proportional to its thrust profile. Section 4 presents the sensor model that was used to generate the clutter images, and the associate measurement model that was used for filtering. • Background Clutter Modeling: Synthetic Scene Generation Model (SSGM)10 software was used to generate high-fidelity background scenes. SSGM uses physics-based models of clouds, atmospheric transmission, and clutter effects to simulate the environment in which ballistic missiles are detected, tracked and engaged. The details of the scene are adjusted according to sensor characteristics such as field of view and distance to the scene.

Figure 2. The end-to-end simulation environment is comprised of phenomenology and tracking toolboxes.

• Target Tracking Toolbox: An object-oriented software toolbox is developed for implementation and evaluation of the tracking algorithms. Each implemented algorithm is an object from a class of algorithms. For example variations of the particle filtering algorithms are objects derived from the class of ParticleFilter. The class properties define the parameters that are used in the algorithm, and the class methods are the operations of the filter such as update, propagation, and plotting of the tracking results. Section 5 presents the theoretical background on Bayesian filtering for target tracking in image sequences. A particle filter implementation of a Bayesian filter is then developed in Section 6. Section 7 develops our feedback approach to particle flow control. In order to evaluate the effect of flow control in particle filtering, the performance of our adaptive particle filter is compared with that of a standard particle filter and an “efficient” estimator, which has performance equal to the Cramer-Rao lower bound (CRLB). The efficient estimator is developed in Section A. Numerical analysis in Section 8 shows the advantage of using flow control in particle filtering in the presence of background clutter.

3. MOTION MODEL OF A BOOSTING TARGET The total acceleration of a boosting target includes thrust, drag, and gravity:11 a = aT + aD + aG .

(1)

The thrust acceleration is aT = (T /m)uT where T is thrust magnitude, m is target mass, and uT is the thrust direction. The gravity acceleration is aG = −(µ/d3 )d, where µ is Earth’s gravitational constant, d is the position vector, from Earth center to the target, d = ∥d∥ is its length. During the boost phase the net thrust, aN = aT + aD , is significantly larger than the gravity. In our studies, it is more convenient to estimate the net acceleration, rather than its components. The thrust profile of Figure (3.a) is used to simulate a boosting target. Target’s acceleration was then computed using the acceleration model (1). Position and velocity of target were computed from integration of the acceleration values. Acceleration and velocity profiles are shown in Figures (3.b) and (3.c), respectively. The profile of the thrust magnitude T depends on the kind of the boosters used in the missile. Solid fuel booster engines burn their fuel with a constant rate, and generate constant thrust during each boost phase. As the fuel burns, the mass of the missile decreases, which results in the nonlinear increase in the acceleration, which can be seen in Figure (3.b). The projection of the target trajectory (in the Earth-centric inertial frame) onto the image plane generates a target trajectory in the image. The state of the target in the image plane is represented by its position p,

5

9

Sample Target Thust Profile

x 10

Norm of Target Acceleration Vector

Norm of Target Velocity Vector

45

3500

40

8

3000

35 2500

6

5

30 Velocity [m/s]

Acceleration [m/s2]

Thrust [Kg.m/s2]

7

25 20

2000

1500

15

4

1000 10

3

2

500

5

0

20

40

60 80 time [second]

100

120

140

0

0

20

40

(a)

60 80 time [second]

100

120

140

(b)

0

0

20

40

60 80 time [second]

100

120

140

(c)

Figure 3. A boosting target with two boost stages: (a) thrust profile [Kg.m/s2 ], (b) acceleration profile [m/s2 ], (c) velocity profile [m/s].

velocity v and acceleration a in units of pixels, pixels/second and pixels/second2 , respectively. In the image plane, we consider a near constant acceleration model for the boosting target. Suppose the image observations arrive at every Ts intervals. Thus the discrete-time equations of motion become

where xk is the 6-dimensional state vector

A is the transition matrix

xk+1 = Axk + wk ,

(2)

  pk xk = vk  , ak

(3)

 I A = 0 0

Ts · I I 0

 Ts2 /2 · I Ts · I  I

(I is the 2 × 2 identity matrix). wk = Bwk is a white Gaussian noise with wk ∼ N (0, q 2 ) where B is  3  Ts /3 · I B = Ts2 /2 · I  Ts · I .

(4)

(5)

Thus the covariance of process noise wk is Qk = E[wk wkT ] = q 2 BB T .

(6)

The process noise power spectral density q 2 is a design parameter. Hence, this model can be used for tracking maneuvering targets by defining an appropriate process noise.

4. SENSOR MODELING Consider an IR sensor that observes a region on the x − y plane, and records measurements at discrete times tk , k = 1, . . . , K. The sensor array consists of N × M pixels. The measurement at pixel (i, j) at time tk is denoted by zk (i, j). Thus, the sensor measurement at time tk is a N × M matrix, denoted by zk (Figure (1) shows a synthetic IR image). Each pixel has dimension ∆x × ∆y. The intensity of a point source is spread according to a point spread function (PSF). The PSF is a 2D Gaussian function with variances σx and σy , which represent the extend of blurring along each dimension. The measurement at pixel (i, j) is a combination of target intensity and the background noise zk (i, j) = Ik (i, j|xk ) + nk (i, j)

(7)

where noise nk is modeled by a Gaussian distribution nk (i, j) ∼ N (µk (i, j), σk2 (i, j)), and xk is target’s state at time tk . For a moving target, one could write Ik (i, j|xk ) = ak (i, j|xk )Ik , where ak (i, j|xk ) is the integral of the PSF over the extent of the pixel (i, j), and over the time between tk−1 and tk ; thus, it represents the contribution of target with state xk to pixel (i, j)2 † . Note that ak (i, j|xk ) will be non-zero only near the pixel in which the target is located, at time tk . Thus, the intensity of pixel (i, j) for a target with intensity Ik and state xk is given by zk (i, j) = ak (i, j|xk )Ik + nk (i, j) . (8) Therefore, sensor measurements can be modeled as zk = h(xk ) + nk .

(9)

As measurements zk arrive, the spatially correlated background noise can be learned over time (see Section 4.1). In this work, the target’s radiance is assumed to be proportional to its thrust profile (Figure (3.a)). Therefore, in our simulations, target brightness drops sharply as it transitions from the first boost stage to the second stage. Figure (4.a) shows the target’s intensity in units of W/(sr.m2 ), and Figure (4.b) shows the SNR over the time span of the simulation.

−3

2.2

Target Radiance

x 10

Signal−to−Noise Ratio for a Boosting Target with Two Boost Stages. 24 22

2

20 18 1.6

16 SNR

Radiance [W/(sr.m2)]

1.8

1.4

14 12

1.2

10 1 8 0.8 0.6

6 0

20

40

60 80 Time [seconds]

100

120

140

4

0

(a) Target Radiance in W/(sr.m2 )

20

40

60

80

100

120

140

time

(b) Signal to noise ratio

Figure 4. (a) Radiance of a boosting target with two stages, (b) the corresponding target SNR computed from the image sequence.

4.1 Signal-to-Noise Ratio and Background Clutter For signal-to-noise ratio (SNR) computation, we make the simplifying assumption that we are dealing with a point target. For point targets, the point spread function tends to a delta function; thus the target only contributes to the pixel in which it falls, and one has { Ik + nk (i, j) pk ∈ (i, j) zk (i, j|xk ) = (10) nk (i, j) otherwise. †

Using a crude approximation of the integral, ak (i, j|xk ) is given by:1 ( ) (i∆x − x)2 (j∆y − y)2 ∆x ∆y exp − − . ak (i, j|xk ) = 2πσx σy 2σx2 2σy2

The background noise at pixel (i, j) is modeled as a Gaussian distribution n(i, j) ∼ N (µ(i, j), σ 2 (i, j)). Suppose at time k target is located in pixel (i, j) (i.e. pk ∈ (i, j)). We define SNR as: SN Rk (i, j|xk ) =

zk (i, j|xk ) − µk (i, j) . σk (i, j)

(11)

From (11), one can see that one must “learn” the features of the background clutter to compute target’s SNR. A running Gaussian average of the image measurements is used to adaptively learn the background, i.e. compute the mean µk and variance σk2 of the background:12 µk σk2

= =

αzk + (1 − α)µk−1 2 α(zk − µk )2 + (1 − α)σk−1

(12) (13)

where α ∈ [0 1] is the learning rate, µk ∈ RN ×M , and σk2 ∈ RN ×M . The learning rate determines how fast the changes in the scene must be associated to changes in the background. When the variations in the background clutter are small, a small value for the learning rate must be used (α ≪ 1). A slow moving target may be “learned” as a change in the background, if large values of the learning rate (α ≃ 1) are used.

5. BAYESIAN FILTERING Let xk = x(tk ) represent the state of a target at time tk . Consider the evolution of the state of the target xk = f (xk−1 , qk−1 ) ,

(14)

where f (·) is possibly a nonlinear function of the state, and qk is the Gaussian process noise. The Markov property holds for system (14), i.e. p(xk |xk−1 , . . . , x0 ) = p(xk |xk−1 ) , where p(xk |xk−1 ) is known as the prior function. The objective of tracking is to recursively estimate the state xk from measurements zk = h(xk , nk ) , (15) where nk is a Gaussian measurement noise. The measurement zk is independent from the past states, i.e. p(zk |zk−1 , . . . , z1 , xk , . . . , x0 ) = p(zk |xk ) , where p(zk |xk ) is known as the likelihood function. We would like to estimate xk using all the available measurements Z 1:k = {z1 , . . . , zk } up to time tk , i.e. we would like to compute the conditional probability density p(xk |Z 1:k ). In Bayesian filtering,5, 13 the posterior p(xk |Z 1:k ) is computed using a two-step recursion. In prediction step the posterior is predicted using the prior function: ∫ 1:k−1 p(xk |Z ) = p(xk |xk−1 )p(xk−1 |Z 1:k−1 )dxk−1 . (16) When measurement zk becomes available, the likelihood function p(zk |xk ) is computed. In the update step the likelihood function is used to update the predicted posterior: p(xk |Z 1:k ) =

p(zk |xk )p(xk |Z 1:k−1 ) , p(zk |Z 1:k−1 )

(17)

where p(zk |Z 1:k−1 ) is the normalization constant. The minimum mean square (MMS) estimate and its covariance are computed as the following: ∫ MMS ˆ k|k x = xk p(xk |Z 1:k )dxk ∫ MMS MS MS T ˆM ˆM Pk|k = (xk − x )(xk − x ) dxk . k|k k|k

The recursive solutions (16) and (17) provide the optimal Bayesian solution for the posterior. Closed form solutions only exist for some special cases such as Kalman filtering. When the analytical solution is intractable, estimation methods such as grid-based Bayesian filtering 5 or particle filtering approximate the optimal Bayesian solution.

Figure 5. Recursive Structure of a Bayesian Filter.

6. PARTICLE FILTERING A technique to implement the recursive Bayesian filter is using sequential Monte Carlo simulations. The key idea of this approach (a.k.a. particle filtering) is to represent the posterior function by a set of random samples t 5 and their associated weights {xi , wi }N i=1 . Formally, the posterior is represented by a discrete approximation p(xk |Z 1:k ) ≈

Nt ∑

wki δ(xk − xik )

(18)

i=1

where δ(·) is dirac delta function. The samples are chosen from an importance density, xi ∼ π(x). If we choose the importance density to be the prior p(xk |xk−1 ), the update equation for the weights becomes the following (recursive) equation: i wki ∝ wk−1 p(zk |xik ) , i = 1, . . . , Nt . ∑ i The weights are later normalized so that i wk = 1. This is equivalent to the measurement update equation in the Bayesian filter. Thus, at each iteration new particles are drawn from the importance density ‡ , their corresponding weights are updated, and the posterior is estimated. The approximated posterior at each iteration is used to obtain the state estimate iteratively. As the number of samples increases, the particle filter’s estimate approaches the optimal Bayesian estimate. To avoid clustering of the particles into a single point (degeneracy), the approximate posterior (18) (which is a discrete density) is resampled whenever the effective number of particles falls below a threshold. During resampling, Nt i.i.d. samples are drawn from the particle set {xik , wki }, and the weights of the new particles i = 1/Nt ; thus the weights become are set to 1/Nt . If resampling is applied at every iteration, we have wk−1 proportional to the likelihood function wki ∝ p(zk |xik ) . (19) The likelihood function p(zk |xk ), and the prior function p(xk |xk−1 ) are chosen based on the measurement model and the dynamical model (respectively) for the specific application in hand. ‡

The prior distribution p(x0 ) could be a Uniform distribution over the state space, e.g. over the image: [0 N ] × [0 M ] pixels.

6.1 Measurement Likelihood Using sensing model (8) the measurement likelihood function given that there is a target in position p becomes: ∏ ∏ ∏ pz (m, n|p) · pz (m, n) (20) p(z|p) = pz (m, n|p) = m,n

Ω(p)

¯ Ω(p)

where pz (m, n|p) is the conditional probability density of target intensity at pixel (m, n). Ω(p) is the set of pixels ¯ around target location p where target intensity is non-zero (in other words a(i, j|p) ̸= 0 in (8)), and Ω(p) is its complement. Note that given target location, target intensity distributions over the pixels are independent, and the measurement likelihood is expressed as the product of individual likelihoods. In data fusion problems it is sometimes easier to work with the log of a probability distribution. Loglikelihoods are computationally more convenient, and numerically more stable. Also, minus log-likelihoods are closely related to the information content of measurements. Therefore, we define the (negative) log-likelihood function: ( ) [ ] ∑ z(n, m|p) − µ(n, m) 2 w = − log p(z|p) = + constant, (21) σ(n, m) Ω(p)

where the constant term depends on the image size and noise properties of the clutter image. In our implementation of the particle filter, we used the (negative) log-likelihood function in the update step, as in (21). Since the particle weights will be normalized, the constant term in (21) is omitted: ∑ ( zk (n, m|xi ) − µk (n, m) )2 k wki = . (22) σ (n, m) k i m,n∈Ω(xk )

The size of the set Ω(xik ) depends on the amount of blurring of the target in the clutter images. Given σx and σy as the amount of blurring along x and y axis in the image, Ω can be chosen such that it contains 3σx and 3σy along the respective directions.

7. PARTICLE FILTERING WITH FLOW CONTROL By constraining the motion of particles using prior knowledge about targets state, one can control the particle flow in a particle filter. Papi et. al 9 show that prior knowledge on target’s state can be incorporated in either the prediction or the update steps of the Bayesian particle filter. Rejection-sampling method uses the prior knowledge in the prediction step by rejecting particles that do not satisfy the constraint, and only propagates those that do. Pseudo-Measurement method uses the prior knowledge in the update step by lowering the weights of the particles that do not satisfy the constraint, resulting in their elimination in the resampling step. The two methods of implementing a constrained Bayesian filter are shown to be theoretically equivalent (for large number of particles).9 Papi et al.9 considered interval “hard” constraints of the form s1 ≤ C(x) ≤ s2 , where s1 and s2 are known boundary values of the constraints. In this work, we consider “soft” constraints of the form C(x) = s + ν, (23) where s is the desired value for C(x) and ν ∼ N (0, Rν ) is the uncertainty in the adherence to the constraint. Tahk and Speyer14 introduced the idea of defining a pseudo-measurement to handle state constraints. Similarly, we define an augmented measurement vector in order to incorporate equality state constraints of the form (23) in the update stage of particle filtering. The augmented measurement vector is given by: [ ] [ ] [ ] z h(x) n ¯= z = + (24) s C(x) ν

where the covariance of the augmented measurement is R = diag(Rn , Rν ). The constraint-based measurement likelihood is now given by: ( ) ( ) 1 T −1 T −1 p(¯ z|x) = exp − 0.5(z − h(x)) Rn (z − h(x)) · exp − 0.5(s − C(x)) Rν (s − C(x)) κ = p(z|x)p(C|x) where κ is a constant normalization factor, and C is the set of states that satisfy the constraint, C = {x : C(x) ∼ N (s, Rν )} , and p(C|x) is the likelihood that the constraint is satisfied. We use the constraint-based likelihood function to control the flow of particles. At time tk the predicted ˆ k|k−1 . By defining an appropriate constraint, we assign higher weights to particles that have target velocity is v velocities close to the target velocity, and penalize particles with much lower or much higher velocities. In other words, the prior knowledge on target velocity is used to control the particle flow. We define C(x) as a time-dependent constraint on target velocity ˆ k|k−1 + µ . vk = v (

Thus we have p(C|xk ) ∝ exp

−1 ˆ k|k−1 )T Rν−1 (vk − v ˆ k|k−1 ) (vk − v 2

(25) ) .

(26)

Therefore, the weight of each particle is updated according to this update equation: wki ∝ p(¯ z|xik ) = p(z|xik )p(C|xik ) , i = 1, . . . , Nt .

(27)

Rν determines how much trustworthy the prior information is, and can be considered as a design parameter. Another way to interpret Rν is to consider it as the degree of “softness” of the constraint. A larger value of Rν means softer constraint, and vice versa. Figure (6) shows the architecture of the “adaptive” particle filter with the constrained update equation.

Figure 6. The adaptive particle filter with constrained update equation. (A patent is pending on this technology.)

An implementation of the above adaptive particle filtering is as the following: • Initialization: Generate Nt particles from the state prior distribution. xi0 ∼ p(x0 ),

w0i = 1/Nt , i = 1, . . . , Nt .

• Repeat the following steps: 1. Constraint Computation: Compute the constraint-based likelihood function p(Ck |xi ), similar to (26), for all particles i = 1, . . . , Nt . 2. Measurement Update: Update the weights of all particles using the likelihood function: wki ∝ p(zk |xik )p(Ck |xik ) , i = 1, . . . , Nt

∑ and normalize to wki := wki / j wkj . 3. State Estimation: Compute the minimum mean square (MMS) estimate and its covariance using: MS ˆM x k|k

ˆk = ≈ x

Nt ∑

wki xik ,

i=1 MMS Pk|k

≈ Pk =

Nt ∑

ˆ k )(xik − x ˆ k )T . wki (xik − x

i=1 t 4. Resampling: Take Nt samples (with replacement) from the set {xik }N i=1 , where the probability of i sampling particle i is wk . 5. Prediction: Take Nt samples qik ∼ p(xk |xk−1 ), and propagate the states:

xik+1 = f (xik , qik ) , i = 1, . . . , Nt .

8. PERFORMANCE OF THE PARTICLE FILTER WITH FLOW CONTROL 8.1 Simulation Scenario To show the effect of using flow control in particle filtering, a scenario is generated where the target’s brightness is dropped to an undetectable level as its trajectory is occluded by clouds. Figure (8.a) shows the trajectory of the target in the clutter image. Target becomes invisible as it gets occluded by the background clutter. The segment with occlusion is highlighted in the image, which corresponds to time interval between t = 75 to t = 90 seconds. Figure (8.b) shows that during that interval the signal-to-noise ratio drops to undetectable values. For performance analysis, we compared the estimation results of our adaptive particle filter with the performance of a benchmark estimator that we call the “efficient” estimator. Before describing the results, let us describe the efficient estimator in more details.

8.2 Efficient Estimator According to Van Trees,15 an estimator that satisfies the Cramer-Rao lower bound (CRLB) with equality is called an efficient estimator. It is shown by Zhang et al.16 that when there is no uncertainty in measurement origin (i.e. true target is always detected, and probability of detection equals one), and the estimator generates only one estimate (i.e. no false alarms), then the performance of the estimator satisfies the CRLB with equality. See Appendix A for details on CRLB. To achieve the CRLB in our target tracking application, we used the true positions of the target as measurements for a 6-state Kalman filter, as shown in Figure (7). The performance of this filter is equivalent to the theoretical Cramer-Rao lower bound, and hence provides the best state estimates possible. The results of this efficient estimator were used as a baseline for performance evaluation of other filters such as the standard and adaptive particle filters.

Figure 7. In the efficient estimator, perfect detections of the target truth are used as measurements for the Kalman filter.

8.3 Tracking Algorithms The performances of the following tracking algorithms are compared: • The standard particle filter, which was presented in Section 6. The 6-state particle filter is tuned to track a single target with SN R > 1, and velocities in the range of [0 2] pixel/frame. A large number of particles (N = 20000) are used because of the large dimension of the state space (D = 6). Since a single target is being tracked, we have assumed that the posterior density is represented by a single Gaussian. Thus, our standard particle filtering algorithm is similar to Gaussian particle filtering (GPF) introduced by Kotecha and Djuri´c.7 • The adaptive particle filter, which is presented in Section 7. A velocity feedback similar to (25) was used to control the flow of the particles. One design parameter is the covariance of the constraint Rv , representing the amount of uncertainty in the velocity feedback value. In our simulations, Rv proportional to the covariance of target velocity estimates. • The efficient estimator, which was explained in Section 8.A, provided the baseline for performance comparison of the particle filtering algorithms. In the above algorithms the background is iteratively learned using a running Gaussian average method. The learning rate was set to α = 0.05, since the target is relatively slow moving in the image sequence. The performances of the estimators are evaluated using the estimation error ˆ k|k − xk , ϵk = x

(28)

and a scalar measure of the estimation covariance matrix √∑ √∑ ck = Pk|k (i, i) = λ2i , i

where λ2i are the singular values of the covariance matrix Pk|k .

i

(29)

Signal−to−Noise Ratio (SNR) 40

1st stage

2nd stage

35

10

30

20

End 25 Signal/noise

30

Occlusion 40 50

20 15

Occlusion

10

60

5

Start 0

70

20

40

60

80

100

−5

120

0

20

40

(a) Occlusion Scenario

60 80 Time (Second)

100

120

140

(b) Target SNR

Figure 8. The effect of clutter on target SNR.

8.4 Summary of Results The estimated trajectories and velocities of the target using standard particle filter and particle filter with flow control are shown in Figure (9). During the occlusion (t ∈ [75 90] sec), the standard particle filter diverges, and its velocity and acceleration estimates diverge substantially from the truth values. This also can be seen in Figures (10.a) and (10.b) where the estimation errors are compared. Figure (10) also displays the estimation errors of the efficient estimator. One can see from Figure (9) that the adaptive particle filter maintains the track by predicting where the target will be after the occlusion, resulting in a better state estimate when the target re-appears. The covariances are compared in Figure (11). One can see that the uncertainty in the estimates grows using the standard particle filter. However, the estimation errors of the adaptive particle filter remain close to that of the efficient estimator. The results show the effectiveness of using velocity feedback to overcome the negative effect of target occlusion due to background clutter.

X velocity V(estimated) [pixel/sec]

Comparison of the Trajectory Estimates 110 Truth PF Adaptive PF

100 90

End

80

0.2 0 −0.2 −0.4 −0.6 −0.8

Truth PF Adaptive PF 20

40

Y

70

50

V(estimated) [pixel/sec]

60 Start

40 30 20 20

30

40

50

60

X

(a) Comparison of the trajectory estimates.

70

60 seconds Y velocity

80

100

120

60 seconds

80

100

120

2 1.5 1

Truth PF Adaptive PF

0.5 0 20

40

(b) Comparison of the velocity estimates.

Figure 9. Estimated position and velocities of the target.

Position Error [pixels]

X Position Error 2 0 −2

PF Adaptive PF Efficient Estimator

−4 −6 −8 30

40

50

60

70 80 sec Y Postion Error

90

100

110

120

Position Error [pixels]

5

0 PF Adaptive PF Efficient Estimator

−5

30

40

50

60

70 80 sec

90

100

110

120

(a) Position Errors

Velocity Error [pixels/sec]

X Velocity Error 1.5 PF Adaptive PF Efficient Estimator

1 0.5 0 −0.5

Velocity Error [pixels/sec]

30

1.5

40

50

60

70 80 sec Y Velocity Error

90

100

90

100

110

120

PF Adaptive PF Efficient Estimator

1 0.5 0 −0.5 30

40

50

60

70 80 sec

110

120

(b) Velocity Errors Figure 10. Comparison of the position and velocity estimation errors.

COV comparison

Covariance Area

50 PF Adaptive PF Efficient Estimator

40 30 20 10 0

30

40

50

60

70 80 sec

90

100

110

120

Figure 11. The (average) radius of the estimation covariance as defined in (29).

9. CONCLUSIONS & FUTURE WORK Using a high-fidelity simulation testbed, we generated a set of images that simulated the observations of an IR sensor. The motion of a boosting target was modeled and projected onto the image plane according the sensor characteristics. The set of clutter image (containing the target) were used as inputs to a number of tracking algorithms. We developed a mechanism to incorporate prior knowledge about target’s state into the update stage of particle filtering. By using the projected target velocity as the prior knowledge, we were able to control the flow of the particles at every iteration. Such an adaptive particle filter performs well in comparison with a standard particle filter, in situations where target is occluded (for a relatively short period of time) or when it has a very small signal-to-noise ratio. The mode changes of a boosting target can be modeled using a Markov process. Application of particle filtering to a jump Markov linear system17, 18 is the subject of the ongoing work.

10. ACKNOWLEDGEMENTS The authors would like to thanks Christopher Damitz (with Lockheed Martin Advanced Technology Center) for his help in generating clutter background images.

APPENDIX A. CRAMER-RAO LOWER BOUND The Cramer-Rao lower bound provides a lower bound on the estimation covariance matrix: Pk|k ≥ Jk−1 ,

(30)

where the matrix inequality means Pk|k − Jk−1 is positive semidefinite, and Jk is the Bayesian version of the Fisher Information Matrix (FIM).16 The Fisher information matrix Jk for a linear system xk+1 zk

= Axk + vk , vk ∼ N (0, Qk ) = Hxk + wk , wk ∼ N (0, Rk )

is computed in a recursion16, 19 ( )−1 −1 Jk+1 = Qk + AJk−1 AT + H T Rk+1 H ,

(31)

with J0 = P0−1 . The first term in (31) represents the prior information at time k, and the second term in (31) represents the new information from the new measurement at time k + 1.

REFERENCES [1] Salmond, D. J. and Birch, H., “A particle filter for track-before-detect,” in [Proceedings of the American Control Conference], 3755 – 3760 (2001). [2] Tonissen, S. and Bar-Shalom, Y., “Maximum likelihood track-before-detect with fluctuating target amplitude,” IEEE Transactiond on Aerospace and Electronic Systems 34 (3), 796 – 809 (1998). [3] Boers, Y. and Drissen, H., “Particle filter based track before detect algorithms,” in [Proceedings of SPIE Vol. 5204 ], Signal and Data Processing of Small Targets, 20–30 (2003). [4] Gordon, N., Salmond, D. J., and Smith, A. F. M., “Novel approach to nonlinear/non-gaussian bayesian state estimation,” IEE Proceedings-F 140, 107 – 113 (1993). [5] Arulampalam, M. S., Maskell, S., Gordon, N., and Clapp, T., “A tutorial on particle filters for online nonlinear/non-gaussian bayesian tracking,” IEEE Transactions of Signal Processing 50, 174 – 188 (2002). [6] Gustafsson, F., Gunnarson, F., Bergman, N., Forssell, U., Jansson, J., Karlsson, R., and Nordlund, P., “Particle filters for positioning, navigation, and tracking,” IEEE Transactions on Signal Processing 50, 425 – 437 (2002). [7] Kotecha, J. and Djuri´c, P. M., “Gaussian particle filtering,” IEEE Transactions on signal processing 51, 2592 – 2601 (2003). [8] Daum, F. and Huang, J., “Smal curvature particle flow for nonlinear filtering,” in [Proceedings of SPIE Defense, Security and Sensing], (2012). [9] Papi, F., Podt, M., Boers, Y., Battistello, G., and Ulmke, M., “Bayes optimal knowledge exploitation for target tracking with hard constraints,” in [Proceedings of 9th IET Data Fusion & Target Tracking Conference: Algorithms & Applications], (2012). [10] Wilcoxen, B. A. and Heckathorn, H. M., “Synthetic scene generation model (SSGM R7.0),” in [Proceedings of SPIE 2742, Targets and Backgrounds: Characterization and Representation II], 57–68 (1996). [11] Li, X. R. and Jilkov, V. P., “A survey of maneuvering target tracking – part II: Ballistic target models,” in [Proceedings of SPIE Conference on Signal and Data Processing of Small Targets], (2001). [12] Wren, C. R., Azarbayejani, A., Darrell, T., and Pentland, A. P., “Pfinder: Real-time tracking of the human body,” IEEE Transations on Pattern Analysis and Machine Intelligence 19, 780 – 786 (1997). [13] Gustafsson, F., “Particle filter theory and practice with positioning applications,” IEEE Aerospace & Electronics Systems Magazine 25, 53 – 81 (2010). [14] Tahk, M. and Speyer, J., “Target tracking problems subject to kinematic constraints,” IEEE Transations of Automatic Control 35, 324 – 326 (1990). [15] Van-Trees, H., [Detection, Estimation, and Modulation Theory, Part I], Wiley, New York (1968). [16] Zhang, X., Willett, P., and Bar-Shalom, Y., “Dynamic cramer-rao bound for target tracking in clutter,” IEEE Transactions on Aerospace and Electronic Systems 41, 1154–1167 (2005). [17] Doucet, A., Gordon, N. J., and Krishnamurthy, V., “Particle filters for state estimation of jump markov linear systems,” IEEE Transactions on Signal Processing 49(3), 613–624 (2001). [18] Driessen, H. and Boers, Y., “Efficient particle filter for jump markov nonlinear systems,” in [IEE ProceedingsRadar, Sonar and Navigation], 152(5), 323–326, IET (2005). [19] Bell, K. L. and Van-Trees, H. L., “Posterior cramer-rao bound for tracking target bearing,” in [13th Annual Workshop on Adaptive Sensor Array Process], (2005).