A NEW TECHNIQUE OF COMPUTER GRAPHICS SYNTHESIS

1.   引言

无线传感器网络(Wireless Sensor Networks, WSN)是由一系列能够对环境做出感知和观测的小型装置通过无线通信组成的自组织网络。由于无线传感器网络成本低廉、易于部署，使得其在监控、安防、交通、定位、医疗等军用和民用领域^[1–3]发挥了越来越重要的作用。

目标跟踪^[4–7]是无线传感器网络的一个重要应用场景。在实际应用中，传感器所需能量通常由电池提供，对整个传感器网络的电池进行充电或更换十分不易。如果所有传感器在所有时刻都开启观测，虽然可以得到较高的跟踪精度，但同时也会带来较大的能量开销。为了延长传感器网络的工作时间，通常需对传感器节点进行自适应的节点调度^[4,5]，或称为节点规划^[7]、节点选择^[8]或传感器分配^[9]等。在每一个决策周期中，节点调度算法选择最优的传感器子集对目标进行观测，以达到跟踪精度和能量功耗的折中。

国内外许多学者对目标跟踪的节点调度算法进行了研究。一些研究将该问题建模为最小化瞬时估计误差或最大化瞬时信息增益的优化问题，建立了基于某个准则下的代价函数，例如熵和相对熵、互信息、先验克拉美罗下界、条件克拉美罗下界、最小均方误差等等^[8,10,11]。文献[9]根据目标预测位置，基于势博弈模型求解纳什均衡，达到了较快的收敛速度。这些贪婪算法仅考虑节点调度算法的瞬时表现，没有考虑对未来的影响，因此可被视为“短视”的策略。

在目标跟踪过程中，当前时刻传感器的调度策略会影响传感器对目标的观测结果，该观测结果会影响当前对目标状态的估计，状态估计又会影响下一时刻传感器的调度策略。显然这是一个序贯的决策过程，当前时刻的决策会影响未来的走向，因此必须综合考虑决策的瞬时表现和长远表现。

为了得到长远考虑的优化决策，可以将该问题建模为一个部分可观测马尔科夫决策过程(Partially Observable Markov Decision Process, POMDP)^[12]。在POMDP中，由于观测噪声或其他原因，目标状态不能根据观测完全确定，文献中称为部分可观测。POMDP研究的是如何根据历史动作和观测序列，在考虑未来代价的情况下求解最优策略的问题。文献[13]提出了POMDP框架下的类卡尔曼估计，并被应用于传感器网络中的主动目标跟踪^[14]。

由于POMDP固有的计算复杂度，只有问题的状态集、动作集和观测集规模很小的情况下，该问题才能被精确求解^[12]。为了解决这个问题，学者们提出了许多近似求解POMDP问题的启发式算法^[14,15]。当代价函数关于置信状态是线性函数时，近似算法有PBVI, PEMA, Perseus, HSVI, FSVI, SARSOP等，但是这些算法不适用于一般的非线性代价函数的情况。QMDP算法^[16]为了处理非线性代价函数，简化未来代价的计算，假设当前状态经过一步转移后，目标状态变为完全可观测，从而基于MDP^[17]计算未来代价，但是该算法不能直接应用于状态集是连续空间的情况。对于连续状态空间的情形，He等人^[18]将Rollout算法^[19]应用到POMDP框架下的节点调度算法中，在计算未来代价时使用基本策略代替最优策略，并采用在线模拟的方法近似Q-value，取得了不错的效果。Li等人^[20,21]做出了和QMDP算法相同的假设，即经过一步转移后，目标状态变为完全可观测的，从而使用CO-Rollout方法^[19]近似Q-value，求解优化策略，简化了未来代价的计算。然而，Rollout和CO-Rollout算法均使用基本策略代替最优策略，对Q-value的近似精度较低；此外，在线模拟的方式计算复杂度过高，难以应用于工程实践。

本文提出一种近似求解算法C-QMDP，该算法通过两个步骤，分别计算节点调度的瞬时代价和未来代价。对于目标跟踪问题，利用MCMC方法推导置信状态的转移，跟踪目标位置，并计算瞬时代价；对于节点调度问题，将目标的状态空间离散化，基于MDP值迭代求解未来代价的近似值，综合考虑瞬时代价和未来代价选择传感器子集。使用C-QMDP算法进行节点调度的优势有以下几点：(1)将节点调度问题建模为POMDP，综合考虑了决策的瞬时表现和长远表现；(2)可以求解MDP最优策略下各个状态的最小损失来近似Q-value，跟踪更加准确；(3)计算未来损失的过程可以离线计算，将各个状态的最小损失存储起来，在线决策时节约了计算量，提高了系统的实时性。

