Constrained Extended Kalman Filter for Target Tracking in Directional Sensor Networks

Abstract

The target tracking problem in directional sensor networks (DSNs) is attracting increasing attention. Unlike the traditional omnidirectional sensor, a directional sensor has a special angle of view. It can offer direction information rather than just the sensing signal measurement with respect to the detected target. The existing tracking approaches in DSNs always separately consider the direction and measurement information; they hardly promise the tracking performance of minimum variance. In this paper, the field of view of directional sensor is approximated to a rectangle; as such the constrained area in which the target is bound to be is constructed. Then, the target tracking problem is formulated as a constrained estimation problem, and a constrained extended Kalman filter (CEKF) tracking algorithm integrating the direction and measurement information is presented; its structural and statistical properties are rigorously derived. It is proved that CEKF is the linear unbiased minimum variance estimator, and CEKF can yield a smaller error covariance than the unconstrained traditional extended Kalman filter using only sensor measurements. Simulation results show that the CEKF has superior tracking performance for directional wireless networks.

1. Introduction

Wireless sensors networks (WSNs) consist of a large number of miniature sensors that integrate physical sensing, data processing, memory, and communication capabilities, promise high-resolution spatial and temporal monitoring of expansive sensing field over extended duration, and hence are an enabling technology of modern cyber physical systems. WSNs have a wide range of civil and military applications such as target tracking, infrastructure monitoring, habitat sensing, and battlefield surveillance [1–3].

Among many tasks of WSNs, tracking moving targets in a sensing field is a challenging yet ubiquitous problem that has attracted much attention. Most of the past work is always based on the assumption of omnidirectional sensors that has an omniangle of sensing range [4–6]. Recently, directional sensors such as camera/video sensors [7], ultrasonic sensors [8], and infrared sensors [9] have been developed and networks consisting of such sensors, directional sensor networks (DSNs), are widely used. Each node in DSNs is composed of multiple directional sensors that have restricted fields of view as angular sectors, and only one orientation can work for a target at each instance. The directional sensor can identify in which sector the target is present or absent, which gives rough direction information of the target. In this work, we consider general directional sensor networks; one node contains an arbitrary number of collocated directional sensors, which is motivated by the fact that the directional sensors have multiple types of field of view angles.

A number of tracking approaches using binary information regarding target's presence or absence within its sensing range have been proposed in recent years. Several binary sensor models are presented to describe whether the object is moving toward or away from the sensor [10], identify targets by the signal strength sensing modalities [11], or explore the fundamental limits of tracking performance with binary proximity sensors [12]. Then, the central unit implements the tracking by some particle filtering style algorithms based on the binary information from each sensor. In [13], each point on the targets' path is computed using the weighted average of the detecting sensors, locations. Then, a line that best fits the newly estimated location and the points on the trajectory established in the recent past is used as the target trajectory. In [14], a tracking method was proposed which first estimates the positions of a target in its most recent past and then fits them with a piece-wise trajectory. However, the abovementioned works discuss only the target tracking problem under the WSNs with omnidirectional sensors. The omnidirectional sensor can only detect the target presence or absence within its sensing range but it cannot get any direction information of the target.

Unlike the traditional omnidirectional sensor, the directional sensor senses the target only in a particular direction. The node could generate sensing direction and measurement information of the target. This information contributes to yielding more accurate results for target tracking compared to omnidirectional sensor networks. In [15], a directional binary sensing model is introduced; each sensor node's sensing region is divided into four angular sectors and each node can identify in which sector the target is present or absent. A real-time and distributed target tracking algorithm with the directional binary sensor networks is given. It combines all angles in the corresponding neighbor list to determine the arc that the target is crossing; the target location is estimated as the middle point of the corresponding arc. In [16], an efficient target tracking solution with the help of beam steering enabled directional sensors is proposed. It adopts the simple centroid algorithm and trilateration localization to compute the location of the target by using the ranges and angles' information with respect to sensor node locations. In [17], the field of vision of each directional sensor is modeled as a line; the data are the times at which sensors detect targets crossing their lines. The estimation problem is formulated as a highly nonconvex optimization problem, and a three-phase algorithm is developed to solve the optimization problem by using Newton's method. However, these methods [15–17] are based on geometric arguments of the sensor locations without taking into account the sensing signal measurements; they hardly obtain nice tracking performance without using prediction and filtering algorithms. In [18], a noncentralized method is put forward for cooperative target tracking in a sensor network where each node contains four 90-degree directional sensors. Each node processes the signal measurements and its own sensors and the binary information about the neighbor nodes that have detecting sensors by way of particle filtering. However, the adopted particle filter is expensive to compute, and this approach does not fully exploit the geometric characteristics of the network itself, ignoring the direction information of the target.

In this work, the direction and measurement information for target tracking is integrated in the Bayesian framework. Firstly, the sector sensing field of directional sensor is approximated to a rectangle; as such the constrained area in which the target is bound to be is constructed. After that, the target tracking problem is formulated as a constrained estimation problem; that is, the target estimated location is restrained within the constrained area. Then, a constrained extended Kalman filter tracking algorithm is presented, and its structural and statistical properties are rigorously derived. It is proved that the proposed constrained extended Kalman filter is the linear unbiased minimum variance estimator, and it can yield a smaller error covariance than the unconstrained traditional extended Kalman filter using only sensor measurement information.

The remainder of the paper is organized as follows. In Section 2, we review the target motion model, establish the directional measurement model, and formulate the state estimation problem in DSNs. A constrained extended Kalman filter for target tracking in DSNs is developed in Section 3. In Section 4, the structural and statistical properties of constrained extended Kalman filter are analyzed. We conduct performance evaluations of the proposed algorithm by simulation comparisons in Section 5. Conclusions are reached in Section 6.

2. Target Tracking Models and Problem Formulation in DSNs

2.1. Target Motion Model

