Optimal Geometries for AOA Localization in the Bayesian Sense

Abstract

This paper considers the optimal sensor placement problem for angle-of-arrival (AOA) target localization in the 2D plane with a Gaussian prior. Optimal sensor locations are analytically determined for a single AOA sensor using the D- and A-optimality criteria and an approximation of the Bayesian Fisher information matrix (BFIM). Optimal sensor placement is shown to align with the minor axis of the prior covariance error ellipse for both optimality criteria. The approximate BFIM is argued to be valid for a sufficiently small prior covariance compared with the target range. Optimal sensor placement results obtained for Bayesian target localization are extended to manoeuvring target tracking. For sensor trajectory optimization subject to turn-rate constraints, numerical search methods based on the D- and A-optimality criteria as well as a new closed-form projection algorithm that aims to achieve alignment with the minor axis of the prior error ellipse are proposed. It is observed that the two optimality criteria generate significantly different optimal sensor trajectories despite having the same optimal sensor placement for the localization of a stationary target. Analysis results and the performance of the sensor trajectory optimization methods are demonstrated with simulation examples. It is observed that the new closed-form projection algorithm achieves superior tracking performance compared with the two numerical search methods.

1. Introduction

In target tracking and localization problems, target–sensor geometries are known to play a significant role in determining the localization and tracking performance. In this paper, we focus on optimal target–sensor geometries for angle-of-arrival (AOA) localization in the 2D plane using a single moving sensor. First, the localization problem is cast as a Bayesian estimation problem, which assumes the availability of prior information in the form of a Gaussian prior for the unknown target location. For this problem, optimal sensor placement results are developed using approximate estimation bounds. Next, the Bayesian estimation problem is extended to target tracking using the Kalman filter, and optimal sensor trajectories are developed to track a manoeuvring target.

Optimal sensor placement has been researched for several decades. Early works included [1,2,3], where the performance of the extended Kalman filter (EKF) [4] and several deterministic (non-Bayesian) estimators was reported for different bearings-only sensor manoeuvres, mostly in sonar applications. In [5], optimal bearings-only sensor manoeuvres for tracking a constant-velocity target were derived using optimal control theory. The sensor trajectory optimization problem was formulated as a partially observable Markov decision problem (POMDP) for a manoeuvring target using the trace of FIM, which is similar to the A-optimality criterion (minimizing trace of inverse FIM), as the reward function in [6]. In [7], the D-optimality criterion, whereby the sensor location is determined to maximize the determinant of the Fisher information matrix (FIM), was adopted to determine optimal sensor trajectories for the localization of a stationary target. Optimal sensor manoeuvres necessary to make a constant-velocity target observable were discussed in [8]. The sensor trajectory optimization problem in the Bayesian sense was considered in [9], where the posterior Cramer–Rao lower bound (PCRLB) [10] was employed to minimize the largest root-mean-square-error (RMSE)-bound approximated by the reciprocal of measurement data contributions to target location information.

A comprehensive analysis of the 2D optimal AOA sensor placement problem for a stationary target was presented in [11,12]. Optimal 3D AOA target–sensor geometries for a stationary target were analysed in [13]. A gradient-descent algorithm for sensor path optimization to minimize the mean-square error of predicted EKF target location estimates was proposed in [14]. A UAV path optimization algorithm that solves a nonlinear programming problem based on the D-optimality criterion was developed in [15] to geolocate a stationary target using a heterogeneous mix of passive payload sensors. In [16], a unified framework was proposed for AOA, range-only, and received signal strength localization when the target was stationary. Optimal target–sensor geometries for maximum a posteriori (MAP) target localization with a Gaussian prior were investigated in [17]. In [18], the optimality criteria for target–sensor geometries in a Kalman filtering setting were analysed for several sensor types using an approximation of the Bayesian FIM. A unified 2D target–sensor geometry optimization framework was proposed in [19] for stationary target localization with a Gaussian prior, reducing the optimization problem to minimization of the modulus of a vector sum, akin to [11]. In [20], optimal sensor placement in 3D space was studied for AOA target localization with a Gaussian prior, employing rotational invariance arguments.

This paper develops optimal sensor placement results for a single AOA sensor at a fixed distance from the mean of the Gaussian prior. To do this, the Bayesian FIM is approximated by replacing the expectation of the contribution form measurement data with its instantaneous value calculated at the mean of the Gaussian prior. It is argued that this approximation is valid when the covariance of the Gaussian prior is relatively small compared with the target range. The optimal sensor placements for the D- and A-optimality criteria are shown to be identical and align with the minor axis of the error ellipse of the prior covariance. In the context of bearings-only manoeuvring target tracking, numerical search methods based on the D- and A-optimality criteria and a new closed-form projection algorithm that attempts to achieve alignment with the minor axis of the prior covariance error ellipse are proposed for sensor trajectory optimization subject to turn-rate constraints. It is observed that the D- and A-optimality criteria yield markedly different optimal sensor trajectories even though they produce identical optimal sensor placement for the localization of a stationary target. The projection algorithm is shown to outperform the other two methods in simulation studies.

This paper is organized as follows. Section 2 investigates the optimal sensor placement problem for a stationary target with a Gaussian prior using the D- and A-optimality criteria. Section 3 extends the results of Section 2 to Kalman filter tracking of a manoeuvring target, proposing two sensor trajectory optimization methods and a new closed-from projection algorithm. Section 4 presents simulation examples to verify the optimal sensor placement results derived in Section 2 and to compare and demonstrate the effectiveness of the sensor trajectory optimization algorithms proposed in Section 3. Concluding remarks are made in Section 5.

2. Optimal Target-Sensor Geometry with Gaussian Prior