2.   问题建模

本文的研究场景描述如下：在一块区域中，部署了M个位置已知的传感器，用于跟踪单个目标。传感器观测结果包含噪声，并且观测范围有限。目标转移符合近似恒定速度(Nearly Constant Velocity, NCV)模型^[22]。存在一个中心处理单元，接收传感器回传的观测数据，进行状态估计和决策，决定下一时刻系统中每个传感器的开关。决策的目标是选择打开一组传感器，综合考虑跟踪精度和传感器功耗，使得跟踪过程累积代价最小。

2.1   场景建模

(1)系统状态：系统状态包含了在2维平面上移动的目标的位置和速度： ${{S}} = {[x,\dot x,y,\dot y]^{\rm{T}}}$ ，其中 $[x,y]$ 表示目标在笛卡尔坐标系下的位置， $\dot x$ 和 $\dot y$ 分别表示目标在对应方向上的速度。此外，采用 ${{{p}}_s} = [x,y]$ 表示目标的位置。

(2)动作动作： ${A}$ 表示对传感器网络中的节点做出开启或关闭的动作： ${{A}} = {[{a_1},{a_2},·\!·\!·,{a_M}]^{\rm{T}}},$ ${a_m} \in \{ 0,1\}$ ，其中M表示传感器网络中传感器的个数， ${a_m} \in \left\{ {0,1} \right\}$ 表示关闭或打开第m个传感器。

(3)状态转移模型：状态转移模型刻画了目标如何从一个状态转移到另一个状态。假设目标的速度变化缓慢，可以使用近似恒定速度模型^[22]：

其中， ${T_s}$ 表示决策周期的时长。 ${v_x}$ 和 ${v_y}$ 是均值为0，方差为 $\sigma _x^2$ 和 $\sigma _y^2$ 的高斯变量，分别表示x和y方向上的加速度的不确定性。

(4)观测模型：传感器观测模型使用电子通信领域常见的RSSI (Received Signal Strength Indicator)观测模型^[23,24]，传感器观测目标发出的能量信号强度与传感器间距离的 $\lambda$ 次方成反比(通常 $\lambda \ge 2$ )。于是观测模型为

其中， ${E_0}$ 表示已知的目标信号功率， ${{{p}}_m}$ 为第m个传感器的位置。 ${w_m}$ 是一个均值为0，方差为 $\sigma _m^2$ 的高斯噪声，表示传感器m的观测噪声。所有传感器的观测组成的集合可以写作 ${{Z}} = {[{z_1},{z_2},·\!·\!·,{z_M}]^{\rm{T}}}$ ，其中 ${z_m}$ 表示第m个传感器采集的观测，k时刻的观测矢量表示为 ${{{Z}}_k}{\rm{ = }}{[{z_{k,1}},{z_{k,2}},·\!·\!·,{z_{k,M}}]^{\rm{T}}}$ 。

此外，定义 ${N_m} = \displaystyle\frac{{{z_m}}}{{\sigma _m^2}}$ 衡量传感器m获得观测 ${z_m}$ 时的相对接收强度，并定义 ${N_{{\rm{th}}}}$ 表示传感器的观测阈值。当 ${N_m} < {N_{{\rm{th}}}}$ 时，可认为此时观测是无效的，目标位置超出了传感器m的观测范围。

事实上，本文算法并不局限于RSSI观测模型，由于后续采用粒子滤波进行状态估计，因此也可运用于其他观测模型。

(5)置信状态：每个时刻可基于历史观测和采取的动作序列得到目标状态的后验分布，该分布被称为目标的置信状态，记作 ${b}$ 。由于目标状态不可直接观测，因此在POMDP中，使用置信状态 ${b}$ 的迁移来刻画目标的运动轨迹。k时刻的置信状态 ${{{b}}_k} = {\rm{P}}({{{S}}_k}|{{{Z}}_{1:k}},{{{A}}_{0:k - 1}})$ ，其中 ${{Z}_{1:k}}$ 表示时刻1到k传感器网络回传的观测序列， ${{{A}}_{0:k - 1}}$ 表示时刻0到k –1传感器网络采取的动作序列。在初始时刻，目标置信状态为 ${{{b}}_0}$ ，表示目标状态的先验。