The motion of a moving target can be described by a 2D state space equation [19, 20]: $$\begin{array}{cc}& X_k=F_{k-\mathrm{1}}{X}_{k-\mathrm{1}}+G_{k-\mathrm{1}}{W}_{k-\mathrm{1}},\end{array}$$ (1)where the state vector $X_k=[x_k,{\dot{x}}_k,y_k,{\dot{y}}_k{]}^T$ specifies the position $(x_k,y_k)$ and velocity $({\dot{x}}_k,{\dot{y}}_k)$ of the target along x-direction and y-direction at $k\text{th}$ time step, respectively. $F_{k-\mathrm{1}}$ and $G_{k-\mathrm{1}}$ are the state transition matrix and the input matrix, respectively. For tracking applications, one has $$\begin{array}{cc}& F_{k-\mathrm{1}}=\left[\begin{array}{cccc}\mathrm{1}& {\Delta }_{k-\mathrm{1}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& {\Delta }_{k-\mathrm{1}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right],\hspace{1em}\hspace{1em}{G}_{k-\mathrm{1}}=\left[\begin{array}{cc}\frac{{\Delta }_{k-\mathrm{1}}^{\mathrm{2}}}{\mathrm{2}}& \mathrm{0}\\ {\Delta }_{k-\mathrm{1}}& \mathrm{0}\\ \mathrm{0}& \frac{{\Delta }_{k-\mathrm{1}}^{\mathrm{2}}}{\mathrm{2}}\\ \mathrm{0}& {\Delta }_{k-\mathrm{1}}\end{array}\right].\end{array}$$ (2)In the above, ${\Delta }_{k-\mathrm{1}}$ is the sampling time interval between two successive time steps k and $k-\mathrm{1}$ . The input noise $W_{k-\mathrm{1}}=[w_x,w_y{]}^T$ is a $\mathrm{2}\times \mathrm{1}$ Gaussian random vector with zero mean and a covariance matrix $Q_{k-\mathrm{1}}$ ; $w_x$ and $w_y$ correspond to the process noises (random accelerations) along the x-axis and y-axis, respectively. Generally, it is assumed that $w_x$ is uncorrelated with $w_y$ ; then $Q_{k-\mathrm{1}}=\text{Diag}\{{\sigma }_x^{\mathrm{2}},{\sigma }_y^{\mathrm{2}}\}$ .

2.2. Directional Measurement Model

Unlike the omnidirectional sensor, the directional sensor has restricted angles of sensing range. The advantage of directional sensor is that the direction of the target could be acquired rather than just the sensing signal measurement.

Let us consider a DSN composed of N identical nodes in a two-dimensional plane; the position of the $n\text{th}$ ( $\mathrm{1}\le n\le N$ ) node is denoted by $(x^n,y^n)$ . Each of the nodes contains M collocated directional sensors providing a coverage of 360 degrees and operates in only one direction with a uniform sensing range r; we assume that a node has no overlap between neighboring directions as shown in Figure 1.

Figure 1

Structure of directional sensor node.

When the node-to-target distance $h_n(k)=\sqrt{(x_k-x^n{)}^{\mathrm{2}}+(y_k-y^n{)}^{\mathrm{2}}}$ is less than r, the target will be detected and one of its own sensors is triggered; the distance $h_n(k)$ will be estimated. We shall denote ${\Omega }_k=\left\{n\hspace{0.17em}|\hspace{0.17em}{h}_n(k)\le r,\hspace{0.17em}\hspace{0.17em}\mathrm{1}\le n\le N\right\}$ to be the set of indices of nodes within the sensing range of the target at time k; the measurement ${\mathbf{z}}_n(k)$ of the $n\text{th}$ node at time k can be expressed as [21, 22] $$\begin{array}{cc}& {\mathbf{z}}_n\left(k\right)=h_n\left(k\right)+v_n\left(k\right)\hspace{1em}n\in {\Omega }_k,\end{array}$$ (3)where $v_n(k)$ is the observation noise at the $n\text{th}$ node, which is assumed to be an independent Gaussian random variable with zero mean and variance ${\sigma }_n^{\mathrm{2}}$ for each node and uncorrelated with $W_{k-\mathrm{1}}$ and $X_k$ .

Denote ${\mathcal{l}}_k$ to be the cardinal number of ${\Omega }_k$ (i.e., the number of elements in ${\Omega }_k$ ); the sensor measurements at the $k\text{th}$ time step may be represented in a matrix form: $$\begin{array}{cc}& Z_k=H_k\left(X_k\right)+V_k,\end{array}$$ (4)where $H_k(X_k)=[h_{\mathrm{1}}(k),\dots ,h_{{\mathcal{l}}_k}(k){]}^T$ and $V_k=[v_{\mathrm{1}}(k),\dots ,v_{{\mathcal{l}}_k}(k){]}^T$ . Moreover, the covariance matrix of $V_k$ is $R_k=\text{Diag}\{{\sigma }_n^{\mathrm{2}};\mathrm{1}\le n\le {\mathcal{l}}_k\}$ .

While the node gets the measurement of the target, its own triggered sensor provides the direction information of the target according to the angle of view in which the target is present. By intersecting the angles of view of triggered sensors that detected the target, we can have preliminary knowledge of the constrained area the target is restrained in. As shown in Figure 2, we illustrate two target detection scenarios in DSNs with $M=\mathrm{4}$ and $M=\mathrm{3}$ . The overlap area of the sensing directions of sensors is a convex set, because there is an unobstructed straight path between any two points in this area. However, it is found that the closed-form algebraic expression of the overlap area is an apparently hitherto unsolved geometrical problem [23]. By dragging a rectangle around the sensing direction of a sensor, we can get a relaxed set that has a closed-form rectangular expression for the overlap area. As shown in Figure 1, the ith sensing direction of the $n\text{th}$ node is represented as $$\begin{array}{cc}& {\mathcal{A}}_{n,i}\triangleq \left\{\left(x,y\right)\hspace{0.17em}|\hspace{0.17em}{l}_{x^n,i}\le x\le u_{x^n,i},\hspace{0.17em}\hspace{0.17em}{l}_{y^n,i}\le y\le u_{y^n,i}\right\},\end{array}$$ (5)where $l_{x^n,i},u_{x^n,i},l_{y^n,i},u_{y^n,i}$ are the lines that are tangential to the boundary edges of the fan-shaped sensing direction of the sensor. Then, the indices set of sensing directions of the nodes where the target locates at time k can be defined as $$\begin{array}{cc}& D_k\triangleq \left\{\left(n,i\right)\hspace{0.17em}|\hspace{0.17em}\left(x_k,y_k\right)\in {\mathcal{A}}_{n,i},\hspace{0.17em}\hspace{0.17em}n\in {\Omega }_k,\hspace{0.17em}\hspace{0.17em}i=\mathrm{1},\dots ,M\right\}.\end{array}$$ (6)

Figure 2

Target detection scenario in a DSN.

By intersection operations, the constrained area of the target at the $k\text{th}$ time step can be described in an approximate rectangular set: $$\begin{array}{l}{\mathcal{B}}_k={\displaystyle \bigcap }_{D_k}‍{\mathcal{A}}_{n,i}\\ =\left\{\left(x,y\right)\hspace{0.17em}|\hspace{0.17em}{x}_{\min }\left(k\right)\le x\le x_{\max }\left(k\right),\right.\\ \hspace{1em}\hspace{1em}\left.y_{\min }\left(k\right)\le y\le y_{\max }\left(k\right)\right\},\end{array}$$ (7)where $$\begin{array}{cc}& x_{\min }\left(k\right)=\underset{\left(n,i\right)\in D_k}{\mathrm{Max}}\left\{l_{x^n,i}\right\},\hspace{1em}\hspace{1em}{x}_{\max }\left(k\right)=\underset{\left(n,i\right)\in D_k}{\mathrm{Min}}\left\{u_{x^n,i}\right\},\\ y_{\min }\left(k\right)=\underset{\left(n,i\right)\in D_k}{\mathrm{Max}}\left\{l_{y^n,i}\right\},\hspace{1em}\hspace{1em}{y}_{\max }\left(k\right)=\underset{\left(n,i\right)\in D_k}{\mathrm{Min}}\left\{u_{y^n,i}\right\}.\end{array}$$ (8)In above equation, $\mathrm{Max}\{·\}$ and $\mathrm{Min}\{·\}$ , respectively, return the biggest and smallest elements in an array. The formula of closed-form rectangle is general enough to cover all possible cases of the overlap of any number of angular sectors. The particular process to model the directional measurement is summarized in Algorithm 1.

Algorithm 1: The directional measurement modeling algorithm.

Require: The sensor positions $$\{(x^n,y^n)\}$$ , and the target state $$X_k$$ .

Ensure: The distance measurements $$Z_k$$ and directional constrained area $${\mathcal{B}}_k$$

(1) Obtain the set of indices of active nodes $${\mathrm{\Omega }}_k$$ ,

(2) Acquire the detecting distance measurements $$Z_k$$ using (4),

(3) Approximate the sensing direction $${\mathcal{A}}_{n,i}$$ by a rectangle using (5),

(4) Define the indices set of triggered sensors $$D_k$$ using (6),

(5) Construct the constrained area $${\mathcal{B}}_k$$ by intersections using (7).

2.3. Problem Formulation

For simplicity, we only consider the problem of tracking a single target, but nevertheless our proposed tracking algorithm still is applicable to multitarget tracking as the targets have been classified. When the target moves through the monitored area, these nodes which have detected the target form a cluster [24, 25]. One of the nodes in the cluster will be selected as the cluster head which also serves as the fusion center of signal and information processing. Cluster members measure their respective node-to-target distances and transmit the distances to the cluster head with the directions and positions' information of their own triggered sensors that detected the target. For convenience, transmission delay or packet loss is ignored. Given the noisy sensor observations $Z_{\mathrm{1}:k}=\{Z_{\mathrm{1}},Z_{\mathrm{2}},\dots ,Z_k\}$ from initial time up to and including time k and directional constrained area ${\mathcal{B}}_k$ , our goal is to obtain a sequential Bayesian estimate of $X_k$ using a KF formulation. This can be formulated as a maximum a posteriori probability (MAP) estimation problem with a constraint on the state vector $X_k$ : $$\begin{array}{l}\mathrm{Max}\hspace{1em}\hspace{1em}\hspace{1em}p\left(X_k\hspace{0.17em}|\hspace{0.17em}{Z}_{\mathrm{1}:k}\right)\\ \text{subject}\mathrm{\hspace{0.17em}\hspace{0.17em}}\text{to}\hspace{1em}\left(x_k,y_k\right)\in {\mathcal{B}}_k,\end{array}$$ (9)where $p(X_k\hspace{0.17em}|\hspace{0.17em}{Z}_{\mathrm{1}:k})$ is the posterior distribution of $X_k$ given the observed value $Z_{\mathrm{1}:k}$ .

3. Constrained Extended Kalman Filter in DSNs

The basic filtering solution to the state estimation problem can be described as a two-stage recursive process of prediction and update. Now, given the distribution $p(X_{k-\mathrm{1}}\hspace{0.17em}|\hspace{0.17em}{Z}_{\mathrm{1}:k-\mathrm{1}})=\mathcal{N}({\widehat{X}}_{k-\mathrm{1}\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}},P_{k-\mathrm{1}\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}})$ , the joint distribution of $X_k$ and $X_{k-\mathrm{1}}$ conditioned on the previous measurements $Z_{\mathrm{1}:k-\mathrm{1}}=\{Z_{\mathrm{1}},Z_{\mathrm{2}},\dots ,Z_{k-\mathrm{1}}\}$ is $$\begin{array}{l}p\left(X_{k-\mathrm{1}},X_k\hspace{0.17em}|\hspace{0.17em}{Z}_{\mathrm{1}:k-\mathrm{1}}\right)\\ =p\left(X_k\hspace{0.17em}|\hspace{0.17em}{X}_{k-\mathrm{1}}\right)p\left(X_{k-\mathrm{1}}\hspace{0.17em}|\hspace{0.17em}{Z}_{\mathrm{1}:k-\mathrm{1}}\right)\\ =\mathcal{N}\left(\left[\begin{array}{l}{\widehat{X}}_{\left.k-\mathrm{1}\right|k-\mathrm{1}}\\ F_{k-\mathrm{1}}{\widehat{X}}_{\left.k-\mathrm{1}\right|k-\mathrm{1}}\end{array}\right],\right.\\ \left.\hspace{1em}\hspace{1em}\hspace{1em}\left[\begin{array}{ll}{P}_{k-\mathrm{1}\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}& P_{k-\mathrm{1}\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}{F}_{k-\mathrm{1}}^T\\ F_{k-\mathrm{1}}{P}_{k-\mathrm{1}\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}& F_{k-\mathrm{1}}{P}_{k-\mathrm{1}\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}{F}_{k-\mathrm{1}}^T+G_{k-\mathrm{1}}{Q}_{k-\mathrm{1}}{G}_{k-\mathrm{1}}^T\end{array}\right]\right).\end{array}$$ (10)The predictive distribution of $X_k$ given the measurement history up to time step $k-\mathrm{1}$ can be calculated by Lemma A.1 in Appendix. Consider $$\begin{array}{l}p\left(X_k\hspace{0.17em}|\hspace{0.17em}{Z}_{\mathrm{1}:k-\mathrm{1}}\right)\\ =\mathcal{N}\left({\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}},P_{\left.k\right|k-\mathrm{1}}\right)\\ =\mathcal{N}\left(F_{k-\mathrm{1}}{\widehat{X}}_{k-\mathrm{1}\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}},F_{k-\mathrm{1}}{P}_{k-\mathrm{1}\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}{F}_{k-\mathrm{1}}^T+G_{k-\mathrm{1}}{Q}_{k-\mathrm{1}}{G}_{k-\mathrm{1}}^T\right).\end{array}$$ (11)

