1 Introduction↩︎

Sensor networks are becoming increasingly important in both military and civilian applications[1], [2]. Multi-target multi-sensor tracking in high-dimensional state spaces is still a challenging computational problem in term of time and resources. When a sensor network attempts to estimate the location of a set of targets at time \(t\), there are two sources of information for the targets’ positions: the sensor measurements at time \(t\); and the prior distribution propagated from time \(t-1\). Classical tracking algorithms (such as the Kalman filter) assume that all probability distributions involved are Gaussian. However, non-directional sensor measurements give rise to probability distributions that are highly non-Gaussian. Past research has employed various particle filters to deal with the non-Gaussian nature of the distributions [3]. However, particle filters typically require high computational complexity to achieve accurate results.

Our research takes a different approach. We take advantage of two facts: (1) typically an analytical expression is available for the measurement distribution; and (2) because of the uncertainty in propagation, the propagated prior is nearly Gaussian. In our tracking filter we use an analytical expression for the log likelihood of the product of measurement distribution \(\times\) Gaussian approximation of the prior, then subsequently refine the approximation by computing the mean and covariance of the product distribution numerically. We also make use of techniques to avoid convergence to a non-global local minimum.

In the following sections we first describe the system model used to develop and test the algorithm; we then outline the steps of the algorithm

2 System model description↩︎

The system model used was taken verbatim from Li and Coates [3] (and their Matlab code was used for model testing [4]). The system consists of a rectangular array of \(S=25\) sensors with a grid spacing of 10 m. Within the region are \(C=4\) targets moving stochastically (and independently) as governed by the equations: \[{\text{\boldsymbol{x}}}_{c,t} := F_{prop}\text{\boldsymbol{x}}_{c,t-1} + \boldsymbol{\nu}_{c,t}\qquad c= 1 \ldots C; ~ t = 1 \ldots T,\] Where \({\text{\boldsymbol{x}}}_{c,t} := [\vec{x}_{c,t},\dot{\vec{x}}_{c,t}]\) is the position+velocity vector of target \(c\) at discrete time index \(t\); \(F_{prop}\) is a transition matrix given by: \[\label{eq:Fprop} F_{prop} := \begin{bmatrix} 1 & 0 & 1 & 0 \\ 0 & 1 & 0 &1 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix},\tag{1}\]

and \(\boldsymbol{\nu}_{c,t} \sim N(0,V)\) is process noise with covariance matrix \(V\) where \[V := \sigma_w^2\begin{bmatrix} 1/3 & 0 & 0.5 & 0 \\ 0 & 1/3 & 0 & 0.5 \\ 0.5 & 0 & 1 & 0 \\ 0 & 0.5 & 0 & 1 \end{bmatrix}, \quad (\sigma_w^2 = 0.05).\] In the following we will often drop the \(t\) index for brevity. The expected measurement signal at sensor \(s\) due to target \(c\) is distance-dependent, and is specified in terms of the quantity \[\label{eq:fsc} f_{s,c} := \frac{A}{||\vec{x}_c - \vec{x}_{[s]}||^p + d_0}, \qquad c= 1 \ldots C, s = 1 \ldots S,\tag{2}\] where \(\vec{x}_{[s]} = (x_{[s]},y_{[s]})\) is the \((x,y)\) position of sensor \(s\), and \(p=1\) was used. The quantity \(f_{s,c}\) expresses the expected signal at sensor \(s\) due to sensor \(c\). Altogether the expected signal at sensor \(s\) given target positions at \(\vec{\boldsymbol{x}}:= [ \vec{x}_1, \ldots \vec{x}_C]\) is: \[\label{eq:Fs} \alpha_s := E[a_s | \vec{\boldsymbol{x}}] = \sum_{c=1}^C f_{s,c} ~~ (s=1 \ldots S).\tag{3}\] The measurement signal \(a_s\) at sensor \(s\) is assumed to have variance \(\sigma_s^2\).