(6)瞬时代价函数：瞬时代价函数表示目标在状态转移过程中每一步所产生的瞬时代价，包含两个方面：跟踪误差以及传感器功耗。使用常数 $\alpha$ 权衡跟踪误差和传感器功耗对最终决策影响的重要性，代价函数可表示为

其中， ${{\hat{ p}}_s} = [\hat x,\hat y]$ 表示目标位置的估计， $c_m^{{\rm{on}}}$ 表示第m个传感器在单位决策时间内开启的功耗。期望是对 ${{{p}}_s} = [x,y]$ 求取的，根据当前置信状态b，可以得到目标真实位置的后验分布 $\,{{b}}({{{p}}_s})$ 。 ${{\hat{ p}}_s}$ 是该后验分布的均值，可视为置信状态b的函数，因此代价函数式(3)相对于置信状态b是非线性的。

2.2   POMDP求解过程

POMDP的迭代过程如下：在k时刻，目标真实状态为 ${{{S}}_k}$ 时，得到观测 ${{{Z}}_k}$ ，此时可计算出目标状态的后验，即置信状态 ${{{b}}_k}$ ，根据置信状态进行决策，采取动作 ${{{A}}_k}$ 并获得代价 ${L_k}$ 。然后，系统状态将转移到 ${{{S}}_{k + 1}}$ ，得到观测 ${{{Z}}_{k{\rm{ + }}1}}$ ，从而置信状态转移到 ${{{b}}_{k{\rm{ + }}1}}$ ，再进行k+1时刻的决策。节点调度问题k 时刻贝叶斯网络如图1所示，灰色节点表示已知变量，白色节点表示未知变量。

图  1  贝叶斯网络

下载: 全尺寸图片 幻灯片

在POMDP框架下，节点调度的目标是，寻找一个策略，即一个置信状态 ${b}$ 到动作 ${A}$ 的映射 ${{π}}$ ，使得能够根据目标的置信状态，动态选取最优的传感器子集对目标进行观测，从而最小化累积损失的期望U。若已给定策略 ${{π}}$ ，则在第k步时，累积损失的期望U为

其中， ${{{A}}_t} = {{π}} ({{{b}}_t})$ ，期望是对所有可能的置信状态和观测求取的。式(4)被称为Bellman方程， $L({{{b}}_k},{{{A}}_k})$ 表示在置信状态为 ${{{b}}_k}$ 时采取动作 ${{{A}}_k}$ 后当前时刻的损失， $\mathbb{E}(U({{{b}}_{k + 1}},{{π}} )|{{{b}}_k},{{{A}}_k})$ 表示采取动作 ${{{A}}_k}$ 后未来累积损失的期望， $0 \le \beta < 1$ 为递减因子，用来权衡瞬时代价和未来代价，并使累积损失有界。

为了更方便地在第k步时求解最优动作 ${{{A}}_k}$ ，定义Q-value为 $Q({{{b}}_k},{{A}})$ 表示在置信状态为 ${{{b}}_k}$ 时采取动作 ${A}$ ，未来采取最优策略的总期望损失，于是，

其中， ${V^*}({{{b}}_k})$ 表示置信状态为 ${{{b}}_k}$ 时，当前及之后每步都采取最优策略所能达到的期望最小累积损失， ${{{π}} ^{\rm{*}}}$ 为要寻找的最优策略。式(5)中，Q-value左半部分表示当前损失，右半部分表示未来损失，其中期望是对所有可能的 ${{{b}}_{k{\rm{ + }}1}}$ 求取的。因此，当置信状态为 ${{{b}}_k}$ 时，如果能够求得每个动作的Q-value，即可采取Q-value最小的动作作为当前置信状态下的最优策略。

3.   节点调度算法

前一节基于POMDP框架对传感器网络进行目标跟踪的场景进行了建模，本节将对该问题进行近似求解。在每个决策周期中，传感器得到目标的观测，并产生下一时刻需要采取的动作。这里，要解决两个方面的问题：(1)如何跟踪目标的置信状态；(2)如何根据置信状态产生下一时刻要采取的动作。

