# Constrained Extended Kalman Filter for Target Tracking in Directional Sensor Networks

Authors: Sha Wen; Zixing Cai; Xiaoqing Hu
Journal: International Journal of Distributed Sensor Networks
 (2015)

---

## 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.


---

## Body

## 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 Section2, 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]:(1)Xk=Fk-1Xk-1+Gk-1Wk-1,where the state vector Xk=[xk,x˙k,yk,y˙k]T specifies the position (xk,yk) and velocity (x˙k,y˙k) of the target along x-direction and y-direction at kth time step, respectively. Fk-1 and Gk-1 are the state transition matrix and the input matrix, respectively. For tracking applications, one has (2)Fk-1=1Δk-1000100001Δk-10001,Gk-1=Δk-1220Δk-100Δk-1220Δk-1.In the above, Δk-1 is the sampling time interval between two successive time steps k and k-1. The input noise Wk-1=[wx,wy]T is a 2×1 Gaussian random vector with zero mean and a covariance matrix Qk-1; wx and wy correspond to the process noises (random accelerations) along the x-axis and y-axis, respectively. Generally, it is assumed that wx is uncorrelated with wy; then Qk-1=Diag{σx2,σy2}.


### 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 ofN identical nodes in a two-dimensional plane; the position of the nth (1≤n≤N) node is denoted by (xn,yn). 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.
(a)
M
=
4
(b)
M
=
3When the node-to-target distancehn(k)=(xk-xn)2+(yk-yn)2 is less than r, the target will be detected and one of its own sensors is triggered; the distance hn(k) will be estimated. We shall denote Ωk=n∣hn(k)≤r,1≤n≤N to be the set of indices of nodes within the sensing range of the target at time k; the measurement zn(k) of the nth node at time k can be expressed as [21, 22](3)znk=hnk+vnkn∈Ωk,where vn(k) is the observation noise at the nth node, which is assumed to be an independent Gaussian random variable with zero mean and variance σn2 for each node and uncorrelated with Wk-1 and Xk.Denotelk to be the cardinal number of Ωk (i.e., the number of elements in Ωk); the sensor measurements at the kth time step may be represented in a matrix form:(4)Zk=HkXk+Vk,where Hk(Xk)=[h1(k),…,hlk(k)]T and Vk=[v1(k),…,vlk(k)]T. Moreover, the covariance matrix of Vk is Rk=Diag{σn2;1≤n≤lk}.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 Figure2, we illustrate two target detection scenarios in DSNs with M=4 and M=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 nth node is represented as(5)An,i≜x,y∣lxn,i≤x≤uxn,i,lyn,i≤y≤uyn,i,where lxn,i,uxn,i,lyn,i,uyn,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(6)Dk≜n,i∣xk,yk∈An,i,n∈Ωk,i=1,…,M.Figure 2
Target detection scenario in a DSN.
(a)
M
=
4
(b)
M
=
3By intersection operations, the constrained area of the target at thekth time step can be described in an approximate rectangular set:(7)Bk=⋂Dk‍An,i=x,y∣xmin⁡k≤x≤xmax⁡k,ymin⁡k≤y≤ymax⁡k,where(8)xmin⁡k=Max⁡n,i∈Dklxn,i,xmax⁡k=Min⁡n,i∈Dkuxn,i,ymin⁡k=Max⁡n,i∈Dklyn,i,ymax⁡k=Min⁡n,i∈Dkuyn,i.In above equation, Max⁡{·} and 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 {(xn,yn)}, and the target state Xk.
Ensure: The distance measurements Zk and directional constrained area Bk
(1) Obtain the set of indices of active nodesΩk,
(2) Acquire the detecting distance measurementsZk using (4),
(3) Approximate the sensing directionAn,i by a rectangle using (5),
(4) Define the indices set of triggered sensorsDk using (6),
(5) Construct the constrained areaBk 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 Z1:k={Z1,Z2,…,Zk} from initial time up to and including time k and directional constrained area Bk, our goal is to obtain a sequential Bayesian estimate of Xk using a KF formulation. This can be formulated as a maximum a posteriori probability (MAP) estimation problem with a constraint on the state vector Xk:(9)Max⁡pXk∣Z1:ksubjecttoxk,yk∈Bk,where p(Xk∣Z1:k) is the posterior distribution of Xk given the observed value Z1:k.