The motion process covariance is not precisely known by the observer, and the filter assumes the following covariance: \[\label{eq:varEst} V' := \begin{bmatrix} 3 & 0 & 0.1 & 0 \\ 0 & 3 & 0 & 0.1 \\ 0.1 & 0 & 0.03 & 0 \\ 0 & 0.1 & 0 & 0.03 \end{bmatrix}.\tag{4}\]

3 Algorithm description↩︎

The steps in our multi-target estimation procedure may be outlined as follows:

(A) Analytically compute the negative log likelihoods (NLL) for the measurement distribution and a Gaussian approximation to the propagated prior together with their gradients and Hessians;

(B) Use the analytical expressions and Hessian minimization of the combined NLL to estimate the maximum likelihood (ML) target vector;

(C) If the NLL evaluated at the estimated target vector is above a threshold, then re-estimate the ML target vector;

(D) If the Hessian of the NLL at the ML vector is not positive definite, then perform modifications to restore positive definiteness;

(E) Generate integration points, using the inverse of the Hessian as a preliminary covariance estimate and making a nonlinear transformation for targets that are close to sensors;

(F) Use the integration points to estimate the mean and covariance for a Gaussian approximation of the new prior;

We fill in the mathematical details for steps (A-F) in the following subsections.

3.1 Calculation of negative log likelihoods, gradients, and Hessians↩︎

3.1.1 Negative log likelihood of measurement distribution↩︎

We use \(a_s\) to denote the received signal at sensor \(s\) at time step \(t\). Based on 2 , the NLL for the measurement distribution \({\mathcal{N}}^{(meas)}(\vec{\boldsymbol{x}})\) is given by \[\label{eq:nu} {\mathcal{N}}^{(meas)}(\vec{\boldsymbol{x}})= \sum_{s=1}^S \frac{(\alpha_s - a_s)^2}{2\sigma_s^2}\tag{5}\] To find the optimum value (minimum) of the total NLL, we must find the gradient of \({\mathcal{N}}^{(meas)}(\vec{\boldsymbol{x}})\) with respect to \([ \vec{x}_1, \ldots \vec{x}_C]\) . Based on 3 5 , the gradient with respect to \(\vec{x}_c\) may be given in terms of the gradients of \(f_{s,c}\): \[\label{eq:gradNmeas} \vec{\nabla}_{\vec{x}_c} {\mathcal{N}}^{(meas)}(\vec{\boldsymbol{x}})= \sum_{s=1}^S \sum_{c=1}^C \left( \frac{\alpha_s - a_s}{\sigma_s^2} \right) \vec{\nabla}_{\vec{x}_c} f_{s,c}.\tag{6}\] Note that \(\vec{\nabla}_{\vec{x}_c} f_{s,c} = \vec{\nabla}_{\vec{r}_{c,s}} f_{s,c}\) where \(\vec{r}_{c,s} \equiv \vec{x}_c - \vec{x}_{[s]}\). It follows that \[\begin{align} \vec{\nabla}_{\vec{x}_c} f_{s,c} &= \frac{-(p/2)A(||\vec{r}_{c,s}||^2)^{p/2-1} }{(||\vec{r}_{c,s}||^p + d_0)^2} \vec{r}_{c,s} \\ & = \frac{-(p/2)A(||\vec{r}_{c,s}||^2)^{p/2-1} }{(||\vec{r}_{c,s}||^p + d_0)^2} \vec{r}_{c,s} \\ & \equiv g(||\vec{r}_{c,s}||^2) \vec{r}_{c,s} \end{align}\] By taking another derivative, we may find the Hessian \(H {\mathcal{N}}^{(meas)}\): \[\label{eq:HNmeas} \begin{align} H {\mathcal{N}}^{(meas)}(\vec{\boldsymbol{x}})= &\sum_{s=1}^S \sum_{c=1}^C \left( \frac{1}{\sigma_s^2} \right) (\vec{\nabla}_{\vec{x}_c} f_{s,c}) \vec{\nabla}_{\vec{x}_c} f_{s,c}^T + \sum_{s=1}^S \sum_{c=1}^C \left( \frac{\alpha_s - a_s}{\sigma_s^2} \right) \vec{\nabla}_{\vec{x}_c} ( \vec{\nabla}_{\vec{x}_c} f_{s,c}), \end{align}\tag{7}\] where (writing \(g(||\vec{r}_{c,s}||^2)\) as \(g\) for brevity) \[\begin{align} \vec{\nabla}_{\vec{x}_c}( \vec{\nabla}_{\vec{x}_c} f_{s,c} ) &= \vec{\nabla}_{\vec{x}_c}( g \cdot \vec{r}_{c,s} ) \\ & = \frac{dg~~~~~~~}{d (||\vec{r}_{c,s}||^2)} (2 \vec{r}_{c,s} \vec{r}_{c,s}^{\,T}) + g \cdot I \\ &= 2 g \cdot \left[ \left(\frac{p/2 - 1}{||\vec{r}_{c,s}||^2} - \frac{2}{||\vec{r}_{c,s}||^p + d_0} \right)(\vec{r}_{c,s}^{\,T} \vec{r}_{c,s} ) + I \right] \\ \end{align}\]

3.1.2 NLL of propagated prior distribution↩︎

The analytical expression for the NLL of the propagated prior may be computed as follows. We assume that the prior distribution is joint Gaussian with mean \((\vec{m}_{\vec{\boldsymbol{x}},t-1} ; \vec{m}_{\vec{\boldsymbol{v}},t-1})\) and covariance \(\Sigma_{t-1}\), so that \[p_{t-1}(\vec{\boldsymbol{x}},\vec{\boldsymbol{v}}) \propto \exp \left(-(\vec{\boldsymbol{x}}- \vec{m}_{\vec{\boldsymbol{x}},t-1};\vec{\boldsymbol{v}}- \vec{m}_{\vec{\boldsymbol{v}},t-1})^T \cdot \Sigma_{t-1}^{-1} (\vec{\boldsymbol{x}}- \vec{m}_{\vec{\boldsymbol{x}},t-1};\vec{\boldsymbol{v}}- \vec{m}_{\vec{\boldsymbol{v}},t-1}) / 2 \right),\] \((\vec{m}_{\vec{\boldsymbol{x}},t-1} ; \vec{m}_{\vec{\boldsymbol{v}},t-1})\) and \(\Sigma_{t-1}\) are the prior’s mean vector and covariance matrix, respectively. This distribution is propagated according to the motion model 1 , except that the actual process noise covariance is replaced with the estimated covariance \(V'\) defined in 4 . It follows that the propagated spatial distribution from time \(t-1\) (i.e., the prior distribution at time \(t\)) is also Gaussian with mean \[\label{eq:propMean} (\vec{m}^{(prop)}_{\vec{\boldsymbol{x}}} ; \vec{m}^{(prop)}_{\vec{\boldsymbol{v}}}) \equiv F_{prop} (\vec{m}_{\vec{\boldsymbol{x}},t-1}; \vec{m}_{\vec{\boldsymbol{v}},t-1})\tag{8}\] and covariance \[\label{eq:propCov} \Sigma_{prop} \equiv F_{prop} \Sigma_{t-1} F_{prop}^T + V'.\tag{9}\] It follows that the NLL of propagated distribution is \({\mathcal{N}}^{(prop)}\), where \[\label{eq:Nprop} {\mathcal{N}}^{(prop)}(\vec{\boldsymbol{x}}) = \frac{1}{2}(\vec{\boldsymbol{x}}- \vec{m}^{(prop)}_{\vec{\boldsymbol{x}}})^T \Sigma_{prop,\vec{\boldsymbol{x}}\vec{\boldsymbol{x}}}^{-1}(\vec{\boldsymbol{x}}- \vec{m}^{(prop)}_{\vec{\boldsymbol{x}}})\tag{10}\] The gradient is \[\label{eq:gradNprop} \vec{\nabla} {\mathcal{N}}^{(prop)}(\vec{\boldsymbol{x}}) = \Sigma_{prop,\vec{\boldsymbol{x}}\vec{\boldsymbol{x}}}^{-1}(\vec{\boldsymbol{x}}- \vec{m}^{(prop)}_{\vec{x}})\tag{11}\] and \[\label{eq:HNprop} H {\mathcal{N}}^{(prop)}(\vec{\boldsymbol{x}}) = \Sigma_{prop,\vec{\boldsymbol{x}}\vec{\boldsymbol{x}}}^{-1}.\tag{12}\]

3.2 Estimation of maximum likelihood target vector (spatial components)↩︎

The minimizer of \({\mathcal{N}} \equiv {\mathcal{N}}^{(meas)} + {\mathcal{N}}^{(prop)}\) corresponds to the ML estimator of the joint location vector for the targets. The Matlab routine “fmincon” was used for the minimization. Constraints were added to require the targets to remain within the grid region. In the minimization, the analytical gradient and Hessian were used, where: \[\vec{\nabla}{\mathcal{N}}(\vec{\boldsymbol{x}}) = \vec{\nabla}{\mathcal{N}}^{(meas)}(\vec{\boldsymbol{x}}) + \vec{\nabla}{\mathcal{N}}^{(prop)}(\vec{\boldsymbol{x}}) \text{ and } H{\mathcal{N}}(\vec{\boldsymbol{x}}) = H{\mathcal{N}}^{(meas)}(\vec{\boldsymbol{x}}) + H{\mathcal{N}}^{(prop)}(\vec{\boldsymbol{x}}),\] and \(\vec{\nabla}{\mathcal{N}}^{(meas)}\), \(\vec{\nabla}{\mathcal{N}}^{(prop)}\), \(H{\mathcal{N}}^{(meas)}\), and \(H{\mathcal{N}}^{(prop)}\) are given by 6 , 11 , 7 , and 12 respectively.

3.3 Consistency check and re-estimation of the target vector↩︎

It is possible that the algorithm may lose track of the targets, or may fail in initial acquisition. Thus we introduce a check for whether or not the ML target location vector is consistent with the measurement distribution. In this section the checking procedure is described, as well as the re-estimation procedure that is applied if the threshold is not met.

According to the signal model, the signal \(a_s\) at sensor \(s\) (\(s=1\ldots S\)) is Gaussian distributed with mean \(\alpha_s\) and variance \(\sigma_s\), and the signals at different sensors are independent. It follows that the the sum of squared normalized signal deviations given by \({\mathcal{N}}^{(meas)}(\vec{\boldsymbol{x}})\) in 5 has a Chi-squared distribution with \(S\) degrees of freedom. This statistic may be computed using signal data and the ML target position vector, and consistency is rejected if the signal exceeds a threshold computed using the known distribution of \(\chi^2_S\) and a user-specified \(p\)-value (in our simulations, we used \(p = 0.0013\), corresponding to 3 standard deviations).

In case the signal does not meet the threshold, a correction procedure is applied to find a target position vector with higher likelihood. Two different correction procedures were developed, as described below.

3.3.1 One-by-one sensor addition correction procedure↩︎

The procedure is based on the observation that when signals from too many sensors are combined, the resulting function has many suboptimal local minima (during optimization, it is easy for the target estimate to get caught in the middle of one of the squares in the grid formed by the sensors). On the other hand, when some sensors are excluded and only sensors that are far from the targets are used, then local minima are avoided but the accuracy of the target position estimates is reduced. So we begin with far-away sensors to steer the target vector to a neighborhood of the correct solution. We then add sensor signals one by one to gradually increase the accuracy without introducing nearby local minima that could prevent the solution from reaching the global optimum. The steps in the procedure may be described as follows:

1. Estimate the target vector, using only sensors on the boundary of the region;

2. Loop over remaining (unused) sensors:

3. Identify the unused sensor which maximizes the minimum distance to current estimated target positions;

4. Add this “maximin” sensor to the list of sensors used to estimate the target vector, and re-estimate the target vector;

5. end loop

3.3.2 Square-hopping correction procedure↩︎

The second correction procedure is somewhat more complicated than the first. This procedure first locates estimated target positions that should excessively large signals at some sensors (compared to the measured signals), then moves these targets near to sensors that are receiving larger signals. For this purpose, we define the “signal excess” due to target \(c\) as \(\epsilon(c)\) where \[\epsilon(c) = \sum_{s=1}^S \max( (\alpha_s - a_s) - Af_{s,c} , 0).\] \(\epsilon(c)\) can be interpreted as the sum over sensors of the reduction in (predicted\(-\)observed) signal that is achieved when the signal from target \(c\) is removed. A large value of \(\epsilon(c)\) is an indicator that target \(c\) is wrongly placed.

On the other hand, we may define the “signal deficit” at sensor \(s\) as

Using these definitions, we may outline the “square-hopping” procedure as follows:

1. Identify the set of \(n_{bad\_tgt}\) of targets associated with the largest values of \(\epsilon(c)\) (our algorithm used \(n_{bad\_tgt}=2\));

2. Remove the targets identified in step (a) and recalculate the expected signal at each sensor

3. Locate the \(n_{bad\_sq}\) grid squares where the sensors at the corners have the largest combined signal deficit (our algorithm used \(n_{bad\_sq} = 12\));

4. List the possible subsets of size \(n_{bad\_tgt}\) of the grid squares computed in step (c);

5. For each subset identified in step (d), place a target at the middle of each square and then run the minimizing algorithm to estimate the MLE target vector

6. If the NLL of the target vector identified in step (e) is below the NLL acceptance threshold, then use the resulting target vector as MLE

7. If the NLL of the target vector identified in step (e) is not below the NLL acceptance threshold, then return to step (d) end loop

3.4 Hessian correction↩︎

Although theoretically the Hessian should be positive definite at a function minimum, in practice the Hessian may fail to be positive definite because of numerical issues. This can happen when at least one target is very close to a sensor. In this case, the target in question is identified by means of the measurement signal, and rows and columns of the Hessian corresponding to that particular target are modified. A description of the Hessian correction procedure is as follows:

1. Initialize an empty list of excluded target-sensor pairs

2. Using the current target vector estimate, locate the closest target-sensor pair that is not in the list, and add to the list

3. Find the MLE target vector estimate with signals from sensors in the excluded list removed, and with targets in the excluded list fixed at their current positions (so they are not free parameters in the optimization)

4. If the Hessian for the non-excluded targets is not positive definite, then go back to step (b)

5. If the Hessian for the non-excluded targets is positive definite, add rows and columns to the Hessian corresponding to the excluded targets. For added rows and columns, all non-diagonal entries are 0, and the diagonal entry is the fixed value \(d_0^{-2}\), which corresponds to a target position variance of \(d_0^2\) (note the inverse of the Hessian is the estimated covariance for the measurement distribution).

3.5 Numerical integration of the combined distribution↩︎

3.5.1 Integral with Gaussian weight function↩︎

Using the method described above, by minimizing \({\mathcal{N}}(\vec{\boldsymbol{x}})\) we obtain an estimate \(\vec{m}_{\vec{\boldsymbol{x}}}\) for the mean of the target vector’s spatial distribution at time \(t\) (denoted by \(p(\vec{\boldsymbol{x}})\)). We may also use the Hessian \(H \equiv H{\mathcal{N}}(\vec{m}_{\vec{\boldsymbol{x}}})\) as a preliminary estimate of the inverse of the spatial covariance. At this point, our goal is to compute better estimates of the spatial mean and covariance via numerical integration. We suppose that the distribution is approximately Gaussian with mean \(\vec{m}_{\vec{\boldsymbol{x}}}\) and spatial covariance \(H^{-1}\). In order to compute moments of the distribution, we will need to compute integrals of the form: \(\int_{{\mathbb{R}}^d} f(\vec{\boldsymbol{x}}) p(\vec{\boldsymbol{x}}) d^d \vec{\boldsymbol{x}}.\) We may change variable as follows. Since \(H\) is the Hessian at a function’s minimum, it is a positive definite symmetric matrix and thus has a Cholesky square root \(H^{1/2}\). We may define: \[\label{eq:ydef} \vec{\boldsymbol{y}}(\vec{\boldsymbol{x}}) \equiv H^{1/2}(\vec{\boldsymbol{x}}- \vec{m}_{\vec{\boldsymbol{x}}}) \implies \vec{\boldsymbol{x}}(\vec{\boldsymbol{y}}) = H^{-1/2}\vec{\boldsymbol{y}}+ \vec{m}_{\vec{\boldsymbol{x}}}.\tag{13}\]

The actual distribution is modeled as a perturbation of the Gaussian distribution with mean \(\vec{m}_{\vec{\boldsymbol{x}}}\) and spatial covariance \(H^{-1}\), so the distribution function can be expressed as:

\[\label{eq:pme} p (\vec{\boldsymbol{x}}) \approx C g(\vec{\boldsymbol{y}}(\vec{\boldsymbol{x}})) \exp \left( -\frac{1}{2} \vec{\boldsymbol{y}}^T \vec{\boldsymbol{y}}\right),\tag{14}\] so that \[\label{eq:pme2} g (\vec{\boldsymbol{y}}(\vec{\boldsymbol{x}})) \approx C^{-1} p(\vec{\boldsymbol{x}}) \exp \left( \frac{1}{2} \vec{\boldsymbol{y}}^T \vec{\boldsymbol{y}}\right),\tag{15}\] where \(g(\vec{\boldsymbol{y}})\) is the perturbation function and \(C\) is a normalizing constant.

This gives us \[\int_{{\mathbb{R}}^d} f(\vec{\boldsymbol{x}}) p(\vec{\boldsymbol{x}}) d^d \vec{\boldsymbol{x}}= \sqrt{|H|} \int_{{\mathbb{R}}^d} h(\vec{\boldsymbol{y}}) \exp \left( -\frac{1}{2} \vec{\boldsymbol{y}}^T \vec{\boldsymbol{y}}\right) d^d \vec{\boldsymbol{y}}\] where \[\label{eq:h} h(\vec{\boldsymbol{y}}) \equiv C f\left(\vec{\boldsymbol{x}}(\vec{\boldsymbol{y}})\right) g(\vec{\boldsymbol{y}}).\tag{16}\] We may write the \(\vec{\boldsymbol{y}}\) integral in polar coordinates as: \[\int_0^{\infty} \left( \int_{S^7} h(r, \Theta) d \Theta \right) r^{d-1} \exp \left( -\frac{r^2}{2}\right) dr,\] and changing variable \(s \equiv r^2/2\) gives: \[\label{eq:GL95int1} 2^{d/2-1} \int_0^{\infty} \left( \int h(\sqrt{2s}, \Theta)d \Theta \right) s^{d/2-1} e^{-s} ds,\tag{17}\]

3.5.2 Application of Gauss-Laguerre quadrature↩︎

The \(s\) integral can be approximated using Gauss-Laguerre quadrature[5]. \(n\)-point Gauss-Laguerre quadrature depends on the generalized Laguerre polynomials \(L^{(d/2-1)}_{n}(z)\) and \(L^{(d/2)}_{n+1}(z)\), and will integrate exactly polynomials in \(s\) of order \(2n-1\). From the recurrence relation in [6], we find for the case \(d=8\): \[\begin{align} L^{(3)}_0(z) &= 1 \\ L^{(3)}_1(z) &= 1 + 3 -z =4-z\\ L^{(3)}_2(z) &= \frac{(3+3-z)L^{(3)}_1(z) - (1+3)L^{(3)}_0(z)}{2} \\ & = \frac{20 - 10x + z^2}{2}\\ L^{(3)}_3(z) &= \frac{(5+3-z)L^{(3)}_2(z) - (2+3)L^{(3)}_1(z)}{3} \\ & = \frac{8-z}{3}L^{(3)}_2(z) - \frac{5}{3}(4 - z) \end{align}\] For two-point Gauss-Laguerre quadrature, the quadrature points are the roots of \(L^{(3)}_2(z)\), which are \(z_\pm = 5 \pm \sqrt{5}\). The weights are \[w_\pm = \frac{(2+3)! \,z_\pm}{2!(2+1)^2 (L_3^{(3)}(z_\pm))^2},\] and noting that \(L_3^{(3)}(x_\pm) = - \frac{5}{3}(4 - x_\pm)\), we have \[w_\pm = \left(\frac{6}{5}\right) \frac{5 \pm \sqrt{5}}{3 \pm \sqrt{5}}.\] It follows that the integral in 17 can be approximated as: \[2^{d/2-1} \left( w_{-} \int h(\sqrt{2z_{-}}, \Theta)d \Theta + w_{+} \int h(\sqrt{2z_{+}}, \Theta)d \Theta \right).\]

3.5.3 Performing angular integration using zero-sum lattice points↩︎

Supposing that we are in \(d\)-dimensional space (\(d=8\) in our particular case), we need to perform the angular integral over a \((d-1)\)-dimensional sphere. This angular integral can be approximated by a sum over the mesh points of a mesh that uniformly covers the sphere. One choice of mesh is the \(d(d+1)\) unit vectors in the root system of the \(d\)-dimesional simplex (zero-sum) lattice. The zero-sum lattice may be defined as the sublattice of \({\mathbb{Z}}^{d+1}\) restricted to the plane \(x_1 + \ldots + x_d +x_{d+1}= 0\)[7]. To obtain a lattice in \({\mathbb{R}}^d\), we need to rotate this plane into the plane \(x_{d+1}=0\). The easiest way to do this is to rotate the unit vector \(\vec{1}\) in \({\mathbb{R}}^{d+1}\) so that it aligns along the elementary basis vector \(\hat{e}_{d+1}\), while leaving all orthogonal vectors fixed. We denote the desired rotation by \(\rho\). Let \(\mathcal{P}\) denote the plane determined by the two vectors \(\vec{1}\) and \(\hat{e}_{d+1}\): note that \(\mathcal{P}\) is fixed by \(\rho\), and all vectors orthogonal to \(\mathcal{P}\) are invariant under \(\rho\). Let \(\hat{\alpha} \equiv \frac{\vec{1} - \hat{e}_{d+1}}{\sqrt{d}}\), so that \(\mathcal{P}\) is spanned by the orthonormal set\(\{\hat{\alpha}, \hat{e}_{d+1}\}\) Then if \(1 \le j \le d\), we may express \(\hat{e}_j\) as the sum of two orthogonal vectors: \[\hat{e}_j = (\hat{e}_j \cdot \hat{\alpha})\hat{\alpha} + (\hat{e}_j - (\hat{e}_j \cdot \hat{\alpha})\hat{\alpha}).\] Since \(\hat{e}_j \cdot \hat{\alpha} = \frac{1}{\sqrt{d}}\) for all \(j \le d\), we have: \[\rho(\hat{e}_j) = \frac{\rho(\hat{\alpha})}{\sqrt{d}} + \rho \left(\hat{e}_j - \frac{\hat{\alpha}}{\sqrt{d}}\right) .\] But the vector \(\hat{e}_j - \frac{\hat{\alpha}}{\sqrt{d}}\) is orthogonal to \(\mathcal{P}\), and is thus fixed under \(\rho\). This implies: \[\label{rho95ej} \rho(\hat{e}_j) = \hat{e}_j + \frac{\rho(\hat{\alpha}) - \hat{\alpha}}{\sqrt{d}}.\tag{18}\] The set of closest neighbors of \(\vec{0}\) in the restricted zero-sum lattice in \({\mathbb{R}}^{d+1}\) consists of all vectors of the form \(\hat{e}_i - \hat{e}_j\) where \(1 \le i,j \le d+1\). When rotating this set into \({\mathbb{R}}^d\) we only need to compute \(\rho(\hat{e}_j - \hat{e}_{d+1})\), where \(1 \le j \le d\), since all other vectors in the set are invariant under \(\rho\). From 18 we have \[\rho(\hat{e}_j - \hat{e}_{d+1}) = \hat{e}_j + \vec{q},\] where \[\vec{q} \equiv \frac{\rho(\hat{\alpha} - \hat{e}_{d+1}) - \hat{\alpha}}{\sqrt{d}}\] Now we have that \(\vec{q} \in \mathcal{P}\) and \(\vec{q} \perp \hat{e}_{d+1}\), which implies that \(\vec{q} \parallel \hat{\alpha}\). The additional condition that \(|| \hat{e}_j + \vec{q} ||^2=2\) implies that. \[\label{eq:q} \vec{q} = \frac{-1 \pm \sqrt{d+1}}{d} (\underbrace{1,1,\ldots , 1}_{d~\textrm{terms}},0).\tag{19}\] There are two possible choices for \(\vec{q}\) (corresponding to rotating \(\vec{1}\) parallel or antiparallel to \(\hat{e}_{d+1}\)), and either may be chosen. Since the final entry of \(\vec{q}\) in 19 is 0, it follows \(\vec{q} \in {\mathbb{R}}^d\), and we may drop the final ‘0’. Consequently, we may generate the \(d(d+1)\) vectors in the root system of the simplex lattice in \({\mathbb{R}}^d\) as follows:

• Take all ordered pairs \((i,j)\) where \(1 \le i,j \le d\)   (\(d(d-1)\) vectors), and form the vectors \(\hat{e}_i - \hat{e}_j\);

• for \(1 \le j \le d\), form the vectors \(\hat{e}_j + \vec{q}\) and \(-\hat{e}_j - \vec{q}\).

We may determine weights on the grid points \(\{\vec{\Theta}_j, j = 1,\ldots J\}\) where \(J \equiv d(d+1)\) as follows. The area of a sphere in \(d\) dimensions (when \(d\) is even) is \(\frac{2\pi^{d/2}}{(d/2-1)!}\). It follows that \(\int d\Theta = \frac{2\pi^{d/2}}{(d/2-1)!}\). In order to give the correct area integral, then each point in the mesh should be assigned a weight \(\frac{2\pi^{d/2}}{(d/2-1)!(d)(d+1)}\). It follows that the integral in 17 can be approximated as: \[\label{eq:GL95int2} A \sum_{j=1}^J \left( w_{-} h(\sqrt{2z_{-}}\vec{\Theta}_j) + w_{+} h(\sqrt{2z_{+}}\vec{\Theta}_j) \right),\tag{20}\] where \[A\equiv \frac{(2\pi)^{d/2}}{(d/2-1)!(d)(d+1)}.\]

3.5.4 Numerical integration summary↩︎

Based on 20 , 14 and 16 the integral \(\int_{{\mathbb{R}}^d} f(\vec{\boldsymbol{x}}) p(\vec{\boldsymbol{x}}) d^d \vec{\boldsymbol{x}}\) can be approximated as: \[\label{eq:int2} \begin{align} \int_{{\mathbb{R}}^d} f(\vec{\boldsymbol{x}}) p(\vec{\boldsymbol{x}}) d^d \vec{\boldsymbol{x}} & \approx \sqrt{|H|} AC \sum_{j=1}^{J} \left( w_{-} h(\sqrt{2z_{-}}\vec{\Theta}_j) + w_{+} h(\sqrt{2z_{+}}\vec{\Theta}_j) \right),\\ & \approx \sqrt{|H|} AC \sum_{j=1}^{J} w_{-} e^{z_{-}} f(\vec{x}_{j}^{\,-}) p(\vec{x}_{j}^{\,-}) + w_{+} e^{z_{+}} f(\vec{\boldsymbol{x}}_{j}^{\,+}) p(\vec{\boldsymbol{x}}_{j}^{\,+}) \end{align}\tag{21}\] where \[\label{eq:ydef2} \vec{\boldsymbol{x}}_{j}^{\,\pm} \equiv \sqrt{2z_{\pm}}H^{-1/2}\vec{\Theta}_j + \vec{m}_{\vec{\boldsymbol{x}}}.\tag{22}\] The relation between p(\(\vec{\boldsymbol{x}}\)) and \({\mathcal{N}}(\vec{\boldsymbol{x}})\) is \[p(\vec{\boldsymbol{x}}) = B \exp(-{\mathcal{N}}(\vec{\boldsymbol{x}})).\] We may avoid calculating normalizing constants by recognizing that \(\int_{{\mathbb{R}}^d} f(\vec{\boldsymbol{x}}) p(\vec{\boldsymbol{x}}) d^d \vec{\boldsymbol{x}}= 1\), which leads to: \[\label{eq:int} \begin{align} 1 &= \sqrt{|H|} ABC \sum_{j=1}^{J} \left( w_{-} e^{z_{-}}\exp(-{\mathcal{N}}(\vec{\boldsymbol{x}}_{j}^{\,-}) + w_{+} e^{z_{+}} \exp(-{\mathcal{N}}(\vec{\boldsymbol{x}}_{j}^{\,+}) \right) \end{align}\tag{23}\] It follows that \[\label{eq:intFinal} \begin{align} &\int_{{\mathbb{R}}^d} f(\vec{\boldsymbol{x}}) p(\vec{\boldsymbol{x}}) d^d \vec{\boldsymbol{x}}\approx C^{\prime} \sum_{j=1}^{J} \left(w_{-} e^{z_{-}} f(\vec{\boldsymbol{x}}_{j}^{\,-}) \exp(-{\mathcal{N}}(\vec{\boldsymbol{x}}_{j}^{\,-}) + w_{+} e^{z_{+}} f(\vec{\boldsymbol{x}}_{j}^{\,+}) \exp(-{\mathcal{N}}(\vec{\boldsymbol{x}}_{j}^{\,+}) \right). \end{align}\tag{24}\] where \[(C^{\prime})^{-1} \equiv \sum_{j=1}^{J} w_{-} e^{z_{-}}\exp(-{\mathcal{N}}(\vec{x}_{j}^{\,-})) + w_{+} e^{z_{+}} \exp(-{\mathcal{N}}(\vec{\boldsymbol{x}}_{j}^{\,+}) ).\]

Summarizing, we have: \[\label{eq:intFinal1} \int_{{\mathbb{R}}^d} f(\vec{\boldsymbol{x}}) p(\vec{\boldsymbol{x}}) d^d \vec{\boldsymbol{x}}\approx C^{\prime} \sum_{k=1}^{2J} p_k f(\vec{\boldsymbol{x}}_{k}),\tag{25}\] where \[\label{eq:defFinal} \begin{align} \vec{x}_{k}&\equiv \begin{cases} \sqrt{2z_{-}}H^{-1/2}\vec{\Theta}_k + \vec{m}_{\vec{\boldsymbol{x}}} & 1 \le k \le J ;\\ \sqrt{2z_{+}}H^{-1/2}\vec{\Theta}_{k-J} + \vec{m}_{\vec{\boldsymbol{x}}} & J+1 \le k \le 2J, \end{cases} \end{align}\tag{26}\] \[\begin{align} p_{k}&\equiv \begin{cases} w_{-} e^{z_{-}} \exp(-{\mathcal{N}}(\vec{\boldsymbol{x}}_{k})) & 1 \le k \le J; \\ w_{+} e^{z_{+}} \exp(-{\mathcal{N}}(\vec{\boldsymbol{x}}_{k})) & J+1 \le k \le 2J, \end{cases} \end{align}\] and \[C^{\prime} \equiv \left( \sum_{j=1}^{2J} p_k \right)^{-1}.\]

3.5.5 Nonlinear transformation of measurement distribution↩︎

When a target is close to a sensor \(s\), the measurement distribution is not Gaussian in Cartesian coordinates. Instead, the measurement distribution is better approximated by a Gaussian in polar coordinates, where the origin of the coordinate system is the closest sensor to the target. The coordinates transformation is given by: \[\label{eq:xy2uv} \begin{align} u& \equiv \sqrt{(x - x_{[s]})^2 + (y - y_{[s]})^2)};\\ v& \equiv u \, \text{atan2}( y - y_{[s]}, x - x_{[s]}). \end{align}\tag{27}\] A Gaussian distribution in these coordinates has the form \[p_{\vec{u}}(\vec{u}) \equiv C^{\prime \prime} \exp( -(\vec{u}-\vec{m}_{\vec{u}})^T \Sigma_{\vec{u},\vec{u}} (\vec{u}-\vec{m}_{\vec{u}}) ),\] where \[\begin{align} \vec{u}&\equiv (u,v);\\ \vec{m}_{\vec{u}} &\equiv \sqrt{x_{c,s}^2 + y_{c,s}^2} \begin{bmatrix} 1 \\ \text{atan2}( y_{c,s}, x_{c,s} ) \end{bmatrix};\\ x_{c,s} &\equiv x_c - x_{[s]}; \qquad y_{c,s} \equiv y_c - y_{[s]}. \end{align}\] It follows that expected values may be computed (in analogy to 25 ) as \(\int_{{\mathbb{R}}^d} f(\,\vec{x}(\vec{u}) \, ) p_{\vec{u}}(\vec{u}) d^d \vec{u}\) By changing variable we may write this as an integral in \(\vec{x}\): \[\int_{{\mathbb{R}}^d} f(\vec{x}) p_{\vec{u}}(\,\vec{u}(\vec{x})\, ) \left| \frac{\partial \vec{u}}{\partial \vec{x}} \right| d^d \vec{x},\] where \[\begin{align} \frac{\partial \vec{u}}{\partial \vec{x}} &= \frac{1}{r}\begin{bmatrix} x & y \\ x \, \text{atan2}(y,x) - y & y \, \text{atan2}(y,x) + x \end{bmatrix}; \\ \left| \frac{\partial \vec{u}}{\partial \vec{x}} \right| &= 1. \end{align}\]

\[\label{eq:defFinal2} \begin{align} \vec{u}_{k}&\equiv \begin{cases} \sqrt{2z_{-}}\Sigma_{\vec{u}\vec{u}}^{1/2}\vec{\Theta}_k + \vec{m}_{\vec{u}} & 1 \le k \le J \\ \sqrt{2z_{+}}\Sigma_{\vec{u}\vec{u}}^{1/2}\vec{\Theta}_{k-J} + \vec{m}_{\vec{u}} & J+1 \le k \le 2J \end{cases} \end{align}\tag{28}\]

The covariance in \(\vec{u}\) coordinates \(\Sigma_{\vec{u}\vec{u}}\) may be computed using the local transformation of \(\Sigma_{\vec{x}\vec{x}}\): \[\Sigma_{\vec{u}\vec{u}} = \left. \frac{\partial \vec{u}}{\partial \vec{x}} \right|_{\vec{x}= \vec{x}_{c,s}} \cdot \Sigma_{\vec{x}\vec{x}} \cdot \left. \frac{\partial \vec{u}}{\partial \vec{x}} \right|_{\vec{x}= \vec{x}_{c,s}} ^{-1},\]

In order to transform the points \(\{\vec{u}_{k}\}\) back to \((x,y)\) coordinates, we use the inverse transformation to 27 : \[\label{eq:uv2xy} x_k = u_k \cos (v_k/u_k) + x_{[s]}; \qquad y_k = u_k \sin(v_k/u_k) + y_{[s]};\tag{29}\]

3.6 Updating mean and covariance↩︎

The updated mean and spatial covariance matrix at time \(t\) are given by: \[\begin{align} \vec{m}_{\vec{x},t} &= \sum_{k=0}^{2J} \vec{x}_{k} p_k \\ \Sigma_{\vec{x}\vec{x}, new} &= \sum_{k=0}^{2J} \vec{x}_{k} \vec{x}_{k}^T p_k - \vec{m}_{\vec{x},t} \vec{m}_{\vec{x},t}^T. \end{align}\]

We must also compute the updated mean and covariance for the velocities for time \(t\). We may use the following facts, which hold for Gaussian distributions[8]: \[\begin{align} \vec{m}_{\vec{v}| \vec{x}} &= \vec{m}_{\vec{v}}^{(prop)} + \Sigma_{\vec{v}\vec{x}}^{(prop)} (\Sigma_{\vec{x}\vec{x}}^{(prop)})^{-1}(\vec{x}- \vec{m}_{\vec{x}}^{(prop)}); \\ \Sigma_{\vec{v}\vec{v}| \vec{x}} &= \Sigma_{\vec{v}\vec{v}}^{(prop)} - \Sigma_{\vec{v}\vec{x}}^{(prop)} (\Sigma_{\vec{x}\vec{x}}^{(prop)})^{-1} \Sigma_{\vec{x}\vec{v}}^{(prop)}. \end{align}\] so we may write: \[\begin{align} \vec{m}_{\vec{v},t} &= \sum_{k=1}^{2J} \vec{m}_{\vec{v}| \vec{x}_k} p_k = \vec{q} + Q \vec{m}_{\vec{x}}, \end{align}\] where \[\begin{align} \vec{q} &\equiv \vec{m}_{\vec{v}}^{(prop)} - \Sigma_{\vec{v}\vec{x}}^{(prop)} (\Sigma_{\vec{x}\vec{x}}^{(prop)})^{-1} \vec{m}_{\vec{x}}^{(prop)} \text{ and }Q \equiv \Sigma_{\vec{v}\vec{x}}^{(prop)} (\Sigma_{\vec{x}\vec{x}}^{(prop)})^{-1}. \end{align}\]

We also have \[\begin{align} &E[\vec{v}\vec{v}^T \, |\, \vec{x}] = \Sigma_{\vec{v}\vec{v}| \vec{x}} + \vec{m}_{\vec{v}| \vec{x}} \vec{m}_{\vec{v}| \vec{x}}^T \implies E[\vec{v}\vec{v}^T] = \Sigma_{\vec{v}\vec{v}| \vec{x}} + \sum_{k=1}^{2J} \vec{m}_{\vec{v}| \vec{x}_k} \vec{m}_{\vec{v}| \vec{x}_k}^T p_k, \end{align}\] so that \[\begin{align} E[\vec{v}\vec{v}^T] = &\Sigma_{\vec{v}\vec{v}| \vec{x}} + \sum_{k=1}^{2J} (\vec{q} + Q \vec{x}_k) (\vec{q} + Q \vec{x}_k)^T p_k\\ = &\Sigma_{\vec{v}\vec{v}| \vec{x}} + \vec{q}\vec{q}^T + \vec{q} (Q \vec{m}_{\vec{x}})^T + (Q \vec{m}_{\vec{x}}) \vec{q}^T + \sum_{k=1}^{2J} Q \vec{x}_k (Q \vec{x}_k)^T p_k. \end{align}\] This gives the following expression for the velocity covariance: \[\begin{align} &\Sigma_{\vec{v}\vec{v},t } = E[\vec{v}\vec{v}^T] - E[\vec{v}]E[\vec{v}^T] \\ & \quad = \Sigma_{\vec{v}\vec{v}| \vec{x}} + \sum_{k=1}^{2J} Q \vec{x}_k (Q \vec{x}_k)^T p_k - Q \vec{m}_{\vec{x}} (Q \vec{m}_{\vec{x}})^T . \end{align}\] The cross-covariance between velocity and position is also needed: \[\begin{align} \Sigma_{\vec{v}\vec{x},t} &=E[\vec{v}\vec{x}^T] - E[\vec{v}]E[\vec{x}]^T \\ &= \sum_{k=1}^{2J} \vec{m}_{\vec{v}| \vec{x}_k} \vec{x}_k^T p_k - \vec{m}_{\vec{v},t} \vec{m}_{\vec{x},t}^T \\ &= \sum_{k}^{2J} p_k (Q \vec{x}_k) (\vec{x}_k - \vec{m}_{\vec{x},t})^T \end{align}\] The complete covariance matrix estimate \(\Sigma_{t}\) is composed of the blocks \(\Sigma_{\vec{x}\vec{x},t}\), \(\Sigma_{\vec{v}\vec{v},t}\), \(\Sigma_{\vec{v}\vec{x},t}\), \(\Sigma_{\vec{v}\vec{x},t}^T\) as follows: \[\Sigma_{t} = \begin{bmatrix} \Sigma_{\vec{x}\vec{x},t} & \Sigma_{\vec{v}\vec{x},t}^T\\ \Sigma_{\vec{v}\vec{x},t} & \Sigma_{\vec{v}\vec{v},t} \end{bmatrix}\] The vector \((\vec{m}_{\vec{x},t}; \vec{m}_{\vec{v},t})\) and matrix \(\Sigma_{t}\) are then used in 8 and@eq:eq:propCov to obtain the propagated Gaussian prior for time \(t+1\), which (when combined with the spatial distribution derived from the measurement) produces position and velocity estimates for time \(t+1\). Tracking of the targets proceeds by iterating this process.

4 Simulations↩︎

4.1 Simulation setup↩︎

The simulation was set up with 4 moving targets in an area of 40 meters by 40 meters with 25 omnidirectional sensors spaced 10 meters apart as shown in 1.

The mathematical model used was described in previous sections. Besides the default parameter values, the following parameter values were used:

For each different set of parameters, 50 random trajectories were generated according to the specifications described in [3]. The initial target vector was chosen from a Gaussian with mean equal to the true target vector, and covariance as a diagonal matrix with spatial entries equal to 100 and velocity entries equal to 0.0005 (as in [3]). For each trajectory, there was one set of measurements on which each filter was ran once. Figure 3 shows position estimates and particles obtained by our filter (designated as the ‘TT filter’ in the following discussion) during a typical simulation.

To compare average errors, we also used the optimal mass transfer (OMAT) metric[9] which matches estimated target positions with the actual target positions so as to minimize total distance error (Li and Coates also used this metric in [3], and it is included in their Matlab testbench). The five different filters tried in the simulations were the PFPF_EDH and PFPF_LEDH filters designed by Li and Coates; , a standard bootstrap particle filter with 1 million points; and the TT filter without and with the nonlinear correction described in Section 3.5.5. The three previous filters used were the best-performing filters investigated by Li and Coates in [3].

For each parameter set, to obtain final OMAT values we averaged over 50 different multitarget tracks over 40 time steps. Simulations were run with MATLAB R2020a on a 64-bit computer with Intel Core i5 dual cores at 2.60 GHz with 6.0 GB total virtual memory.

4.2 Filter performance comparisons↩︎

Figure 2 shows the OMAT metric (averaged over 50 multitarget tracks) at each time step for the five tracking algorithms, for \(\sigma_w^2 = 0.1\) and other parameters at baseline.

Table ¿tbl:table:baseline? compares the performance for the different algorithms for \(\sigma_w^2 = 0.1\) and other parameters at baseline. The TT filter gives 36% lower error than the PFPF_LEDH, with less than 1/10 of the execution time.

Figures 3-5 compare the accuracy of the different algorithms for ranges of values for \(\sigma_w^2\), measurement variance, filter’s estimate of spatial covariance, and process noise covariance respectively. The TT algorithms consistently gave the smallest error across all parameter values. Figures 4 and 5 which shows that only the TT algorithm shows consistent improvements when the filter’s spatial covariance estimate or the process noise spatial covariance is decreased.

4.3 Variants of the TT filter↩︎

Figures 6-8 compare the performance of several variants of the TT algorithm. The following variants are included:

• TT baseline: TT algorithm with the nonlinear correction described in Section 3.5.5, the hopping target correction described in Section 3.3.2, and the 1 by 1 correction method described in Section 3.3.1. This option also uses the randomized initialization used by Li and Coates in [3].

• TT linear: same as previous, but without the nonlinear correction described in Section 3.5.5;

• TT noHopping: same as TT nonLinear, except without the hopping target correction described in Section 3.3.2;

• TT fixedInit: TT algorithm without any prior estimate of the target position vector at \(t=0\). Targets were initialized to lie within 5 m of the central sensor, and the ML702LE estimation procedure described in Section 3.1.1 was used with the first measurement to obtain the target vector at the first time step \(t=1\).

• TT hopping no 1 by1: Same as TT nonLinear above, except without the 1 by 1 correction method described in Section 3.3.1;

• TT no hopping no 1 by1: Same as TT nonLinear above, except without neither the 1 by 1 correction method described in Section 3.3.1 nor the hopping target correction described in Section 3.3.2;

The figures show that all variants perform nearly the same except when none of the corrections described in Section 3.3 are applied. We conclude that either correction method performs equally well, and there is no need for secondary correction. We note also that the filters of Li and Coates require an initial guess target vector that is not too far from the true target vector to achieve convergence, while the TT filter converges as long as targets are placed near the center of the sensing grid. This represents another distinct advantage of the TT filter.

5 Conclusion↩︎

We have documented the design,validation, and verification of a higher-performing, lower complexity multi-target tracking algorithm for a sensor network. Our TT filter achieves over 30% reduction in error over the most accurate filter described in [3] (as measured by the OMAT metric) while reducing the complexity by over 90%. The TT algorithm also works equally well in scenarios where there is no prior knowledge of targets before tracking begins.

For future work, we may examine larger systems with more targets and sensors. We will also improve the uncertainty estimation for targets. Currently, only the covariance of the targets is computed at time \(t\). Since the target position distribution is not Gaussian, this is not an accurate indication of the target’s uncertainty. By making use of nonlinear transformations, we hope to improve the uncertainty region.

Ackowledgements↩︎

This research was supported by the Air Force Research Laboratory, and has been approved for public release: distribution unlimited (Case Number: 88ABW-2020-2798).

References↩︎