首先，考虑如何描述状态迁移的问题。在状态转移模型和观测模型皆为线性方程、噪声均为高斯噪声的场景下，可以使用卡尔曼滤波器求得目标置信状态迁移过程的解析解。在更一般的情况下，线性高斯假设不一定成立，例如本文中的观测模型，此时无法获得解析解。在非线性系统中，粒子滤波^[25]是一种常用的近似计算目标置信状态迁移的方法。它使用蒙特卡洛采样近似后验概率的积分，可以解决转移方程或观测方程非线性的问题。

接下来，考虑节点调度问题。根据式(6)，只需找到让 $Q({{{b}}_k},{{A}})$ 最小的动作，即可保证当前时刻到未来的累积损失的期望最小。因此，决策过程的重点是估计Q-value。精确地求解Q-value是十分困难的，尤其是当状态集、观测集较大的时候，计算复杂度十分庞大。不过，Q-value并不需要被精确计算，因为Q-value的意义在于度量动作的好坏，寻找最优的动作。即使Q-value近似有误差，只要依然能够通过Q-value的大小对动作的优劣进行排序，选出最优动作即可。本文基于C-QMDP的算法框架进行Q-value近似，离线计算未来损失，在线计算当前损失，可以大大降低计算复杂度，提高决策的实时性。

3.1   求解当前时刻瞬时代价

在每个决策周期中，需要根据现有观测，更新对目标状态的估计。在第1节中，已经定义了置信状态 ${{{b}}_k} = {\rm{P}}({{{S}}_k}|{{{Z}}_{1:k}},{{{A}}_{0:k - 1}})$ 来表示k时刻的目标状态的后验分布。将 $({{{Z}}_{1:k{\rm{ - }}1}},{{{A}}_{0:k - 1}})$ 记作 ${{{I}}_k}$ ，表示k时刻之前的历史观测和动作。由贝叶斯迭代^[26]公式，当没有获得k时刻的观测 ${{{Z}}_k}$ 时，根据上一时刻的置信状态 ${{{b}}_{k - 1}}$ ，可以得到k时刻目标状态 ${{{S}}_k}$ 的预测分布：

得到k时刻观测 ${{{Z}}_k}$ 后，可以求得目标状态 ${{{S}}_k}$ 的后验分布，即

(8)

其中， $\gamma = \displaystyle\int\nolimits_{{{{S}}_k}} {{\rm{P}}({{{Z}}_k}|{{{S}}_k},{{{A}}_{k - 1}}){\rm{P}}({{{S}}_k}|{{{I}}_k}){\rm{d}}{{{S}}_k}}$ 。

根据式(8)，可以获得置信状态 ${{{b}}_{k{\rm{ - }}1}}$ 到 ${{{b}}_k}$ 的状态迁移。由于积分难以获得解析解，粒子滤波使用蒙特卡洛模拟的方式，使用N个粒子去近似置信状态 ${{{b}}_k}$ ，以数值方法求解上述积分：

其中， ${{S}}_k^{(i)}$ 表示第k步时第i个粒子的状态值， $w_k^{(i)}$ 表示第k步时第i个粒子的权重， ${δ}$ 为狄拉克函数。

粒子滤波的实现有很多种，本文使用常见的SIR滤波器^[25]。权值更新公式为 $w_k^{(i)} = w_{k - 1}^{(i)}{\rm{P}}\left({{{Z}}_k}|{{S}}_k^{(i)}\right)$ 。假设不同传感器之间的观测噪声相互条件独立，因此有

其中， ${a_{k,m}}$ 表示k时刻传感器m的开关状态，当 ${a_{k,m}}{\rm{ = }} 1$ 时，传感器m的观测参与计算。

重采样后粒子权重均为 ${1}/{N}$ ，即 $w_{k - 1}^{(i)} = {1}/{N}$ 。因此，权值更新公式可以写作

之后对权值进行归一化：

根据当前N个粒子的状态和权重，可以求得当前目标状态的MMSE估计，即后验均值：

基于粒子滤波算法，可以更新目标置信状态的状态迁移。于是，在粒子滤波的时刻k，可以将此时的粒子集 ${{{b}}_k}$ 代入式(3)，求解当前时刻的瞬时代价：