To update the state estimate using latest observation $Z_k$ , consider a first order Taylor series expansion of $H_k(X_k)$ around ${\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}$ ; thus $$\begin{array}{cc}& Z_k\approx H_k\left({\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}\right)+{\dot{H}}_k\left(X_k-{\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}\right)+V_k,\end{array}$$ (12)where ${\dot{H}}_k=\left[{\left.\left(\partial H_k(X_k)/\partial X_k^T\right)\right|}_{X_k={\widehat{X}}_{\left.k\right|k-\mathrm{1}}}\right]$ is a ${\mathcal{l}}_k\times \mathrm{4}$ matrix.

According to (12), one has the parameters of the posterior distribution by Bayes rule: $$\begin{array}{cc}& p\left(X_k\hspace{0.17em}|\hspace{0.17em}{Z}_{\mathrm{1}:k}\right)\propto p\left(Z_k\hspace{0.17em}|\hspace{0.17em}{X}_k\right)p\left(\left.X_k\right|Z_{\mathrm{1}:k-\mathrm{1}}\right),\end{array}$$ (13)where $$\begin{array}{l}p\left(Z_k\hspace{0.17em}|\hspace{0.17em}{X}_k\right)\\ =\frac{\mathrm{1}}{{\left(\mathrm{2}\pi \right)}^{{\mathcal{l}}_k/\mathrm{2}}{\left(\det \left\{R_k\right\}\right)}^{\mathrm{1}/\mathrm{2}}}\\ \hspace{1em}·\exp \left\{-\frac{\mathrm{1}}{\mathrm{2}}{\left(Z_k-{\dot{H}}_k{X}_k-H_k\left({\widehat{X}}_{k|k-\mathrm{1}}\right)+{\dot{H}}_k{\widehat{X}}_{k|k-\mathrm{1}}\right)}^T\right.\\ \left.\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\times \hspace{0.17em}{R}_k^{-\mathrm{1}}\left(Z_k-{\dot{H}}_k{X}_k-H_k\left({\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}\right)+{\dot{H}}_k{\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}\right)\phantom{\frac{\mathrm{1}}{\mathrm{2}}}\right\},\\ p\left(X_k\hspace{0.17em}|\hspace{0.17em}{Z}_{\mathrm{1}:k-\mathrm{1}}\right)\\ =\frac{\mathrm{1}}{{\left(\mathrm{2}\pi \right)}^{\mathrm{2}}{\left(\det \left\{P_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}\right\}\right)}^{\mathrm{1}/\mathrm{2}}}\\ \hspace{1em}·\exp \left\{-\frac{\mathrm{1}}{\mathrm{2}}{\left(X_k-{\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}\right)}^T{P}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}^{-\mathrm{1}}\left(X_k-{\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}\right)\right\}.\end{array}$$ (14)