In tracking problems with a state space that can be modelled or approximated as a Gauss–Markov process, the Kalman filter has been extensively used to compute the Gaussian prior for state estimates from the noisy sensor measurements available in each recursion in the form of the predicted state estimate and predicted state covariance. Starting with the Gaussian prior $\mathcal{N}({\mathit{x}}_0,{\mathit{P}}_0)$, where ${\mathit{x}}_0$ is the mean and ${\mathit{P}}_0$ is the covariance of the prior, the objective of optimal sensor placement is to determine a sensor location $\mathit{s}$ at a fixed distance of $d=\parallel \mathit{d}\parallel$, where $\mathit{d}={\mathit{x}}_0-\mathit{s}$, from the mean of the Gaussian prior (or predicted target location estimate) ${\mathit{x}}_0$ so that sensor measurements collected at the new sensor location will optimize a well-defined objective function that is related to the Bayesian localization performance. The new measurements together with the Gaussian prior are then used to compute the filtered state estimate and covariance, which are optimized in terms of target–sensor geometry. We consider two optimality criteria for sensor placement: namely, the D-optimality and A-optimality criteria [21,22,23], which are commonly used in practice. Referring to Figure 1, the task of geometry optimization is reduced to finding a range vector $\mathit{d}$ or bearing angle $\theta \left({\mathit{x}}_0\right)$ pivoted at the mean of the Gaussian prior with $d=\parallel \mathit{d}\parallel$ fixed, which gives the location of the sensor to satisfy the chosen optimization criterion.

The bearing measurements collected by the sensor located at $\mathit{s}$ are given by

$$\theta =\theta \left(\mathit{x}\right)+w$$

(1)

where $\mathit{x}\sim \mathcal{N}({\mathit{x}}_{\mathit{0}},{\mathit{P}}_{\mathit{0}})$ is the Bayesian prior, $w\sim \mathcal{N}(0,{\sigma }^2)$ is the bearing angle noise, and

$$\theta \left(\mathit{x}\right)={tan}^{-1}(x_2-s_2,x_1-s_1),-\pi \le h\left(\mathit{x}\right)\le \pi$$

(2)

is the true bearing angle with $\mathit{x}={[x_1,x_2]}^T$ and $\mathit{s}={[s_1,s_2]}^T$. In (2), ${tan}^{-1}(·)$ is the four-quadrant arc-tangent. The objective of Bayesian estimation is to determine an estimate for the unknown random target location $\mathit{x}$ from the bearing measurement $\theta$ and the knowledge of the Gaussian prior with mean ${\mathit{x}}_0$ and covariance ${\mathit{P}}_0$.

The optimality criteria considered in this paper employ estimation bounds obtained from the FIM or CRLB. In Bayesian estimation problems that involve random unknown parameters, these bounds are replaced by the Bayesian FIM (BFIM) or Bayesian CRLB (BCRLB). The BFIM for the single-sensor AOA localization problem is defined as [24]

$$\mathsf{\Phi }={\mathit{K}}_0+E_{\mathit{x}}\left\{\frac{1}{{\sigma }^2{d}^2\left(\mathit{x}\right)}\mathit{u}\left(\mathit{x}\right){\mathit{u}}^T\left(\mathit{x}\right)\right\}$$

(3)

where $E_{\mathit{x}}\{·\}$ denotes the expectation over $\mathit{x}$,

$$\mathit{u}\left(\mathit{x}\right)=\left[\begin{array}{c}-sin\theta \left(\mathit{x}\right)\\ cos\theta \left(\mathit{x}\right)\end{array}\right]$$

(4)

is the unit vector orthogonal to the range vector, and $d\left(\mathit{x}\right)=\parallel \mathit{x}-\mathit{s}\parallel$ is the target range from the sensor positioned at $\mathit{s}$. The matrix ${\mathit{K}}_0={\mathit{P}}_0^{-1}$ represents the contribution of the a priori information, and

$$E_{\mathit{x}}\left\{\frac{1}{{\sigma }^2{d}^2\left(\mathit{x}\right)}\mathit{u}\left(\mathit{x}\right){\mathit{u}}^T\left(\mathit{x}\right)\right\}$$

(5)

is the contribution of data. The inverse of the BFIM gives the BCRLB.

Equation (5) can be rewritten as

$$\begin{array}{cc}\hfill E_{\mathit{x}}\left\{\frac{1}{{\sigma }^2{d}^2\left(\mathit{x}\right)}\mathit{u}\left(\mathit{x}\right){\mathit{u}}^T\left(\mathit{x}\right)\right\}& =\frac{1}{{\sigma }^2}{E}_{\mathit{x}}\left\{\frac{1}{d^2\left(\mathit{x}\right)}\mathit{u}\left(\mathit{x}\right){\mathit{u}}^T\left(\mathit{x}\right)\right\}\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill & =\frac{1}{{\sigma }^2}\left[\begin{array}{cc}E\left\{\frac{{sin}^2\theta \left(\mathit{x}\right)}{{\parallel \mathit{x}-\mathit{s}\parallel }^2}\right\}& -E\left\{\frac{sin\theta \left(\mathit{x}\right)cos\theta \left(\mathit{x}\right)}{{\parallel \mathit{x}-\mathit{s}\parallel }^2}\right\}\\ -E\left\{\frac{sin\theta \left(\mathit{x}\right)cos\theta \left(\mathit{x}\right)}{{\parallel \mathit{x}-\mathit{s}\parallel }^2}\right\}& E\left\{\frac{{cos}^2\theta \left(\mathit{x}\right)}{{\parallel \mathit{x}-\mathit{s}\parallel }^2}\right\}\end{array}\right]\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill & \approx \frac{1}{{\sigma }^2{d}^2}\mathit{u}\left({\mathit{x}}_0\right){\mathit{u}}^T\left({\mathit{x}}_0\right),\mathit{u}\left({\mathit{x}}_0\right)=\left[\begin{array}{c}-sin\theta \left({\mathit{x}}_0\right)\\ cos\theta \left({\mathit{x}}_0\right)\end{array}\right]\hfill \end{array}$$

(8)

which gives

