无线传感器网络分布式量化卡尔曼曼滤滤波

第 28卷第 12期 2011年 12月 控 制 理 论 与 应 用 Control Theory & Applications Vol. 28 No. 12 Dec. 2011 无无无线线线传传传感感感器器器网网网络络络分分分布布布式式式量量量化化化卡卡卡尔尔尔曼曼曼滤滤滤波波波 文文文章章章编编编号号号: 1000−8152(2011)12−1729−11 陈军勇1, 邬依林1,2, 祁 恬1 (1. 华南理工大学自动化科学与工程学院,广东广州 510640; 2. 广东第二师范学院计算机科学系,广东广州 510310) 摘要: 本文针对无线传感器网络中的目标跟踪问, 研究了分布式量化卡尔曼滤波问题. 由于网络中存在能 量和带宽限制, 传感器传输的数据必须经过量化处理. 考虑一个线性离散随机动态系统, 首先提出了一种动 态Lloyd-Max量化器并设计了其在线更新,然后基于贝叶斯原理导出了递归形式的最优量化卡尔曼滤波器,同 时给出了一种渐近等价的迭代算法,并进一步分析了量化卡尔曼滤波器的稳定性. 最后,仿真结果验证了所设计算 法的可行性与有效性. 关键词: 无线传感器网络;分布式量化卡尔曼滤波;动态Lloyd-Max量化器;稳定性 中图分类号: TP393 文献标识码: A Distributed quantized Kalman filtering for wireless sensor networks CHEN Jun-yong1, WU Yi-lin1,2, QI Tian1 (1. College of Automation Science and Technology, South China University of Technology, Guangzhou Guangdong 510640, China; 2. Faculty of Computer Science, Guangdong University of Education, Guangzhou Guangdong 510310, China) Abstract: We study the distributed quantized Kalman filtering for the target-tracking in wireless sensor networks (WSNs). Because of the constraints on power and bandwidth in WSNs, sensor data have to be quantized before transmis- sion. A linear discrete-time stochastic dynamic system is employed for this purpose. First, a dynamic Lloyd-Max quantizer is adopted and the corresponding online update scheme is designed. Then, the optimal recursive quantized Kalman filter is derived based on the Bayesian principles, and an asymptotically equivalent iterative algorithm is developed. The sta- bility of the quantized Kalman filter is analyzed. Simulation results show the feasibility and effectiveness of the designed algorithms. Key words: wireless sensor networks; distributed quantized Kalman filtering; dynamic Lloyd-Max quantizer; stability 1 引引引言言言(Introduction) 无线传感器网络(WSNs)近年来成为研究热点, 其由大量低功耗的分布式传感器节点组成, 通过节 点之间的协同式工作可以完成许多复杂的任务,如 目标定位与跟踪、网络控制、环境监测等[1∼5]. 带宽 受限下的分布式估计问题是无线传感器网络研究 的一个重要方向,由于传感器节点能量和带宽的限 制,为了延长网络使用寿命,节点一般需要对观测进 行量化编码后才发送, 接收端对量化信号解码并进 行融合处理. 因此, 网络中分布式估计问题的量化 器、估计器以及节点通信设计等是研究重点. 多年来,国内外一些研究者提出了多种能量和带 宽受限下的分布式估计算法. Gubner[1]针对随机性 参数的分布式量化估计问题,设计了给定线性形式 下的最优估计器,并得到了局部最优量化器. Zhang 等[3]基于文献[1]的思想首先研究了静态参数的线性 最优分布式估计问题,然后扩展到动态系统的量化 卡尔曼滤波问题,其设计了递归形式的估计器和动 态更新的量化器, 但文中只讨论了单传感器观测情 形, 且量化器的更新需要融合中心在每个采样周期 反馈状态预测信息到节点, 这增加了通信代价并带 来一定局限性. Ribeiro等[6]研究了量化比特率仅为 一位的量化卡尔曼滤波器(SOI–KF),然后Msechu等 在文献 [7]中扩展到了多位情形,他们虽然考虑了多 传感器网络环境, 但在每个采样周期只有一个传感 器激活, 且量化器的更新亦需要融合中心广播状态 预测信息到各节点. 此外, Sun等[8]研究了线性量化 卡尔曼滤波问题,其采用均匀量化器,但量化器的更 新亦需要融合中心广播状态预测信息到各节点. 针对上述问题,本文研究了多传感器多位的分布 式量化卡尔曼滤波问题.首先提出了一种新的动态 Lloyd–Max量化器并设计了其在线自动更新方案, 然后导出了递归形式的最优量化卡尔曼滤波器, 同 时给出了一种渐近等价的迭代算法(迭代量化卡尔 曼滤波器(IQKF)),最后分析了量化卡尔曼滤波器的 稳定性. 本文提出的动态Lloyd-Max量化器的更新无 收稿日期: 2010−12−11;收修改稿日期: 2011−07−05. 基金项目: 国家自然科学基金重点资助项目(60834003);国家自然科学基金资助项目(60774057);国家科技部“973”计划资助项目 (2010CB731802). 1730 控 制 理 论 与 应 用 第 28卷 需融合中心的反馈信息,这节省了通信代价,且适用 性增强. 迭代量化卡尔曼滤波器算法在存在丢包或 延迟的网络环境下亦适用,这增强了算法的鲁棒性. 对于不稳定系统,量化卡尔曼滤波器稳定的临界比 特率在给定量化器下由系统矩阵的不稳定特征值决 定,这与文献[9, 10]关于量化控制问题的结论保持一 致.最后,仿真结果验证了所得结论的正确性和算法 的高效性. 2 问问问题题题描描描述述述(Problem formulation) 考虑一个存在融合中心的传感器网络,其中包含 有N个分布式传感器节点,且节点之间没有通信,它 们各自只与融合中心通信. 在每个采样周期,各节点对同一个动态目标进行 独立观测, 由于网络中的带宽限制以及数字信号传 输的要求, 传感器的观测需要经过量化并编码成数 字信号,然后经过无噪声信道发送到融合中心,融合 中心基于接收到的量化观测对目标状态进行实时估 计.无线传感器网络分布式估计系统如图1所示. 图 1 无线传感器网络分布式估计示意图 Fig. 1 Distributed estimation scheme for WSNs 上述目标的动态模型由下列线性离散随机动态 系统描述: x(k) = Ax(k − 1) + u(k), k = 1, 2, · · · ,∞. (1) 传感器n的观测为 yn(k) = hx(k) + vn(k), n = 1, 2, · · · , N, (2) 其中: 状态变量x(k) ∈ Rp×1, 初始状态x(0) ∼ N (x0,M0), 转移矩阵A ∈ Rp×p, 驱动噪声u(k) ∈ Rp×1为零均值高斯白噪声, 其协方差阵E{u(k) · uT(l)} = Cuδkl, Cu ∈ Rp×p,传感器观测yn(k) ∈ R, 其中h ∈ R1×p为观测向量,观测噪声vn(k) ∈ R为零 均值高斯白噪声, 其方差E{vn(k)vTm(l)} = σ2nδnm · δkl,假定{u(k), vn(k), x(0)}相互独立. 传感器n的量 化比特率给定为qn位,其中qn > 0, qn ∈ Z, qn = 0表 示传感器未激活. 2.1 动动动态态态量量量化化化器器器(Dynamic quantizer) 本节首先介绍标量量化器的概念,然后针对本文 的动态目标跟踪问题，提出动态量化器的设计及其 在线更新问题[3, 11, 12]. 量化是一个把连续信号转换成离散信号的处 理过程. 考虑一维标量信源y ∈ R, 设量化比特率 为q位, q > 0, q ∈ Z,定义量化分区Ri := [τi, τi+1), i = 0, 1, · · · , 2q − 1, 其中{τi}2qi=0为量化阈值集, 且τ0 = −∞, τ2q =∞. 量化准则如图2所示,即 b = i, iff y ∈ Ri, (3) 其中: b为量化输出, i ∈ {0, 1, · · · , 2q − 1}, 其需要 转换为二进制数字信号然后才发送. 图 2 量化器 Fig. 2 Quantizer 在接收端,解码器根据量化规则对数字信号进行 解码,将信息从数字信号的形式恢复到其原来形式, 解码规则如下: y ∈ Ri or y = yˆi, iff b = i, (4) 其中: 量化值yˆi ∈ Ri, i = 0, 1, · · · , 2q − 1, yˆi的具体 值根据量化准则而定. 对于一个具有给定分布的随机信源y,设其概率 密度函数为p(y),定义平均失真 D := 2q−1∑ i=0 w τi+1 τi (y − yˆi)2p(y)dy, (5) 则存在一种最优量化器使其平均失真达到最小, 即著名的Lloyd-Max量化器[11, 12], 其在给定量化比 特率q下, 量化分区可由信源分布p(y)唯一确定, 即 {τi}2q−1i=1 = F LM(p(y)). 由于本文考虑的是动态目标的跟踪问题,在不同 的采样时刻,传感器的观测信号yn(k)的分布不一样, 为了得到最优的量化效果,本文采用动态Lloyd-Max 量化方式, 即在每个采样周期, 各传感器的Lloyd- Max量化分区根据采样信号的分布进行动态更新. 假定各传感器以及融合中心存有系统模型的信息, 因而关于动态量化器的一个基本问题是: 传感器的 编码器和融合中心的解码器如何根据量化观测和系 统模型等信息实现在线自动更新? 2.2 基基基于于于量量量化化化数数数据据据的的的极极极小小小均均均方方方误误误差差差估估估计计计 (Minimum mean square error estimation based on quantized data) 本节根据贝叶斯估计原理提出基于量化数据的 最优卡尔曼滤波问题. 设bn(k)为传感器n在k时刻对观测yn(k)的量化 输出, 其中bn(k) ∈ {0, 1, · · · , 2qn − 1}. 定义k时刻 所有传感器的量化观测集为Bk := {b1(k), b2(k), · · · , bN(k)},整个传感器网络从初始观测时刻1到时 刻k的量化观测集为B1:k := {B1, B2, · · · , Bk}. 本 文用xˆ(k|B1:k)表示基于所有接收到的量化观测B1:k 第 12期 陈军勇等: 无线传感器网络分布式量化卡尔曼滤波 1731 得到的对当前状态x(k)的估计. 分布式量化估计结 构如图3所示,其中{Sn}Nn=1表示传感器集合. 图 3 分布式量化估计结构图 Fig. 3 Distributed quantized-estimation scheme 本文采用均方误差(MSE)作为估计性能准则,定 义k时刻估计器的误差协方差阵为 M(k|B1:k) := E{[x(k)−xˆ(k|B1:k)][x(k)−xˆ(k|B1:k)]T|B1:k−1}, (6) 则xˆ(k|B1:k)的均方误差(MSE)可用该误差协方差阵 的迹表示,即MSE := tr{M(k|B1:k)}. 根据贝叶斯估计原理[13],融合中心在k时刻基于 B1:k对状态x(k)的极小均方误差估计(MMSE)即为 下面的条件均值: xˆ(k|B1:k) = E{x(k)|B1:k} =w Rp x(k)p(x(k)|B1:k)dx(k), (7) 其中p(·|·)表示条件概率密度函数. 在任意k时刻,状态的后验分布p(x(k)|B1:k)可以 根据系统模型(1)和(2)以及B1:k递归确定. 基于量化 观测的卡尔曼滤波器([Pq]–[Cq])可以表示为: [Pq]预测. 假设k时刻p(x(k − 1)|B1:k−1)已知,则 状态x(k)在k时刻的先验分布p(x(k)|B1:k−1)为 p(x(k)|B1:k−1) =w Rp p(x(k)|x(k − 1))p(x(k − 1)|B1:k−1)dx(k − 1), (8) 其中p(x(k)|x(k − 1)) = N [x(k);Ax(k − 1), Cu]. [Cq]校正. 当融合中心接收到新的量化观测Bk 时,根据贝叶斯原理[13]，即有 p(x(k)|B1:k)=p(x(k)|B1:k−1)P{Bk|x(k), B1:k−1}P{Bk|B1:k−1} . (9) 由于量化引入的非线性影响, 基于量化观测的 贝叶斯最优估计算法计算复杂, 因此对状态后验分 布p(x(k)|B1:k)作高斯近似,即 假假假设设设 1 在k时刻,基于量化观测B1:k的状态后 验分布p(x(k)|B1:k)假定为高斯分布,即 p(x(k)|B1:k) = N [x(k); xˆ(k|B1:k),M(k|B1:k)]. (10) 注注注 1 类似的假设见文献 [7, 14],当量化比特率比较 大的时候,模拟结果显示其基本成立. 量化卡尔曼滤波器的预测部分[Pq]和校正部分 [Cq]是本文研究的核心问题, 即融合中心如何根据 历史量化观测对当前状态进行预测? 当获得当前的 量化观测后,如何对状态预测进行校正? 进一步,当 得到状态的最优估计后, 如何分析量化卡尔曼滤波 器的稳定性与量化器以及量化比特率之间的关系? 3 动动动态态态Lloyd-Max量量量化化化器器器(Dynamic Lloyd- Max quantizer) 本节首先简介Lloyd-Max量化器的一些基本特 征[11, 12], 然后提出一种动态Lloyd-Max量化器更新 方案. 3.1 Lloyd-Max量量量化化化器器器(Lloyd-Max quantizer) 若信源 y服从高斯分布, 设 y ∼ N (µy, σ2y), 则 Lloyd-Max量化分区可表示为 {τi}2q−1i=1 = F LM(µy, σ2y), (11) 其具体实现见文献[12],且其平均失真可表示为 D = K(q)σ2y 4q , (12) 其中K(q)为一随q缓慢递增的收敛序列, 1 6 K(q) 6 √ 3pi 2 ,由于其值变化不大,亦可近似视其为常数. 3.2 动动动态态态Lloyd-Max量量量化化化器器器及及及其其其更更更新新新(Dynamic Lloyd-Max quantizer and update) 下面讨论动态Lloyd-Max量化器的具体设计问 题. 设{τ (n)i (k)}2 qn i=0表示传感器n在k时刻的量化阈 值集,其中τ (n)0 (k)=−∞, τ (n)2qn (k)=∞. 记R(n)i (k) := [τ (n)i (k), τ (n) i+1(k)),则动态量化规则如下: bn(k) = i, iff yn(k) ∈ R(n)i (k), (13) 其中量化输出i ∈ {0, 1, · · · , 2qn − 1}. 在解码端, 解码器n根据接收到的bn(k)进行解 码,解码规则如下: yn(k) ∈ R(n)i (k), iff bn(k) = i. (14) 为了得到动态更新的Lloyd-Max量化器, 采用局 部量化卡尔曼滤波器的方式. 定义b(n)1:k := {bn(1), bn(2), · · · , bn(k)}, 则传感器n在k时刻基于局部量 化观测集b(n)1:k得到的最优状态估计为 xˆ(k|b(n)1:k) := E(x(k)|b(n)1:k) =w Rp x(k)p(x(k)|b(n)1:k)dx(k), (15) 且局部估计误差协方差阵为 1732 控 制 理 论 与 应 用 第 28卷 M(k|b(n)1:k) := E{[x(k)− xˆ(k|b(n)1:k)][x(k)− xˆ(k|b(n)1:k)]T|b(n)1:k−1}. (16) 根据贝叶斯原理,传感器n的局部量化卡尔曼滤 波器([Plq]–[Clq])可以表述如下: [Plq]预测.在k时刻,假定xˆ(k−1|b(n)1:k−1)和M(k− 1|b(n)1:k−1)已经得到, 同假设1类似, 对状态后验分布 p(x(k − 1)|b(n)1:k−1)作高斯近似,即 p(x(k − 1)|b(n)1:k−1) = N [x(k − 1); xˆ(k − 1|b(n)1:k−1),M(k − 1|b(n)1:k−1)]. (17) 根据系统模型(1)则有 p(x(k)|b(n)1:k−1) = N [x(k); xˆ(k|b(n)1:k−1),M(k|b(n)1:k−1)], (18) 其中: xˆ(k|b(n)1:k−1) = Axˆ(k − 1|b(n)1:k−1), M(k|b(n)1:k−1) = AM(k − 1|b(n)1:k−1)AT + Cu. (19) [Clq]校正. 当传感器的量化观测为bn(k)得到后, 有 p(x(k)|b(n)1:k) =N [x(k); xˆ(k|b(n)1:k−1),M(k|b(n)1:k−1)]· P{bn(k)|x(k), b(n)1:k−1} P{bn(k)|b(n)1:k−1} . (20) 传感器n的动态Lloyd-Max量化器更新规则如下 所述: 在系统初始时刻 k = 0时, 由于已知x(0) ∼ N (x0,M0),则根据系统模型(1)和(2),对下一时刻的 观测yn(1)有先验预测 p(yn(1)) = N [yn(1);µyn(1), σ2yn(1)], (21) 其中: µyn(1)=hAx0, σ2yn(1)=h(AM0AT+Cu)hT+ σ2n. 因此, 在编码端, 传感器n可根据式(11)和式(21) 初始化编码器n在观测时刻k = 1的Lloyd-Max量化 分区,即 {τ (n)i (1)}2 qn−1 i=1 = F LM n (µyn(1), σ 2 yn(1) ), (22) 其中F LMn (·)表示传感器n的Lloyd-Max量化分区生 成函数. 在解码端,融合中心亦可根据式(22)初始化 对应的解码器n. 当k > 1时,在编码端,传感器n首先对观测yn(k) 量化得到bn(k)并发送出去, 然后根据bn(k)运行局 部量化卡尔曼滤波的校正算法[Clq]得到xˆ(k|b(n)1:k)和 M(k|b(n)1:k),其具体校正算法将在下一节详细讨论. 当传感器获得xˆ(k|b(n)1:k)和M(k|b(n)1:k)后, 根据局 部量化卡尔曼滤波器的预测[Plq]即可得 p(x(k + 1)|b(n)1:k) = N [x(k + 1); xˆ(k + 1|b(n)1:k),M(k + 1|b(n)1:k)], (23) 其中: xˆ(k + 1|b(n)1:k) = Axˆ(k|b(n)1:k), M(k + 1|b(n)1:k) = AM(k|b(n)1:k)AT + Cu. (24) 又根据式(23)以及观测模型(2)即有 p(yn(k + 1)|b(n)1:k) = N [yn(k + 1);hxˆ(k + 1|b(n)1:k), hM(k + 1|b(n)1:k)hT + σ2n]. (25) 因此, 传感器n可根据式(11)和式(25)对编码器n 进行更新, 确定其下一个时刻k + 1的Lloyd-Max量 化分区,即 {τ (n)i (k + 1)}2 qn−1 i=1 = F LMn (hxˆ(k + 1|b(n)1:k), hM(k + 1|b(n)1:k)hT + σ2n). (26) 在解码端,融合中心亦可基于接收到的新的量化 观测bn(k)运行局部量化卡尔曼滤波器算法,并对解 码器n进行更新,为下一时刻的量化作准备,具体步 骤同编码器一致. 传感器n的动态Lloyd-Max量化器更新过程可用 图4表示. 图 4 传感器n的动态Lloyd-Max量化器更新示意图 Fig. 4 Update of dynamic Lloyd-Max quantizer for sensor-n 为了简化分析,对于局部量化卡尔曼滤波器的校 正部分[Clq]的解, 下面将结合融合中心的全局量化 卡尔曼滤波问题一起进行讨论. 4 基基基于于于量量量化化化观观观测测测的的的卡卡卡尔尔尔曼曼曼滤滤滤波波波器器器(Kalman filtering based on quantized observations) 首先回顾一下无量化效应时的卡尔曼滤波问题. 4.1 基基基于于于分分分布布布式式式观观观测测测的的的卡卡卡尔尔尔曼曼曼滤滤滤波波波器器器(Kalman filtering based on distributed observations) 假设融合中心接收到的为未量化的观测, 即Bk = {y1(k), y2(k), · · · , yN(k)},则量化卡尔曼滤波器 即简化为卡尔曼滤波器[13, 15]. 定义: Y (k) := [y1(k) y2(k) · · · yN(k)]T, 第 12期 陈军勇等: 无线传感器网络分布式量化卡尔曼滤波 1733 Y1:k := {Y (1), Y (2), · · · , Y (k)}, xˆ(k|Y1:k) := E{x(k)|Y1:k}, M(k|Y1:k) := E{[x(k)− xˆ(k|Y1:k)][x(k)− xˆ(k|Y1:k)]T}. (27) 基于分布式观测Y (k)的卡尔曼滤波器([Pkf ]– [Ckf ])可表示为: [Pkf ]预预预测测测. 假设k时刻xˆ(k − 1|Y1:k−1)和M(k − 1|Y1:k−1)已获得,则有: xˆ(k|Y1:k−1) = Axˆ(k − 1|Y1:k−1), M(k|Y1:k−1) = AM(k − 1|Y1:k−1)AT + Cu. (28) [Ckf ]校校校正正正. 记 ϕ := N∑ n=1 1 σ2n , φ := 1 ϕ [ 1 σ21 1 σ22 · · · 1 σ2N ], (29) 则基于观测Y1:k对状态x(k)的贝叶斯最优估计可表 示为 xˆ(k|Y1:k) = xˆ(k|Y1:k−1)+ g(k)[φY (k)− hxˆ(k|Y1:k−1)], M(k|Y1:k) =M(k|Y1:k−1)− g(k)hM(k|Y1:k−1), (30) 其中增益 g(k) = M(k|Y1:k−1)hT hM(k|Y1:k−1)hT + 1/ϕ. (31) [Ckf ]的证明见附录1. 4.2 量量量化化化卡卡卡尔尔尔曼曼曼滤滤滤波波波器器器(Quantized Kalman filter- ing) 下面研究融合中心的量化卡尔曼滤波问题. 在k时刻, 假设融合中心已经得到k − 1时刻状 态的最优估计xˆ(k − 1|B1:k−1)和估计误差协方差阵 M(k − 1|B1:k−1),则根据假设1和系统模型(1),有 p(x(k)|B1:k−1) = N [x(k); xˆ(k|B1:k−1),M(k|B1:k−1)], (32) 其中: xˆ(k|B1:k−1) = Axˆ(k − 1|B1:k−1), M(k|B1:k−1) = AM(k − 1|B1:k−1)AT + Cu. (33) 式(32)和式(33)给出了量化卡尔曼滤波器预测部 分[Pq]的解,下面研究其校正部分的问题. 根据Y (k)的定义以及式(32),可得 p(Y (k)|B1:k−1) = N [Y (k);Hxˆ(k|B1:k−1),HM(k|B1:k−1)HT+R], (34) 其中: H = [hT · · · hT]T,H ∈ RN×p,R = diag{σ21, · · · , σ2N}. 下面的定理1给出了量化卡尔曼滤波器校正部分 [Cq]的解. 定定定理理理 1 考虑系统模型(1)和(2)以及式(13)中定 义的bn(k),在k时刻,设融合中心接收到的量化观测 集为Bk,则有 xˆ(k|B1:k)=xˆ(k|B1:k−1)+ g¯(k)[φYˆ (k|B1:k−1, Bk)−hxˆ(k|B1:k−1)], (35) 且误差协方差阵为 M(k|B1:k) =M(k|B1:k−1)− g¯(k)hM(k|B1:k−1)+ λ(k)g¯(k)g¯T(k), (36) 其中: g¯(k) = M(k|B1:k−1)hT hM(k|B1:k−1)hT + 1/ϕ, Yˆ (k|B1:k−1, Bk) = E(Y (k)|B1:k−1, Bk), λ(k) = φ∆Y (k)φT, ∆Y (k) = E{Y˜ (k)Y˜ T(k)|B1:k−1}, Y˜ (k) = Y (k)− Yˆ (k|B1:k−1, Bk). (37) 证证证 见附录2. 注注注 2 定理1中量化卡尔曼滤波器的误差协方差阵 (36)可分为两部分: 前两项为基于分布式观测的卡尔曼滤 波器的误差协方差阵,最后一项为量化引起的增量部分,其 中标量系数λ(k)反映了量化效果, ∆Y (k)为量化误差协方差 阵. 易知,当{qn → ∞}Nn=1时,则Yˆ (k|B1:k−1, Bk) → Y (k), 且∆Y (k) → 0, λ(k) → 0,因此量化卡尔曼滤波器简化为分 布式卡尔曼滤波器. 定理1给出了量化卡尔曼滤波器的最优解, 但 是根据Y (k)的预测分布式(34)可知, 对于式(35)中 的Yˆ (k|B1:k−1, Bk)以及式(36)中λ(k)的∆Y (k), 其分 别需要计算一N维高斯先验分布的随机变量Y (k) 在给定量化观测Bk下的条件均值和量化误差协方 差阵,当N比较大时,解这类问题需要蒙特卡罗积分 方法[16],运算负担较重. 下面给出一种渐近等价于定理1的IQKF,其简单 易实现,运算代价也因此降低. 定义量化观测集B1:nk :={b1(k), b2(k),· · ·, bn(k)}, 且 xˆ(n)(k|B1:k−1) := E { x(k)|B1:k−1, B1:nk } , M (n)(k|B1:k−1) := E{[x(k)− xˆ(n)(k|B1:k−1)] [x(k)− xˆ(n)(k|B1:k−1)]T|B1:k−1, B1:n−1k }, (38) 其中n = 0, 1, 2, · · · , N ,且 xˆ(0)(k|B1:k−1) := xˆ(k|B1:k−1), M (0)(k|B1:k−1) :=M(k|B1:k−1). (39) 1734 控 制 理 论 与 应 用 第 28卷 同假设1类似, 对状态后验分布p(x(k)|B1:k−1, B1:nk )作高斯近似,即 p(x(k)|B1:k−1, B1:nk ) = N [x(k); xˆ(n)(k|B1:k−1),M (n)(k|B1:k−1)]. (40) 根据上式即可得 p(yn+1(k)|B1:k−1, B1:nk ) = N [yn+1(k);hxˆ(n)(k|B1:k−1), hM (n)(k|B1:k−1)hT + σ2n+1]. (41) 迭代量化卡尔曼滤波器的预测部分同式(32)和 式(33),其校正部分由下面的定理2给出. 定定定理理理 2 对于定理1中的量化卡尔曼滤波器的 校正部分,有如下渐近等价的迭代算法: 对迭代序列n = 1, 2, · · · , N ,有 xˆ(n)(k|B1:k−1) = xˆ(n−1)(k|B1:k−1)+ g(n)(k)[yˆn(k|B1:k−1, B1:n−1k , bn(k))− hxˆ(n−1)(k|B1:k−1)], (42) 且误差协方差阵为 M (n)(k|B1:k−1) = M (n−1)(k|B1:k−1)− g(n)(k)hM (n−1)(k|B1:k−1)+ ∆ (n) yn(k) g(n)(k)g(n)T(k), (43) 其中: g(n)(k) = M (n−1)(k|B1:k−1)hT hM (n−1)(k|B1:k−1)hT + σ2n , yˆn(k|B1:k−1, B1:n−1k , bn(k)) = E[yn(k)|B1:k−1, B1:n−1k , bn(k)], ∆ (n) yn(k) = E{y˜2n(k)|B1:k−1, B1:n−1k }, y˜n(k) = yn(k)− yˆn(k|B1:k−1, B1:n−1k , bn(k)). (44) 在每个采样周期, 融合中心迭代计算式(42)和 式(43)N次, 即可得 xˆ(k|B1:k) = xˆ(N)(k|B1:k−1), M(k|B1:k) =M (N)(k|B1:k−1). 证证证 见附录3. 注注注 3 对于IQKF算法,在每个采样周期,融合中心每 接收到一个传感器的量化观测,即可对状态估计进行更新, 最终得到全局最优估计.因此,其在存在丢包或延迟的网络 环境下亦适用,融合中心在每个采样周期可以根据当前接 收到的量化观测形成对状态当前的最优估计. 4.3 局局局部部部量量量化化化卡卡卡尔尔尔曼曼曼滤滤滤波波波器器器(Local quantized Kalman filtering) 传感器的编码器和融合中心的解码器的在线更 新都需要通过局部量化卡尔曼滤波来实现. 式(18) 和式(19)给出了局部量化卡尔曼滤波器预测部分 [Plq]的解, 而对于其校正部分[Clq], 可根据定理1得 到下面的推论. 推推推论论论 1 考虑状态空间模型(1)和(2)以及式(13) 中定义的bn(k), 在k时刻, 假定传感器n和融合中心 已经得到预测估计xˆ(k|b(n)1:k−1)和预测误差协方差阵 M(k|b(n)1:k−1), 则当传感器n和融合中心得到局部量 化观测bn(k)后,即有 xˆ(k|b(n)1:k)=xˆ(k|b(n)1:k−1)+gn(k)[yˆn(k|b(n)1:k−1, bn(k)]− hxˆ(k|b(n)1:k−1)], (45) 且其局部量化估计误差协方差阵为 M(k|b(n)1:k) =M(k|b(n)1:k−1)− (1− K(qn) 4qn )gn(k)hM(k|b(n)1:k−1), (46) 其中: gn(k) = M(k|b(n)1:k−1)hT hM(k|b(n)1:k−1)hT + σ2n , yˆn(k|b(n)1:k−1, bn(k)] = E[yn(k)|b(n)1:k−1, bn(k)]. (47) 证证证 根据定理1,只需考虑单个传感器观测的特 殊情形,此时φ = 1, ϕ = 1/σ2n,则由式(35)(36)即可 得式(45)以及 M(k|b(n)1:k) =M(k|b(n)1:k−1)− gn(k)hM(k|b(n)1:k−1)+ ∆yn(k)gn(k)g T n (k), (48) 其中 ∆yn(k) = E{[yn(k)− yˆn(k|b(n)1:k−1, bn(k)]2|b(n)1:k−1}. (49) 又由式(18)可得 p(yn(k)|b(n)1:k−1) = N [yn(k);hxˆ(k|b(n)1:k−1), hM(k|b(n)1:k−1)hT+ σ2n]. (50) 因此, 根据高斯信源的Lloyd-Max量化器平均失 真表达式(12)即可得 ∆yn(k) = K(qn)(hM(k|b(n)1:k−1)hT + σ2n) 4qn . (51) 把式(51)代入式(48),化简即可得式(46). 证毕. 5 量量量化化化卡卡卡尔尔尔曼曼曼滤滤滤波波波器器器的的的稳稳稳定定定性性性(Stability of quantized Kalman filtering) 本节分析定理1中量化卡尔曼滤波器的稳定性问 题.定理1适用于任何采用标量量化器的量化卡尔曼 滤波问题,研究预测误差协方差阵M(k|B1:k−1)的收 敛性与量化器以及量化比特率的关系. 简记M qk := M(k|B1:k−1), 根据式(33)(36)可得 第 12期 陈军勇等: 无线传感器网络分布式量化卡尔曼滤波 1735 黎卡提方程 M qk+1 =AM q kA T + Cu −Ag¯(k)hM qkAT+ λ(k)Ag¯(k)g¯T(k)AT, (52) 其中: k = 1, 2, · · · ,∞,初始值M q1 = AM0AT+Cu. 对于M qk ,下面给出一个定义: 定定定义义义 1 若∀M0 º 0,有 lim k→∞ M qk = M q ∞,其中 M q∞ Â 0, 且与M0无关, 则称M qk绝对收敛; 若∀M0 º 0, 有 lim k→∞ M qk ≺ +∞, 则称M qk绝对有界; 若∃M0 º 0,有 lim k→∞ M qk = +∞,则称M qk条件发散. 黎卡提方程(52)中的λ(k)反映了量化效果,其与 量化器以及各传感器的量化比特率有关. 众所周 知, 若系统矩阵A稳定, 则M qk必绝对收敛或绝对有 界[15, 17], 若系统矩阵A不稳定, 则M qk的收敛性与量 化器以及量化比特率{qn}Nn=1有关. 下面重点考虑A不稳定的情形,分别对单传感器 系统(N = 1)和多传感器系统(N > 1)的M qk的收敛 性进行讨论. 5.1 单单单传传传感感感器器器系系系统统统(One-sensor system) 首先考虑网络中只存在一个传感器观测的情形, 即N = 1. 不失一般性,设此时的观测为传感器n的 观测,量化比特率为qn位. 记标量系数α := 1 − K(qn) 4qn , 函数gn(M) := MhT hMhT + σ2n ,则有下面的定理. 定定定理理理 3 假定(A,C 12u )能控, (A, h)能观, A不稳 定, N = 1, 则在动态Lloyd-Max量化器下, M qk绝对 收敛的充要条件为 K(qn) 4qn < 1∏ i |λui |2 , (53) 其中{λui }为A的不稳定特征值集合, 则其绝对收敛 的临界比特率可表示为 qc = min{q|K(q)4q < 1∏ i |λui |2 , q > 0, q ∈ Z}, (54) 即若qn>qc,则M qk绝对收敛,且唯一正定收敛值M q∞ 可由M = AMAT+Cu−αAgn(M)hMAT解得;若 qn < qc,则M qk条件发散. 证证证 在单传感器观测下,传感器的局部量化卡尔 曼滤波即为融合中心的量化卡尔曼滤波,因此根据 式(19)以及推论1中的式(46),式(52)可简化为 M qk+1 = AM q kA T +Cu − αAgn(M)hM qkAT, (55) 这里M qk表示单传感器观测下的预测误差协方差阵, 即M qk =M(k|b(n)1:k−1). 对于黎卡提方程(55), 其收敛性已有相关的结 论[17∼19]. 根据Elia[19], 由于(A,C 12u )能控, (A, h)能 观, A不稳定,且rank(h) = 1,因此其绝对收敛的充 要条件为 α > 1− 1∏ i |λui |2 , (56) 即为式(53). 定理3中其他相关结论亦可根据文献 [17∼19]得到. 证毕. 注注注 4 在单传感器观测下, Lloyd-Max量化器在 极小化均方误差意义上实际上是最优量化器[7], 这 亦可从单传感器观测下的误差协方差阵(48)可看出, 极小化均方误差等价于极小化观测的平均失真. 注注注 5 若视K(q) = K为常数, 则M qk绝对收敛 的临界比特率为1 2 log2K ∏ i |λui |2, 其与系统矩阵的 不稳定特征值以及量化器性能指标K有关. 从信 息论[20]的角度看, K的理论最小值为1,则可简化为 qc = log2 ∏ i |λui |, 这与文献 [9, 10]的结论是一致的, 但是K = 1需对信号分块编码且块长为无穷大的时 候才能达到,对于本文的实时目标跟踪系统,这是难 以实现的, 因此本文的临界比特率不仅与系统矩阵 的不稳定特征值有关,亦与量化器性能指标K有关. 5.2 多多多传传传感感感器器器系系系统统统(Multi-sensor system) 现在分析多传感器下(N > 1)黎卡提方程(52)的 收敛性问题.设网络中的总比特率为Q, Q = N∑ n=1 qn, 其中Q > 0, Q ∈ Z. 首先不限定量化器的具体形式, 则有如下定理. 定定定理理理 4 假定(A,C 12u )能控, (A, h)能观, A不稳 定, N > 1,则对任何标量量化器, M qk绝对有界的必 要条件为 K(Q) 4Q < 1∏ i |λui |2 . (57) 上述必要条件亦等价于Q > qc. 证证证 见附录4. 注注注 6 定理4说明在多传感器系统中, 对任何标量量 化器, Mqk绝对有界所需的最小总比特率亦不会低于qc, 亦 即传感器数目增加不能降低黎卡提方程(52)稳定的临界总 比特率,其仍然只与系统矩阵的不稳定特征值以及K(q)有 关. 注注注 7 由定理4的证明可知, 由于Y ′(k) = φY (k), 若 对Y ′(k)采用Lloyd-Max量化器, 则式(A25)取等号, 因此∀k, 有Mqk = Zk,此时定理4中的必要条件则为充要条件,且此 时的量化器在均方误差意义上为最优量化器,但由于Y ′(k) 为各传感器观测的线性组合, 若N = 1, 则Y ′(k)为一维观 测, 定理4退化为定理3, 若N > 1, 则需要对传感器的观测 进行联合编码才能实现. 1736 控 制 理 论 与 应 用 第 28卷 在分布式量化卡尔曼滤波问题中, 传感器需要 独立编码各自的观测. 对于本文提出的动态Lloyd- Max量化器,根据定理3易得下面的推论. 推推推论论论 2 假定(A,C 12u )能控, (A, h)能观, A不稳 定, N > 1, 则在动态Lloyd-Max量化器下, M qk绝对 有界的充要条件为∃n ∈ {1, 2, · · · , N},有qn > qc. 证证证 充分性. 若∃n ∈ {1, 2, · · · , N},有qn > qc, 则根据定理3,由于(A,C 12u )能控, (A, h)能观, A不稳 定, 则至少一个传感器的量化观测可以使M qk绝对 有界,又根据激活的传感器数目增加M qk必降低的事 实,可知M qk必绝对有界. 必要性(反证法). 若∀n, 有qn < qc, 则根据定 理3和文献 [18], 由于(A,C 12u )能控, (A, h)能观, A不 稳定,可知∃M0,有 lim k→∞ M(k|b(n)1:k−1) = +∞, ∀n,而 根据动态Lloyd-Max量化器更新规则,可知此时所有 传感器的量化分区将随k增加而趋于无穷大,量化失 真亦将趋于无穷大, 则 lim k→∞ λ(k) = +∞, 又根据黎 卡提方程(52),可知 lim k→∞ M qk = +∞,这与M qk绝对有 界矛盾. 证毕. 6 仿仿仿真真真结结结果果果(Simulation results) 本节将通过MATLAB仿真验证文中所设计算法 的可行性和有效性. 考虑下面的二维目标跟踪系统: x(k) = [ 1 T 0 1 ] x(k − 1) + [ T 2/2 T ] u(k). (58) 上述系统实际为一常速运动模型, 其中: 状态向量 x(k) := [s(k) s˙(k)]T, s(k)和s˙(k)分别为目标位置 和速度, T为采样周期,驱动噪声为u(k). 传感器的观测模型为 yn(k) = [1 0]x(k) + vn(k), n = 1, 2, (59) 其中: v1(k) ∼ N (0, 1), v2(k) ∼ N (0, 2). 在全部仿真试验中, 系统初始条件设为x(0) ∼ N (x0,M0), 其中x0 = (0, 5)T, M0 = 0.3I2, 驱动噪 声u(k) ∼ N (0, 1),采样周期T = 0.1 s. 对于均方误 差(MSE), 本文采用5000次蒙特卡罗仿真结果取平 均值得到. 本文分析基于动态Lloyd-Max量化器的IQKF在 不同带宽限制下的估计性能.根据目标运动模型(58) 和推论2可知M qk绝对有界的临界比特率qc = 1位. 图5给出 IQKF–1位和 IQKF–2位以及基于未量化观 测的分布式卡尔曼滤波器(KF)的理论和实际MSE曲 线比较. IQKF–1位和 IQKF–2位分别指两个传感器 的量化比特率均为1位和均为2位的情形. 图5(a)为 理论MSE比较, 图5(b)为实际MSE比较. IQKF和 KF的理论MSE分别通过文中M (N)(k|B1:k−1)和 M(k|Y1:k)的显示解得到, 实际MSE指状态估计和 状态真实值之间的均方误差. (a)理论MSE比较 (b)实际MSE比较 图 5 各种滤波器的理论MSE和实际MSE比较 Fig. 5 Comparison of theoretical MSE and empirical MSE for various filters 从图5(a)和图5(b)可以看出, 分布式卡尔曼滤波 器(KF)的MSE为IQKF的MSE的理论不可达下界,且 IQKF–2位的MSE小于IQKF–1位的MSE, 这说明量 化比特率愈大, IQKF估计性能愈好.另外,由于都存 在某个传感器的量化比特率不小于qc,因此随着k增 加, IQKF–1位和IQKF–2位的MSE都趋于稳定. IQKF 在量化比特率都为最小的1位时都能很好的跟踪目 标,这说明了IQKF的高效性. 图 6给出了 IQKF的理论MSE和实际MSE曲线 比较. 从图 6可以看出,不同量化比特率下 IQKF的 理论MSE和实际MSE都很接近, 这说明文中所得 M (N)(k|B1:k−1)的显示解的有效性. 另外, IQKF–2 位的MSE差异比IQKF-1位的MSE差异更小, 这是因 为量化比特率愈大,文中对状态后验分布的高斯近 似愈与实际相符,因此理论MSE和实际MSE愈接近. 图 6 IQKF的理论MSE和实际MSE比较 Fig. 6 Comparison of theoretical MSE and empirical MSE for IQKF 第 12期 陈军勇等: 无线传感器网络分布式量化卡尔曼滤波 1737 为了更好地分析本文提出的 IQKF算法的性能, 本 文 对 IQKF和Sun等[8]提 出 的QKFQI(quantized Kalman filter based on quantized innovations)算法的 估计性能进行了比较,其中QKFQI采用均匀量化器, 且量化器的更新需要融合中心反馈预测信息.图7给 出了IQKF–2位和QKFQI–2位算法的实际MSE曲线 比较, 从图7可以看出, 本文的IQKF估计器性能优 于QKFQI,且由于本文的动态Lloyd-Max量化器更新 过程无需融合中心的反馈信息,因此通信代价亦得 到节省. 图 7 IQKF和QKFQI的实际MSE比较 Fig. 7 Comparison of empirical MSE for IQKF and QKFQI 7 结结结论论论(Conclusions) 本文基于无线传感器网络研究了分布式量化卡 尔曼滤波问题.提出了一种新的动态Lloyd-Max量化 器并设计了其在线更新方案,基于贝叶斯原理得到 了递归形式的最优量化卡尔曼滤波器, 并与标准卡 尔曼滤波器进行了比较, 同时给出了一种渐近等价 的迭代算法(IQKF).对于不稳定系统,量化卡尔曼滤 波器稳定的临界比特率由量化器以及系统矩阵的不 稳定特征值决定,与传感器数目和观测噪声无关.动 态Lloyd-Max量化器更新无需融合中心的反馈信息, 通信代价小、适用性强. IQKF算法运算代价低,且在 存在丢包或延迟的网络环境中亦适用. 因此,本文所 提出的设计方法能得到较广泛的应用. 参参参考考考文文文献献献(References): [1] GUBNER J A. Distributed estimation and quantization[J]. IEEE Transactions on Information Theory, 1993, 39(4): 1456 – 1459. [2] LI J L, ALREGIB G. Distributed estimation in energy-constrained wireless sensor networks[J]. IEEE Transactions on Signal Process- ing, 2009, 57(10): 3746 – 3758. [3] ZHANG K S, LI X R. Optimal sensor data quantization for best lin- ear unbiased estimation fusion[C] //The 43rd IEEE Conference on Decision and Control. New York: IEEE, 2004: 2656 – 2661. [4] MARIC I, YATES R D. Bandwidth and power allocation for cooper- ative strategies in Gaussian relay networks[J]. IEEE Transactions on Information Theory, 2010, 56(4): 1880 – 1889. [5] KHAN U A, MOURA J M F. Distributing the Kalman filter for large- scale systems[J]. IEEE Transactions on Signal Processing, 2008, 56(10): 4919