Distributed Iteratively Quantized Kalman Filtering for ... - CiteSeerX
Distributed Iteratively Quantized Kalman Filtering for Wireless Sensor Networks Eric J. Msechu∗, Stergios I. Roumeliotis† , Alejandro Ribeiro∗ , and Georgios B. Giannakis∗ (contact author)
Abstract— Estimation and tracking of generally nonstationary Markov processes is of paramount importance for applications such as localization and navigation. In this context, ad hoc wireless sensor networks (WSNs) offer distributed Kalman ¿ltering (KF) based algorithms with documented merits over centralized alternatives. Adhering to the limited power and bandwidth resources WSNs must operate with, this paper introduces a novel distributed KF estimator based on quantized measurement innovations. The quantized observations and the distributed nature of the iteratively quantized KF algorithm are amenable to the resource constraints of the ad hoc WSNs. Analysis and simulations show that KF-like tracking based on m bits of iteratively quantized innovations communicated among sensors exhibits MSE performance identical to a KF based on analog-amplitude observations applied to an observation model with noise variance increased by a factor of [1 − (1 − 2/π)m ]−1 . With minimal communication overhead, the mean-square error (MSE) of the distributed KF-like tracker based on 2-3 bits is almost indistinguishable from that of the clairvoyant KF.
quantized observations are generated as the sign of the innovation (SoI) sequence, a ¿lter with complexity and performance very close to the clairvoyant KF based on the analog-amplitude observations can be derived. Even though promising, the approach of [8] is limited to a particular 1-bit per observation quantizer. This paper builds on and considerably broadens the scope of [8] by addressing the middle ground between estimators based on severely quantized (1-bit) data and those based on unquantized data. The end result is a quantizer-estimator that offers desirable trade-offs between bandwidth requirements (dictating the number of quantization bits used for inter-sensor communications) and overall tracking performance (assessed by the meansquare state estimation error).
I. I NTRODUCTION
x˙ c (t) = Ac (t)xc (t) + uc (t)
Consider an ad-hoc wireless sensor network (WSN) deployed to track a Markov stochastic process. Each sensor node acquires observations which are noisy linear transformation of a common state. The sensors then transmit observations to each other in order to form a state estimate. If observations were available at a common location, minimum mean-square error (MMSE) estimates could be obtained using a Kalman ¿lter (KF). However, since observations are distributed in space and there is limited communication bandwidth, the observations have to be quantized before transmission. Thus, the original estimation problem is transformed into distributed state estimation based on quantized observations. Quantizing observations to estimate a parameter of interest, is not the same as quantizing a signal for later reconstruction [3]. Instead of a reconstruction algorithm, the objective is ¿nding, e.g., MMSE optimal, estimators using quantized observations [7]. Furthermore, optimal quantizers for reconstruction are, generally, different from optimal quantizers for estimation. State estimation using quantized observations is a non-linear estimation problem that can be solved using e.g., unscented (U)KFs [4] or particle ¿lters [2]. It was shown in [8] that if
where Ac (t) ∈ Rp×p denotes the state transition matrix and uc (t) is assumed zero-mean white Gaussian process with covariance matrix E{uc (t) uTc (τ )} = Cuc (t)δc (t − τ ). The k-th sensor Sk records scalar observations
II. M ODELS AND PROBLEM STATEMENT Consider an ad-hoc WSN whose K sensor nodes {Sk }K k=1 Keywords: wireless sensor networks, distributed state estimation, estimate a multivariate real-valued random process xc (t) ∈ Rp , Kalman ¿ltering, quantized observations, limited-rate communication. where c denotes continuous-time. The state equation is given as
Prepared through collaborative participation in the Communication and Networks Consortium sponsored by the U. S. Army Research Lab under the Collaborative Technology Alliance Program, Cooperative Agreement DAAD1901-2-0011. The U. S. Government is authorized to reproduce and distribute reprints for Government purposes notwithstanding any copyright notation thereon. ∗ Dept. of ECE, Univ. of Minnesota, 200 Union Street SE, Minneapolis, MN 55455 , Email:{emsechu,aribeiro,georgios}@ece.umn.edu, Tel/fax:(612) 6267781/ 625-4583; † Dept. of CSE, Univ. of Minnesota, 200 Union Street SE, Minneapolis, MN 55455, Email:[email protected], Tel/fax:(612) 626-7507/ (612) 625-0572.
978-1-4244-2110-7/08/$25.00 ©2007 IEEE
646
yc,k (t) = hTc,k (t)xc (t) + vc,k (t)
(1)
(2)
p
where hc,k (t) ∈ R denotes the regression vector, and vc,k (t) is a temporally and spatially white zero-mean Gaussian noise process with covariance E{vc,k (t) vc,l (τ )} = cvc (t)δ(t − τ )δkl . It is further assumed that uc (t) is independent of both vc,k (t) and xc (t0 ), where t0 is an arbitrary initial reference time. The discrete-time counterpart of (1) is obtained using the def t initions Φ(t2 , t1 ) := exp t12 Ac (t)dt , x(n) := xc (nTs ), and nTs u(n) := (n−1)T Φ(nTs , τ )uc (τ )dτ where Ts is the sampling s period. From [6, Section 4.9] the discrete-time state is x(n) = A(n) x(n − 1) + u(n)
(3)
where u(n) is zero-mean white Gaussian with covariance ma nTs T trix Cu (n) = (n−1)T Φ(nT s , τ )Cuc (τ )Φ (nTs , τ )dτ ; and s A(n) := Φ(nTs , (n − 1)Ts ). The discrete-time counterpart of the observation equation (2) is obtained as yk (n) = hTk (n) x(n) + vk (n)
(4)
where yk (n) := yc,k (nTs ) is obtained by uniform sampling of (2) followed by low- or band-pass ¿ltering with bandwidth Ts , leading to zero-mean white Gaussian discrete-time noise vk (n) with variance cv (n) = cvc (nTs )/Ts [6, Section 4.9]. Supposing that A(n), Cu (n), hk (n) and cv (n) are available ∀ n, k, the goal of the WSN is for each sensor Sk to form
an estimate of x(n). Estimating x(n) necessitates each sensor Sk to communicate yk (n) to all other sensors {Sl }K l=1,l=k . Communication takes place over the shared wireless channel that we assume can afford transmission of a single packet of m bits per time slot n. This leads to a one-to-one correspondence between time n and sensor index k and allows us to drop the sensor argument k in (4).
A. MMSE estimation with quantized observations Let the quantization at time index n be de¿ned as b(n) := qn [y(n)], then given current and past messages b1:n := ˆ (n|b1:n ) {b(1), b(2), . . . , b(n)} we are interested in estimates x of the state x(n) using the information in b1:n . The mean-square error (MSE) of the estimator is obtained as the trace of the x(n|b1:n ) − error covariance matrix (ECM) M(n|b1:n ) := E{[ˆ x(n)][ˆ x(n|b1:n ) − x(n)]T }, i.e., tr{M(n|b1:n )}. The MMSE estimator is the conditional mean, see e.g., [6, Chapter 5], ˆ (n|b1:n ) = E{x(n)|b1:n } := x(n) p[x(n)|b1:n ] dx(n). (5) x Rp
ˆ (n|y1:n ) in (5) and its ECM M(n|y1:n ) are given as x ˆ (n|y1:n ) = x ˆ (n|y1:n−1 )+ x M(n|y1:n−1 )h(n) y˜(n|y1:n−1 ) hT (n)M(n|y1:n−1 )h(n) + cv (n) M(n|y1:n ) = M(n|y1:n−1 )− M(n|y1:n−1 )h(n)hT (n)M(n|y1:n−1 ) . hT (n)M(n|y1:n−1 )h(n) + cv (n)
(10)
(11)
Computations for the KF iteration in [P1]-[C1] require a few algebraic operations per time-step n whereas (6) - (7) require: (i) multi-dimensional numerical integration to obtain the predicted pdf p[x(n)|b1:n−1 ] in (7) and to evaluate the expectation in (5); and (ii) numerical update of the posterior pdf p[x(n)|b1:n ] in (6). This high computational cost is inherent to non-linear models (non-linearity in this paper is due to quantization of the observations) and motivates lower complexity approximations. By using a Gaussian approximation for the prior pdf p[x(n)|b1:n−1 ], see e.g., [5], tracking of the potentially intractable pdf p[x(n)|b1:n−1 ] simpli¿es to keeping track of its mean and covariance matrix. This leads to iterations similar to [P1]-[C1] as detailed in the next section.
ˆ (n|b1:n ), the posterior III. D ISTRIBUTED ITERATIVELY QUANTIZED K ALMAN FILTER To obtain a closed-form expression for x distribution p[x(n)|b1:n ] has to be known and the integral in (5) To design the iteratively quantized Kalman ¿lter (IQKF), let needs to be computed. In principle, p[x(n)|b1:n ] can be obtained sensors rely on m-bit binary messages b(n) := b(1:m) (n) := recursively using Bayes’ rule as follows: [b(1) (n), . . . , b(m) (n)], with the i-th bit b(i) (n) de¿ned as the sign of innovations [cf. (13)] conditioned on the previous messages Pr{b(n)|x(n), b1:n−1 } (6) b1:n−1 := [b(1), . . . , b(n − 1)] and previous bits b(1:i−1) (n) := p[x(n)|b1:n ] = p[x(n)|b1:n−1 ] Pr{b(n)|b1:n−1 } [b(1) (n), . . . , b(i−1) (n)] of the current (n-th) message. Speci¿cally, let where Pr{b(n)|x(n), b1:n−1 } and Pr{b(n)|b1:n−1 } depend on the quantization rule qn [y(n)]. If p[x(n − 1)|b1:n−1 ] is known, yˆ(0) (n|b1:n−1 ) :=ˆ y(n|b1:n−1 ) = E {y(n)|b1:n−1 } the prior pdf p[x(n)|b1:n−1 ] can be obtained as yˆ(i) (n|b1:n−1 ) :=E y(n)|b1:n−1 , b(1:i) (n) , for i ≥ 1 (12) p[x(n)|b1:n−1 ] (7) stand for MMSE estimates of y(n) using past messages b1:n−1 and the ¿rst i bits of the current message denoted as b(1:i) (n). = p[x(n)|x(n−1), b1:n−1 ]p[x(n−1)|b1:n−1]dx(n−1) The i-th bit of the current message, b(i) (n), is de¿ned as Rp where due to Gaussian and Markov properties of the random process x(n) in (3), p[x(n)|x(n−1), b1:n−1 ] = p[x(n)|x(n−1)] = N [x(n); A(n)x(n−1), Cu (n)]. If b(n) = y(n), i.e., un-quantized observations are used, both conditional pdfs in (7) and (6) are Gaussian and it suf¿ces to propagate their means and covariance matrices only. This leads to the following Kalman ¿lter recursions [6, Chapter 5.3]: [P1] KF prediction step. Given the previous estimate ˆ (n − 1|y1:n−1 ) and its ECM M(n − 1|y1:n−1 ), x ˆ (n|y1:n−1 ) the predicted estimate x := } and its ECM M(n|y E{x(n)|y1:n−1 1:n−1 ) := E [ˆ x(n|y1:n−1 ) − x(n)][ˆ x(n|y1:n−1 ) − x(n)]T are obtained as
y (i−1) (n|b1:n−1 )]. b(i) (n) := sign[y(n) − yˆ(i−1)(n|b1:n−1 )] := sign[˜ (13) Our goal is to iteratively compute the estimates ˆ (i) (n|b1:n−1 ) := E x(n)|b1:n−1 , b(1:i) (n) x (14) ˆ (i−1) (n|b1:n−1 ) and b(i) (n). When using m bits, based on x b(1:m) (n), we will refer to the resulting algorithm as m-IQKF. A possible option for de¿ning the quantizer would be to x(i) (n|b1:n−1 ). However, set yˆ(i) (n|b1:n−1 ) equal to hT (n)ˆ (i) T (i) yˆ (n|b1:n−1 ) = h (n)ˆ x (n|b1:n−1 ) if i ≥ 1 as explained next. Using the observations (4), the de¿nitions of yˆ(i) (n|b1:n−1 ) ˆ (i) (n|b1:n−1 ) in (14), we obtain in (12), and x yˆ(i) (n|b1:n−1 ) = E hT (n)x(n) + v(n)|b1:n−1 , b(1:i) (n)
x(i) (n|b1:n−1 ) + E{v(n)|b1:n−1 , b(1:i) (n)}. (15) = hT (n)ˆ ˆ (n|y1:n−1 )=A(n)ˆ x x(n−1|y1:n−1 ) (8) M(n|y1:n−1 )=A(n)M(n−1|y1:n−1 )AT(n)+Cu(n). (9) The noise estimate E{v(n)|b1:n−1 , b(1:i) (n)} is not necessarily zero - see also (22). Therefore, in orderto obtain yˆ(i) (n|b1:n−1) ˆ (i) (n|b1:n−1 ) as well as E v(n)|b1:n−1 , b(1:i) (n) [C1] KF correction step. Consider the predicted data we need x E{y(n)|y1:n−1 } := yˆ(n|y1:n−1 ) = hT (n)ˆ x(n|y1:n−1 ) and which is achieved by augmenting the state vector x(n) with the their innovation y˜(n|y1:n−1 ) := y(n) − yˆ(n|y1:n−1 ). Then noise term v(n) as described in the next section.
647
A. State augmentation state the estimates augmentation Through E v(n)|b1:n−1 , b(1:i) (n) can be obtained so that the predicted observation yˆ(i) (n|b1:n−1 ) in (15) can be evaluated. Speci¿cally, ˘ (n) := [xT (n), v(n)]T , u ˘ (n) := [uT (n), v(n)]T , let x ˘ ˘ h(n) := [hT (n), 1]T , and A(n) de¿ned with A(n) as the leading p × p submatrix and with zeros in all other entries; then model in (3)-(4) can be rewritten as ˘ (n) = x
˘ ˘ (n − 1) + u ˘ (n) A(n) x T ˘ ˘ (n) + v˘(n) h (n) x
ˆ ˘ (0) (n|b1:n−1 ) := where we used the de¿nitions x (0) ˆ ˘ ˘ ˘ x(n|b1:n−1 ) and M (n|b1:n−1 ) := M(n|b 1:n−1 ). For time index n, (20) and (21) are repeated m-times. ˆ ˘ (n) ˘ (n|b x The MMSE estimate of x given b1:n is 1:n ) = (1:m) ˘ (n)|b1:n−1 , b E {˘ x(n)|b1:n } = E x (n) = ˆ ˘ ˘ (m) (n|b1:n−1 ). The corresponding ECM is M(n|b x 1:n ) := ˘ (m) (n|b1:n−1 ). M
(16) The state estimates x ˆ (i) (n|b1:n−1 ) in (14) are the ¿rst p comˆ ˘ (i) (n|b1:n−1 ), i.e., y˘(n) = (17) ponents of the augmented state estimate x ˆ ˆ (i) (n|b1:n−1 ) = [x ˘ (i) (n|b1:n−1 )]1:p and the noise estimate is x where the new observation noise v˘(n) = 0, ∀n (by construc- the (p + 1)-st component. tion). Note that the covariance matrix of the augmented driving noise is a block-diagonal matrix Cu˘ (n) ∈ R(p+1)×(p+1) with Corollary 1 For i = 1, E v(n)|b1:n−1 , b(1:i) (n) = 0. [Cu˘ (n)]1:p,1:p = Cu (n) and [Cu˘ (n)]p+1,p+1 = cv (n). The augmented state formulation (16)-(17) increases the diˆ mension of the state vector but is otherwise equivalent to (3)-(4). ˘ (i) (n|b Proof: The (p+1)-st entry of x 1:n−1 ) in (20) is the However, it has the appealing property that MMSE estimates noise estimate E v(n)|b1:n−1 , b(1:i) (n) . Thus, for i = 1 ˘ (n) contain MMSE estimates of the of the augmented state x original state x(n) and of the observation noise v(n) which E{v(n)|b1:n−1 , b(n)}
allows computation of yˆ(i) (n|b1:n−1 ) in (12) as follows. Since ˘ ˘ (0) (n|b1:n−1 )h(n)] 2 [M p+1 T T ˘ b(n) = E{v(n)|b1:n−1 } + y˘(n) = h (n)˘ x(n) = h (n)x(n) + v(n) = y(n), it follows π (i) (i) T (i) T (0) ˘ ˘ ˘ ˆ ˆ ˘ (n) M (n|b ) h(n) h ˘ ) = h (n) x (n|b ) that yˆ (n|b1:n−1 ) = y˘ (n|b1:n−1 1:n−1 1:n−1
ˆ˘ (i) (n|b1:n−1 ) = E x ˘ (n)|b1:n−1 , b(1:i) (n) and where x cv (n) 2 b(n) = 0. (22) = y˘ˆ(i) (n|b1:n−1 ) = E y˘(n)|b1:n−1 , b(1:i) (n) . Using the augπ hT (n)M(n|b1:n−1 )h(n) + cv (n) mented state model, we obtain the MMSE estimates in (14) using the algorithm detailed in Proposition 1. The last equality follows since E {v(n)|b1:n−1 } = E{v(n)} = 0 ˘ ˘ (0) (n|b1:n−1 )h(n) and M = cv (n). Proposition 1 Consider the augmented state space p+1 model in (16)-(17). De¿ne the augmented state The similarity of the m-IQKF in Proposition 1 and the clairˆ˘ (n|b1:n−1 ) estimates x := E{˘ x(n)|b1:n−1 } and voyant KF based on un-quantized observations (cf. [P1]-[C1]) ˆ˘ (i−1) (n|b1:n−1 ) ˘ (n)|b1:n−1 , b(1:i−1) (n) . Let is quite remarkable. The ECM updates in (21) are identical to := E x x ˆ ˆ˘ (n|b1:n−1 ) − x ˘ ˘ (n|b1:n−1 ) − x ˘ (n)][x ˘ (n)]T } the ECM updates of the KF in (11) except for the scale factor M(n|b1:n−1 ) := E{[x (i−1) (i−1) ˆ ˘ ˘ (n|b1:n−1 ) := E{[x (n|b1:n−1 ) − 2/π. The variance penalties associated with the m-IQKF are and M ˆ˘ (i−1) (n|b1:n−1 ) − x ˘ (n)]T } denote the corresponding investigated in the following corollaries. ˘ (n)][x x ECMs for i = 1, . . . , m. Construct the messages b(1:m) (n) ˆ ˘ T (n)x ˘ (i−1) (n|b1:n−1 ) . Corollary 2 Consider the m-IQKF algorithm in Proposition 1 in (13) as b(i) (n) := sign y(n) − h ˆ˘ (i) (n|b1:n−1 ) is obtained from the recursion: and de¿ne the ECM reduction after i = 1, . . . , m iteraThe estimate x ˘ (0) (n|b1:n−1 ) − M ˘ (i) (n|b1:n−1 ), where ˘ i (n) := M tions as ΔM ˆ˘ (n − 1|b1:n−1 ) and its ECM ˘ (0) 2 i ˘ [P2] Given the previous estimate x M (n|b1:n−1 ) := M(n|b ). With c 1:n−1 i := 1 − (1 − π ) , the ˘ M(n−1|b ), form 1:n−1 error covariance reduction after i iterations is given as ˘ ˘ˆ (n|b1:n−1 )=A(n) ˘ˆ (n−1|b1:n−1) x x (18) ˘ ˘T ˘ (0) ˘ (0) ˘ i (n) = ci M (n|b1:n−1 )h(n)h (n)M (n|b1:n−1 ) . ˘ ˘ ˘ ˘T ΔM M(n|b 1:n−1 )=A(n)M(n−1|b1:n−1 )A (n)+Cu ˘ (n)(19) ˘ ˘ T (n)M ˘ (0) (n|b1:n−1 )h(n) h = [C2] Assuming that p[˘ x(n)|b1:n−1 , b(1:i−1) (n)] (23) ˆ˘ (i−1) (n|b1:n−1 ), M ˘ (i−1) (n|b1:n−1 )], the MMSE N [˘ x(n); x ˆ ˘ i (n) as a summation of differences ˘ (i) (n|b1:n−1 ) := E{˘ estimate x x(n)|b1:n−1 , b(1:i) (n)} Proof: We ¿rst write ΔM ˘ (i) (n|b1:n ) are obtained between successive ECM matrices and the corresponding ECM M iteratively from ˘ (0) (n|b1:n−1 ) − M ˘ (i) (n|b1:n−1 ) ˘ i (n) := M ΔM ˆ ˆ˘ (i) (n|b1:n−1 ) = x ˘ (i−1) (n|b1:n−1 ) + x i
˘ (i−1) (n|b1:n−1 )h(n) ˘ M 2 ˘ (j−1) (n|b1:n−1 ) − M ˘ (j) (n|b1:n−1 ) = M b(i)(n) (20) π ˘T j=1 ˘ ˘ (i−1) (n|b1:n−1 )h(n) h (n)M i ˘ h ˘ T (n)M ˘ (j−1) (n|b1:n−1 )h(n) ˘ (j−1) (n|b1:n−1 ) 2 M ˘ (i−1) (n|b1:n−1 ) − ˘ (i) (n|b1:n−1 ) = M M = T (j−1) ˘ ˘ ˘ π j=1 h (n)M (n|b1:n−1 )h(n) ˘ (i−1) (n|b1:n−1 )h(n) ˘ h ˘ T (n)M ˘ (i−1) (n|b1:n−1 ) 2M (24) ˘ ˘ T (n)M ˘ (i−1) (n|b1:n−1 )h(n) π h (21) where the last equality follows from (21). Next, we multiply
648
TABLE I P ER STEP FACTOR , cm , AND NOISE PENALTY FACTOR FOR IQKF
˘ ˘ (i) (n|b1:n−1 ) from (21) by h(n) M to obtain ˘ ˘ (i−1) (n|b1:n−1 )h(n) ˘ ˘ (i) (n|b1:n−1 )h(n) =M − M (i−1) T (i−1) ˘ ˘ ˘ ˘ ˘ (n|b1:n−1 )h(n)h (n)M (n|b1:n−1 )h(n) 2M ˘ ˘ T (n)M ˘ (i−1) (n|b1:n−1 )h(n) h 2 ˘ (i−1) ˘ M (n|b1:n−1 )h(n) . = 1− π
π
bits, m cm (1/cm − 1)100%
(25)
1 0.637=2/π 57.08%
2 0.868 15.21%
3 0.952 5.04%
4 0.983 1.77%
(26) B. Performance analysis of the m-IQKF In the previous section, we quanti¿ed the per correction step ˘ m (n) in the m-IQKF as a function of the ˘ ˘ (i−1) (n|b1:n−1 )h(n) ECM reduction ΔM Repeating steps (25)-(26) for M ˘ ˘ (i) (n|b1:n−1 )h(n) with decreasing index i yields M = number of bits m used for iteratively quantizing the observations i (0) ˘ ˘ (n|b1:n−1 )h(n) which when substituted in (24) y(n). We next compare MSE performance of m-IQKF with that 1− π2 M of the clairvoyant KF, when both prediction and correction steps with i = j − 1, results in are considered, by deriving the continuous-time Algebraic Riccati ˘ i (n) = ΔM Equations (ARE) for both ¿lters. j−1 ˘ (0) i Consider ¿rst the discrete time ARE for the m-IQKF [P2]T (0) ˘ ˘
˘ M (n|b1:n−1 )h(n)h (n)M (n|b1:n−1 ) 2 2 ˘ ˘ + 1) := M(n + 1|b1:n ), . [C2]. To simplify notation let M(n 1− ˘ ˘ T (n)M ˘ (0) (n|b1:n−1 )h(n) π j=1 π h ˘ ˘ M(n) := M(n|b1:n−1 ), M(n + 1) := M(n + 1|b1:n ), M(n) := M(n|b1:n−1 ), and substitute (27) in (19) to obtain the m-IQKF Upon invoking the geometric sum identity ARE for the ECM of x ˘ (n) as i 2 2 j−1 = 1 − (1 − π2 )i , (23) follows. j=1 1 − π π ˘ ˘ ˘ ˘T (29) From the KF’s ECM in (11), it follows readily that the ECM re- M(n + 1) =A(n + 1)M(n)A (n + 1) + Cu˘ (n + 1) T T ˘ ˘ ˘ ˘ ˘ ˘ duction for un-quantized observations M(n|b1:n−1 )−M(n|b1:n ) A(n + 1)M(n)h(n)h (n)M(n)A (n + 1) − cm . corresponds with cm = 1 in (23). Therefore, the ECM reduction ˘ ˘T (n)M(n) ˘ h(n) h achieved by quantizing to m bits is cm times smaller than the one ˘ ˘ ˘ + 1), Cu˘ (n + 1), h(n) in (29), and M(n) achieved with analog-amplitude observations. Values of cm are Substituting for A(n shown in Table I. With only 4 bits of quantization, the value of and writing only the leading p × p sub-matrices of the resulting cm is just 2% less than the value for the clairvoyant KF (cm = 1). expression, yields the ARE for the ECM of x(n) as From Corollary 2 we can relate the predicted MSE M(n + 1) =A(n + 1)M(n)AT(n + 1) + C (n + 1) (30) u ˘ ˘ tr{M(n|b 1:n−1 )} and the corrected MSE tr{M(n|b1:n )} afT T A(n + 1)M(n)h(n)h (n)M(n)A (n + 1) (1:m) (n), i.e., ter observing the m-bits of n-th observation, b − cm . T (n)M(n)h(n) + c (n) h ˘ ˘ ˘ v M(n|b1:n ) = M(n|b1:n−1 ) − ΔMm (n). Substituting for ˘ m (n) we obtain Interestingly, the resulting ARE for the ECM of the m-IQKF ΔM [cf. (30)] becomes identical to the ARE of the clairvoyant KF ˘ ˘ M(n|b 1:n ) = M(n|b1:n−1 )− [6, Chapter 5.11] as limm→∞ cm = 1. The implication of these ˘ ˘T ˘ ˘ relations follows from the continuous-time ARE for the m-IQKF M(n|b 1:n−1 )h(n)h (n)M(n|b1:n−1 ) − cm (27) given in Proposition 2. T ˘ ˘ ˘ h (n)M(n|b1:n−1 )h(n) Equation (27) is instructive for understanding the MSE perfor- Proposition 2 The continuous-time ECM Mc (t) = Mc (nTs ) := mance of m-IQKF since for cm = 1, (27) coincides with the cor- limTs →0 M(n; Ts ) for the m-IQKF of (18)-(21) is the solution rection ECM of the KF in (11). Furthermore, limm→∞ cm → 1, of the following differential (Riccati) equation: implying that for in¿nite number of quantization bits we recover ˙ c (t) =Ac (t)Mc (t) + Mc (t)AT (t) + Cu (t) M c the clairvoyant KF. This important consistency result is summa −1 cvc (t) rized in the following corollary. − Mc (t)hc (t) hTc (t)Mc (t) . (31) cm Proof: The derivation steps from (30) to (31) are the same Corollary 3 As the number of quantization bits m → ∞, the as those followed for deriving the KF Riccati equation from the correction step ECM at time n of the m-IQKF converges to the corresponding ECMs [6, pp. 259]. The only difference is that the ECM of the clairvoyant KF given that b1:n−1 = y1:n−1 , i.e., scale factor, cm , in (30) equals 1 in the case of the KF. ˘ ˘ (m) (n|b1:n−1 ) = M(n|y (28) lim M From the clairvoyant KF ARE [6, pp. 259] 1:n ) . m→∞
Proof: Since limm→∞ cm = limm→∞ 1 − (1 − π2 )m = 1, ˘ (m) (n|b1:n−1 ) [cf. (27)] and M(n|y ˘ then for cm = 1, M 1:n ) [cf. (11)] become identical, which results in (28). Corollary 3 establishes that m-IQKF asymptotically (as the number of quantization bits increases) achieves the per correction step MSE performance of the clairvoyant KF. As demonstrated by simulations in Section IV, even a small number of quantization bits (m = 2 or m = 3) renders the MSE performance of m-IQKF indistinguishable from that of the clairvoyant KF.
˙ K (t) =Ac (t)MK (t) + MK (t)AT (t) + Cu (t) M c c c c
−1 T hc (t)MK c (t)
− MK c (t)hc (t) [cvc (t)] MK c (t)
(32)
denotes the continuous-time ECM of KF, it is where evident that the variance of the continuous-time observation noise of the m-IQKF is ampli¿ed by 1/cm compared to that of the clairvoyant KF. This is the price paid for using quantized observations, b(n), instead of the analog-amplitude ones, y(n). In Table I the percentage observation noise variance increase (1/cm −1)100% versus number of quantization bits, m is shown.
649
IV. S IMULATIONS
distributed estimation, it was shown that quantization using 2 to 3 bits improves the performance of the Kalman-like algorithm to virtually coincide with the optimal state estimator (Kalman ¿lter), with only minimal increase in computation. This result was corroborated by simulation results. It was further shown that true ¿lter covariances are consistent with Noisy measurements at sensor Sk are given as yk (t) = xc1 (t) + through simulations 1 . analytical ones θk xc2 (t) + vk (t) where θk is a parameter at sensor k. A discreteR EFERENCES time equivalent model is used with sampling time Ts = 1. [1] Y. Bar-Shalom and X. Li, Estimation and Tracking: Principles, Techniques, WSN data corresponding to K = 2 sensors with [θ1 θ2 ] = and Software. Artech House, 1993. [0.1 0.2] are simulated. The measurement and state driving noise [2] P. Djuric, J. Kotecha, J. Zhang, Y. Huang, T. Ghirmai, M. Bugallo, and J. Miguez, “Particle ¿ltering,” IEEE Sig. Proc. Mag., vol. 20, pp. 19–38, variances are cv (t) = 1 and cu (t) = 1. The MSE plots in Fig. 1 Sep. 2003. and Fig. 2 illustrate the evolution of the MSEs, obtained from [3] R. M. Gray, “Quantization in task-driven sensing and distributed processing,” the trace of the respective ECMs, against the time index n. in Proc. Int. Conf. Acoustics, Speech, Signal Processing, vol. 5, pp. V–1049– V–1052, Toulouse, France, May 14-19, 2006. In the ¿rst simulation setup, Fig. 1, the MSE of the IQKF, given [4] S. Julier and J. Uhlmann, “Unscented ¿ltering and nonlinear estimation,” by tr{M(n|b1:n )}, is compared to the MSE of the clairvoyant Proc. IEEE, vol. 92, pp. 401–422, Mar. 2004. KF. The simulations for IQKF are performed for 1, 2 and 3 bits. [5] J. Kotecha and P. Djuric, “Gaussian particle ¿ltering,” IEEE Trans. on Sig. Processing, vol. 51, pp. 2602–2612, Oct. 2003. The plots demonstrate that there is a substantial MSE reduction [6] P. S. Maybeck, Stochastic Models, Estimation and Control – Vol.1. Academic when going from 1 bit of quantization to 2 bits and little MSE Press, 1979. reduction for higher number of quantization bits. With 2 bits of [7] H. Papadopoulos, G. Wornell, and A. Oppenheim, “Sequential signal encoding from noisy measurements using quantizers with dynamic bias control,” quantization the MSE of IQKF is virtually coincident with that IEEE Trans. on Infor. Theory, vol. 47, pp. 978–1002, 2001. of the clairvoyant KF as was postulated by the analytical values [8] A. Ribeiro, G. B. Giannakis, and S. I. Roumeliotis, “SOI-KF: Distributed in Table I whereby c2 = 86.8%. Kalman ¿ltering with low-cost communications using the sign of innovations,” IEEE Trans. on Sig. Proc., vol. 54, no. 12, pp. 4782–4795, Dec 2006. Model consistency checks comparing the empirically obtained MSEs with analytical MSEs are shown in Fig. 2. Note that KF 1íbit IQKF the analytical MSE is obtained from the trace of the ECM 2íbit IQKF 3íbit IQKF ˆ (n|b1:n )][x(n) − x ˆ (n|b1:n )]T |b1:n }. M(n|b1:n ) := E{[x(n) − x The empirical MSE is the sample estimator of tr{M(n|b1:n )} obtained as a sample average of the squared estimation errors. The plots depict both predictor MSE tr{M(n|b1:n−1 )} and estimator MSE tr{M(n|b1:n )}. The consistency check reveals that the empirical and analytical MSEs are nearly identical. Fig. 3 shows alternative model consistency tests for time index, n the IQKF using the normalized estimation error squared Fig. 1. IQKF estimator MSE tr{M(n|b1:n )} from Monte Carlo data for 1-3 (NEES) tests of [1, Ch. 5.4]. NEES r(n) := [x(n) − bits of quantization; MSE for clairvoyant KF is included for comparison. ˆ (n|b1:n )]T [M(n|b1:n )]−1 [x(n) − x ˆ (n|b1:n )] is postulated to x IQKF predictor MSEí analytical IQKF predictor MSEí empirical have a χ2 pdf with p degrees of freedom (since the p-dimensional IQKF estimator MSEí analytical IQKF estimator MSE í empirical ˜ (n) := x(n)− x ˆ (n|b1:n ) is assumed zero-mean Gaussian with x covariance M(n|b1:n ) if the estimator is consistent with the model). Under the hypothesis that the estimator is consistent, L realizations of the NEES statistics {ri (n)}L i=1 each χ2 distributed with p degrees of freedom, lead to a χ2 distribution with Lp degrees of freedom. This is checked by running Monte Carlo simulations and computing the sample average NEES rˆ(n) := L time index, n 1 r (n) and then de¿ning an acceptance (con¿dence) Fig. 2. Comparing empirical MSE from Monte Carlo data and analytical MSE, i i=1 L region (for the consistent hypothesis). If l ≤ rˆ(n) ≤ u, then for predictor ECM M(n|b1:n−1 ) and estimator ECM M(n|b1:n ). the estimator is consistent; lower and upper bounds l and u are obtained from Pr{ˆ rn ∈ [l, u]} = 1 − α, where 1 − α is the probability of acceptance region. Using L = 50 realizations, 200 time samples, a state space of p = 2 dimensions, and α = 0.05 (i.e., 95% region), we observed that indeed only 7% of the 100 time samples lie outside the 95% acceptance region. The two-dimensional state space model simulated is x˙ c1 (t) 0 1 0 xc1 (t) x˙ c (t) := = + uc (t) . x˙ c2 (t) xc2 (t) 0 0 1
5.5
meanísquare error (MSE)
5
4.5 4
3.5 3
2.5 2
1.5 1
5
10
15
20
25
30
35
40
45
50
20
25
30
35
40
45
50
12
meanísquare error (MSE)
10
8
6
4
2
0
5
10
15
Normalized Est. Error Squared, NEES
4
3.5
NEES, r(n)
3
2.5 2
1.5 1
V. C ONCLUDING REMARKS
0.5
Recursive state estimation based on quantized observations was considered. Multi-bit quantization was done using iterative binary quantizing of the measurement innovations. Joint quantization and estimation was used to develop a Kalman-like algorithm. Motivated by the need to ¿nd quanti¿able trade-offs between estimation performance and number of quantization bits for
0
50
100
150
200
time index, n
Fig. 3. Consistency test using normalized estimation error squared (NEES), r(n), compared with the 95% probability region using 2 bits of quantization. 1 The views and conclusions contained in this document are those of the authors and should not be interpreted as representing the of¿cial policies of the Army Research Laboratory or the U. S. Government.
650
Abstract— Estimation and tracking of generally nonstationary Markov processes is of paramount importance for applications such as localization and navigation. In this context, ad hoc wireless sensor networks (WSNs) offer distributed Kalman ¿ltering (KF) based algorithms with documented merits over centralized alternatives. Adhering to the limited power and bandwidth resources WSNs must operate with, this paper introduces a novel distributed KF estimator based on quantized measurement innovations. The quantized observations and the distributed nature of the iteratively quantized KF algorithm are amenable to the resource constraints of the ad hoc WSNs. Analysis and simulations show that KF-like tracking based on m bits of iteratively quantized innovations communicated among sensors exhibits MSE performance identical to a KF based on analog-amplitude observations applied to an observation model with noise variance increased by a factor of [1 − (1 − 2/π)m ]−1 . With minimal communication overhead, the mean-square error (MSE) of the distributed KF-like tracker based on 2-3 bits is almost indistinguishable from that of the clairvoyant KF.
quantized observations are generated as the sign of the innovation (SoI) sequence, a ¿lter with complexity and performance very close to the clairvoyant KF based on the analog-amplitude observations can be derived. Even though promising, the approach of [8] is limited to a particular 1-bit per observation quantizer. This paper builds on and considerably broadens the scope of [8] by addressing the middle ground between estimators based on severely quantized (1-bit) data and those based on unquantized data. The end result is a quantizer-estimator that offers desirable trade-offs between bandwidth requirements (dictating the number of quantization bits used for inter-sensor communications) and overall tracking performance (assessed by the meansquare state estimation error).
I. I NTRODUCTION
x˙ c (t) = Ac (t)xc (t) + uc (t)
Consider an ad-hoc wireless sensor network (WSN) deployed to track a Markov stochastic process. Each sensor node acquires observations which are noisy linear transformation of a common state. The sensors then transmit observations to each other in order to form a state estimate. If observations were available at a common location, minimum mean-square error (MMSE) estimates could be obtained using a Kalman ¿lter (KF). However, since observations are distributed in space and there is limited communication bandwidth, the observations have to be quantized before transmission. Thus, the original estimation problem is transformed into distributed state estimation based on quantized observations. Quantizing observations to estimate a parameter of interest, is not the same as quantizing a signal for later reconstruction [3]. Instead of a reconstruction algorithm, the objective is ¿nding, e.g., MMSE optimal, estimators using quantized observations [7]. Furthermore, optimal quantizers for reconstruction are, generally, different from optimal quantizers for estimation. State estimation using quantized observations is a non-linear estimation problem that can be solved using e.g., unscented (U)KFs [4] or particle ¿lters [2]. It was shown in [8] that if
where Ac (t) ∈ Rp×p denotes the state transition matrix and uc (t) is assumed zero-mean white Gaussian process with covariance matrix E{uc (t) uTc (τ )} = Cuc (t)δc (t − τ ). The k-th sensor Sk records scalar observations
II. M ODELS AND PROBLEM STATEMENT Consider an ad-hoc WSN whose K sensor nodes {Sk }K k=1 Keywords: wireless sensor networks, distributed state estimation, estimate a multivariate real-valued random process xc (t) ∈ Rp , Kalman ¿ltering, quantized observations, limited-rate communication. where c denotes continuous-time. The state equation is given as
Prepared through collaborative participation in the Communication and Networks Consortium sponsored by the U. S. Army Research Lab under the Collaborative Technology Alliance Program, Cooperative Agreement DAAD1901-2-0011. The U. S. Government is authorized to reproduce and distribute reprints for Government purposes notwithstanding any copyright notation thereon. ∗ Dept. of ECE, Univ. of Minnesota, 200 Union Street SE, Minneapolis, MN 55455 , Email:{emsechu,aribeiro,georgios}@ece.umn.edu, Tel/fax:(612) 6267781/ 625-4583; † Dept. of CSE, Univ. of Minnesota, 200 Union Street SE, Minneapolis, MN 55455, Email:[email protected], Tel/fax:(612) 626-7507/ (612) 625-0572.
978-1-4244-2110-7/08/$25.00 ©2007 IEEE
646
yc,k (t) = hTc,k (t)xc (t) + vc,k (t)
(1)
(2)
p
where hc,k (t) ∈ R denotes the regression vector, and vc,k (t) is a temporally and spatially white zero-mean Gaussian noise process with covariance E{vc,k (t) vc,l (τ )} = cvc (t)δ(t − τ )δkl . It is further assumed that uc (t) is independent of both vc,k (t) and xc (t0 ), where t0 is an arbitrary initial reference time. The discrete-time counterpart of (1) is obtained using the def t initions Φ(t2 , t1 ) := exp t12 Ac (t)dt , x(n) := xc (nTs ), and nTs u(n) := (n−1)T Φ(nTs , τ )uc (τ )dτ where Ts is the sampling s period. From [6, Section 4.9] the discrete-time state is x(n) = A(n) x(n − 1) + u(n)
(3)
where u(n) is zero-mean white Gaussian with covariance ma nTs T trix Cu (n) = (n−1)T Φ(nT s , τ )Cuc (τ )Φ (nTs , τ )dτ ; and s A(n) := Φ(nTs , (n − 1)Ts ). The discrete-time counterpart of the observation equation (2) is obtained as yk (n) = hTk (n) x(n) + vk (n)
(4)
where yk (n) := yc,k (nTs ) is obtained by uniform sampling of (2) followed by low- or band-pass ¿ltering with bandwidth Ts , leading to zero-mean white Gaussian discrete-time noise vk (n) with variance cv (n) = cvc (nTs )/Ts [6, Section 4.9]. Supposing that A(n), Cu (n), hk (n) and cv (n) are available ∀ n, k, the goal of the WSN is for each sensor Sk to form
an estimate of x(n). Estimating x(n) necessitates each sensor Sk to communicate yk (n) to all other sensors {Sl }K l=1,l=k . Communication takes place over the shared wireless channel that we assume can afford transmission of a single packet of m bits per time slot n. This leads to a one-to-one correspondence between time n and sensor index k and allows us to drop the sensor argument k in (4).
A. MMSE estimation with quantized observations Let the quantization at time index n be de¿ned as b(n) := qn [y(n)], then given current and past messages b1:n := ˆ (n|b1:n ) {b(1), b(2), . . . , b(n)} we are interested in estimates x of the state x(n) using the information in b1:n . The mean-square error (MSE) of the estimator is obtained as the trace of the x(n|b1:n ) − error covariance matrix (ECM) M(n|b1:n ) := E{[ˆ x(n)][ˆ x(n|b1:n ) − x(n)]T }, i.e., tr{M(n|b1:n )}. The MMSE estimator is the conditional mean, see e.g., [6, Chapter 5], ˆ (n|b1:n ) = E{x(n)|b1:n } := x(n) p[x(n)|b1:n ] dx(n). (5) x Rp
ˆ (n|y1:n ) in (5) and its ECM M(n|y1:n ) are given as x ˆ (n|y1:n ) = x ˆ (n|y1:n−1 )+ x M(n|y1:n−1 )h(n) y˜(n|y1:n−1 ) hT (n)M(n|y1:n−1 )h(n) + cv (n) M(n|y1:n ) = M(n|y1:n−1 )− M(n|y1:n−1 )h(n)hT (n)M(n|y1:n−1 ) . hT (n)M(n|y1:n−1 )h(n) + cv (n)
(10)
(11)
Computations for the KF iteration in [P1]-[C1] require a few algebraic operations per time-step n whereas (6) - (7) require: (i) multi-dimensional numerical integration to obtain the predicted pdf p[x(n)|b1:n−1 ] in (7) and to evaluate the expectation in (5); and (ii) numerical update of the posterior pdf p[x(n)|b1:n ] in (6). This high computational cost is inherent to non-linear models (non-linearity in this paper is due to quantization of the observations) and motivates lower complexity approximations. By using a Gaussian approximation for the prior pdf p[x(n)|b1:n−1 ], see e.g., [5], tracking of the potentially intractable pdf p[x(n)|b1:n−1 ] simpli¿es to keeping track of its mean and covariance matrix. This leads to iterations similar to [P1]-[C1] as detailed in the next section.
ˆ (n|b1:n ), the posterior III. D ISTRIBUTED ITERATIVELY QUANTIZED K ALMAN FILTER To obtain a closed-form expression for x distribution p[x(n)|b1:n ] has to be known and the integral in (5) To design the iteratively quantized Kalman ¿lter (IQKF), let needs to be computed. In principle, p[x(n)|b1:n ] can be obtained sensors rely on m-bit binary messages b(n) := b(1:m) (n) := recursively using Bayes’ rule as follows: [b(1) (n), . . . , b(m) (n)], with the i-th bit b(i) (n) de¿ned as the sign of innovations [cf. (13)] conditioned on the previous messages Pr{b(n)|x(n), b1:n−1 } (6) b1:n−1 := [b(1), . . . , b(n − 1)] and previous bits b(1:i−1) (n) := p[x(n)|b1:n ] = p[x(n)|b1:n−1 ] Pr{b(n)|b1:n−1 } [b(1) (n), . . . , b(i−1) (n)] of the current (n-th) message. Speci¿cally, let where Pr{b(n)|x(n), b1:n−1 } and Pr{b(n)|b1:n−1 } depend on the quantization rule qn [y(n)]. If p[x(n − 1)|b1:n−1 ] is known, yˆ(0) (n|b1:n−1 ) :=ˆ y(n|b1:n−1 ) = E {y(n)|b1:n−1 } the prior pdf p[x(n)|b1:n−1 ] can be obtained as yˆ(i) (n|b1:n−1 ) :=E y(n)|b1:n−1 , b(1:i) (n) , for i ≥ 1 (12) p[x(n)|b1:n−1 ] (7) stand for MMSE estimates of y(n) using past messages b1:n−1 and the ¿rst i bits of the current message denoted as b(1:i) (n). = p[x(n)|x(n−1), b1:n−1 ]p[x(n−1)|b1:n−1]dx(n−1) The i-th bit of the current message, b(i) (n), is de¿ned as Rp where due to Gaussian and Markov properties of the random process x(n) in (3), p[x(n)|x(n−1), b1:n−1 ] = p[x(n)|x(n−1)] = N [x(n); A(n)x(n−1), Cu (n)]. If b(n) = y(n), i.e., un-quantized observations are used, both conditional pdfs in (7) and (6) are Gaussian and it suf¿ces to propagate their means and covariance matrices only. This leads to the following Kalman ¿lter recursions [6, Chapter 5.3]: [P1] KF prediction step. Given the previous estimate ˆ (n − 1|y1:n−1 ) and its ECM M(n − 1|y1:n−1 ), x ˆ (n|y1:n−1 ) the predicted estimate x := } and its ECM M(n|y E{x(n)|y1:n−1 1:n−1 ) := E [ˆ x(n|y1:n−1 ) − x(n)][ˆ x(n|y1:n−1 ) − x(n)]T are obtained as
y (i−1) (n|b1:n−1 )]. b(i) (n) := sign[y(n) − yˆ(i−1)(n|b1:n−1 )] := sign[˜ (13) Our goal is to iteratively compute the estimates ˆ (i) (n|b1:n−1 ) := E x(n)|b1:n−1 , b(1:i) (n) x (14) ˆ (i−1) (n|b1:n−1 ) and b(i) (n). When using m bits, based on x b(1:m) (n), we will refer to the resulting algorithm as m-IQKF. A possible option for de¿ning the quantizer would be to x(i) (n|b1:n−1 ). However, set yˆ(i) (n|b1:n−1 ) equal to hT (n)ˆ (i) T (i) yˆ (n|b1:n−1 ) = h (n)ˆ x (n|b1:n−1 ) if i ≥ 1 as explained next. Using the observations (4), the de¿nitions of yˆ(i) (n|b1:n−1 ) ˆ (i) (n|b1:n−1 ) in (14), we obtain in (12), and x yˆ(i) (n|b1:n−1 ) = E hT (n)x(n) + v(n)|b1:n−1 , b(1:i) (n)
x(i) (n|b1:n−1 ) + E{v(n)|b1:n−1 , b(1:i) (n)}. (15) = hT (n)ˆ ˆ (n|y1:n−1 )=A(n)ˆ x x(n−1|y1:n−1 ) (8) M(n|y1:n−1 )=A(n)M(n−1|y1:n−1 )AT(n)+Cu(n). (9) The noise estimate E{v(n)|b1:n−1 , b(1:i) (n)} is not necessarily zero - see also (22). Therefore, in orderto obtain yˆ(i) (n|b1:n−1) ˆ (i) (n|b1:n−1 ) as well as E v(n)|b1:n−1 , b(1:i) (n) [C1] KF correction step. Consider the predicted data we need x E{y(n)|y1:n−1 } := yˆ(n|y1:n−1 ) = hT (n)ˆ x(n|y1:n−1 ) and which is achieved by augmenting the state vector x(n) with the their innovation y˜(n|y1:n−1 ) := y(n) − yˆ(n|y1:n−1 ). Then noise term v(n) as described in the next section.
647
A. State augmentation state the estimates augmentation Through E v(n)|b1:n−1 , b(1:i) (n) can be obtained so that the predicted observation yˆ(i) (n|b1:n−1 ) in (15) can be evaluated. Speci¿cally, ˘ (n) := [xT (n), v(n)]T , u ˘ (n) := [uT (n), v(n)]T , let x ˘ ˘ h(n) := [hT (n), 1]T , and A(n) de¿ned with A(n) as the leading p × p submatrix and with zeros in all other entries; then model in (3)-(4) can be rewritten as ˘ (n) = x
˘ ˘ (n − 1) + u ˘ (n) A(n) x T ˘ ˘ (n) + v˘(n) h (n) x
ˆ ˘ (0) (n|b1:n−1 ) := where we used the de¿nitions x (0) ˆ ˘ ˘ ˘ x(n|b1:n−1 ) and M (n|b1:n−1 ) := M(n|b 1:n−1 ). For time index n, (20) and (21) are repeated m-times. ˆ ˘ (n) ˘ (n|b x The MMSE estimate of x given b1:n is 1:n ) = (1:m) ˘ (n)|b1:n−1 , b E {˘ x(n)|b1:n } = E x (n) = ˆ ˘ ˘ (m) (n|b1:n−1 ). The corresponding ECM is M(n|b x 1:n ) := ˘ (m) (n|b1:n−1 ). M
(16) The state estimates x ˆ (i) (n|b1:n−1 ) in (14) are the ¿rst p comˆ ˘ (i) (n|b1:n−1 ), i.e., y˘(n) = (17) ponents of the augmented state estimate x ˆ ˆ (i) (n|b1:n−1 ) = [x ˘ (i) (n|b1:n−1 )]1:p and the noise estimate is x where the new observation noise v˘(n) = 0, ∀n (by construc- the (p + 1)-st component. tion). Note that the covariance matrix of the augmented driving noise is a block-diagonal matrix Cu˘ (n) ∈ R(p+1)×(p+1) with Corollary 1 For i = 1, E v(n)|b1:n−1 , b(1:i) (n) = 0. [Cu˘ (n)]1:p,1:p = Cu (n) and [Cu˘ (n)]p+1,p+1 = cv (n). The augmented state formulation (16)-(17) increases the diˆ mension of the state vector but is otherwise equivalent to (3)-(4). ˘ (i) (n|b Proof: The (p+1)-st entry of x 1:n−1 ) in (20) is the However, it has the appealing property that MMSE estimates noise estimate E v(n)|b1:n−1 , b(1:i) (n) . Thus, for i = 1 ˘ (n) contain MMSE estimates of the of the augmented state x original state x(n) and of the observation noise v(n) which E{v(n)|b1:n−1 , b(n)}
allows computation of yˆ(i) (n|b1:n−1 ) in (12) as follows. Since ˘ ˘ (0) (n|b1:n−1 )h(n)] 2 [M p+1 T T ˘ b(n) = E{v(n)|b1:n−1 } + y˘(n) = h (n)˘ x(n) = h (n)x(n) + v(n) = y(n), it follows π (i) (i) T (i) T (0) ˘ ˘ ˘ ˆ ˆ ˘ (n) M (n|b ) h(n) h ˘ ) = h (n) x (n|b ) that yˆ (n|b1:n−1 ) = y˘ (n|b1:n−1 1:n−1 1:n−1
ˆ˘ (i) (n|b1:n−1 ) = E x ˘ (n)|b1:n−1 , b(1:i) (n) and where x cv (n) 2 b(n) = 0. (22) = y˘ˆ(i) (n|b1:n−1 ) = E y˘(n)|b1:n−1 , b(1:i) (n) . Using the augπ hT (n)M(n|b1:n−1 )h(n) + cv (n) mented state model, we obtain the MMSE estimates in (14) using the algorithm detailed in Proposition 1. The last equality follows since E {v(n)|b1:n−1 } = E{v(n)} = 0 ˘ ˘ (0) (n|b1:n−1 )h(n) and M = cv (n). Proposition 1 Consider the augmented state space p+1 model in (16)-(17). De¿ne the augmented state The similarity of the m-IQKF in Proposition 1 and the clairˆ˘ (n|b1:n−1 ) estimates x := E{˘ x(n)|b1:n−1 } and voyant KF based on un-quantized observations (cf. [P1]-[C1]) ˆ˘ (i−1) (n|b1:n−1 ) ˘ (n)|b1:n−1 , b(1:i−1) (n) . Let is quite remarkable. The ECM updates in (21) are identical to := E x x ˆ ˆ˘ (n|b1:n−1 ) − x ˘ ˘ (n|b1:n−1 ) − x ˘ (n)][x ˘ (n)]T } the ECM updates of the KF in (11) except for the scale factor M(n|b1:n−1 ) := E{[x (i−1) (i−1) ˆ ˘ ˘ (n|b1:n−1 ) := E{[x (n|b1:n−1 ) − 2/π. The variance penalties associated with the m-IQKF are and M ˆ˘ (i−1) (n|b1:n−1 ) − x ˘ (n)]T } denote the corresponding investigated in the following corollaries. ˘ (n)][x x ECMs for i = 1, . . . , m. Construct the messages b(1:m) (n) ˆ ˘ T (n)x ˘ (i−1) (n|b1:n−1 ) . Corollary 2 Consider the m-IQKF algorithm in Proposition 1 in (13) as b(i) (n) := sign y(n) − h ˆ˘ (i) (n|b1:n−1 ) is obtained from the recursion: and de¿ne the ECM reduction after i = 1, . . . , m iteraThe estimate x ˘ (0) (n|b1:n−1 ) − M ˘ (i) (n|b1:n−1 ), where ˘ i (n) := M tions as ΔM ˆ˘ (n − 1|b1:n−1 ) and its ECM ˘ (0) 2 i ˘ [P2] Given the previous estimate x M (n|b1:n−1 ) := M(n|b ). With c 1:n−1 i := 1 − (1 − π ) , the ˘ M(n−1|b ), form 1:n−1 error covariance reduction after i iterations is given as ˘ ˘ˆ (n|b1:n−1 )=A(n) ˘ˆ (n−1|b1:n−1) x x (18) ˘ ˘T ˘ (0) ˘ (0) ˘ i (n) = ci M (n|b1:n−1 )h(n)h (n)M (n|b1:n−1 ) . ˘ ˘ ˘ ˘T ΔM M(n|b 1:n−1 )=A(n)M(n−1|b1:n−1 )A (n)+Cu ˘ (n)(19) ˘ ˘ T (n)M ˘ (0) (n|b1:n−1 )h(n) h = [C2] Assuming that p[˘ x(n)|b1:n−1 , b(1:i−1) (n)] (23) ˆ˘ (i−1) (n|b1:n−1 ), M ˘ (i−1) (n|b1:n−1 )], the MMSE N [˘ x(n); x ˆ ˘ i (n) as a summation of differences ˘ (i) (n|b1:n−1 ) := E{˘ estimate x x(n)|b1:n−1 , b(1:i) (n)} Proof: We ¿rst write ΔM ˘ (i) (n|b1:n ) are obtained between successive ECM matrices and the corresponding ECM M iteratively from ˘ (0) (n|b1:n−1 ) − M ˘ (i) (n|b1:n−1 ) ˘ i (n) := M ΔM ˆ ˆ˘ (i) (n|b1:n−1 ) = x ˘ (i−1) (n|b1:n−1 ) + x i
˘ (i−1) (n|b1:n−1 )h(n) ˘ M 2 ˘ (j−1) (n|b1:n−1 ) − M ˘ (j) (n|b1:n−1 ) = M b(i)(n) (20) π ˘T j=1 ˘ ˘ (i−1) (n|b1:n−1 )h(n) h (n)M i ˘ h ˘ T (n)M ˘ (j−1) (n|b1:n−1 )h(n) ˘ (j−1) (n|b1:n−1 ) 2 M ˘ (i−1) (n|b1:n−1 ) − ˘ (i) (n|b1:n−1 ) = M M = T (j−1) ˘ ˘ ˘ π j=1 h (n)M (n|b1:n−1 )h(n) ˘ (i−1) (n|b1:n−1 )h(n) ˘ h ˘ T (n)M ˘ (i−1) (n|b1:n−1 ) 2M (24) ˘ ˘ T (n)M ˘ (i−1) (n|b1:n−1 )h(n) π h (21) where the last equality follows from (21). Next, we multiply
648
TABLE I P ER STEP FACTOR , cm , AND NOISE PENALTY FACTOR FOR IQKF
˘ ˘ (i) (n|b1:n−1 ) from (21) by h(n) M to obtain ˘ ˘ (i−1) (n|b1:n−1 )h(n) ˘ ˘ (i) (n|b1:n−1 )h(n) =M − M (i−1) T (i−1) ˘ ˘ ˘ ˘ ˘ (n|b1:n−1 )h(n)h (n)M (n|b1:n−1 )h(n) 2M ˘ ˘ T (n)M ˘ (i−1) (n|b1:n−1 )h(n) h 2 ˘ (i−1) ˘ M (n|b1:n−1 )h(n) . = 1− π
π
bits, m cm (1/cm − 1)100%
(25)
1 0.637=2/π 57.08%
2 0.868 15.21%
3 0.952 5.04%
4 0.983 1.77%
(26) B. Performance analysis of the m-IQKF In the previous section, we quanti¿ed the per correction step ˘ m (n) in the m-IQKF as a function of the ˘ ˘ (i−1) (n|b1:n−1 )h(n) ECM reduction ΔM Repeating steps (25)-(26) for M ˘ ˘ (i) (n|b1:n−1 )h(n) with decreasing index i yields M = number of bits m used for iteratively quantizing the observations i (0) ˘ ˘ (n|b1:n−1 )h(n) which when substituted in (24) y(n). We next compare MSE performance of m-IQKF with that 1− π2 M of the clairvoyant KF, when both prediction and correction steps with i = j − 1, results in are considered, by deriving the continuous-time Algebraic Riccati ˘ i (n) = ΔM Equations (ARE) for both ¿lters. j−1 ˘ (0) i Consider ¿rst the discrete time ARE for the m-IQKF [P2]T (0) ˘ ˘
˘ M (n|b1:n−1 )h(n)h (n)M (n|b1:n−1 ) 2 2 ˘ ˘ + 1) := M(n + 1|b1:n ), . [C2]. To simplify notation let M(n 1− ˘ ˘ T (n)M ˘ (0) (n|b1:n−1 )h(n) π j=1 π h ˘ ˘ M(n) := M(n|b1:n−1 ), M(n + 1) := M(n + 1|b1:n ), M(n) := M(n|b1:n−1 ), and substitute (27) in (19) to obtain the m-IQKF Upon invoking the geometric sum identity ARE for the ECM of x ˘ (n) as i 2 2 j−1 = 1 − (1 − π2 )i , (23) follows. j=1 1 − π π ˘ ˘ ˘ ˘T (29) From the KF’s ECM in (11), it follows readily that the ECM re- M(n + 1) =A(n + 1)M(n)A (n + 1) + Cu˘ (n + 1) T T ˘ ˘ ˘ ˘ ˘ ˘ duction for un-quantized observations M(n|b1:n−1 )−M(n|b1:n ) A(n + 1)M(n)h(n)h (n)M(n)A (n + 1) − cm . corresponds with cm = 1 in (23). Therefore, the ECM reduction ˘ ˘T (n)M(n) ˘ h(n) h achieved by quantizing to m bits is cm times smaller than the one ˘ ˘ ˘ + 1), Cu˘ (n + 1), h(n) in (29), and M(n) achieved with analog-amplitude observations. Values of cm are Substituting for A(n shown in Table I. With only 4 bits of quantization, the value of and writing only the leading p × p sub-matrices of the resulting cm is just 2% less than the value for the clairvoyant KF (cm = 1). expression, yields the ARE for the ECM of x(n) as From Corollary 2 we can relate the predicted MSE M(n + 1) =A(n + 1)M(n)AT(n + 1) + C (n + 1) (30) u ˘ ˘ tr{M(n|b 1:n−1 )} and the corrected MSE tr{M(n|b1:n )} afT T A(n + 1)M(n)h(n)h (n)M(n)A (n + 1) (1:m) (n), i.e., ter observing the m-bits of n-th observation, b − cm . T (n)M(n)h(n) + c (n) h ˘ ˘ ˘ v M(n|b1:n ) = M(n|b1:n−1 ) − ΔMm (n). Substituting for ˘ m (n) we obtain Interestingly, the resulting ARE for the ECM of the m-IQKF ΔM [cf. (30)] becomes identical to the ARE of the clairvoyant KF ˘ ˘ M(n|b 1:n ) = M(n|b1:n−1 )− [6, Chapter 5.11] as limm→∞ cm = 1. The implication of these ˘ ˘T ˘ ˘ relations follows from the continuous-time ARE for the m-IQKF M(n|b 1:n−1 )h(n)h (n)M(n|b1:n−1 ) − cm (27) given in Proposition 2. T ˘ ˘ ˘ h (n)M(n|b1:n−1 )h(n) Equation (27) is instructive for understanding the MSE perfor- Proposition 2 The continuous-time ECM Mc (t) = Mc (nTs ) := mance of m-IQKF since for cm = 1, (27) coincides with the cor- limTs →0 M(n; Ts ) for the m-IQKF of (18)-(21) is the solution rection ECM of the KF in (11). Furthermore, limm→∞ cm → 1, of the following differential (Riccati) equation: implying that for in¿nite number of quantization bits we recover ˙ c (t) =Ac (t)Mc (t) + Mc (t)AT (t) + Cu (t) M c the clairvoyant KF. This important consistency result is summa −1 cvc (t) rized in the following corollary. − Mc (t)hc (t) hTc (t)Mc (t) . (31) cm Proof: The derivation steps from (30) to (31) are the same Corollary 3 As the number of quantization bits m → ∞, the as those followed for deriving the KF Riccati equation from the correction step ECM at time n of the m-IQKF converges to the corresponding ECMs [6, pp. 259]. The only difference is that the ECM of the clairvoyant KF given that b1:n−1 = y1:n−1 , i.e., scale factor, cm , in (30) equals 1 in the case of the KF. ˘ ˘ (m) (n|b1:n−1 ) = M(n|y (28) lim M From the clairvoyant KF ARE [6, pp. 259] 1:n ) . m→∞
Proof: Since limm→∞ cm = limm→∞ 1 − (1 − π2 )m = 1, ˘ (m) (n|b1:n−1 ) [cf. (27)] and M(n|y ˘ then for cm = 1, M 1:n ) [cf. (11)] become identical, which results in (28). Corollary 3 establishes that m-IQKF asymptotically (as the number of quantization bits increases) achieves the per correction step MSE performance of the clairvoyant KF. As demonstrated by simulations in Section IV, even a small number of quantization bits (m = 2 or m = 3) renders the MSE performance of m-IQKF indistinguishable from that of the clairvoyant KF.
˙ K (t) =Ac (t)MK (t) + MK (t)AT (t) + Cu (t) M c c c c
−1 T hc (t)MK c (t)
− MK c (t)hc (t) [cvc (t)] MK c (t)
(32)
denotes the continuous-time ECM of KF, it is where evident that the variance of the continuous-time observation noise of the m-IQKF is ampli¿ed by 1/cm compared to that of the clairvoyant KF. This is the price paid for using quantized observations, b(n), instead of the analog-amplitude ones, y(n). In Table I the percentage observation noise variance increase (1/cm −1)100% versus number of quantization bits, m is shown.
649
IV. S IMULATIONS
distributed estimation, it was shown that quantization using 2 to 3 bits improves the performance of the Kalman-like algorithm to virtually coincide with the optimal state estimator (Kalman ¿lter), with only minimal increase in computation. This result was corroborated by simulation results. It was further shown that true ¿lter covariances are consistent with Noisy measurements at sensor Sk are given as yk (t) = xc1 (t) + through simulations 1 . analytical ones θk xc2 (t) + vk (t) where θk is a parameter at sensor k. A discreteR EFERENCES time equivalent model is used with sampling time Ts = 1. [1] Y. Bar-Shalom and X. Li, Estimation and Tracking: Principles, Techniques, WSN data corresponding to K = 2 sensors with [θ1 θ2 ] = and Software. Artech House, 1993. [0.1 0.2] are simulated. The measurement and state driving noise [2] P. Djuric, J. Kotecha, J. Zhang, Y. Huang, T. Ghirmai, M. Bugallo, and J. Miguez, “Particle ¿ltering,” IEEE Sig. Proc. Mag., vol. 20, pp. 19–38, variances are cv (t) = 1 and cu (t) = 1. The MSE plots in Fig. 1 Sep. 2003. and Fig. 2 illustrate the evolution of the MSEs, obtained from [3] R. M. Gray, “Quantization in task-driven sensing and distributed processing,” the trace of the respective ECMs, against the time index n. in Proc. Int. Conf. Acoustics, Speech, Signal Processing, vol. 5, pp. V–1049– V–1052, Toulouse, France, May 14-19, 2006. In the ¿rst simulation setup, Fig. 1, the MSE of the IQKF, given [4] S. Julier and J. Uhlmann, “Unscented ¿ltering and nonlinear estimation,” by tr{M(n|b1:n )}, is compared to the MSE of the clairvoyant Proc. IEEE, vol. 92, pp. 401–422, Mar. 2004. KF. The simulations for IQKF are performed for 1, 2 and 3 bits. [5] J. Kotecha and P. Djuric, “Gaussian particle ¿ltering,” IEEE Trans. on Sig. Processing, vol. 51, pp. 2602–2612, Oct. 2003. The plots demonstrate that there is a substantial MSE reduction [6] P. S. Maybeck, Stochastic Models, Estimation and Control – Vol.1. Academic when going from 1 bit of quantization to 2 bits and little MSE Press, 1979. reduction for higher number of quantization bits. With 2 bits of [7] H. Papadopoulos, G. Wornell, and A. Oppenheim, “Sequential signal encoding from noisy measurements using quantizers with dynamic bias control,” quantization the MSE of IQKF is virtually coincident with that IEEE Trans. on Infor. Theory, vol. 47, pp. 978–1002, 2001. of the clairvoyant KF as was postulated by the analytical values [8] A. Ribeiro, G. B. Giannakis, and S. I. Roumeliotis, “SOI-KF: Distributed in Table I whereby c2 = 86.8%. Kalman ¿ltering with low-cost communications using the sign of innovations,” IEEE Trans. on Sig. Proc., vol. 54, no. 12, pp. 4782–4795, Dec 2006. Model consistency checks comparing the empirically obtained MSEs with analytical MSEs are shown in Fig. 2. Note that KF 1íbit IQKF the analytical MSE is obtained from the trace of the ECM 2íbit IQKF 3íbit IQKF ˆ (n|b1:n )][x(n) − x ˆ (n|b1:n )]T |b1:n }. M(n|b1:n ) := E{[x(n) − x The empirical MSE is the sample estimator of tr{M(n|b1:n )} obtained as a sample average of the squared estimation errors. The plots depict both predictor MSE tr{M(n|b1:n−1 )} and estimator MSE tr{M(n|b1:n )}. The consistency check reveals that the empirical and analytical MSEs are nearly identical. Fig. 3 shows alternative model consistency tests for time index, n the IQKF using the normalized estimation error squared Fig. 1. IQKF estimator MSE tr{M(n|b1:n )} from Monte Carlo data for 1-3 (NEES) tests of [1, Ch. 5.4]. NEES r(n) := [x(n) − bits of quantization; MSE for clairvoyant KF is included for comparison. ˆ (n|b1:n )]T [M(n|b1:n )]−1 [x(n) − x ˆ (n|b1:n )] is postulated to x IQKF predictor MSEí analytical IQKF predictor MSEí empirical have a χ2 pdf with p degrees of freedom (since the p-dimensional IQKF estimator MSEí analytical IQKF estimator MSE í empirical ˜ (n) := x(n)− x ˆ (n|b1:n ) is assumed zero-mean Gaussian with x covariance M(n|b1:n ) if the estimator is consistent with the model). Under the hypothesis that the estimator is consistent, L realizations of the NEES statistics {ri (n)}L i=1 each χ2 distributed with p degrees of freedom, lead to a χ2 distribution with Lp degrees of freedom. This is checked by running Monte Carlo simulations and computing the sample average NEES rˆ(n) := L time index, n 1 r (n) and then de¿ning an acceptance (con¿dence) Fig. 2. Comparing empirical MSE from Monte Carlo data and analytical MSE, i i=1 L region (for the consistent hypothesis). If l ≤ rˆ(n) ≤ u, then for predictor ECM M(n|b1:n−1 ) and estimator ECM M(n|b1:n ). the estimator is consistent; lower and upper bounds l and u are obtained from Pr{ˆ rn ∈ [l, u]} = 1 − α, where 1 − α is the probability of acceptance region. Using L = 50 realizations, 200 time samples, a state space of p = 2 dimensions, and α = 0.05 (i.e., 95% region), we observed that indeed only 7% of the 100 time samples lie outside the 95% acceptance region. The two-dimensional state space model simulated is x˙ c1 (t) 0 1 0 xc1 (t) x˙ c (t) := = + uc (t) . x˙ c2 (t) xc2 (t) 0 0 1
5.5
meanísquare error (MSE)
5
4.5 4
3.5 3
2.5 2
1.5 1
5
10
15
20
25
30
35
40
45
50
20
25
30
35
40
45
50
12
meanísquare error (MSE)
10
8
6
4
2
0
5
10
15
Normalized Est. Error Squared, NEES
4
3.5
NEES, r(n)
3
2.5 2
1.5 1
V. C ONCLUDING REMARKS
0.5
Recursive state estimation based on quantized observations was considered. Multi-bit quantization was done using iterative binary quantizing of the measurement innovations. Joint quantization and estimation was used to develop a Kalman-like algorithm. Motivated by the need to ¿nd quanti¿able trade-offs between estimation performance and number of quantization bits for
0
50
100
150
200
time index, n
Fig. 3. Consistency test using normalized estimation error squared (NEES), r(n), compared with the 95% probability region using 2 bits of quantization. 1 The views and conclusions contained in this document are those of the authors and should not be interpreted as representing the of¿cial policies of the Army Research Laboratory or the U. S. Government.
650