$$\mathsf{\Phi }\approx {\mathit{K}}_0+\frac{1}{{\sigma }^2{d}^2}\mathit{u}\left({\mathit{x}}_0\right){\mathit{u}}^T\left({\mathit{x}}_0\right).$$

(9)

The approximation in (9) is valid for a sufficiently “small” prior covariance compared with the target range. Here, we measure the size of prior covariance by its trace (see (16)). Referring to Figure 2, we have

$$\begin{array}{cc}\hfill sin\theta \left(\mathit{x}\right)& =sin\theta ({\mathit{x}}_0+\eta )\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill & =\frac{x_{0,2}-s_2+{\eta }_2}{\sqrt{{(x_{0,1}-s_1+{\eta }_1)}^2+{(x_{0,2}-s_2+{\eta }_2)}^2}}\hfill \end{array}$$

(11)

where $\eta ={[{\eta }_1,{\eta }_2]}^T\sim \mathcal{N}(\mathit{0},{\mathit{P}}_0)$, and ${\mathit{x}}_0={[x_{0,1},x_{0,2}]}^T$. Taking the expectation of the squared sine function yields

$$\begin{array}{cc}\hfill E\left\{{sin}^2\theta \left(\mathit{x}\right)\right\}& \approx \frac{E\left\{{(x_{0,2}-s_2+{\eta }_2)}^2\right\}}{E\{{(x_{0,1}-s_1+{\eta }_1)}^2+{(x_{0,2}-s_2+{\eta }_2)}^2\}}\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill & \approx \frac{{(x_{0,2}-s_2)}^2+E\left\{{\eta }_2^2\right\}}{d^2+E\left\{{\eta }_1^2\right\}+E\left\{{\eta }_2^2\right\}}\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill & \approx {\left(\frac{x_{0,2}-s_2}{d}\right)}^2\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill & \approx {sin}^2\theta \left({\mathit{x}}_0\right)\hfill \end{array}$$

(15)

if

$$d^2\gg E\left\{{\eta }_1^2\right\}+E\left\{{\eta }_2^2\right\}=tr{\mathit{P}}_0.$$

(16)

Here, tr denotes trace. Applying the same line of reasoning to $E\{sin\theta (\mathit{x})cos\theta (\mathit{x}\left)\right\}$ and $E\left\{{cos}^2\theta \left(\mathit{x}\right)\right\}$, we conclude that (16) is the condition that must be met to justify (9).

2.1. D-Optimality Criterion

The D-optimality criterion aims to maximize the determinant of the Fisher information matrix (FIM). For the Bayesian estimation problem considered here, the FIM is replaced by the Bayesian FIM (BFIM). Considering the optimization problem described in Figure 1, the optimal placement for the sensor is obtained from the solution of

$$\underset{\mathit{s}}{max}|\mathsf{\Phi }|$$

(17)

where $|·|$ denotes the determinant. The optimization problem in (17) determines the sensor location $\mathit{s}$ that maximizes the determinant of BFIM for a given Gaussian prior.

Noting that in (9) ${\mathit{K}}_0$ is a square matrix and $\mathit{u}\left({\mathit{x}}_0\right)$ is a column vector, the determinant of the BFIM can be rewritten as a sum of two terms [25]

$$\begin{array}{cc}\hfill |\mathsf{\Phi }|& \approx \left|{\mathit{K}}_0+\frac{1}{{\sigma }^2{d}^2}\mathit{u}\left({\mathit{x}}_0\right){\mathit{u}}^T\left({\mathit{x}}_0\right)\right|\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill & \approx |{\mathit{K}}_0|+\frac{1}{{\sigma }^2{d}^2}{\mathit{u}}^T\left({\mathit{x}}_0\right){\mathit{K}}_0^{*}\mathit{u}\left({\mathit{x}}_0\right)\hfill \end{array}$$

(19)

where ${\mathit{K}}_0^{*}$ is the adjoint of ${\mathit{K}}_0$, defined by

$${\mathit{K}}_0^{*}=|{\mathit{K}}_0|{\mathit{K}}_0^{-1}=\left|{\mathit{K}}_0\right|{\mathit{P}}_0.$$

(20)

Thus, for a given Gaussian prior ${\mathit{P}}_0$ and fixed d, the optimization problem in (17) reduces to

$$\underset{\theta \left({\mathit{x}}_0\right)}{max}{\mathit{u}}^T\left({\mathit{x}}_0\right){\mathit{P}}_0\mathit{u}\left({\mathit{x}}_0\right).$$

(21)

Since $\mathit{u}\left({\mathit{x}}_0\right)$ is a unit vector, (21) is a problem of quadratic form maximization over the unit circle. Using the eigenvalues of ${\mathit{P}}_0$, denoted by ${\lambda }_1,{\lambda }_2$ with ${\lambda }_1\ge {\lambda }_2$, the solution of (21) is given by [26]:

$$\begin{array}{c}{\lambda }_1=\underset{\theta \left({\mathit{x}}_0\right)}{max}{\mathit{u}}^T\left({\mathit{x}}_0\right){\mathit{P}}_0\mathit{u}\left({\mathit{x}}_0\right)\end{array}$$

(22)

$$\begin{array}{c}{\nu }_1=\underset{\mathit{u}\left({\mathit{x}}_0\right)}{arg\; max}{\mathit{u}}^T\left({\mathit{x}}_0\right){\mathit{P}}_0\mathit{u}\left({\mathit{x}}_0\right)\end{array}$$

(23)

where ${\nu }_1$ is an orthonormal eigenvector of ${\mathit{P}}_0$ corresponding to ${\lambda }_1$. The optimal bearing angle for the sensor, ${\theta }_{\mathrm{opt}}$, is easily obtained from the optimal unit vector ${\nu }_1$ by noting that it is orthogonal to the range vector (see Figure 1). In other words, the optimal range vector ${\mathit{d}}_{\mathrm{opt}}$ must be aligned with the minor axis of the error ellipsoid of the Gaussian prior, as shown in Figure 3.

