Multiple Target Tracking with Quantized Measurements: A Standard Bayesian Approach Lawrence D. Stone Stephen L. Anderson Metron Inc. Reston, USA
[email protected] /
[email protected] Abstract—This paper shows how the use of standard Bayesian likelihood functions and particle filters provides a simple and general approach to tracking with quantized measurements. We demonstrate this result with two examples. The first involves a maneuvering target and the second multiple targets. In both cases the motion model for the target(s) is non-Gaussian and the quantized measurements correspond to a set of possible locations composed of disjoint regions in 2-space. For the multiple target example, we apply a nonlinear particle filter version of JPDA to perform the data association and tracking. Keywords—quantized measurements, tracking, maneuvering targets, multiple targets, particle filters, JPDA
I.
INTRODUCTION
There have been a number of papers, [1]-[8], describing how to perform estimation and tracking with quantized measurements and imprecise measurements or likelihood functions. References [5] and [6] claim that extensions of standard Bayesian inference are required in these situations and give examples to justify their claims. Reference [7] demonstrated that the examples used in [5] and [6] can be easily and correctly handled using standard Bayesian inference and likelihood functions. References [2]-[4] tackle the problem of tracking a single target using quantized measurements. In these papers the motion model is Gaussian and the measurements are quantized linear functions of target state with Gaussian errors. Reference [2] proposes a particle filter approach to tracking a single target with quantized measurements and provides a number of examples. References [3]-[4] base their approach on [1] which involves a rather complex set of discrete approximations to truncated Gaussian distributions or mixtures of Gaussian distributions. The purpose of this paper is two-fold. The first is to show by example that a particle filter naturally and simply tracks a single maneuvering target in a situation where the measurements are complicated and quantized and the motion model of the target is not Gaussian. The second is to show how particle filters with standard Bayesian likelihood functions can perform multiple target tracking even with complicated quantized measurements. This is done through the use of a particle filter implementation of the non-linear JPDA algorithm described in Chapter 4 of [9].
II.
STANDARD BAYESIAN INFERENCE FORMULATION
In order to clarify what we mean by a standard Bayesian approach, we give the formulation of the basic Bayesian inference problem that is presented in [9] and is consistent with [10]. There is an unknown parameter Θ that we wish to estimate. There is a prior probability distribution p0 on Θ such that
p0 (θ ) = Pr {Θ = θ } for θ ∈ S
(1)
where S is the set of possible values of and Pr indicates either probability or probability density as appropriate. We obtain a measurement Z which takes values in a set Z . The measurement Z is a random variable whose distribution depends on θ . For each value of θ ,
Pr {Z = z | Θ = θ } for z ∈ Z is a probability density (distribution) on Z . However, to define the likelihood function, we view the measurement Z = z as data and fix it so that the likelihood function l ( z | θ ) ≡ Pr {Z = z | Θ = θ } for θ ∈ S
(2)
becomes a function of θ . When we use the term likelihood function, we mean we are holding the measurement z fixed and letting the parameter θ vary. We write this function as l ( z | ⋅) . Note that l ( z | ⋅) is not, in general, a probability distribution on S . That is,
∫
S
l ( z | θ) d θ
may or may not equal 1. We use the notation l ( z | ⋅) rather than a notation like p ( z | θ ) to emphasize this possibility. If we receive a measurement Z = z , we compute the posterior distribution on Θ by
p1 (θ | Z = z ) =
l ( z | θ ) p0 (θ )
∫
S
l ( z | θ ′) p0 (θ ′)d θ ′
for θ ∈ S
(3)
In [7] we demonstrated how to characterize quantized measurements in terms of standard Bayesian likelihood functions. Briefly the method is as follows.
incorporation of the measurement at time t ′ . The posterior in (7) is typically resampled before the motion update is performed. If desired, other proposal distributions can used to improve particle filter performance as discussed in [11].
For convenience we represent the quantized measurements by a set I of integers, i.e., Z = I . A measurement Z = i for i ∈ I means that that the parameter θ is in a subset Si ⊂ S . Each set Si can be composed of the union of disjoint subsets, for example disjoint intervals of the real numbers. The likelihood function for this quantized measurement is simply
In the following example we consider an underwater acoustic detection situation that goes beyond the quantized measurement examples considered in [5] and [6]. In particular, the bins are unions of disjoint intervals. In this example we employ a particle filter as described above to track a maneuvering target.
⎧1 if θ ∈ Si ⎪ l (i | θ ) = ⎪ ⎨ ⎪ ⎪ ⎩0 otherwise.
A. Likelihood Function for Acoustic Detection When a passive acoustic sensor is located in a deep water region of the ocean, the sound propagation conditions often produce detection regions that are disjoint. For example, there may be good detection conditions from the sensor’s location out to range 5 nm. This is typically called the direct path region. In addition there are often convergence zone regions at ranges of roughly 30 nm, 60 nm, and even farther out. A convergence zone is a region where the acoustic rays converge and produce low propagation loss and increased detection probability for the sensor. Suppose the convergence zones are 5 nm wide. It is often the case that the uncertainty about the source level of a potential target means that although one cannot calculate the detection probability as a function of range, one does know that if a target is detected it is in one of these regions. In this case, a detection means that the target is in one of these regions, i.e., its range is in the set S d where
III.
QUANTIZED MEASUREMENTS
(4)
As mentioned above, l (i | ⋅) is not a probability density, and it may not integrate to 1. In fact
∫
S
l (i | θ ) d θ = ∫ 1 d θ .
(5)
Si
It is possible to account for noise in (4) as was done in [7]. IV.
TRACKING A MANUEVERING TARGET USING QUANTIZED MEASURMENTS
One can track moving targets using sensors with quantized measurements. Since the measurements do not satisfy the linear Gaussian assumptions required for a Kalman filter, we perform the tracking using a particle filter as described in [11] or Chapter 3 of [9]. The general procedure is straightforward. The particles are motion updated to the time of a measurement. The likelihood function for a quantized measurement is applied to the weight of each particle to produce the posterior distribution on target state. The particles are resampled and then motion updated to the time of the next measurement. Specifically, suppose the target state distribution at time t is represented by the set of particles
{( x j (t ), w j (t )} for
j = 1,…, J
(6)
where x j (t ) is the state of the j th particle at time t and w j (t ) is its probability. This set of particles represents a discrete probability approximation to the distribution on target state at the time t . Suppose we obtain a quantized measurement Z = z . Let l ( z | ⋅) be the likelihood function for this measurement. The posterior distribution on target state at time t is given by
{( x j (t ), w j (t )} for
j = 1,…, J
(7)
where
l ( z | x j (t )) w j (t )
j
w (t ) =
∑
J
(
)
l z | x j ′ (t ) w j ′ (t ) j ′=1
.
(8)
If the next measurement is received at time t ′ > t , the posterior particle filter representation in (7) can be motionupdated to the time t ′ to act as a proposal distribution for the
S d = [0,5] ∪ [27.5,32.5] ∪ [57.5, 62.5] .
(9)
Generally one does not know the edges of the regions exactly. Depending on the source level of the target and the level of ambient noise in the ocean, these regions can be a bit larger or smaller than the nominal numbers in (9). Let εi be mutually independent normally distributed random variables with mean 0 and variance σi2 for i = 1, 2,3 . These are added to the nominal edges of the regions to account for this uncertainty. We denote a measurement or detection by Z = 1. Let r denote the range of the target from the sensor. The likelihood function ld for a detection is
ld (1| r ) = Pr {Z = 1| target at range r } ⎧r ≤ 5 + ε1 or 27.5 − ε2 ≤ r ≤ 32.5 + ε2 ⎪ ⎫ ⎪ ⎪ = Pr ⎪ (10) ⎨ ⎬ ⎪ ⎪ or 57.5 62.5 r ε ε − ≤ ≤ + 3 3 ⎪ ⎪ ⎩ ⎭ ⎧ ⎪ Pr {r ≤ 5 + ε1 } for 0 ≤ r ≤ 20 ⎪ ⎪ ⎪ ≈ ⎨Pr {27.5 − ε2 ≤ r ≤ 32.5 + ε2 } for 20 < r ≤ 45 ⎪ ⎪ ⎪ ⎪ ⎩Pr {57.5 − ε3 ≤ r ≤ 62.5 + ε3 } for 45 < r < ∞ The approximation in the last line of (10) is essentially an equality if σi2 < 4 for i = 1, 2,3. Let η (⋅, 0, σ 2 ) be the probability density function for a Gaussian distribution with mean 0 and variance σ 2 . The notation η (⋅, 0, σ 2 ) is used to indicate the function of one variable obtained by fixing the values of the 2nd and 3rd variables at 0 and σ 2 . Let
Φ ( z, σ 2 ) = ∫
z
−∞
η ( y, 0, σ 2 ) dy for −∞ < z < ∞ . (11)
Then the likelihood function in (10) becomes
⎧⎪1−Φ(r − 5, σ12 ) for 0 ≤ r ≤ 20 ⎪⎪ ⎪⎪ 2 2 ⎪⎪min {1−Φ(27.5 − r , σ2 ),1−Φ(r − 32.5, σ2 )} ⎪ (12) ld (1| r ) = ⎪⎨ for 20 < r ≤ 45 ⎪⎪ 2 2 ⎪⎪min {1−Φ(57.5 − r , σ3 ),1−Φ(r − 62.5, σ3 )} ⎪⎪ ⎪⎪ for 45 < r < ∞. ⎪⎩ Figure 1 shows the likelihood function ld (1| ⋅) when σi2 = 0.5 for i = 1, 2,3. 1
0.8
Probability
0.6
0.4
0.2
0 0
10
20
30 40 Range (NM)
50
60
The particles change velocity at exponentially distributed times with mean 1 hour. When a velocity change takes place, a new velocity is chosen from a distribution that produces a mean change of 30 degrees in course and 2 kn in speed. More details on this motion model are given in Section 1.3.3 of [9]. This a reasonable motion model for a submarine. Note that the velocity distribution is not Gaussian and that the motion model allows for maneuvers. This motion model is a generalization of the random tour motion model given in [12]. C. Sensor Field and Target Track For this example we placed four passive acoustic sensors in the corners of the square region in Figure 2. The red line shows the track of the target used to generate sensor measurements. The green areas indicate regions where only one sensor can detect the target. In the red areas two sensors can detect the target. Sensor measurements are received every 5 minutes. If the target is in one of the green regions at the time of a measurement, it produces a bearing with error according to the Gaussian density in (13). The likelihood function in (14) is used to update the weights on the particles. If the target is in a red region at the time of a measurement, a bearing is produced by two sensors. The update in (8) is performed by setting the likelihood function l ( z | ⋅) equal to the product of the likelihood functions corresponding to the two detections.
70
Fig. 1. Likelihood function for accoustic detection in a convergence zone environment.
When a detection occurs in the example considered below, the sensor provides a bearing estimate β . The error in this estimate is assumed to be Gaussian with mean 0 and standard The likelihood function for this deviation σb = 5º . measurement is ⎛ ⎞ ⎛ x − x2s ⎟⎞ ⎟ , 0, σb2 ⎟⎟⎟ lb (β | ( x1 , x2 )) = η ⎜⎜⎜β − arctan ⎜⎜⎜ 2 s ⎟ ⎜⎝ x1 − x1 ⎟⎠ ⎟⎠ ⎜⎝
(13)
where ( x1 , x2 ) is the position of a particle and ( x1s , x2s ) is the position of the sensor making the detection. The combined likelihood function for the detection and bearing estimate is lc ((1, θ ) | ( x1 , x2 )) = lb (θ | ( x1 , x2 )) ld (1| r ( x1 , x2 ))
(14)
where r ( x1 , x2 ) is the distance from the sensor to ( x1 , x2 ) . B. Target Motion Model The target’s initial position is uniform over the 100 nm by 100 nm square region shown in Figure 2. We used 25,000 particles to represent possible target paths. Since the initial location distribution is uniform, we initialized each particle’s position by making a draw from the position distribution resulting from the likelihood function for the first detection. Each particle’s initial speed was drawn from a uniform distribution on the interval 2 to 20 kn. Its initial course was drawn from a uniform distribution over the interval 0 to 360°.
Fig. 2. Sensor field and target track. In green areas, one sensor can detect the target. In red areas, two sensors can detect. The sensors are located in the corners of the 100 nm by 100 nm square.
D. Results Figures 3-6 show the distribution on target position at times 0.50, 1.25, 3.75, and 4.50 hrs. A sample of 100 of the 25,000 particles is shown. In Figures 3 and 4, the target has been detected only by the sensor in the lower left corner of the sensor field. Figure 3 shows three modes corresponding to the three possible detection regions for this sensor. In Figure 4 at 1.25 hrs, the target is again detected by the original sensor, but because of the speed constraints, the detection zone farthest
70 Target 60
50
40
nm
from the sensors is almost incompatible with past detections and there are few particles in this zone. Thus there are essentially two modes in the distribution. Between Figure 4 and Figure 5, the target maneuvers and remains undetected. When detection occurs (on the original sensor) at 3.75 hrs, we see three modes emerge once again. This happens because it is possible for particles (and the target) to reach all detection zones during this interval of non-detection. At 4.5 hours the target is detected by both sensors at the top of the sensor field which resolves the target distribution to a single mode.
30
20
This example shows that one can use standard Bayesian likelihood functions and a particle filter to track a maneuvering target with quantized measurements. In addition it shows that one can combine quantized and continuous measurements in an easy and correct fashion using likelihood functions. In the simulation used to produce the measurements for this example, detection occurs with probability 1 in the detection regions. In the multiple target example given in Section V, we will use detection probabilities less than 1 in the simulation that produces the measurements.
10
0
-10 -10
0
10
20
30 nm
40
50
60
70
60
70
Fig. 5. Position distribution at 3.75 hrs. 70 Target
70 Target
60
60 50
50 40
nm
nm
40 30
30 20
20 10
10 0
0 -10 -10
-10 -10
0
10
20
30 nm
40
50
60
70
0
10
20
30 nm
40
50
Fig. 6. Position distribution at 4.5 hrs.
Fig. 3. Postion distribution at 0.5 hrs.
V.
70 Target
We consider a second example where there are two targets moving at constant velocity and crossing paths as shown in Figure 7. The area of interest is a 140 nm by 140 nm square with six sensors located in the centers of the six light blue disks in Figure 7. The color scheme indicates the number of sensors that can detect a target in the colored regions with dark blue = 0, light blue = 1, green = 2, and orange = 3. The sensor characteristics and likelihood functions are the same as those in the previous example except that assume detection probability is less than one. We use the same motion model as in the maneuvering target example (with exception that mean time to velocity change is 2 hours) to generate the 25,000 particle paths used for this example.
60
50
nm
40
30
20
10
0
-10 -10
TRACKING MULTIPLE TARGETS WITH QUANTIZED MEASUREMENTS
0
10
20
30 nm
40
50
Fig. 4. Postion distribution at 1.25 hrs.
60
70
i for i = 1, 2. In the particle filter algorithm for JPDA given in Section 4.5.5 of [9], we set Pd ( n | xnj (tk )) = pd ld (1| ρ ( xnj (tk )) where ρ ( x) is the range from the state x to the sensor reporting the detection and set f (β | n, xnj (tk )) = lb (β | xnj (tk )) .
(16)
Since there are no false alarms, we take g 0 ( yk ) = 1 . In the case where yk = (β1 , β2 ) , there are two possible association hypotheses γ1 and γ 2 . For γ1 , β1 is associated to target 1 and β2 is associated to target 2. For γ 2 , the reverse holds. Finally we set the prior association probabilities Pr{γ i } = 0.5 for i = 1, 2 in both cases. With the above identifications, we employed the JPDA particle filter algorithm in Section 4.5.5 to compute the results shown below. Fig. 7. Sensor field for multiple target example
Each of the six sensors sends a report every 6 minutes, but the reports are staggered so that only one sensor reports each minute. To generate detections, the simulation uses the probability of detection versus range function obtained by multiplying the likelihood function in Figure 1 by 0.7. Detections are independent from sensor to sensor and target to target. Bearing measurements and errors are simulated as in the first example. There are no false alarms. A. JPDA Particle Filter Algorithm for Two Targets For this example, we employed a particle filter implementation of the non-linear Joint Probabilistic Data Association (JPDA) algorithm described in Section 4.5 of [9]. We assumed the number of targets (two) is known to the tracker. At each minute the report consists of 0, 1, or 2 detections. With each detection, a bearing is reported. The distribution on the state of target n is represented by J = 25, 000 particles
{( xnj (t ), w nj (t )} for
j = 1,…, J
(15)
where xnj (t ) is the state of the j th particle at time t and w nj (t ) is its probability. When detections are reported there are one or two detections in the scan. Thus the report at time tk is either yk = β1 or yk = (β1 , β2 ) where βi is the bearing resulting from the ith detection. The tracker assumes that both targets have the same unknown detection probability Pd (r ) = pd ld (1| r ) as a function of range r where pd = 0.5 . This value was chosen to represent ignorance of the detection probability. It would be better to replace this with a distribution, perhaps uniform on (0,1), or to perform a sensitivity check by assuming a number of different detection probabilities and comparing the results. In the case where yk = β1 , there are two scan association hypotheses, γ1 and γ 2 where γ i assigns the detection to target
B. Results Figures 8-19 show the results of the simulation. Figure 8 shows the probability that β1 is associated to target 1 for each detection as calculated in the JPDA particle filter algorithm in Section 4.5.5 of [9]. The bearings produced by target 1 are labeled with x’s and those by target 2 with o’s. There are no detections in the period starting after 3 hours to just before 5 hours. At 4.5 hours, the targets cross. At 5 hours the association probabilities are close to 0.5 and remain so until about 6.5 hours. Before 3 hours, when the targets are well separated, the bearings from target 1 have probability 1 of being associated with target 1 while those from target 2 have probability 0. As the targets come closer, the association probabilities come close to 0.5. After the targets cross and separate, the association probabilities return to 1 or 0. Figures 9, 11, 12, 14, 16, and 18 show the marginal distribution on the position of targets 1 and 2 at 0.4, 1.7, 3.3, 4.8, 5.2, and 6.6 hours. The track of target 1 is shown in red, target 2 in blue. A sample of the 100 of the particle positions is plotted for each target. Red circles are from target 1, blue from target 2. Figures 10, 13, 15, 17, and 19 show the velocity distributions at the above times. At 0.4 hours (Fig. 9), the targets are well separated but the distributions are multimodal. Figure 10 shows there is little velocity information. In Figure 11, the position distributions are unimodal and separated. In Figure 12 the targets are coming close together. Their position distributions are over lapping and the association probabilities (Figure 8) are close to 0.5. Figure 13 shows the velocity distributions at 3.3 hours. The distribution for target 1 has velocities predominately heading downward with a few scattered points heading upwards. For target 2, the reverse holds. This separation in velocity will allow the tracker to maintain track through the crossing. Figures 14 and 15 show the position and velocity distributions at 4.8 hours. They show a lot of mixing of the two targets and substantial spread because of the lack of detections.
1
80
Tgt1 Obs Tgt2 Obs
Tgt1 Tgt2 Tgt1 Tgt2
0.9
70 0.8
60 0.7
nm
40
0.5
30
0.4
20
0.3
10
0.2
0.1
0
0 -20
0
1
2
3
4
5
6
7
8
-10
0
10
20
30 nm
40
50
60
70
80
Fig. 11. Postion distribution at 1.7 hrs. Detection on target 1
9
Time (hrs)
Fig. 8. Association probability for target 1. 80 Tgt1 Tgt2 Tgt1 Tgt2
70
80 Tgt1 Tgt2 Tgt1 Tgt2
70
60
60
50
40
40
30 30
20 20
10 10
0 -20
0 -30
-20
-10
0
10
20 nm
30
40
50
60
70
-10
0
10
20
30 nm
40
50
60
70
80
80
Fig. 12. Position distributions at 3.3 hrs. Detection on target 2.
Fig. 9. Position distributions at 0.4 hrs. Detection on taraget 2.
30 30
Tgt1 Tgt2
Tgt1 Tgt2
20 20
10 10
kn
nm
nm
50
kn
Association probability
50 0.6
0
-10
-10
-20
-20
-30 -30
0
-20
-10
0 kn
10
20
30
Fig. 10. Velocity distributions at 0.4 hrs. Detection on target 2.
-30 -30
-20
-10
0 kn
10
20
30
Fig. 13. Velocity distributions at 3.3 hrs. Detection on target 2.
hours, Figures 18 and 19 show that the position and velocity distributions have separated. Figure 8 shows the association probabilities at this time are again close to either 0 or 1.
80 Tgt1 Tgt2 Tgt1 Tgt2
70
30
60
Tgt1 Tgt2
50
nm
20
40 10
30
kn
20
0
10 -10
0 -20
-10
0
10
20
30 nm
40
50
60
70
80
Fig. 14. Position distributions at 4.8 hrs. No detection
-20
30 Tgt1 Tgt2
-30 -30
-20
-10
0 kn
10
20
30
Fig. 17. Velocity distributions at 5.2 hours. Detection on target 1
20
80 Tgt1 Tgt2 Tgt1 Tgt2
10
kn
70
0
60
50
nm
-10
-20
40
30
20
-30 -30
-20
-10
0 kn
10
20
30 10
Fig. 15. Velocity distirbutions at 4.8 hrs. No detection 0 -20
-10
0
10
20
80 Tgt1 Tgt2 Tgt1 Tgt2
70
30 nm
40
50
60
70
80
Fig. 18. Postion distributions at 6.6 hrs. Detection on target 1. 30
60
Tgt1 Tgt2
nm
50
20
40
10 30
kn
20
0
10
0 -20
-10 -10
0
10
20
30 nm
40
50
60
70
80
Fig. 16. Position distributions at 5.2 hrs. Detection on target 1.
Figures 16 and 17 show the position and velocity distributions at 5.2 hours. The targets are still close, so the association probabilities are near 0.5, and the position distributions overlap. The velocity distributions have maintained their separation but with more scatter. At 6.6
-20
-30 -30
-20
-10
0 kn
10
20
30
Fig. 19. Velocity distributions at 6.6 hrs. Dectection of target 1
C. Extension to Feature Information We now add feature (non-kinematic) information to the problem. In addition to measuring the bearing to the target when a detection occurs, the sensor also measures the frequency of the acoustic signal. Specifically, suppose that each target radiates acoustic energy at a single frequency ω1 for target 1 and ω2 for target 2.
estimates so there is continuous bearing information as well as quantized range information. There are multiple sensors to deal with and the targets not only move, but they cross paths. In spite of these complications, the use of standard Bayesian likelihood functions to represent the measurements and a nonlinear particle filter version of JPDA successfully tracked the targets.
To incorporate this information, we expand the target state space to include the (constant) frequency ω radiated by the target and assume a prior distribution on ω that is uniform over [990hz,1010hz] . Frequency measurements are assumed to have independent Gaussian errors with mean 0 and standard deviation σ f = 2hz . The likelihood function for a frequency measurement ϕ is
The extension of the example, shows how conceptually simple it is to add non-kinematic information. This illustrates the power of likelihood functions. Once one has converted the measurements into likelihood functions and combined them as above, the tracking process become conceptually easy and straight-forward. Special assumptions and probability models are not required.
l f (ϕ | ω ) = n (ϕ − ω , 0, σ f ) ,
VI.
and the function f in (16) becomes f ((β , ϕ ) | n, xnj (tk )) = lb (β | xnj (tk )) l f (ϕ | ω ) .
(17)
The above example was rerun using this additional information with ω1 = 995hz and ω2 = 1005hz . For this run, we used the JPDA algorithm described above with the definition of f given in (17) and the extended state space that includes frequency. Figure 20 shows the association probabilities for target 1 for each detection. One can see that with the additional information, the association probabilities are almost always equal to 1 or 0. This means that even if the two targets were to be “reflected” off one another when they meet, we would still be able maintain tracks on the targets and not confuse them. 1 Tgt1 Obs Tgt2 Obs
0.9
0.8
Association probability
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0
1
2
3
4 5 Time (hrs)
6
7
8
9
Fig. 20. Association probabilitiy for target 1 with frequency information
D. Summary of Example This example shows that multiple targets can be successfully tracked in situations with complex quantized measurements. In the example, detections produce bearing
CONCLUSIONS
Reference [7] showed that the quantized and imprecise measurement cases considered by [5] and [6] can be handled with standard Bayesian likelihood functions. This paper continues by showing that standard Bayesian likelihood functions coupled with particle filters can easily and naturally handle maneuvering and multiple target tracking problems even when the measurements are complex quantized measurements coupled with continuous measurements. REFERENCES [1]
R. E. Curry, W. E. Vander Velde, and J. E. Potter, “Nonlinear estimation with quantized measurements – PCM, predictive quantization, and data compression” IEEE Trans Information Theory, Vol IT-16, No 2, March 1970, pp. 152-161. [2] Y. Ruan, P. Willet, and Alan Marrs, “Fusion of quantized measurements via particle filtering” in Proceedings IEEE Aerospace Conference, Big Sky Montana, March 2003, Vol 4 pp 1967 – 1978. [3] Z. Duan, V. P. Jilkov, and X. R. Li, “State estimation with quantized measurements: approximate MMSE approach” in Proceedings of the 11th International Conference on Information Fusion, Cologne Germany, June 30 – July 3, 2008. [4] Z. Duan, X. R. Li, and V. P. Jilkov, “State estimation with point and set measurements” Proceedings of the 13th International Conference on Information Fusion, Ediburgh, UK, July 26 – 29, 2010. [5] R. Mahler, “General Bayes Filtering of Quantized Measurements” Proceedings of the 14th International Conference on Information Fusion, Chicago, USA, July 5-8, 2011 [6] B. Ristic, Bayesian estimation with imprecise likelihoods: random set approach” IEEE Signal Processing Letters Vol 18, No 7 July 2011, pp. 395-398. [7] L. D. Stone, “Standard Bayesian approach to quantized measurements and imprecise likelihoods” Proceedings of the 16th International Conference on Information Fusion, Istanbul, Turkey, July 9-12, 2013 [8] A. Gning, S. Julier, and L. Mihaylova, “Non-linear state estimtation using imprecise samples” Proceedings of the 16th International Conference on Information Fusion, Istanbul, Turkey, July 9-12, 2013. [9] L. D. Stone, R. L. Streit, T. L. Corwin, and K. L. Bell, Bayesian Multiple Target Tracking 2ed, Boston: Artech House 2014. [10] J. Berger, Statistical Decsion Theory and Bayesian Analysis, 2nd ed, New York, NY: Springer-Verlag, 1985. [11] B. Ristic, S. Arulampalam, and N. Gordon, Beyond the Kalman Filter. Boston: Artech House, 2004. [12] A. Washburn, “Probabilty density of a moving particle,” Operations Research Vol 17, September/October 1969, p. 861.