The MAP estimate associated with the posterior $p(X_k\hspace{0.17em}|\hspace{0.17em}{Z}_{\mathrm{1}:k})$ in (13) is the vector $X_k$ that maximizes the likelihood $p(X_k\hspace{0.17em}|\hspace{0.17em}{Z}_{\mathrm{1}:k})$ , which is identical to minimizing its negative logarithm: $$\begin{array}{l}L\left(X_k\right)={\left[Z_k-{\dot{H}}_k{X}_k-H_k\left({\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}\right)+{\dot{H}}_k{\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}\right]}^T\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\times R_k^{-\mathrm{1}}\left[Z_k-{\dot{H}}_k{X}_k-H_k\left({\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}\right)+{\dot{H}}_k{\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}\right]\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\left(X_k-{\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}\right)}^T{P}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}^{-\mathrm{1}}\left(X_k-{\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}\right)+c_{\mathrm{1}}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=X_k^T\left({\dot{H}}_k^T{R}_k^{-\mathrm{1}}{\dot{H}}_k+P_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}^{-\mathrm{1}}\right)X_k\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}-\mathrm{2}\left[{\left(Z_k+{\dot{H}}_k{\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}-H_k\left({\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}\right)\right)}^T\right.\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\left.\hspace{1em}\phantom{{\left(Z_k+{\dot{H}}_k{\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}-H_k\left({\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}\right)\right)}^T}\times R_k^{-\mathrm{1}}{\dot{H}}_k+{\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}^T{P}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}^{-\mathrm{1}}\right]X_k+c_{\mathrm{2}},\end{array}$$ (15)where $c_{\mathrm{1}}$ and $c_{\mathrm{2}}$ are constant terms not dependent on $X_k$ .

Let $$\begin{array}{cc}& B_k={\dot{H}}_k^T{R}_k^{-\mathrm{1}}{\dot{H}}_k+P_{k|k-\mathrm{1}}^{-\mathrm{1}}\\ d_k=P_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}^{-\mathrm{1}}{\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}+{\dot{H}}_k^T{R}_k^{-\mathrm{1}}\left[Z_k+{\dot{H}}_k{\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}-H_k\left({\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}\right)\right].\end{array}$$ (16)

Equation (15) can be simplified into $$\begin{array}{cc}& L\left(X_k\right)=X_k^T{B}_k{X}_k-\mathrm{2}{d}_k^T{X}_k+c_{\mathrm{2}}.\end{array}$$ (17)

The Kalman filter estimate for target tracking in DSNs can be derived by finding an $X_k$ such that the logarithm function is minimized and $X_k$ locates in the directional constrained area ${\mathcal{B}}_k$ . So the MAP estimation problem can be expressed as $$\begin{array}{l}\mathrm{Min}\hspace{1em}\hspace{1em}\hspace{1em}{X}_k^T{B}_k{X}_k-\mathrm{2}{d}_k^T{X}_k\\ \text{subject}\mathrm{\hspace{0.17em}\hspace{0.17em}}\text{to}\hspace{1em}\left(x_k,y_k\right)\in {\mathcal{B}}_k.\end{array}$$ (18)

Note that this problem is called a relaxation of the target tracking problem in DSNs, since we have replaced the original overlap area constraints with a looser rectangular area constraint. It is obvious that the optimal value of (18) is less than or equal to the optimal value of original tracking problem, since we minimize the same objective function over a larger set ${\mathcal{B}}_k$ .

Theorem 1.

The estimation problem defined in (18) is convex.

Proof.

Since ${\mathcal{B}}_k$ is a convex set, noting that $B_k$ is a positive definite matrix, then the objective function $X_k^T{B}_k{X}_k-\mathrm{2}{d}_k^T{X}_k$ is strictly convex; as such we immediately claim that the estimation problem is a convex optimization problem.

Thus, the KKT necessary conditions in (18) at the minimum $X_k^{\mathrm{\ast }}$ , ${\lambda }_k^{\mathrm{\ast }}$ , and ${\mu }_k^{\mathrm{\ast }}$ require that $$\begin{array}{cc}& B_k{X}_k^{\mathrm{\ast }}-d_k-{\lambda }_k^{\mathrm{\ast }}+{\mu }_k^{\mathrm{\ast }}=\mathrm{0}\\ {\lambda }_k^{\mathrm{\ast }}\ge \mathrm{0},\hspace{1em}\hspace{1em}{\mu }_k^{\mathrm{\ast }}\ge \mathrm{0}\\ {\lambda }_k^{\mathrm{\ast }T}\left(l_k-X_k^{\mathrm{\ast }}\right)=\mathrm{0}\\ {\mu }_k^{\mathrm{\ast }T}\left(X_k^{\mathrm{\ast }}-u_k\right)=\mathrm{0}\\ l_k\le X_k^{\mathrm{\ast }}\le u_k,\end{array}$$ (19)where $l_k=[x_{\min }(k),-\infty ,y_{\min }(k),-\infty {]}^T$ and $u_k=[x_{\max }(k),\infty ,y_{\max }(k),\infty {]}^T$ . The problem (18) is a box-constrained quadratic programming problem of optimization problems, and there are two major strategies existing in the literatures. The first one is the active set strategy [26], which generates iterates on a face of the feasible box until either a minimizer of the objective function on that face or a point on the boundary of that face is reached. The other solution strategy is the interior point algorithm [27], which consists of a series of parameterized barrier functions that are minimized using Newton's method.

By an interior reflective Newton method [28] that exhibits strong convergence properties and effectively deals with infinite bounds, called Quadprog (see the details in [28]), the solution to the above problem can be briefly described in Algorithm 2.

Algorithm 2: Interior reflective Newton algorithm.

Initially set: the iteration step $$t=\mathrm{0}$$ , $${\mathrm{\lambda }}_k^{(\mathrm{0})}=\mathrm{0}$$ , $${\mathrm{\mu }}_k^{(\mathrm{0})}=\mathrm{0}$$ and $$X_k^{(\mathrm{0})}=-B_k^{-\mathrm{1}}{d}_k$$ . If $$X_k^{(\mathrm{0})}$$ is feasible, then the solution

is: $$X_k^{\mathrm{*}}=X_k^{(\mathrm{0})}$$ , $${\mathrm{\lambda }}_k^{\mathrm{*}}={\mathrm{\lambda }}_k^{(\mathrm{0})}$$ , $${\mathrm{\mu }}_k^{\mathrm{*}}={\mathrm{\mu }}_k^{(\mathrm{0})}$$ .

Looping: If $$X_k^{(\mathrm{0})}$$ is not feasible, the components of $$X_k^{(\mathrm{0})}$$ that violate the bound constraints are reset to the midpoint of the box