Some remarks are in order here:

• If the Gaussian prior has a circular error ellipse with ${\mathit{P}}_0$ given by a scaled identity matrix, (21) becomes

  $$\underset{\theta \left({\mathit{x}}_0\right)}{max}{\mathit{u}}^T\left({\mathit{x}}_0\right)\mathit{u}\left({\mathit{x}}_0\right)\mathrm{or}\underset{\theta \left({\mathit{x}}_0\right)}{max}1$$

  (24)

  which means that optimality is achieved by any bearing angle $\theta \left({\mathit{x}}_0\right)$.
• In all other cases, there are two optimal bearing angles aligned with the minor axis of the error ellipse, producing two possible optimal sensor locations with range vectors $\pm {\mathit{d}}_{\mathrm{opt}}$ symmetric about ${\mathit{x}}_0$.

2.2. A-Optimality Criterion

The objective of the A-optimality criterion is to minimize the trace of the BCRLB or the inverse BFIM. In this case, the optimal sensor placement is obtained from

$$\underset{\mathit{s}}{min}tr\left({\mathsf{\Phi }}^{-1}\right).$$

(25)

Applying the matrix inversion lemma [27] to the approximate BFIM in (9), we get

$${\mathsf{\Phi }}^{-1}\approx {\mathit{P}}_0\left(\mathit{I}-\frac{\mathit{u}\left({\mathit{x}}_0\right){\mathit{u}}^T\left({\mathit{x}}_0\right){\mathit{P}}_0}{{\sigma }^2{d}^2+{\mathit{u}}^T\left({\mathit{x}}_0\right){\mathit{P}}_0\mathit{u}\left({\mathit{x}}_0\right)}\right)$$

(26)

It is clear that to solve (25), we need to maximize the trace of the second term on the right-hand side of (26); i.e.,

$$\begin{array}{c}\underset{\mathit{s}}{max}tr\left(\frac{{\mathit{P}}_0\mathit{u}\left({\mathit{x}}_0\right){\mathit{u}}^T\left({\mathit{x}}_0\right){\mathit{P}}_0}{{\sigma }^2{d}^2+{\mathit{u}}^T\left({\mathit{x}}_0\right){\mathit{P}}_0\mathit{u}\left({\mathit{x}}_0\right)}\right)\end{array}$$

(27)

or

$$\begin{array}{c}\hfill \underset{\mathit{u}\left({\mathit{x}}_0\right)}{max}\frac{{\mathit{u}}^T\left({\mathit{x}}_0\right){\mathit{P}}_0^2\mathit{u}\left({\mathit{x}}_0\right)}{{\sigma }^2{d}^2+{\mathit{u}}^T\left({\mathit{x}}_0\right){\mathit{P}}_0\mathit{u}\left({\mathit{x}}_0\right)}.\end{array}$$

(28)

To solve (28) for $\mathit{u}\left({\mathit{x}}_0\right)$, let $R=\frac{{\mathit{u}}^T\left({\mathit{x}}_0\right){\mathit{P}}_0^2\mathit{u}\left({\mathit{x}}_0\right)}{{\sigma }^2{d}^2+{\mathit{u}}^T\left({\mathit{x}}_0\right){\mathit{P}}_0\mathit{u}\left({\mathit{x}}_0\right)}$ and set its gradient equal to zero,

$$\frac{\partial }{\partial \mathit{u}\left({\mathit{x}}_0\right)}R=\mathit{0}$$

(29)

which results in

$$\begin{array}{c}\frac{2{\mathit{P}}_0^2\mathit{u}\left({\mathit{x}}_0\right)({\sigma }^2{d}^2+{\mathit{u}}^T\left({\mathit{x}}_0\right){\mathit{P}}_0\mathit{u}\left({\mathit{x}}_0\right))-2\left({\mathit{u}}^T\left({\mathit{x}}_0\right){\mathit{P}}_0^2\mathit{u}\left({\mathit{x}}_0\right)\right){\mathit{P}}_0\mathit{u}\left({\mathit{x}}_0\right)}{{({\sigma }^2{d}^2+{\mathit{u}}^T\left({\mathit{x}}_0\right){\mathit{P}}_0\mathit{u}\left({\mathit{x}}_0\right))}^2}=\mathit{0}\end{array}$$

(30)

$$\begin{array}{c}{\mathit{P}}_0^2\mathit{u}\left({\mathit{x}}_0\right)({\sigma }^2{d}^2+{\mathit{u}}^T\left({\mathit{x}}_0\right){\mathit{P}}_0\mathit{u}\left({\mathit{x}}_0\right))=\left({\mathit{u}}^T\left({\mathit{x}}_0\right){\mathit{P}}_0^2\mathit{u}\left({\mathit{x}}_0\right)\right){\mathit{P}}_0\mathit{u}\left({\mathit{x}}_0\right)\end{array}$$

(31)

$$\begin{array}{c}{\mathit{P}}_0^2\mathit{u}\left({\mathit{x}}_0\right)=R{\mathit{P}}_0\mathit{u}\left({\mathit{x}}_0\right).\end{array}$$

(32)

Left-multiplying both sides of the above equation with ${\mathit{P}}_0^{-1}$ finally gives

$${\mathit{P}}_0\mathit{u}\left({\mathit{x}}_0\right)=R\mathit{u}\left({\mathit{x}}_0\right)$$

(33)

where the nonzero scalar R is an eigenvalue of ${\mathit{P}}_0$ and the unit vector $\mathit{u}\left({\mathit{x}}_0\right)$ is the corresponding eigenvector. We therefore conclude that R is maximized when $\mathit{u}\left({\mathit{x}}_0\right)={\nu }_1$, which is the eigenvector of ${\mathit{P}}_0$ associated with its largest eigenvalue ${\lambda }_1$. Note that this optimality result is identical to that for the A-optimality criterion derived in Section 2.1.