## 2.1. Target Motion Model

The motion of a moving target can be described by a 2D state space equation [19, 20]:(1)Xk=Fk-1Xk-1+Gk-1Wk-1,where the state vector Xk=[xk,x˙k,yk,y˙k]T specifies the position (xk,yk) and velocity (x˙k,y˙k) of the target along x-direction and y-direction at kth time step, respectively. Fk-1 and Gk-1 are the state transition matrix and the input matrix, respectively. For tracking applications, one has (2)Fk-1=1Δk-1000100001Δk-10001,Gk-1=Δk-1220Δk-100Δk-1220Δk-1.In the above, Δk-1 is the sampling time interval between two successive time steps k and k-1. The input noise Wk-1=[wx,wy]T is a 2×1 Gaussian random vector with zero mean and a covariance matrix Qk-1; wx and wy correspond to the process noises (random accelerations) along the x-axis and y-axis, respectively. Generally, it is assumed that wx is uncorrelated with wy; then Qk-1=Diag{σx2,σy2}.


## 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 ofN identical nodes in a two-dimensional plane; the position of the nth (1≤n≤N) node is denoted by (xn,yn). 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.
(a)
M
=
4
(b)
M
=
3When the node-to-target distancehn(k)=(xk-xn)2+(yk-yn)2 is less than r, the target will be detected and one of its own sensors is triggered; the distance hn(k) will be estimated. We shall denote Ωk=n∣hn(k)≤r,1≤n≤N to be the set of indices of nodes within the sensing range of the target at time k; the measurement zn(k) of the nth node at time k can be expressed as [21, 22](3)znk=hnk+vnkn∈Ωk,where vn(k) is the observation noise at the nth node, which is assumed to be an independent Gaussian random variable with zero mean and variance σn2 for each node and uncorrelated with Wk-1 and Xk.Denotelk to be the cardinal number of Ωk (i.e., the number of elements in Ωk); the sensor measurements at the kth time step may be represented in a matrix form:(4)Zk=HkXk+Vk,where Hk(Xk)=[h1(k),…,hlk(k)]T and Vk=[v1(k),…,vlk(k)]T. Moreover, the covariance matrix of Vk is Rk=Diag{σn2;1≤n≤lk}.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 Figure2, we illustrate two target detection scenarios in DSNs with M=4 and M=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 nth node is represented as(5)An,i≜x,y∣lxn,i≤x≤uxn,i,lyn,i≤y≤uyn,i,where lxn,i,uxn,i,lyn,i,uyn,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(6)Dk≜n,i∣xk,yk∈An,i,n∈Ωk,i=1,…,M.Figure 2
Target detection scenario in a DSN.
(a)
M
=
4
(b)
M
=
3By intersection operations, the constrained area of the target at thekth time step can be described in an approximate rectangular set:(7)Bk=⋂Dk‍An,i=x,y∣xmin⁡k≤x≤xmax⁡k,ymin⁡k≤y≤ymax⁡k,where(8)xmin⁡k=Max⁡n,i∈Dklxn,i,xmax⁡k=Min⁡n,i∈Dkuxn,i,ymin⁡k=Max⁡n,i∈Dklyn,i,ymax⁡k=Min⁡n,i∈Dkuyn,i.In above equation, Max⁡{·} and 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 {(xn,yn)}, and the target state Xk.
Ensure: The distance measurements Zk and directional constrained area Bk
(1) Obtain the set of indices of active nodesΩk,
(2) Acquire the detecting distance measurementsZk using (4),
(3) Approximate the sensing directionAn,i by a rectangle using (5),
(4) Define the indices set of triggered sensorsDk using (6),
(5) Construct the constrained areaBk 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 Z1:k={Z1,Z2,…,Zk} from initial time up to and including time k and directional constrained area Bk, our goal is to obtain a sequential Bayesian estimate of Xk using a KF formulation. This can be formulated as a maximum a posteriori probability (MAP) estimation problem with a constraint on the state vector Xk:(9)Max⁡pXk∣Z1:ksubjecttoxk,yk∈Bk,where p(Xk∣Z1:k) is the posterior distribution of Xk given the observed value Z1: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 distributionp(Xk-1∣Z1:k-1)=N(X^k-1∣k-1,Pk-1∣k-1), the joint distribution of Xk and Xk-1 conditioned on the previous measurements Z1:k-1={Z1,Z2,…,Zk-1} is(10)pXk-1,Xk∣Z1:k-1=pXk∣Xk-1pXk-1∣Z1:k-1=NX^k-1k-1Fk-1X^k-1k-1,Pk-1∣k-1Pk-1∣k-1Fk-1TFk-1Pk-1∣k-1Fk-1Pk-1∣k-1Fk-1T+Gk-1Qk-1Gk-1T.The predictive distribution of Xk given the measurement history up to time step k-1 can be calculated by Lemma A.1 in Appendix. Consider (11)pXk∣Z1:k-1=NX^k∣k-1,Pkk-1=NFk-1X^k-1∣k-1,Fk-1Pk-1∣k-1Fk-1T+Gk-1Qk-1Gk-1T.To update the state estimate using latest observationZk, consider a first order Taylor series expansion of Hk(Xk) around X^k∣k-1; thus(12)Zk≈HkX^k∣k-1+H˙kXk-X^k∣k-1+Vk,where H˙k=∂Hk(Xk)/∂XkTXk=X^kk-1 is a lk×4 matrix.According to (12), one has the parameters of the posterior distribution by Bayes rule:(13)pXk∣Z1:k∝pZk∣XkpXkZ1:k-1,where(14)pZk∣Xk=12πlk/2det⁡⁡Rk1/2·exp⁡⁡-12Zk-H˙kXk-HkX^k|k-1+H˙kX^k|k-1T×Rk-1Zk-H˙kXk-HkX^k∣k-1+H˙kX^k∣k-112,pXk∣Z1:k-1=12π2det⁡⁡Pk∣k-11/2·exp⁡-12Xk-X^k∣k-1TPk∣k-1-1Xk-X^k∣k-1.The MAP estimate associated with the posteriorp(Xk∣Z1:k) in (13) is the vector Xk that maximizes the likelihood p(Xk∣Z1:k), which is identical to minimizing its negative logarithm:(15)LXk=Zk-H˙kXk-HkX^k∣k-1+H˙kX^k∣k-1T×Rk-1Zk-H˙kXk-HkX^k∣k-1+H˙kX^k∣k-1+Xk-X^k∣k-1TPk∣k-1-1Xk-X^k∣k-1+c1=XkTH˙kTRk-1H˙k+Pk∣k-1-1Xk-2Zk+H˙kX^k∣k-1-HkX^k∣k-1TZk+H˙kX^k∣k-1-HkX^k∣k-1T×Rk-1H˙k+X^k∣k-1TPk∣k-1-1Xk+c2,where c1 and c2 are constant terms not dependent on Xk.Let(16)Bk=H˙kTRk-1H˙k+Pk|k-1-1dk=Pk∣k-1-1X^k∣k-1+H˙kTRk-1Zk+H˙kX^k∣k-1-HkX^k∣k-1.Equation (15) can be simplified into(17)LXk=XkTBkXk-2dkTXk+c2.The Kalman filter estimate for target tracking in DSNs can be derived by finding anXk such that the logarithm function is minimized and Xk locates in the directional constrained area Bk. So the MAP estimation problem can be expressed as(18)Min⁡XkTBkXk-2dkTXksubjecttoxk,yk∈Bk.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 Bk.Theorem 1.
The estimation problem defined in (18) is convex.Proof.
SinceBk is a convex set, noting that Bk is a positive definite matrix, then the objective function XkTBkXk-2dkTXk 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 Xk∗, λk∗, and μk∗ require that(19)BkXk∗-dk-λk∗+μk∗=0λk∗≥0,μk∗≥0λk∗Tlk-Xk∗=0μk∗TXk∗-uk=0lk≤Xk∗≤uk,where lk=[xmin⁡(k),-∞,ymin⁡(k),-∞]T and uk=[xmax⁡(k),∞,ymax⁡(k),∞]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=0, λk(0)=0, μk(0)=0 and Xk(0)=-Bk-1dk. If Xk(0) is feasible, then the solution
is:Xk*=Xk(0), λk*=λk(0), μk*=μk(0).
Looping: If Xk(0) is not feasible, the components of Xk(0) that violate the bound constraints are reset to the midpoint of the box
defined by the bounds, and then uses Quadprog method
[Xk(t),λk(t),μk(t)]=Quadprog{Bk,dk,lk,uk,Xk(0)}
Stopping and output: if L(Xk(t-1))-L(Xk(t)) is smaller than the specified tolerance, then the solution
is:Xk*=Xk(t), λk*=λk(t), μk*=μk(t).Remark 2.
Nonzero elements of the Lagrange multipliersλk∗ and μk∗ indicate active constraints at the solution Xk∗; 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, ifA 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 Xk∗ and λk∗, μk∗ satisfy the KKT conditions in (19).Proof.
SinceBk>0 and L(Xk) is convex, there exists only an optimal point for minimizing L(Xk) with unconstrained conditions; that is, there is a unique solution Xk⋆ such that(20)LXk⋆=Min⁡Xk∈R4LXk=Min⁡Xk∈Ψk⋃‍R4-ΨkLXk,=Min⁡Min⁡Xk∈ΨkLXk,Min⁡Xk∈R4-ΨkLXk,where Ψk={Xk∣lk≤Xk≤uk,Xk∈R4}. The above equation indicates that there must be a solution Xk∗ such that L(Xk) is minimal for Xk∈Ψk; it means that the problem (18) is strictly feasible. By the Slaters constraint qualification [29], there exist dual optimal λk∗, μk∗, which together with Xk∗ satisfy the KKT conditions in (19).Letbk be the column vector composed of active constrained boundaries whose Lagrange multiplier elements are nonzero, and Λk is the column vector that contains the indices of elements of λk∗ and μk∗ corresponding to the active constrained boundaries in bk; then, one has(21)Aki,j=0ifbki,1∈bk,j∉Λk1ifbki,1∈bk,j∈Λkforj=1,2,…,4.It is noted thatAk has the same number of rows with bk. From Remark 2, these active boundary conditions may be expressed in terms of the constraints that are satisfied with equality at Xk∗; the solution of (18) would be a solution of the equivalent equality-constrained problem:(22)Min⁡XkTBkXk-2dkTXksubjecttoAkXk=bk.Here, the KKT necessary conditions for a minimum are(23)BkXk∗-dk+AkTνk∗=0,AkXk∗=bk,which we can write as(24)BkAkTAk0Xk∗νk∗=dkbk.Since each component of the minimumXk∗ cannot be at the lower and upper bound simultaneously, the matrix Ak is full row rank; as such one can get the solution by the Schur complement:(25)Xk∗=Bk-1dk-Bk-1AkTAkBk-1AkT-1AkBk-1dk-bk.From (16), according to the matrix inversion lemma, it has(26)Bk-1dk=H˙kTRk-1H˙k+Pk∣k-1-1-1×Pk|k-1-1X^k|k-1+H˙kTRk-1Zk+H˙kX^k|k-1-HkX^k|k-1=X^k∣k-1+Pk|k-1H˙kTH˙kPk∣k-1H˙kT+Rk-1×Zk+H˙kX^k∣k-1-HkX^k∣k-1-H˙kX^k∣k-1=X^k|k-1+Pk∣k-1H˙kTH˙kPk∣k-1H˙kT+Rk-1×Zk-HkX^k∣k-1≜X^k∣k.Theorem 4.
Given the dynamic model (1) and directional measurement models (7) and (4), the Kalman filter estimate is(27)Xk∗=X^k∣k-Bk-1AkTAkBk-1AkT-1AkX^k∣k-bk.Proof.
Substituting (26) into (25), one can easily conclude the theorem.Remark 5.
We can find thatX^k∣k is just the state estimate of the unconstrained traditional extended Kalman filter for the system (1) and (4) [30]. The constrained state estimation Xk∗ depends on the state estimate X^k∣k which only considers the noisy distance measurements Zk, and the constrained boundaries recorrect X^k∣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 Xk∗ as given by (27) is unbiased.Proof.
Consider(28)Xk-Xk∗=Xk-X^k∣k+Bk-1AkTAkBk-1AkT-1×AkX^k∣k-bk=I-Bk-1AkTAkBk-1AkT-1AkXk-X^k∣k,where I is a unit matrix; we take the expectation of both sides of the above equation to obtain(29)EXk-Xk∗=I-Bk-1AkTAkBk-1AkT-1AkEXk-X^k∣k=0.The right side of the above equation is zero; this is because E{Xk}=E{X^k∣k} in the unconstrained traditional extended Kalman filter. Therefore, the constrained extended Kalman filter estimate Xk∗ is an unbiased estimate of Xk.Theorem 7.
Given the dynamic model (1) and directional measurement models (7) and (4), the Kalman filter estimate Xk∗ as given by (27) has a smaller error covariance than the estimate X^k∣k of the unconstrained traditional extended Kalman filter which only considers the noisy sensor measurements Zk.Proof.
From (28), the covariance of the estimate error of Xk∗ is given by(30)EXk-Xk∗Xk-Xk∗T=I-Bk-1AkTAkBk-1AkT-1Ak·EXk-X^k∣kXk-X^k∣kT×I-Bk-1AkTAkBk-1AkT-1AkT.
From (26), using the matrix inversion lemma, it can obtain(31)EXk-X^k∣kXk-X^k∣kT=Pk∣k-1-Pk∣k-1H˙kTH˙kPk∣k-1H˙kT+Rk-1H˙kPk|k-1=Pk∣k-1-1+H˙kTRk-1H˙k-1=Bk-1≜Pk∣k.Substituting (31) into (30) yields(32)EXk-Xk∗Xk-Xk∗T=Bk-1-Bk-1Bk-1AkTAkBk-1AkT-1AkT-Bk-1AkTAkBk-1AkT-1AkBk-1+Bk-1AkTAkBk-1AkT-1Ak×Bk-1Bk-1AkTAkBk-1AkT-1AkT=Bk-1-Bk-1AkTAkBk-1AkT-1AkBk-1≜Pk∣k∗.Since Ak is full rank, Bk-1AkT(AkBk-1AkT)-1AkBk-1 is positive definite [29]; we can conclude that(33)Pk∣k∗<Bk-1=Pk∣k.Theorem 8.
Given the dynamic model (1) and directional measurement models (7) and (4), the Kalman filter estimate Xk∗ as given by (27) is the linear minimum variance estimator.Proof.
From (27) and (26), it is obvious that Xk∗ is a linear estimator that is linear in the measurement Zk. Suppose a linear estimate of Xk for any pair of a∈R4 and C∈R4×lk, denoted by X¯k=a+CZk; let ρ=a-Xk∗+CHk(X^k∣k-1), and then the error covariance of X¯k is(34)EXk-X¯kXk-X¯kT=EXk-a-CZkXk-a-CZkT=EXk-Xk∗-ρ-CZk-HkX^k∣k-1Xk-Xk∗-ρ-CZk-HkX^k|k-1T×Xk-Xk∗-ρ-CZk-HkX^k|k-1T=EXk-Xk∗Xk-Xk∗T+ρρT+CEZk-HkX^k∣k-1Zk-HkX^k∣k-1TCT-EXk-Xk∗Zk-HkX^k|k-1TCT-CEZk-HkX^k∣k-1Xk-Xk∗T.
From (12) and (26), define Kk=Pk∣k-1H˙kT(H˙kPk∣k-1H˙kT+Rk)-1; one can get(35)Zk-HkX^k∣k-1=H˙kXk-X^k∣k-1+VkXk-X^k∣k=I-KkH˙kXk-X^k∣k-1-KkVkKkRk=I-KkH˙kPk∣k-1H˙kT.
Combining (35) with (28), using Lemma 6, we obtain(36)EXk-Xk∗Zk-HkX^k∣k-1T=EI-Bk-1AkTAkBk-1AkT-1Ak×Xk-X^k∣kH˙kXk-X^k∣k-1+VkT=I-Bk-1AkTAkBk-1AkT-1Ak×I-KkH˙kPk∣k-1H˙kT-KkRk=0.
Equation (34) can be transformed into(37)EXk-X¯kXk-X¯kT=Pk∣k∗+ρρT+CH˙kPk∣k-1H˙kT+RkCT.
SinceρρT≥0 and C(H˙kPk∣k-1H˙kT+Rk)CT≥0, one can get(38)EXk-X¯kXk-X¯kT≥Pk∣k∗.Note that if and only if ρ=0 and C=0, that is, X¯k=Xk∗, the error covariance is the minimum, equal to Pk∣k∗. This states that Xk∗ is the linear minimum variance estimator.Remark 9.
From Theorems7 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 X^k|k-1∈Ψk; that is, λk∗=0, μk∗=0, and then Ak and bk are empty sets; they are set to 0. In this case, it follows from (27) and (32) that Xk∗=X^k∣k and Pk∣k∗=Pk∣k; the directional constraint Bk has no effect on the state update; as such the proposed constrained extended Kalman filter simplifies to the traditional unconstrained extended Kalman filter.