其中， ${{\hat{ p}}_{s,k}}$ 表示目标状态估计 ${{\hat{ S}}_k}$ 中的目标位置部分， ${{p}}_{s,k}^{(i)}$ 表示粒子状态 ${{S}}_k^{(i)}$ 中的目标位置部分。

3.2   基于C-QMDP算法求解未来代价

根据式(4)，由于 ${{{b}}_k}$ 和 ${A}$ 都是已知的，可直接计算当前损失。然而，求解未来损失，即 $\mathbb{E}({V^*}({{{b}}_{k + 1}})|$ ${{{b}}_k},{{A}})$ ，是十分困难的。

C-QMDP算法中将Q-value近似计算分为两部分，在POMDP的框架下，基于3.1节中粒子滤波算法更新置信状态，将 ${{{b}}_k}$ 和 ${A}$ 直接代入式(3)计算当前损失。求解未来损失时，假设经过一步转移后，目标状态变为完全可观测的，此时使用目标状态 ${S}$ 而不是置信状态 ${b}$ 的迁移来描述目标的移动，未来损失由求解 $\mathbb{E}({V^*}({{{b}}_{k + 1}})|{{{b}}_k},{{A}})$ 转变为求解 $\mathbb{E}({V^*}({{{S}}_{k + 1}})|$ ${{{b}}_k},{{A}})$ ，其中目标状态用粒子的值表示。

为了求解 $\mathbb{E}({V^*}({{{S}}_{k + 1}})|{{{b}}_k},{{A}})$ ，需要计算置信状态为 ${{{b}}_k}$ 时执行动作 ${A}$ 后的置信状态 ${{{b}}_{k{\rm{ + }}1}}$ ，即k+1时刻 ${{{S}}_{k + 1}}$ 可能的位置和概率。由于k+1时刻的观测 ${{{Z}}_{k + 1}}$ 在当前时刻是未知的，因此需要根据 ${{{b}}_k}$ 和 ${A}$ 预测 ${{{Z}}_{k + 1}}$ 的值，从而计算 ${{{b}}_{k{\rm{ + }}1}}$ 的预测分布。 ${{\hat{ S}}_k}$ 经过式(1)一步转移后，作为k+1时刻 ${{{S}}_{k + 1}}$ 位置的估计。采取动作 ${A}$ 后，根据式(2)计算 ${{{Z}}_{k + 1}}$ 的预测值，代入3.1节中即可计算k+1时刻粒子的位置和权重。由于假设计算未来代价时目标完全可观测，因此k+1时刻开始可用粒子集合中的粒子代表目标状态 ${S}$ ，粒子权重w表示目标在这个状态的概率。在MDP的框架下，使用值迭代离线计算各状态 ${S}$ 采取最优策略时对应的损失 ${V^*}({{S}})$ 。于是，置信状态 ${{{b}}_k}$ 采取动作 ${A}$ 对应的未来损失的期望为

其中，概率 ${\rm{P}}({{{S}}_{k + 1}}|{{{b}}_k},{{A}})$ 是通过蒙特卡洛采样方式计算的，使用粒子的权重表示概率的值。然而，MDP只能直接应用于状态空间离散的情况。为了在连续状态空间中应用MDP算法计算式(15)中的 ${V^{\rm{*}}}({{S}}_{k + 1}^{(i)})$ ，可采取将连续状态空间离散化的方法。

基于C-QMDP求解最优策略的步骤如下：

(1)首先，将目标的状态空间进行离散化，以便应用MDP求解。使用网格划分的方法，将目标状态的各个分量在取值范围内等间隔划分。设离散化后的目标状态一共有D个，由 ${\tilde{ S}}$ 表示，则 ${{\tilde{ S}}_d}$ 表示目标处于第d个状态。