3. Application to Tracking

In this section, sensor waypoint optimization algorithms are devised to embed the optimal sensor placement results derived in Section 2 into the Kalman filter. As a specific application, bearings-only manoeuvring target tracking is considered. When the target is moving, it is often the case that the target dynamics and constraints on the motion of a single sensor, such as turn-rates and distances between successive waypoints, do not allow strictly optimal sensor placement geometries to be achieved from one Kalman filter recursion to the next. We develop sensor trajectory optimization methods that respect dynamic sensor constraints.

The principle we follow is based on the treatment of each Kalman filter recursion as solving a Bayesian target localization problem with a Gaussian prior available from the previous recursion and measurements taken at a new optimized sensor location to compute filtered state estimates and an updated prior for the next Kalman filter recursion. Figure 4 captures the computational steps of a Kalman filter recursion with sensor waypoint optimization embedded into it. The details of how optimal sensor waypoints are computed from Kalman filter parameters are discussed later in the section.

The single moving sensor collects AOA measurements from a manoeuvring target at time instants $k=0,1,2,\dots$. The sensor location at time k is denoted by ${\mathit{s}}_k$. The process equation for the target is

$${\mathit{x}}_{k+1}=\mathit{F}{\mathit{x}}_k+{\mathit{n}}_k,k=0,1,\dots$$

(34)

where ${\mathit{x}}_k={[x_k,{\dot{x}}_k,y_k,{\dot{y}}_k]}^T$ is the target state vector with ${[x_k,y_k]}^T$ and ${[{\dot{x}}_k,{\dot{y}}_k]}^T$ denoting the target location and velocity, respectively, at time k. In (34) the dynamical constraint (the state transition matrix) is given by

$$\mathit{F}=\left[\begin{array}{cc}\mathit{A}& \mathit{0}\\ \mathit{0}& \mathit{A}\end{array}\right],\mathit{A}=\left[\begin{array}{cc}1& T\\ 0& 1\end{array}\right]$$

(35)

where T denotes the time interval between discrete-time instants k. The process noise ${\mathit{n}}_k$ accounts for unknown target manoeuvres and is zero-mean white Gaussian with covariance

$$\mathit{Q}=\left[\begin{array}{cc}{q}_x\mathit{B}& \mathit{0}\\ \mathit{0}& q_y\mathit{B}\end{array}\right],\mathit{B}=\left[\begin{array}{cc}{T}^4/4& T^3/2\\ T^3/2& T^2\end{array}\right]$$

(36)

where $q_x$ and $q_y$ are often determined from maximum target acceleration [28].

The AOA measurement equation is

$$z_k=h\left({\mathit{x}}_k\right)+w_k$$

(37)

where $h(·)$ is the bearing angle of the target from the sensor location [see (2)], and $w_k\sim \mathcal{N}(\mathit{0},{\sigma }^2)$ is the bearing measurement noise. As the measurement Equation (37) is nonlinear, the extended Kalman filter (EKF) is often used to estimate the target state vector, which is given by the recursion:

State Prediction:

$$\begin{array}{cccc}\hfill & {\mathit{x}}_{k|k-1}\hfill & =\mathit{F}{\mathit{x}}_{k-1|k-1}\hfill & \mathrm{mean}\mathrm{of}\mathrm{prior}\hfill \end{array}$$

(38)

$$\begin{array}{ccc}{\mathit{P}}_{k|k-1}\hfill & =\mathit{F}{\mathit{P}}_{k-1|k-1}{\mathit{F}}^T+\mathit{Q}\hfill & \mathrm{covariance}\mathrm{of}\mathrm{prior}\hfill \end{array}$$

(39)

State Update:

$$\begin{array}{ccc}\hfill & \hfill {\tilde{z}}_k& =z_k-h\left({\mathit{x}}_{k|k-1}\right)\hfill \end{array}$$

(40)

$$\begin{array}{ccc}\hfill & {\mathit{G}}_k\hfill & ={\mathit{P}}_{k|k-1}{\mathit{h}}_k^T{({\mathit{h}}_k{\mathit{P}}_{k|k-1}{\mathit{h}}_k^T+{\sigma }^2)}^{-1}\hfill \end{array}$$

(41)

$$\begin{array}{cccc}\hfill & {\mathit{x}}_{k|k}\hfill & ={\mathit{x}}_{k|k-1}+{\mathit{G}}_k{\tilde{z}}_k\hfill & \mathrm{state}\mathrm{estimate}\hfill \end{array}$$

(42)

$$\begin{array}{cccc}\hfill & {\mathit{P}}_{k|k}\hfill & =(\mathit{I}-{\mathit{G}}_k{\mathit{h}}_k){\mathit{P}}_{k|k-1}\hfill & \mathrm{covariance}\mathrm{of}\mathrm{state}\mathrm{estimate}\hfill \end{array}$$

(43)

where ${\mathit{x}}_{k|k-1}$ is the state prediction at time k given all measurements up to time $k-1$, and ${\mathit{x}}_{k|k}$ is the filtered state estimate at time k. The EKF replaces the nonlinear measurement Equation (37) with

$$z_k={\mathit{h}}_k{\mathit{x}}_k+w_k$$

(44)

where ${\mathit{h}}_k$ is the Jacobian of $h\left({\mathit{x}}_k\right)$ evaluated at ${\mathit{x}}_{k|k-1}$:

$${\mathit{h}}_k=\frac{1}{d_{k|k-1}}\left[\begin{array}{cccc}{u}_1\left({\mathit{x}}_{k|k-1}\right)& 0& u_2\left({\mathit{x}}_{k|k-1}\right)& 0\end{array}\right].$$

(45)

Here, $d_{k|k-1}$ is the target–sensor range estimate computed from ${\mathit{x}}_{k|k-1}$ and ${\mathit{s}}_k$, and $\mathit{u}(·)={[u_1(·),u_2(·)]}^T$ is the unit vector defined in (4).