defined by the bounds, and then uses Quadprog method

           $$\left[X_k^{(t)},{\mathrm{\lambda }}_k^{(t)},{\mathrm{\mu }}_k^{(t)}\right]=\mathrm{\hspace{0.17em}\hspace{0.17em}}\text{Quadprog}\mathrm{\hspace{0.17em}\hspace{0.17em}}\left\{B_k,d_k,l_k,u_k,X_k^{(\mathrm{0})}\right\}$$

Stopping and output: if $$L(X_k^{(t-\mathrm{1})})-L(X_k^{(t)})$$ is smaller than the specified tolerance, then the solution

is: $$X_k^{\mathrm{*}}=X_k^{(t)}$$ , $${\mathrm{\lambda }}_k^{\mathrm{*}}={\mathrm{\lambda }}_k^{(t)}$$ , $${\mathrm{\mu }}_k^{\mathrm{*}}={\mathrm{\mu }}_k^{(t)}$$ .

Remark 2.

Nonzero elements of the Lagrange multipliers ${\lambda }_k^{\mathrm{\ast }}$ and ${\mu }_k^{\mathrm{\ast }}$ indicate active constraints at the solution $X_k^{\mathrm{\ast }}$ ; it means that the corresponding components of the solution are on these constraint boundaries.

4. Properties of the Constrained Extended Kalman Filter in DSNs

In this section, we study the structural properties of CEKF to give the closed-form algebraic expression of the state estimate. Based on the structural property, some statistical properties of CEKF are analysed. In the development below, if A and B are square matrices with the same dimension, we use the notation $A>B$ to indicate that $A-B$ is positive definite. The entry in the ith row and jth column of a matrix A is denoted as $A[i,j]$ .

4.1. Structural Properties of CEKF

Lemma 3.

Given the quadratic optimization problem (18), there exists a unique pair of primal and dual optimal $X_k^{\mathrm{\ast }}$ and ${\lambda }_k^{\mathrm{\ast }}$ , ${\mu }_k^{\mathrm{\ast }}$ satisfy the KKT conditions in (19).

Proof.

Since $B_k>\mathrm{0}$ and $L(X_k)$ is convex, there exists only an optimal point for minimizing $L(X_k)$ with unconstrained conditions; that is, there is a unique solution $X_k^{\star }$ such that $$\begin{array}{l}L\left(X_k^{\star }\right)=\underset{X_k\in {\mathbb{R}}^{\mathrm{4}}}{\mathrm{Min}}L\left(X_k\right)=\underset{X_k\in \left\{{\Psi }_k{\displaystyle \bigcup }‍{\mathbb{R}}^{\mathrm{4}}-{\Psi }_k\right\}}{\mathrm{Min}}L\left(X_k\right),\\ =\mathrm{Min}\left\{\underset{X_k\in {\Psi }_k}{\mathrm{Min}}L\left(X_k\right),\underset{X_k\in {\mathbb{R}}^{\mathrm{4}}-{\Psi }_k}{\mathrm{Min}}L\left(X_k\right)\right\},\end{array}$$ (20)where ${\Psi }_k=\{X_k\hspace{0.17em}|\hspace{0.17em}{l}_k\le X_k\le u_k,\hspace{0.17em}\hspace{0.17em}{X}_k\in {\mathbb{R}}^{\mathrm{4}}\}$ . The above equation indicates that there must be a solution $X_k^{\mathrm{\ast }}$ such that $L(X_k)$ is minimal for $X_k\in {\Psi }_k$ ; it means that the problem (18) is strictly feasible. By the Slaters constraint qualification [29], there exist dual optimal ${\lambda }_k^{\mathrm{\ast }}$ , ${\mu }_k^{\mathrm{\ast }}$ , which together with $X_k^{\mathrm{\ast }}$ satisfy the KKT conditions in (19).

Let $b_k$ be the column vector composed of active constrained boundaries whose Lagrange multiplier elements are nonzero, and ${\Lambda }_k$ is the column vector that contains the indices of elements of ${\lambda }_k^{\mathrm{\ast }}$ and ${\mu }_k^{\mathrm{\ast }}$ corresponding to the active constrained boundaries in $b_k$ ; then, one has $$\begin{array}{cc}& A_k\left[i,j\right]=\left\{\begin{array}{cc}\mathrm{0}\hfill & \text{if}\hspace{0.17em}\hspace{0.17em}{b}_k\left[i,\mathrm{1}\right]\in b_k,\hspace{0.17em}\hspace{0.17em}j\notin {\Lambda }_k\hfill \\ \mathrm{1}\hfill & \text{if}\hspace{0.17em}\hspace{0.17em}{b}_k\left[i,\mathrm{1}\right]\in b_k,\hspace{0.17em}\hspace{0.17em}j\in {\Lambda }_k\hfill \end{array}\right.\hspace{1em}\text{for}\hspace{0.17em}\hspace{0.17em}j=\mathrm{1,2},\dots ,\mathrm{4}.\end{array}$$ (21)

It is noted that $A_k$ has the same number of rows with $b_k$ . From Remark 2, these active boundary conditions may be expressed in terms of the constraints that are satisfied with equality at $X_k^{\mathrm{\ast }}$ ; the solution of (18) would be a solution of the equivalent equality-constrained problem: $$\begin{array}{l}\mathrm{Min}\hspace{1em}\hspace{1em}\hspace{1em}{X}_k^T{B}_k{X}_k-\mathrm{2}{d}_k^T{X}_k\\ \text{subject}\mathrm{\hspace{0.17em}\hspace{0.17em}}\text{to}\hspace{1em}{A}_k{X}_k=b_k.\end{array}$$ (22)

Here, the KKT necessary conditions for a minimum are $$\begin{array}{cc}& B_k{X}_k^{\mathrm{\ast }}-d_k+A_k^T{\nu }_k^{\mathrm{\ast }}=\mathrm{0},\\ A_k{X}_k^{\mathrm{\ast }}=b_k,\end{array}$$ (23)which we can write as $$\begin{array}{cc}& \left[\begin{array}{cc}{B}_k& A_k^T\\ A_k& \mathrm{0}\end{array}\right]\left[\begin{array}{c}{X}_k^{\mathrm{\ast }}\\ {\nu }_k^{\mathrm{\ast }}\end{array}\right]=\left[\begin{array}{c}{d}_k\\ b_k\end{array}\right].\end{array}$$ (24)

Since each component of the minimum $X_k^{\mathrm{\ast }}$ cannot be at the lower and upper bound simultaneously, the matrix $A_k$ is full row rank; as such one can get the solution by the Schur complement: $$\begin{array}{cc}& X_k^{\mathrm{\ast }}=B_k^{-\mathrm{1}}{d}_k-B_k^{-\mathrm{1}}{A}_k^T{\left(A_k{B}_k^{-\mathrm{1}}{A}_k^T\right)}^{-\mathrm{1}}\left[A_k{B}_k^{-\mathrm{1}}{d}_k-b_k\right].\end{array}$$ (25)

From (16), according to the matrix inversion lemma, it has $$\begin{array}{l}{B}_k^{-\mathrm{1}}{d}_k={\left[{\dot{H}}_k^T{R}_k^{-\mathrm{1}}{\dot{H}}_k+P_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}^{-\mathrm{1}}\right]}^{-\mathrm{1}}\\ \times \left[P_{k|k-\mathrm{1}}^{-\mathrm{1}}{\widehat{X}}_{k|k-\mathrm{1}}\right.\\ \left.\hspace{0.5em}+{\dot{H}}_k^T{R}_k^{-\mathrm{1}}\left(Z_k+{\dot{H}}_k{\widehat{X}}_{k|k-\mathrm{1}}-H_k\left({\widehat{X}}_{k|k-\mathrm{1}}\right)\right)\right]\\ ={\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}+P_{k|k-\mathrm{1}}{\dot{H}}_k^T{\left({\dot{H}}_k{P}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}{\dot{H}}_k^T+R_k\right)}^{-\mathrm{1}}\\ \times \left[Z_k+{\dot{H}}_k{\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}-H_k\left({\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}\right)-{\dot{H}}_k{\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}\right]\\ ={\widehat{X}}_{k|k-\mathrm{1}}+P_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}{\dot{H}}_k^T{\left({\dot{H}}_k{P}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}{\dot{H}}_k^T+R_k\right)}^{-\mathrm{1}}\\ \hspace{1em}\times \left[Z_k-H_k\left({\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}\right)\right]\\ \triangleq {\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k}.\end{array}$$ (26)