(2)然后，根据状态转移方程，计算离散状态两两之间的转移概率 ${\rm{P}}({\tilde{ S'}}|{\tilde{ S}})$ 。离散状态共有D个，故共需计算 ${{D(D - 1)}}/{2}$ 个概率值。

求解转移概率的解析解是困难的，转移概率多重积分的解析解无法获得，可以使用蒙特卡洛采样代替积分运算。当计算转移概率 ${\rm{P}}\left({\tilde{ S'}}|{\tilde{ S}}\right)$ 时，在状态 ${\tilde{ S}}$ 中等概率采样 $\tilde N$ 个点，记作 $\left[{{\tilde{ S}}^{(1)}},{{\tilde{ S}}^{(2)}},·\!·\!·,{{\tilde{ S}}^{(\tilde N)}}\right]$ 。状态 ${\tilde{ S'}}$ 是根据网格等间隔划分出的一个离散状态区域，假设该离散区域所代表的x, y方向上坐标和速度的界限为 $[{x_{\min }},{x_{\max }},{\dot x_{\min }},{\dot x_{\max }},{y_{\min }},{y_{\max }},{\dot y_{\min }},{\dot y_{\max }}]$ 。以x方向为例，对于采样点 ${{\tilde{ S}}^{(i)}}$ ，根据状态转移模型式(1)，有

其中， ${x^{(i)}}$ , ${\dot x^{(i)}}$ , ${T_s}$ 均为已知量，根据式(16)，可以计算出 ${v_x}$ 的上下界。由于 ${v_x}$ 是高斯分布，因此可以计算出 ${\rm{P}}_x^{(i)}\left({\tilde{ S'}}|{\tilde{ S}}\right)$ 的值。如果下界大于上界，则该概率为0。进一步，由x和y方向状态转移概率独立，可得 ${{\rm{P}}^{(i)}}\left({\tilde{ S'}}|{\tilde{ S}}\right) = {\rm{P}}_x^{(i)}\left({\tilde{ S'}}|{\tilde{ S}}\right){\rm{P}}_y^{(i)}\left({\tilde{ S'}}|{\tilde{ S}}\right)$ 。最后，离散状态的状态转移概率 ${\rm{P}}\left({\tilde{ S'}}|{\tilde{ S}}\right)$ 就是所有的采样点求得的转移概率的均值。

(3)接下来，根据代理代价函数^[19]，计算每个状态对应的瞬时代价 $\tilde L$ 。在POMDP框架中，代价函数是置信状态和动作到代价的映射；在MDP框架中，代价函数是离散的目标状态和动作到代价的映射。因此，不能沿用式(3)中定义的代价函数。

代理代价函数选取的标准是可以有效反映真实代价，计算出的Q-value可用于评判动作的优劣。本文使用累积观测质量的倒数作为代理代价函数，即

其中， ${N_m}({{S}})$ 表示传感器m在状态 ${S}$ 处相对接收强度，

当计算状态 ${{\tilde{ S}}_d}$ 的代价时，只需对式(18)作二重积分或使用蒙特卡洛算法即可。

(4)此时，已经求得了MDP框架下所有要素的值。可以使用MDP值迭代^[17]的方法，求解各个离散状态对应的期望最小损失。离散状态 ${{\tilde{ S}}_d}$ 对应的期望最小损失记作 $V\!\, \raisebox{4pt}{\hat} \ \left({{\tilde{ S}}_d}\right)$ 。

(5)以上步骤均为离线计算。在线决策时，如当前置信状态为 ${{{b}}_k}$ ，根据粒子滤波置信状态迁移得到 ${{{b}}_{k{\rm{ + }}1}}$ 的粒子集合。根据连续状态离散化的结果，将每一个粒子 ${{S}}_{k + 1}^{(i)}$ 对应到相应的离散状态 ${\tilde{ S}}_{k + 1}^{(i)}$ 中。然后可基于 $V\!\, \raisebox{4pt}{\hat} \ \left({{\tilde{ S}}_d}\right)$ ，计算出对应的未来损失的期望：

求得未来损失后，代入式(5)，由算法1第8行即可计算置信状态 ${{{b}}_k}$ 对应的Q-value的近似值，然后选取使得Q-value最小的动作 ${{{A}}_k}$ 作为当前最优策略。

整体算法过程如表1的算法1所示。

表  1  C-QMDP算法

算法 1　C-QMDP算法 　输入：置信状态(包括粒子状态和粒子权重)： 　 ${{{b}}_k} = \left({{S}}_k^{(1)},\tilde w_k^{(1)},{{S}}_k^{(2)},\tilde w_k^{(2)},·\!·\!·,{{S}}_k^{(N)},\tilde w_k^{(N)}\right)$ 　输出：最优动作A 　(1) function C-QMDP( ${{{b}}_k}$) 　(2) $V\!\, \raisebox{4pt}{\hat} \leftarrow {\rm{MDP\_discrete\_value\_iteration()}}$ 　(3)　 for all control actions A do 　(4) ${{{b}}_{k + 1}} \leftarrow {\rm{Particle\_filter(}}{{{b}}_k},{{A}}{\rm{)}}$ 　(5) for i = 1 : N do 　(6) ${\tilde{ S}}_{k + 1}^{(i)} \leftarrow {{S}}_{k + 1}^{(i)}$ 　(7) end for 　(8) $Q({{{b}}_k},{{A}}) = L({{{b}}_k},{{A}}) + \beta \sum\limits_{i = 1}^N {\tilde w_k^{(i)}V\!\, \raisebox{4pt}{\hat}\left({\tilde{ S}}_{k + 1}^{(i)}\right)}$ 　(9) end for 　(10) return $\arg \mathop {\max }\limits_{{A}} Q({{{b}}_k},{{A}})$ 　(11) end function

下载: 导出CSV 

| 显示表格

4.   仿真结果

为了验证基于C-QMDP近似Q-value的节点调度算法的性能，本文与传感器全部开启、文献[27]的最近点方法(Closest Point Approach, CPA)算法、文献[21]的CO-Rollout算法进行仿真比较。选用CO-Rollout算法比较的原因是该算法适用于本文中观测方程与代价函数均为非线性的场景，且同样是将本文场景建模为POMDP问题并求解。传感器全部开启的方法指的是在每个时刻将所有的传感器开启，对目标状态进行观测。CPA是一种“贪婪”算法，在每个决策时间，根据对目标状态的估计，选择距离目标最近的m个传感器对目标进行观测。由于观测能量信号强度与目标和传感器间距离的 $\lambda$ 次方成反比，因此选择距离目标较近的传感器进行观测是合理的。但是，CPA也有其不合理性，一方面它没有考虑传感器间的差异，不同传感器具有不同的观测精度和使用功耗，另一方面没有考虑长期损失，只考虑了当前时刻的观测效果。

仿真中，将M=20个观测精度不同的传感器随机分布在160 m×120 m的平面区域内，对移动目标进行跟踪。x轴坐标范围为[–80, 80], y轴坐标范围为[0, 160]。目标根据式(1)定义的转移模型移动，初始位置为(2,7)，初始速度为(1,2)，方差 $\sigma _x^2$ =0.3, $\sigma _y^2$ =0.5。决策周期 ${T_s}$ =2 s。目标初始能量强度 ${E_0}$ =1000, $\lambda$ =2。在粒子滤波中，使用N=10000个粒子近似目标的后验分布，更新目标的置信状态。粒子滤波初始化时，从均值为0，方差为3的高斯分布中采样N个粒子，每个粒子是一个包含位置和速度的4维向量。

通过仿真验证C-QMDP算法的跟踪性能，跟踪结果如图2所示。可以看出，当目标在平面区域随机移动时，C-QMDP算法能够在误差允许的情况下，实现对目标位置的实时跟踪。

图  2  目标真实轨迹与估计轨迹比较

下载: 全尺寸图片 幻灯片

接下来，比较不同算法下目标跟踪的累积误差，计算50次后取平均，结果如图3所示。跟踪误差的计算公式为 $\mathbb{E}\left[{\left\| {{{{\hat{ p}}}_s} - {{{p}}_s}} \right\|^2}\right]$ , k时刻累积误差 $\varDelta$ 为时刻1到k跟踪误差的和。可见，传感器全部开启时，拥有较高的跟踪精度。C-QMDP算法近似Q-value进行节点调度，使用较少的传感器组合对目标进行观测，跟踪效果接近传感器全部开启的跟踪精度。CO-Rollout算法计算未来代价时使用基本策略代替最优策略进行估计，跟踪效果略差于C-QMDP算法的跟踪精度。CPA算法只是贪婪地选取距离目标最近的传感器，没有考虑每个传感器的特性和长期损失，因此跟踪精度较其余算法更差。不过该算法的实现和计算都十分简单，易于使用。

节点调度算法的目的是对跟踪误差和能量消耗进行权衡。在式(3)中，使用 $\alpha$ 对两者的量纲进行统一，并对重要性进行权衡。 $\alpha$ 是一个预先给定的常数，取值应在合理的区间内。仿真中，取 $\alpha$ =20，使得跟踪误差和传感器功耗对最终决策的影响大致相当，从而可以较好地体现节点调度算法的作用。传感器每一步的代价为0.1～0.8，不同传感器具有不同的使用代价。比较不同算法目标跟踪过程中的累积总代价L，时刻k的累积总代价表示时刻1到时刻k总代价的和，计算50次取平均，结果如图4所示。

图  3  目标跟踪过程中的累积跟踪误差

下载: 全尺寸图片 幻灯片

图  4  目标跟踪过程中的累积总代价

下载: 全尺寸图片 幻灯片

对文中各算法的性能进行定量分析，跟踪误差的计算公式为 $\mathbb{E}\left[{\left\| {{{{\hat{ p}}}_s} - {{{p}}_s}} \right\|^2}\right]$ ，传感器功耗的计算公式为 $\displaystyle\sum\nolimits_{m = 1}^M {c_m^{{\rm{on}}}{a_m}}$ ，总代价的计算见式(3)，仿真中取 $\alpha$ =20，计算50次后取平均，各项指标在整个跟踪过程的累积值的比较如表2所示。

表  2  C-QMDP与其余算法性能比较

算法	跟踪误差	传感器功耗	总代价
全部开启	1.3647	70.0722	97.3662
CPA	2.9132	24.0217	82.2857
CO-Rollout	1.8371	20.3493	57.0913
C-QMDP	1.5542	20.0906	51.1746

下载: 导出CSV 

| 显示表格

由图4和表2可知，虽然将所有传感器全部开启可以获得较高的跟踪精度，但是会带来很大的能量消耗，导致整体的代价很高，累积总代价上涨速度较快。当迭代到第5次后，全部开启的累积总代价超过了CPA算法。通过CO-Rollout和C-QMDP近似Q-value的节点调度算法每次选择了较少的传感器组合进行目标跟踪，能量消耗明显小于全部开启的方法，因此保持了较小的累积总代价。由于CO-Rollout算法计算未来代价时使用基本策略代替最优策略进行估计，估计误差较C-QMDP算法更大，因此总代价高于基于C-QMDP近似Q-value的节点调度算法。

图  5  不同时刻节点调度策略比较

下载: 全尺寸图片 幻灯片

比较不同算法各时刻的节点调度策略，如图5所示。可以看出，在每个时刻，CPA算法根据对目标位置的估计，选择最近的m个传感器。在各个时刻不同的算法选择了不同数量和位置的传感器子集对目标进行观测，其中C-QMDP算法往往可以选择更加合理的传感器子集。同时，C-QMDP算法可以根据目标位置估计的精确程度，灵活调整传感器子集的数量，当目标位置估计较为精确时，将选用较少数量的传感器以减少功耗；反之，则选用较多数量的传感器以降低跟踪误差。

此外，由于C-QMDP算法将未来代价的计算在线下进行，节约了在线计算时间。因此在相同的粒子数目、传感器数目的情况下，C-QMDP算法相比CO-Rollout算法需要更少的计算资源。

5.   结束语

本文将无线传感器网络中的目标跟踪问题建模为部分可观测马尔科夫决策过程(POMDP)，并提出了C-QMDP算法对非线性代价函数的POMDP问题进行近似求解。C-QMDP算法基于粒子滤波推导当前置信状态的转移，然后根据置信状态计算当前损失；计算未来损失时，将目标状态离散化，使用MDP框架的值迭代方法求解最小期望损失。C-QMDP算法一方面提高了计算准确度，得到了更好的传感器调度策略，另一方面将未来损失的计算转移到线下进行，减小了在线决策时的计算量。仿真结果表明，C-QMDP算法在计算准确度和在线计算时间方面均比CPA和CO-Rollout算法有所改进。

值得注意的是，C-QMDP算法也可适用于其他类似的观测方程和代价函数非线性的序贯决策问题中。如果观测方程线性且噪声是高斯噪声，可以使用卡尔曼滤波代替文中的粒子滤波推导置信状态的迁移，进一步降低计算量和准确性。对于连续状态的离散化，本文使用的是等间隔划分的方法，后续可以研究是否有更好的离散化方法，或采用近似方法直接求解连续状态空间的MDP问题。此外，可以对代理代价函数做更多的尝试，选择更能准确近似真实Q-value的代理代价函数。