The moving AOA sensor is assumed to travel with a constant velocity, which means that the distance between successive waypoints is constant, i.e., $\parallel {\mathit{s}}_k-{\mathit{s}}_{k-1}\parallel =s$. Assuming a maximum turn-rate of $\pm {\vartheta }_{max}$ in azimuth, the next waypoint is constrained to lie on an arc defined by

$${\mathit{s}}_k={\mathit{s}}_{k-1}+s\upsilon \left({\vartheta }_k\right),|{\vartheta }_k-{\vartheta }_{k-1}|<{\vartheta }_{max}$$

(46)

where $\upsilon \left({\vartheta }_k\right)={[cos{\vartheta }_k,sin{\vartheta }_k]}^T$ is the sensor heading vector with heading angle ${\vartheta }_k$ at time instant k.

The recursive BFIM for the Kalman filter tracking problem is given by [24]

$${\mathsf{\Phi }}_k={(\mathit{F}{\mathsf{\Phi }}_{k-1}{\mathit{F}}^T+\mathit{Q})}^{-1}+\frac{1}{{\sigma }^2}{E}_{{\mathit{x}}_k}\left\{{\tilde{\mathit{h}}}_k^T{\tilde{\mathit{h}}}_k\right\}$$

(47)

where ${\tilde{\mathit{h}}}_k$ is the Jacobian matrix in (45) calculated at target state ${\mathit{x}}_k$. Equation (47) has the same structure as (3) in that it is the sum of prior information ${(\mathit{F}{\mathsf{\Phi }}_{k-1}{\mathit{F}}^T+\mathit{Q})}^{-1}$ and contribution from measurements $\frac{1}{{\sigma }^2}{E}_{{\mathit{x}}_k}\left\{{\tilde{\mathit{h}}}_k^T{\tilde{\mathit{h}}}_k\right\}$. It is necessary to simplify (47) so that readily available Kalman filter estimates can be used rather than resorting to computationally expensive Monte Carlo simulations to calculate the expectation.

The contribution from measurements is approximated by

$$\frac{1}{{\sigma }^2}{E}_{{\mathit{x}}_k}\left\{{\tilde{\mathit{h}}}_k^T{\tilde{\mathit{h}}}_k\right\}\approx \frac{1}{{\sigma }^2}{\mathit{h}}_k^T{\mathit{h}}_k$$

(48)

where ${\mathit{h}}_k$ is the Jacobian in (45). Thus, using ${\mathit{P}}_{k|k-1}^{-1}$ as the prior information and (48) as the contribution from measurements, we have

$${\mathsf{\Phi }}_k\approx {\mathit{P}}_{k|k-1}^{-1}+\frac{1}{{\sigma }^2}{\mathit{h}}_k^T{\mathit{h}}_k$$

(49)

where both ${\mathit{P}}_{k|k-1}$ and ${\mathit{h}}_k$ are calculated by the EKF. In the following subsections, we show how to apply the D- and A-optimality criteria to the approximate recursive BFIM expression in (49) to derive trajectory optimization algorithms.

3.1. Sensor Trajectory Optimization Using D-Optimality

Referring to (21) and (49), the optimal waypoint for the sensor at time k using the D-optimality for the EKF takes the following form:

$$\underset{{\mathit{s}}_k\in {\mathcal{S}}_k}{max}{\mathit{h}}_k{\mathit{P}}_{k|k-1}{\mathit{h}}_k^T$$

(50)

which can be rewritten as

$$\underset{{\mathit{s}}_k\in {\mathcal{S}}_k}{max}\frac{1}{d_{k|k-1}^2}{\mathit{u}}^T\left({\mathit{x}}_{k|k-1}\right){\mathit{P}}_{\mathrm{loc},k|k-1}\mathit{u}\left({\mathit{x}}_{k|k-1}\right)$$

(51)

where ${\mathcal{S}}_k$ is the set of permissible waypoints compliant with the turn-rate

$${\mathcal{S}}_k=\{{\mathit{s}}_k|{\mathit{s}}_k-{\mathit{s}}_{k-1}=s\upsilon \left({\vartheta }_k\right),|{\vartheta }_k-{\vartheta }_{k-1}|<{\vartheta }_{max}\}$$

(52)

and

$${\mathit{P}}_{\mathrm{loc},k|k-1}=\left[\begin{array}{cc}{p}_{\mathrm{loc},k|k-1}(1,1)& p_{\mathrm{loc},k|k-1}(1,2)\\ p_{\mathrm{loc},k|k-1}(2,1)& p_{\mathrm{loc},k|k-1}(2,2)\end{array}\right]$$

(53)

is the $2\times 2$ covariance matrix for predicted target location, which is extracted from ${\mathit{P}}_{k|k-1}$ as shown below:

$${\mathit{P}}_{k|k-1}=\left[\begin{array}{cccc}{p}_{\mathrm{loc},k|k-1}(1,1)& \ast & p_{\mathrm{loc},k|k-1}(1,2)& \ast \\ \ast & \ast & \ast & \ast \\ p_{\mathrm{loc},k|k-1}(2,1)& \ast & p_{\mathrm{loc},k|k-1}(2,2)& \ast \\ \ast & \ast & \ast & \ast \end{array}\right].$$

(54)

Note that as $d_{k|k-1}$ also depends on ${\mathit{s}}_k$ (i.e., the fixed range constraint does not apply), (51) does not have a simple closed-form solution. This means that it must be solved by a numerical search over a finite number of permissible waypoints contained in the set ${\mathcal{S}}_k$.

3.2. Sensor Trajectory Optimization Using A-Optimality

Using (28), the A-optimality criterion for sensor waypoints is given by