Theorem 4.

Given the dynamic model (1) and directional measurement models (7) and (4), the Kalman filter estimate is $$\begin{array}{cc}& X_k^{\mathrm{\ast }}={\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k}-B_k^{-\mathrm{1}}{A}_k^T{\left(A_k{B}_k^{-\mathrm{1}}{A}_k^T\right)}^{-\mathrm{1}}\left[A_k{\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k}-b_k\right].\end{array}$$ (27)

Proof.

Substituting (26) into (25), one can easily conclude the theorem.

Remark 5.

We can find that ${\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k}$ is just the state estimate of the unconstrained traditional extended Kalman filter for the system (1) and (4) [30]. The constrained state estimation $X_k^{\mathrm{\ast }}$ depends on the state estimate ${\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k}$ which only considers the noisy distance measurements $Z_k$ , and the constrained boundaries recorrect ${\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k}$ to restrict the state estimate within the boundaries.

4.2. Statistical Properties of CEKF

Lemma 6.

Given the dynamic model (1) and directional measurement models (7) and (4), the Kalman filter estimate $X_k^{\mathrm{\ast }}$ as given by (27) is unbiased.

Proof.

Consider $$\begin{array}{l}{X}_k-X_k^{\mathrm{\ast }}=X_k-{\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k}+B_k^{-\mathrm{1}}{A}_k^T{\left(A_k{B}_k^{-\mathrm{1}}{A}_k^T\right)}^{-\mathrm{1}}\\ \times \left[A_k{\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k}-b_k\right]\\ =\left[I-B_k^{-\mathrm{1}}{A}_k^T{\left(A_k{B}_k^{-\mathrm{1}}{A}_k^T\right)}^{-\mathrm{1}}{A}_k\right]\left(X_k-{\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k}\right),\end{array}$$ (28)where I is a unit matrix; we take the expectation of both sides of the above equation to obtain $$\begin{array}{l}E\left\{X_k-X_k^{\mathrm{\ast }}\right\}\\ \hspace{1em}\hspace{1em}\hspace{0.17em}=\left[I-B_k^{-\mathrm{1}}{A}_k^T{\left(A_k{B}_k^{-\mathrm{1}}{A}_k^T\right)}^{-\mathrm{1}}{A}_k\right]E\left\{\left(X_k-{\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k}\right)\right\}\\ \hspace{1em}\hspace{1em}\hspace{0.17em}=\mathrm{0}.\end{array}$$ (29)The right side of the above equation is zero; this is because $E\{X_k\}=E\{{\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k}\}$ in the unconstrained traditional extended Kalman filter. Therefore, the constrained extended Kalman filter estimate $X_k^{\mathrm{\ast }}$ is an unbiased estimate of $X_k$ .

Theorem 7.

Given the dynamic model (1) and directional measurement models (7) and (4), the Kalman filter estimate $X_k^{\mathrm{\ast }}$ as given by (27) has a smaller error covariance than the estimate ${\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k}$ of the unconstrained traditional extended Kalman filter which only considers the noisy sensor measurements $Z_k$ .

Proof.

From (28), the covariance of the estimate error of $X_k^{\mathrm{\ast }}$ is given by $$\begin{array}{l}E\left\{\left(X_k-X_k^{\mathrm{\ast }}\right){\left(X_k-X_k^{\mathrm{\ast }}\right)}^T\right\}\\ \hspace{1em}=\left[I-B_k^{-\mathrm{1}}{A}_k^T{\left(A_k{B}_k^{-\mathrm{1}}{A}_k^T\right)}^{-\mathrm{1}}{A}_k\right]\\ \hspace{1em}\hspace{1em}·E\left\{\left(X_k-{\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k}\right){\left(X_k-{\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k}\right)}^T\right\}\\ \hspace{1em}\hspace{1em}\times {\left[I-B_k^{-\mathrm{1}}{A}_k^T{\left(A_k{B}_k^{-\mathrm{1}}{A}_k^T\right)}^{-\mathrm{1}}{A}_k\right]}^T.\end{array}$$ (30)

From (26), using the matrix inversion lemma, it can obtain $$\begin{array}{l}E\left\{\left(X_k-{\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k}\right){\left(X_k-{\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k}\right)}^T\right\}\\ =P_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}-P_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}{\dot{H}}_k^T{\left({\dot{H}}_k{P}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}{\dot{H}}_k^T+R_k\right)}^{-\mathrm{1}}{\dot{H}}_k{P}_{k|k-\mathrm{1}}\\ ={\left[P_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}^{-\mathrm{1}}+{\dot{H}}_k^T{R}_k^{-\mathrm{1}}{\dot{H}}_k\right]}^{-\mathrm{1}}\\ =B_k^{-\mathrm{1}}\\ \triangleq P_{k\hspace{0.17em}|\hspace{0.17em}k}.\end{array}$$ (31)Substituting (31) into (30) yields $$\begin{array}{l}E\left\{\left(X_k-X_k^{\mathrm{\ast }}\right){\left(X_k-X_k^{\mathrm{\ast }}\right)}^T\right\}\\ \hspace{1em}=B_k^{-\mathrm{1}}-B_k^{-\mathrm{1}}{\left[B_k^{-\mathrm{1}}{A}_k^T{\left(A_k{B}_k^{-\mathrm{1}}{A}_k^T\right)}^{-\mathrm{1}}{A}_k\right]}^T\\ \hspace{1em}\hspace{1em}-B_k^{-\mathrm{1}}{A}_k^T{\left(A_k{B}_k^{-\mathrm{1}}{A}_k^T\right)}^{-\mathrm{1}}{A}_k{B}_k^{-\mathrm{1}}\\ \hspace{1em}\hspace{1em}+\left[B_k^{-\mathrm{1}}{A}_k^T{\left(A_k{B}_k^{-\mathrm{1}}{A}_k^T\right)}^{-\mathrm{1}}{A}_k\right]\\ \hspace{1em}\hspace{1em}\times B_k^{-\mathrm{1}}{\left[B_k^{-\mathrm{1}}{A}_k^T{\left(A_k{B}_k^{-\mathrm{1}}{A}_k^T\right)}^{-\mathrm{1}}{A}_k\right]}^T\\ \hspace{1em}=B_k^{-\mathrm{1}}-B_k^{-\mathrm{1}}{A}_k^T{\left(A_k{B}_k^{-\mathrm{1}}{A}_k^T\right)}^{-\mathrm{1}}{A}_k{B}_k^{-\mathrm{1}}\\ \hspace{1em}\triangleq P_{k\hspace{0.17em}|\hspace{0.17em}k}^{\mathrm{\ast }}.\end{array}$$ (32)Since $A_k$ is full rank, $B_k^{-\mathrm{1}}{A}_k^T(A_k{B}_k^{-\mathrm{1}}{A}_k^T{)}^{-\mathrm{1}}{A}_k{B}_k^{-\mathrm{1}}$ is positive definite [29]; we can conclude that $$\begin{array}{cc}& P_{k\hspace{0.17em}|\hspace{0.17em}k}^{\mathrm{\ast }}<B_k^{-\mathrm{1}}=P_{k\hspace{0.17em}|\hspace{0.17em}k}.\end{array}$$ (33)