## 4.1. Structural Properties of CEKF

Lemma 3.
Given the quadratic optimization problem (18), there exists a unique pair of primal and dual optimal Xk∗ and λk∗, μk∗ satisfy the KKT conditions in (19).Proof.
SinceBk>0 and L(Xk) is convex, there exists only an optimal point for minimizing L(Xk) with unconstrained conditions; that is, there is a unique solution Xk⋆ such that(20)LXk⋆=Min⁡Xk∈R4LXk=Min⁡Xk∈Ψk⋃‍R4-ΨkLXk,=Min⁡Min⁡Xk∈ΨkLXk,Min⁡Xk∈R4-ΨkLXk,where Ψk={Xk∣lk≤Xk≤uk,Xk∈R4}. The above equation indicates that there must be a solution Xk∗ such that L(Xk) is minimal for Xk∈Ψk; it means that the problem (18) is strictly feasible. By the Slaters constraint qualification [29], there exist dual optimal λk∗, μk∗, which together with Xk∗ satisfy the KKT conditions in (19).Letbk be the column vector composed of active constrained boundaries whose Lagrange multiplier elements are nonzero, and Λk is the column vector that contains the indices of elements of λk∗ and μk∗ corresponding to the active constrained boundaries in bk; then, one has(21)Aki,j=0ifbki,1∈bk,j∉Λk1ifbki,1∈bk,j∈Λkforj=1,2,…,4.It is noted thatAk has the same number of rows with bk. From Remark 2, these active boundary conditions may be expressed in terms of the constraints that are satisfied with equality at Xk∗; the solution of (18) would be a solution of the equivalent equality-constrained problem:(22)Min⁡XkTBkXk-2dkTXksubjecttoAkXk=bk.Here, the KKT necessary conditions for a minimum are(23)BkXk∗-dk+AkTνk∗=0,AkXk∗=bk,which we can write as(24)BkAkTAk0Xk∗νk∗=dkbk.Since each component of the minimumXk∗ cannot be at the lower and upper bound simultaneously, the matrix Ak is full row rank; as such one can get the solution by the Schur complement:(25)Xk∗=Bk-1dk-Bk-1AkTAkBk-1AkT-1AkBk-1dk-bk.From (16), according to the matrix inversion lemma, it has(26)Bk-1dk=H˙kTRk-1H˙k+Pk∣k-1-1-1×Pk|k-1-1X^k|k-1+H˙kTRk-1Zk+H˙kX^k|k-1-HkX^k|k-1=X^k∣k-1+Pk|k-1H˙kTH˙kPk∣k-1H˙kT+Rk-1×Zk+H˙kX^k∣k-1-HkX^k∣k-1-H˙kX^k∣k-1=X^k|k-1+Pk∣k-1H˙kTH˙kPk∣k-1H˙kT+Rk-1×Zk-HkX^k∣k-1≜X^k∣k.Theorem 4.
Given the dynamic model (1) and directional measurement models (7) and (4), the Kalman filter estimate is(27)Xk∗=X^k∣k-Bk-1AkTAkBk-1AkT-1AkX^k∣k-bk.Proof.
Substituting (26) into (25), one can easily conclude the theorem.Remark 5.
We can find thatX^k∣k is just the state estimate of the unconstrained traditional extended Kalman filter for the system (1) and (4) [30]. The constrained state estimation Xk∗ depends on the state estimate X^k∣k which only considers the noisy distance measurements Zk, and the constrained boundaries recorrect X^k∣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 Xk∗ as given by (27) is unbiased.Proof.
Consider(28)Xk-Xk∗=Xk-X^k∣k+Bk-1AkTAkBk-1AkT-1×AkX^k∣k-bk=I-Bk-1AkTAkBk-1AkT-1AkXk-X^k∣k,where I is a unit matrix; we take the expectation of both sides of the above equation to obtain(29)EXk-Xk∗=I-Bk-1AkTAkBk-1AkT-1AkEXk-X^k∣k=0.The right side of the above equation is zero; this is because E{Xk}=E{X^k∣k} in the unconstrained traditional extended Kalman filter. Therefore, the constrained extended Kalman filter estimate Xk∗ is an unbiased estimate of Xk.Theorem 7.
Given the dynamic model (1) and directional measurement models (7) and (4), the Kalman filter estimate Xk∗ as given by (27) has a smaller error covariance than the estimate X^k∣k of the unconstrained traditional extended Kalman filter which only considers the noisy sensor measurements Zk.Proof.
From (28), the covariance of the estimate error of Xk∗ is given by(30)EXk-Xk∗Xk-Xk∗T=I-Bk-1AkTAkBk-1AkT-1Ak·EXk-X^k∣kXk-X^k∣kT×I-Bk-1AkTAkBk-1AkT-1AkT.
From (26), using the matrix inversion lemma, it can obtain(31)EXk-X^k∣kXk-X^k∣kT=Pk∣k-1-Pk∣k-1H˙kTH˙kPk∣k-1H˙kT+Rk-1H˙kPk|k-1=Pk∣k-1-1+H˙kTRk-1H˙k-1=Bk-1≜Pk∣k.Substituting (31) into (30) yields(32)EXk-Xk∗Xk-Xk∗T=Bk-1-Bk-1Bk-1AkTAkBk-1AkT-1AkT-Bk-1AkTAkBk-1AkT-1AkBk-1+Bk-1AkTAkBk-1AkT-1Ak×Bk-1Bk-1AkTAkBk-1AkT-1AkT=Bk-1-Bk-1AkTAkBk-1AkT-1AkBk-1≜Pk∣k∗.Since Ak is full rank, Bk-1AkT(AkBk-1AkT)-1AkBk-1 is positive definite [29]; we can conclude that(33)Pk∣k∗<Bk-1=Pk∣k.Theorem 8.
Given the dynamic model (1) and directional measurement models (7) and (4), the Kalman filter estimate Xk∗ as given by (27) is the linear minimum variance estimator.Proof.
From (27) and (26), it is obvious that Xk∗ is a linear estimator that is linear in the measurement Zk. Suppose a linear estimate of Xk for any pair of a∈R4 and C∈R4×lk, denoted by X¯k=a+CZk; let ρ=a-Xk∗+CHk(X^k∣k-1), and then the error covariance of X¯k is(34)EXk-X¯kXk-X¯kT=EXk-a-CZkXk-a-CZkT=EXk-Xk∗-ρ-CZk-HkX^k∣k-1Xk-Xk∗-ρ-CZk-HkX^k|k-1T×Xk-Xk∗-ρ-CZk-HkX^k|k-1T=EXk-Xk∗Xk-Xk∗T+ρρT+CEZk-HkX^k∣k-1Zk-HkX^k∣k-1TCT-EXk-Xk∗Zk-HkX^k|k-1TCT-CEZk-HkX^k∣k-1Xk-Xk∗T.
From (12) and (26), define Kk=Pk∣k-1H˙kT(H˙kPk∣k-1H˙kT+Rk)-1; one can get(35)Zk-HkX^k∣k-1=H˙kXk-X^k∣k-1+VkXk-X^k∣k=I-KkH˙kXk-X^k∣k-1-KkVkKkRk=I-KkH˙kPk∣k-1H˙kT.
Combining (35) with (28), using Lemma 6, we obtain(36)EXk-Xk∗Zk-HkX^k∣k-1T=EI-Bk-1AkTAkBk-1AkT-1Ak×Xk-X^k∣kH˙kXk-X^k∣k-1+VkT=I-Bk-1AkTAkBk-1AkT-1Ak×I-KkH˙kPk∣k-1H˙kT-KkRk=0.
Equation (34) can be transformed into(37)EXk-X¯kXk-X¯kT=Pk∣k∗+ρρT+CH˙kPk∣k-1H˙kT+RkCT.
SinceρρT≥0 and C(H˙kPk∣k-1H˙kT+Rk)CT≥0, one can get(38)EXk-X¯kXk-X¯kT≥Pk∣k∗.Note that if and only if ρ=0 and C=0, that is, X¯k=Xk∗, the error covariance is the minimum, equal to Pk∣k∗. This states that Xk∗ is the linear minimum variance estimator.Remark 9.
From Theorems7 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 X^k|k-1∈Ψk; that is, λk∗=0, μk∗=0, and then Ak and bk are empty sets; they are set to 0. In this case, it follows from (27) and (32) that Xk∗=X^k∣k and Pk∣k∗=Pk∣k; the directional constraint Bk 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 Bk in (7) as the location estimation of the target.The experiment scenario is a square sensing field of the size8×8 m2; each sensor has M=4 sensing directions with a detection radius r=3m. 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 Wk-1 can be approximated by a variable acceleration with Qk-1=Diag{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=50 and σn2=5, (ii) high density and low noise condition with N=50 and σn2=0.05, (iii) low density and high noise condition with N=10 and σn2=5, and (iv) low density and low noise condition with N=10 and σn2=0.05.The initial state estimate and the corresponding covariance matrix for these tracking algorithms are chosen to be(39)X^0∣0=12.00.73209.00.0T,P0∣0=0.01×Diag1111.100 Monte Carlo runs are carried out and the results are summarized in the following. Figure3 shows the tracking trajectories of the four tracking algorithms under different simulation scenarios. Figure 4 describes the estimated location biases (xk,yk)-(x^k∣k,y^k∣k) of target at each step for different tracking algorithms under four simulation scenarios, where (x^k∣k,y^k∣k) comes from the estimated value X^k∣k. Figure 5 represents the mean squared errors (MSEs) of the state estimate E{Xk-X^k2} at each step for these three Kalman filters (CEKF, TEKF, and TUKF).Figure 3
Tracking trajectories under different scenarios.
(a)
High density and high noise
(b)
High density and low noise
(c)
Low density and high noise
(d)
Low density and low noiseFigure 4
The estimated location biases of target under different scenarios.
(a)
High density and high noise
(b)
High density and low noise
(c)
Low density and high noise
(d)
Low density and low noiseFigure 5
The mean squared errors of the Kalman filter estimate under different scenarios.
(a)
High density and high noise
(b)
High density and low noise
(c)
Low density and high noise
(d)
Low density and low noiseIn 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 Figures4 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 Figure4, 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 Figure5, 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) E{(xk,yk)-(x^k∣k,y^k∣k)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 densityand 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.



---

*Source: 158570-2015-05-17.xml*