$$\underset{{\mathit{s}}_k\in {\mathcal{S}}_k}{max}\frac{{\mathit{h}}_k^T{\mathit{P}}_{k|k-1}^2{\mathit{h}}_k}{{\sigma }^2+{\mathit{h}}_k^T{\mathit{P}}_{k|k-1}{\mathit{h}}_k}$$

(55)

which is obtained by substituting ${\mathit{h}}_k$ for ${\mathit{u}}^T\left({\mathit{x}}_0\right)/d$ and ${\mathit{P}}_{k|k-1}$ for ${\mathit{P}}_0$ into (28).

Different from the D-optimality criterion in (51), (55) depends not only on the target range $d_{k|k-1}$ through the Jacobian ${\mathit{h}}_k$, but also the bearing noise variance ${\sigma }^2$. Again, it is not straightforward to solve (55) for the optimal ${\mathit{s}}_k$. A numerical search over the members of the set ${\mathcal{S}}_k$ is necessary.

3.3. Projection Algorithm: A Closed-Form Solution

The D- and A-optimality solutions for determining optimal sensor waypoints described above can be computationally expensive, especially if the numerical search must be carried out over a large number of candidate waypoints in the set ${\mathcal{S}}_k$. In this subsection, we present an alternative closed-form solution, called the projection algorithm, inspired by the ultimate objective of aligning the sensor with the minor axis of the prior covariance error ellipse.

The idea behind the projection algorithm is illustrated in Figure 5. The next waypoint ${\mathit{s}}_k$ is chosen to guide the sensor towards the closest point ${\mathit{\psi }}_k$ that is aligned with the minor axis of the target location prior ${\mathit{P}}_{\mathrm{loc},k|k-1}$ and is at the same distance from the mean of the prior ${\mathit{x}}_{\mathrm{loc},k|k-1}={[x_{k|k-1},y_{k|k-1}]}^T$ as the estimated target–sensor range $d_{k|k-1}$. The next waypoint is found by projecting the waypoint vector ${\mathit{s}}_k-{\mathit{s}}_{k-1}$ to ${\mathit{\psi }}_k-{\mathit{s}}_{k-1}$ subject to the turn-rate constraint. If the projection causes the sensor heading angle to exceed the turn-rate, the next waypoint is chosen to have the maximum turn-rate. This projection also brings the sensor closer to the target with $d_k<d_{k|k-1}$. The reduction in $d_k$ is proportional to how far ${\mathit{\psi }}_k$ is from ${\mathit{s}}_{k-1}$. If ${\mathit{x}}_{\mathrm{loc},k|k-1}-{\mathit{s}}_{k-1}$ is aligned with the major axis of the error ellipse, which represents the worst geometry, the distance between ${\mathit{s}}_{k-1}$ and ${\mathit{\psi }}_k$ is maximized and $d_k$ will have the maximum reduction.

This behaviour makes sense because the optimality of a target–sensor geometry is improved at the maximum rate by moving the sensor directly towards the optimal sensor location for the given Gaussian prior with mean ${\mathit{x}}_{\mathrm{loc},k|k-1}$ and covariance ${\mathit{P}}_{\mathrm{loc},k|k-1}$.

A detailed description of the projection algorithm is provided in Algorithm 1.

Algorithm 1 Projection algorithm.

Input: ${\mathit{s}}_{k-1}$, ${\mathit{x}}_{\mathrm{loc},k|k-1}$, ${\mathit{P}}_{\mathrm{loc},k|k-1}$, ${\vartheta }_{max}$, ${\vartheta }_{k-1}$, s

Output:  ${\mathit{s}}_k$

Calculate ${\nu }_2$, eigenvector of ${\mathit{P}}_{\mathrm{loc},k|k-1}$ associated with its smallest eigenvalue

Calculate $d_{k|k-1}=\parallel {\mathit{x}}_{\mathrm{loc},k|k-1}-{\mathit{s}}_{k-1}\parallel$

Calculate ${\psi }_k^{\pm }={\mathit{x}}_{\mathrm{loc},k|k-1}\pm d_{k|k-1}{\nu }_2$

If $\parallel {\psi }_k^{+}-{\mathit{s}}_{k-1}\parallel <\parallel {\psi }_k^{-}-{\mathit{s}}_{k-1}\parallel$

${\psi }_k={\psi }_k^{+}$

Else

${\psi }_k={\psi }_k^{-}$

End

If ${\mathit{s}}_{k-1}={\psi }_k$ (sensor is already aligned with prior covariance minor axis)

Return ${\mathit{s}}_k={\mathit{s}}_{k-1}+s{[cos{\vartheta }_{k-1},sin{\vartheta }_{k-1}]}^T$ (no change in sensor heading)

Else

Calculate $\Delta \mathit{s}=s\frac{{\psi }_k-{\mathit{s}}_{k-1}}{\parallel {\psi }_k-{\mathit{s}}_{k-1}\parallel }$

Calculate $\Delta \vartheta =\angle \Delta \mathit{s}-{\vartheta }_{k-1}$, $-\pi \le \Delta \vartheta \le \pi$ (∠ denotes bearing angle)

If $|\Delta \vartheta |<{\vartheta }_{max}$

Return ${\mathit{s}}_k={\mathit{s}}_{k-1}+\Delta \mathit{s}$

Else

Calculate ${\vartheta }_k={\vartheta }_{k-1}+sign(\Delta \vartheta ){\vartheta }_{max}$

Return ${\mathit{s}}_k={\mathit{s}}_{k-1}+s{[cos{\vartheta }_k,sin{\vartheta }_k]}^T$

End

End

4. Simulation Examples

This section presents simulation examples to verify the optimization results and to demonstrate the performance of the sensor trajectory optimization algorithms developed in Section 2 and Section 3. In the first set of simulations, the focus is on optimal target–sensor geometries using the D- and A-optimality criteria for Bayesian target localization. The Gaussian prior has zero mean ${\mathit{x}}_0={[0,0]}^T$ and covariance