Theorem 8.

Given the dynamic model (1) and directional measurement models (7) and (4), the Kalman filter estimate $X_k^{\mathrm{\ast }}$ as given by (27) is the linear minimum variance estimator.

Proof.

From (27) and (26), it is obvious that $X_k^{\mathrm{\ast }}$ is a linear estimator that is linear in the measurement $Z_k$ . Suppose a linear estimate of $X_k$ for any pair of $a\in {\mathbb{R}}^{\mathrm{4}}$ and $C\in {\mathbb{R}}^{\mathrm{4}\times {\mathcal{l}}_k}$ , denoted by ${\overline{X}}_k=a+C{Z}_k$ ; let $\rho =a-X_k^{\mathrm{\ast }}+C{H}_k({\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}})$ , and then the error covariance of ${\overline{X}}_k$ is $$\begin{array}{l}E\left\{\left(X_k-{\overline{X}}_k\right){\left(X_k-{\overline{X}}_k\right)}^T\right\}\\ =E\left\{\left(X_k-a-C{Z}_k\right){\left(X_k-a-C{Z}_k\right)}^T\right\}\\ =E\left\{\left[X_k-X_k^{\mathrm{\ast }}-\rho -C\left(Z_k-H_k\left({\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}\right)\right)\right]\phantom{{\left.X_k-X_k^{\mathrm{\ast }}-\rho -C\left(Z_k-H_k\left({\widehat{X}}_{k|k-\mathrm{1}}\right)\right)\right]}^T}\right.\\ \hspace{1em}\hspace{1em}\hspace{1em}\left.\times {\left[X_k-X_k^{\mathrm{\ast }}-\rho -C\left(Z_k-H_k\left({\widehat{X}}_{k|k-\mathrm{1}}\right)\right)\right]}^T\right\}\\ =E\left\{\left(X_k-X_k^{\mathrm{\ast }}\right){\left(X_k-X_k^{\mathrm{\ast }}\right)}^T\right\}+\rho {\rho }^T\\ \hspace{1em}+CE\left\{{\left[Z_k-H_k\left({\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}\right)\right]\left[Z_k-H_k\left({\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}\right)\right]}^T\right\}{C}^T\\ \hspace{1em}-E\left\{{\left[X_k-X_k^{\mathrm{\ast }}\right]\left[Z_k-H_k\left({\widehat{X}}_{k|k-\mathrm{1}}\right)\right]}^T\right\}{C}^T\\ \hspace{1em}-CE\left\{\left[Z_k-H_k\left({\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}\right)\right]{\left[X_k-X_k^{\mathrm{\ast }}\right]}^T\right\}.\end{array}$$ (34)

From (12) and (26), define $K_k=P_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}{\dot{H}}_k^T({\dot{H}}_k{P}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}{\dot{H}}_k^T+R_k{)}^{-\mathrm{1}}$ ; one can get $$\begin{array}{cc}& Z_k-H_k\left({\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}\right)={\dot{H}}_k\left(X_k-{\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}\right)+V_k\\ X_k-{\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k}=\left(I-K_k{\dot{H}}_k\right)\left(X_k-{\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}\right)-K_k{V}_k\\ K_k{R}_k=\left[I-K_k{\dot{H}}_k\right]P_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}{\dot{H}}_k^T.\end{array}$$ (35)

Combining (35) with (28), using Lemma 6, we obtain $$\begin{array}{l}E\left\{\left[X_k-X_k^{\mathrm{\ast }}\right]{\left[Z_k-H_k\left({\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}\right)\right]}^T\right\}\\ =E\left\{\left[I-B_k^{-\mathrm{1}}{A}_k^T{\left(A_k{B}_k^{-\mathrm{1}}{A}_k^T\right)}^{-\mathrm{1}}{A}_k\right]\right.\\ \hspace{1em}\hspace{1em}\hspace{1em}\left.\times \left(X_k-{\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k}\right){\left[{\dot{H}}_k\left(X_k-{\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}\right)+V_k\right]}^T\right\}\\ =\left[I-B_k^{-\mathrm{1}}{A}_k^T{\left(A_k{B}_k^{-\mathrm{1}}{A}_k^T\right)}^{-\mathrm{1}}{A}_k\right]\\ \hspace{1em}\times \left[\left(I-K_k{\dot{H}}_k\right)P_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}{\dot{H}}_k^T-K_k{R}_k\right]=\mathrm{0}.\end{array}$$ (36)

Equation (34) can be transformed into $$\begin{array}{l}E\left\{\left(X_k-{\overline{X}}_k\right){\left(X_k-{\overline{X}}_k\right)}^T\right\}\\ \hspace{1em}\hspace{1em}=P_{k\hspace{0.17em}|\hspace{0.17em}k}^{\mathrm{\ast }}+\rho {\rho }^T+C\left({\dot{H}}_k{P}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}{\dot{H}}_k^T+R_k\right)C^T.\end{array}$$ (37)

Since $\rho {\rho }^T\ge \mathrm{0}$ and $C({\dot{H}}_k{P}_{k\hspace{0.17em}|\hspace{0.17em}k-\mathrm{1}}{\dot{H}}_k^T+R_k)C^T\ge \mathrm{0}$ , one can get $$\begin{array}{cc}& E\left\{\left(X_k-{\overline{X}}_k\right){\left(X_k-{\overline{X}}_k\right)}^T\right\}\ge P_{k\hspace{0.17em}|\hspace{0.17em}k}^{\mathrm{\ast }}.\end{array}$$ (38)Note that if and only if $\rho =\mathrm{0}$ and $C=\mathrm{0}$ , that is, ${\overline{X}}_k=X_k^{\mathrm{\ast }}$ , the error covariance is the minimum, equal to $P_{k\hspace{0.17em}|\hspace{0.17em}k}^{\mathrm{\ast }}$ . This states that $X_k^{\mathrm{\ast }}$ is the linear minimum variance estimator.

Remark 9.

From Theorems 7 and 8, one can see that, under the boundary constrained measurement model, the traditional extended Kalman filter cannot provide the minimum variance estimate, and our proposed constrained extended Kalman filter is the linear unbiased minimum variance estimator. The only time that there are no active boundary constraints at the solution of (18) is when ${\widehat{X}}_{k|k-\mathrm{1}}\in {\Psi }_k$ ; that is, ${\lambda }_k^{\mathrm{\ast }}=\mathrm{0}$ , ${\mu }_k^{\mathrm{\ast }}=\mathrm{0}$ , and then $A_k$ and $b_k$ are empty sets; they are set to $\mathrm{0}$ . In this case, it follows from (27) and (32) that $X_k^{\mathrm{\ast }}={\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k}$ and $P_{k\hspace{0.17em}|\hspace{0.17em}k}^{\mathrm{\ast }}=P_{k\hspace{0.17em}|\hspace{0.17em}k}$ ; the directional constraint ${\mathcal{B}}_k$ has no effect on the state update; as such the proposed constrained extended Kalman filter simplifies to the traditional unconstrained extended Kalman filter.