$${\mathit{P}}_0=\left[\begin{array}{cc}27.6047& -14.7721\\ -14.7721& 22.3953\end{array}\right]$$

(56)

with eigenvalues ${\lambda }_1=40$ and ${\lambda }_2=10$. The minor axis of the error ellipse make an angle of ${50}^{\circ }$ with the positive x-axis. The AOA sensor is allowed to be located on a circle of radius $d=50$ km centred at the mean of the Gaussian prior ${\mathit{x}}_0$, and the bearing angle noise standard deviation is $\sigma =5^{\circ }$.

Figure 6 shows the D- and A-optimality measures, $|\mathsf{\Phi }|$ and $tr{\mathsf{\Phi }}^{-1}$, respectively, versus the bearing angle $\theta$ in the range $0\le \theta <\pi$ for the approximate and exact BFIM given by (9) and (3), respectively. The exact BFIM was calculated using 50,000 Monte Carlo runs for each bearing angle. As evident from Figure 6, the simulated optimal bearing angles are not significantly different for the approximate and exact BFIM. This is also backed up by the close proximity of the two curves in Figure 6. The optimal bearing angle obtained from the approximate BFIM aligns with the minor axis of the error ellipse at ${\theta }_{\mathrm{opt}}={50}^{\circ }$ for both D- and A-optimality criteria, which is in agreement with the analytical results derived in Section 2.

To confirm that (9) is a valid approximation only for sufficiently small ${\mathit{P}}_0$ compared with the target range, we repeated the previous simulations for ${\mathit{P}}_0$ increased by a factor of two. Figure 7 shows the resulting D- and A-optimality measures for the approximate and exact BFIM. While the approximate and exact BFIM still yield almost the same optimal bearing angles, the curves corresponding to them exhibit significant discrepancy, in particular at bearing angles away from ${\theta }_{\mathrm{opt}}$.

The optimal target–sensor geometry for the simulated scenario is depicted in Figure 8. As expected, at the optimal bearing angle, the target–sensor range vector is perfectly aligned with the minor axis of the error ellipse of the Gaussian prior. Note that there are, in fact, two optimal sensor locations. The other one has the bearing angle ${\theta }_{\mathrm{opt}}-180=-130$ degrees.

In the next set of simulations, we consider sensor trajectory optimization in a bearings-only manoeuvring target tracking problem. The algorithms in (51), (55) and Algorithm 1 are simulated for a single realization of target manoeuvres. The process noise parameter for the target is $q_x=q_y={10}^{-4}$ m${}^2$/s${}^4$. The bearing angle noise is assumed to be $\sigma =4^{\circ }$. The initial target dynamics are

$${\mathit{x}}_{0|-1}={[0,0,0,0]}^T$$

(57)

and

$${\mathit{P}}_{0|-1}=\left[\begin{array}{cccc}19.7500& 0& 9.0933& 0\\ 0& 0& 0& 0\\ 9.0933& 0& 9.2500& 0\\ 0& 0& 0& 0\end{array}\right].$$

(58)

The time interval between successive EKF recursions is set equal to $T=10$ s. The sensor is initially located at ${\mathit{s}}_0=[-26.8116,-22.4976]$ km and moves with a constant speed of 25 m/s (90 km/h). The maximum turn-rate for the sensor is ${30}^{\circ }$ per 10 s. The separation between successive sensor waypoints is $s=0.25$ km. For the sensor trajectory optimization methods based on the D- and A-optimality criteria, 10 uniformly spaced points are used for numerical search over sensor waypoints in the set ${\mathcal{S}}_k$.

The simulated sensor trajectory for the D-optimality criterion in (51) is depicted in Figure 9. The 2-$\sigma$ error ellipses for predicted target location estimates are plotted every 50 time instants. The initial error ellipse is drawn in black and all others are in grey. The optimal sensor trajectory achieves a rapid reduction in the size of error ellipses. Following the initial approach, the sensor chases the target by circling around it.

Figure 10 shows the simulated sensor trajectory for the A-optimality criterion in (55). The optimal sensor trajectory has a markedly different behaviour to that observed for the D-optimality criterion (see Figure 9) in that it seems to favour circling the target more than getting close to it initially. As a consequence, it takes longer to achieve a significant reduction in error ellipses than the A-optimality criterion.

The simulation results for the projection algorithm in Algorithm 1 are shown in Figure 11. The sensor follows a more direct route towards the target than both the D- and A-optimality methods, followed by circling manoeuvres. This is expected to produce faster convergence to the minimum estimation error than the D- and A-optimality methods at the expense of somewhat larger estimation error initially. This observation is confirmed by the root-mean-square error (RMSE) of target location estimates shown in Figure 12, which were computed from 5000 Monte Carlo simulation runs.

5. Conclusions

Optimal sensor placement for the Bayesian AOA target localization problem was considered. The concept of updating prior information with data from measurements was extended to the Kalman filter tracking problem. Using an approximation of the BFIM, the optimal sensor placement for a given Gaussian prior was shown to be aligned with the minor axis of the prior error ellipse for both D- and A-optimality criteria. By way of simulations, this result was shown to match the optimal sensor placement result for the exact BFIM, which was numerically calculated using Monte Carlo simulations. The D- and A-optimality criteria were adopted for optimal sensor guidance in target tracking applications. Simple methods requiring a numerical search over sensor heading were developed and demonstrated in simulations. A new method, called the projection algorithm, was also developed predicated on the optimal sensor placement result for the D- and A-optimality criteria. The efficacy of the projection algorithm in achieving fast minimization of tracking error was confirmed by numerical simulations. Even though the D- and A-optimality criteria share the same optimal sensor placement result for a stationary target, they generate quite different sensor trajectories in a target tracking setting where optimal target–sensor geometries cannot be realized instantaneously because of constrained sensor dynamics. The A-optimality method was observed to favour circular motion around the target initially. This results in delays in minimizing the tracking error to a desirable level.