5. Numerical Simulations

In this section, we use Monte Carlo simulations to analyse the efficacy of the proposed constrained extended Kalman filter (CEKF) target tracking algorithm in DSNs, compared with the traditional EKF (TEKF) as well as traditional UKF (TUKF) based on the measurement model (4) and the geometric centroid method that adopts the centroid of the directional constrained area ${\mathcal{B}}_k$ in (7) as the location estimation of the target.

The experiment scenario is a square sensing field of the size $\mathrm{8}\times \mathrm{8}$  m²; each sensor has $M=\mathrm{4}$ sensing directions with a detection radius $r=\mathrm{3}$  m. A moving target travels within the sensing field at a constant angular velocity of 0.244 rad/s along a circle centered at (12 m, 12 m) with a radius of 0.3 meters. The observation interval is 0.2 seconds (5 Hz). The process noise $W_{k-\mathrm{1}}$ can be approximated by a variable acceleration with $Q_{k-\mathrm{1}}=\text{Diag}\{\mathrm{0.1138,0.1138}\}$ . Consider four directional sensor network scenarios with two different node densities as well as two noise conditions: (i) high density and high noise condition with $N=\mathrm{50}$ and ${\sigma }_n^{\mathrm{2}}=\mathrm{5}$ , (ii) high density and low noise condition with $N=\mathrm{50}$ and ${\sigma }_n^{\mathrm{2}}=\mathrm{0.05}$ , (iii) low density and high noise condition with $N=\mathrm{10}$ and ${\sigma }_n^{\mathrm{2}}=\mathrm{5}$ , and (iv) low density and low noise condition with $N=\mathrm{10}$ and ${\sigma }_n^{\mathrm{2}}=\mathrm{0.05}$ .

The initial state estimate and the corresponding covariance matrix for these tracking algorithms are chosen to be $$\begin{array}{cc}& {\widehat{X}}_{\mathrm{0}\hspace{0.17em}|\hspace{0.17em}\mathrm{0}}={\left[\begin{array}{cccc}\mathrm{12.0}& \mathrm{0.7320}& \mathrm{9.0}& \mathrm{0.0}\end{array}\right]}^T,\\ P_{\mathrm{0}\hspace{0.17em}|\hspace{0.17em}\mathrm{0}}=\mathrm{0.01}\times \text{Diag}\left\{\begin{array}{cccc}1& 1& 1& 1\end{array}\right\}.\end{array}$$ (39)

100 Monte Carlo runs are carried out and the results are summarized in the following. Figure 3 shows the tracking trajectories of the four tracking algorithms under different simulation scenarios. Figure 4 describes the estimated location biases $‖(x_k,y_k)-({\widehat{x}}_{k\hspace{0.17em}|\hspace{0.17em}k},{\widehat{y}}_{k\hspace{0.17em}|\hspace{0.17em}k})‖$ of target at each step for different tracking algorithms under four simulation scenarios, where $({\widehat{x}}_{k\hspace{0.17em}|\hspace{0.17em}k},{\widehat{y}}_{k\hspace{0.17em}|\hspace{0.17em}k})$ comes from the estimated value ${\widehat{X}}_{k\hspace{0.17em}|\hspace{0.17em}k}$ . Figure 5 represents the mean squared errors (MSEs) of the state estimate $E\{{‖X_k-{\widehat{X}}_k‖}^{\mathrm{2}}\}$ at each step for these three Kalman filters (CEKF, TEKF, and TUKF).

Figure 3

Tracking trajectories under different scenarios.

Figure 4

The estimated location biases of target under different scenarios.

Figure 5

The mean squared errors of the Kalman filter estimate under different scenarios.

In the low noise condition, the impact of the noise has negligible magnitude; then the distance measurements are very accurate, all the Kalman filtering algorithms achieve similar performances, and the estimated location biases are smaller than those of the CENT method. Comparing the two subpictures (b) and (d) in Figures 4 and 5, we can find that because the number of nodes in the high density condition is more than that of the low density condition, the observations' information with respect to the target including distance measurements and sensing directions in the high density condition is more than that of the low density condition; then the four tracking algorithms under the high density condition can obtain better tracking accuracy than the low density condition.

However, in the high noise condition, the distance measurement errors become bigger; the location estimate of TEKF and TUKF cannot guarantee better estimation accuracy; the tracking trajectories in TEKF and TUKF greatly diverge from the true target trajectory, even exceeding the boundaries in which the target locates; then the CENT method can perform better than the traditional Kalman filters. In contrast, CEKF integrates both distance measurements and sensing directions information; it effectively restrains the target estimated position within the boundary; as such CEKF yields the best tracking performance. Comparing the two subpictures (a) and (c) in Figure 4, we can further see that, in the high density condition, the directional constrained area of the target would be much smaller than that of the low density condition; hence CENT also performs better in the high density condition.

As shown in Figure 5, the MSE in the CEKF is the lowest among the Kalman filters under all the simulation conditions; it evidently verifies the validity of Theorems 7 and 8. In addition, the averaged root-mean square of the target estimated position errors (RMSEs) $\sqrt{E\{{‖(x_k,y_k)-({\widehat{x}}_{k\hspace{0.17em}|\hspace{0.17em}k},{\widehat{y}}_{k\hspace{0.17em}|\hspace{0.17em}k})‖}^{\mathrm{2}}\}}$ for the four tracking algorithms is summarised in Table 1. Results are also consistent with those in Figures 4 and 5. In Table 1, CEKF can acquire the smallest averaged RMSE under the four scenarios due to merging both distance measurements and sensing directions information about the target. The averaged RMSEs in TEKF and TUKF are lower than those of CENT under the low noise condition; this is because the target position estimated from the precise distance measurements in TEKF and TUKF will almost be in the constrained area. In addition, it shows that TUKF slightly outweighs TEKF for highly manoeuvred targets. The above analysis results indicate that the proposed constrained extended Kalman filter target tracking algorithm promises superior tracking performance for directional wireless networks.

Table 1

Averaged RMSEs of the four tracking algorithms under different scenarios.

Algorithms	High density and high noise	High density and low noise	Low density and high noise	Low density and low noise
TEKF	0.4665	0.0790	0.9640	0.1903
TUKF	0.4613	0.0789	0.9192	0.1783
CENT	0.1676	0.1676	0.5925	0.5925
CEKF	0.1373	0.0618	0.3587	0.1472

6. Conclusion

In this paper, we propose a constrained extended Kalman filter for target tracking in directional sensor networks, which effectively integrates both directions and measurements' information of the directional sensors. Its structural and statistical properties are also derived; it is found that, by incorporating the directional measurement model, the proposed constrained extended Kalman filter is the linear unbiased minimum variance estimator and promises smaller tracking error covariance compared with the unconstrained traditional extended Kalman filter only using sensor measurements' information. In the future work, we will extend our study of target tracking in DSNs from the ideal disk sensing range detection model to realistic irregular sensing range detection model.