0 引言

在基于无线传感器网络 (Wireless Sensor Network, WSN) 的应用中，传感器节点检测到的信息若没有准确的位置信息将变得毫无价值^[1]。但考虑到成本和能量限制，一般只有少数参考节点的位置是已知的，其他大部分节点 (称为待定位节点) 通过邻近的参考节点的位置和与其之间的距离等观测信息确定自己的位置，参考节点较少时定位误差较大。近年来出现的协作定位技术中，待定位节点之间也进行通信，从而可有效提高定位精度，引起国内外学者的广泛关注，并提出了许多协作定位算法，如最小二乘 (Least Square, LS)^[2]、最大似然 (Maximum Likelihood, ML)^[3]、多维尺度MAP (MultiDimensional Scaling MAP, MDS-MAP)^[4]以及卡尔曼滤波^[5]和粒子滤波^[6]等。

在近似贝叶斯推理的实现方法中，基于因子图的消息传递算法特别适用于大规模概率模型和分布式运算^[7]，其中置信传播 (Belief Propagation, BP) 方法在树图下可得到精确的边缘概率密度分布函数，在有环图下也常常能获得较好的性能。近年有学者将BP应用于无线网络协作定位中，但非线性测距模型会导致消息计算非常复杂。一种解决方法是使用非参数化消息。Ihler等^[8]将非参数化置信传播 (Nonparametric Belief Propagation, NBP) 用于静态无线传感器网络节点定位，Wymeersch等^[9]进一步提出了适用于动态网络的SPAWN (Sum-Product Algorithm over a Wireless Network) 算法。基于NBP的定位算法使用大量粒子表示消息，能灵活表示多种形式的消息，当粒子数足够多时可获得很高的定位精度，但计算复杂度和通信开销非常高^[10]。Li等^[11]提出一种基于序贯粒子的BP算法，利用吉布斯采样对消息进行粒子化，有效降低了计算复杂度，但通信负载依然较高。另一种解决方法是将测距模型线性化并采用参数化消息。北京理工大学武楠等利用一阶泰勒展开将测距模型线性化，提出了高斯BP节点定位算法^[12-13]。相对于BP方法，平均场 (Mean Field, MF) 方法更适用于指数型消息，消息更新规则更加简单。Pedersen等^[14]提出一种基于MF的静态网络节点定位算法，该算法采用高斯型消息，并使用最小化KL散度 (Kullback-Leibler Divergence, KLD) 将非高斯消息近似为高斯消息，但近似过程中出现了第一类合流超几何函数，致使计算复杂度非常高，且估计性能与基于BP的算法相比有所下降。

针对上述问题，本文将高性能的BP方法和低复杂度的MF方法结合起来，提出一种低复杂度低通信开销的协作定位算法，并考虑参考节点位置可能存在误差，以避免误差累积引起的定位精度下降问题。

1 系统模型与因子图

考虑二维平面上包含参考节点和待定位节点的无线传感器网络，定义符号A和M分别表示参考节点和待定位节点的集合。所有节点在区域内可独立随机地移动，k时刻，节点i的位置变量记作x_i^k$ \buildrel \Delta \over = $[x_i1^k, x_i2^k]^T。k=0时，假设x_i⁰的先验概率密度函数 (Probability Density Function, PDF) p(x_i⁰) 是高斯的，其均值为μ_i⁰，协方差矩阵为V_i⁰$ \buildrel \Delta \over = $σ_{i, 0}²·I_2×2(其中I_2×2是单位矩阵)，则可记作p(x_i⁰)=N (x_i⁰; μ_i⁰, V_i⁰)。

假设节点i∈M的状态转移模型为：

$ {\boldsymbol{x}}_i^k = {\boldsymbol{x}}_i^{k - 1} + {\boldsymbol{s}}_i^k \cdot \Delta T + {\boldsymbol{n}}_i^k $

(1)

其中:s_i^k是速度矢量，ΔT是时间间隔，n_i^k是均值为0、方差为σ²_ni^k·I_2×2的系统噪声，记作n_i^k~N (n_i^k; 0, σ²_ni^k·I_2×2)。由此可得状态转移函数为：

$ p({\boldsymbol{x}}_i^k|{\boldsymbol{x}}_i^{k - 1}) = {\rm{N}}({\boldsymbol{x}}_i^k;{\kern 1pt} {\kern 1pt} {\kern 1pt} {\boldsymbol{x}}_i^{k - 1} + {\boldsymbol{s}}_i^k \cdot \Delta T,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \sigma _{{\boldsymbol{n}}_i^k}^2 \cdot {{\boldsymbol{I}}_{2 \times 2}}) $

(2)

节点间的通信和测距半径记为L。k时刻，若待定位节点i和节点j(j∈A∪M) 间的距离小于L，称两者互为邻居，并记作j∈N_i^k(N_i^k表示k时刻节点i邻居节点标号的集合)。为便于讨论，假设节点间的通信和测距链路是对称的，若节点j时节点i的邻居，则节点i也是节点j的邻居，即j∈N_i^k$ \Leftrightarrow $i∈N_j^k。k时刻，节点i测得的与节点j之间的距离为：

$ d_{i \leftarrow j}^k = \left\| {{\boldsymbol{x}}_j^k - {\boldsymbol{x}}_i^k} \right\| + \omega _{ij}^k $

(3)

其中ω_ij^k~N (ω_ij^k; 0, σ_{ij, k}²) 是测量噪声。则已知x_i^k和x_j^k时d_i←j^k的PDF (也称为似然函数) 为：

$ \begin{array}{l} p(d_{i \leftarrow j}^k|{\boldsymbol{x}}_i^k,{\boldsymbol{x}}_j^k) = {\rm{N}}\left( {d_{i \leftarrow j}^k;{\kern 1pt} {\kern 1pt} {\kern 1pt} \left\| {{\boldsymbol{x}}_j^k - {\boldsymbol{x}}_i^k} \right\|,{\kern 1pt} {\kern 1pt} {\kern 1pt} \sigma _{ij,k}^2} \right)\\ \;\;\;\; \propto \exp \left\{ {{{(d_{i \leftarrow j}^k - \left\| {{\boldsymbol{x}}_j^k - {\boldsymbol{x}}_i^k} \right\|)}^2}/2\sigma _{ij,k}^2} \right\} \end{array} $

(4)

k时刻，所有节点位置变量和测量距离的集合分别记作X^k$ \buildrel \Delta \over = ${x_i^k:i∈A∪M}和D^k$ \buildrel \Delta \over = ${d_i←j^k:∀i∈M, j∈N_i^k}，并将0时刻到K时刻所有节点位置变量和测量距离的集合分别记作X^0:K$ \buildrel \Delta \over = ${X⁰, X¹, …, X^K}和D^1:K$ \buildrel \Delta \over = ${D¹, D², …, D^K}。根据贝叶斯准则，从0时刻到K时刻节点位置变量的联合后验PDF可表示为：

$ \begin{array}{l} p({X^{0:K}}|{D^{1:K}}) \propto p({X^{0:K}})p({D^{1:K}}|{X^{0:K}}) = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \prod\limits_{k = 1}^K {\prod\limits_{a \in A} {p({\boldsymbol{x}}_a^k)} } \prod\limits_{i \in M} {p({\boldsymbol{x}}_i^0)} p({\boldsymbol{x}}_i^k|{\boldsymbol{x}}_i^{k - 1}) \times \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \prod\limits_{j \in N_i^k} {p(d_{i \leftarrow j}^k|{\boldsymbol{x}}_i^k,{\boldsymbol{x}}_j^k)} \end{array} $

(5)

其中p(x_a^k) 是k时刻参考节点a位置的先验PDF，其均值为已知的测量位置${\boldsymbol{\tilde x}}_a^k$，协方差矩阵为${\boldsymbol{\tilde V}}_a^k = \sigma _{a,k}^2 \cdot {{\boldsymbol{I}}_{2 \times 2}}$，则有$p({\boldsymbol{x}}_a^k) = {\rm{N}}({\boldsymbol{x}}_a^k;{\boldsymbol{\tilde x}}_a^k,{\boldsymbol{\tilde V}}_a^k)$。

因子图作为图模型的一种，可根据全局函数的因子分解形象地描述变量与因子之间的关系。由式 (5) 的因子分解可得到图 1所示的因子图，图中包含两类节点，即函数节点 (矩形) 和变量节点 (圆圈)。每个函数节点和与之有关的变量节点相连，其中$f_i^{k|k - 1} \buildrel \Delta \over = p({\boldsymbol{x}}_i^k|{\boldsymbol{x}}_i^{k - 1})$，$f_{ij}^k \buildrel \Delta \over = p(d_{i \leftarrow j}^k|{\boldsymbol{x}}_i^k,{\boldsymbol{x}}_j^k)$。消息传递算法在因子图上定义了两类消息：变量节点到函数节点的消息和函数节点到变量节点的消息。不同的消息传递算法中，两类消息的计算式不同。

通常情况下因子图是无向图，但由于网络中节点间的连通关系会随时间变化，因此不计算从当前时刻到过去时刻的消息。此外，本算法考虑参考节点的位置可能有一定误差的情况，将其看作变量。在各时刻的迭代计算中，参考节点的位置也可以更新。

2 基于BP-MF的协作定位算法

消息传递算法按照一定准则在因子图上以迭代的方式计算边缘后验的近似，即变量的置信。某一变量的置信定义为所有相关函数节点到该变量消息的乘积。由图 1可知，x_i^k的置信b(x_i^k) 由两类消息构成：一是来自状态转移函数f_i^k|k-1的消息，这里称为预测消息，二是来自似然函数f_ij^k和f_ji^k(其中j∈N_i^k) 的协作消息。由式 (1) 可知，节点的状态转移模型是线性的，且只与自身的运动速度和前一时刻的位置有关；由式 (3) 可知，节点间的测距模型是非线性的，与自身和邻居节点的位置有关。基于此，将定位过程分为基于BP方法的预测消息计算和基于MF方法的协作消息计算。为降低通信开销，所有消息均约束为高斯函数，对非线性测距模型引起的非高斯置信，利用泰勒展开将其近似为高斯消息。

2.1 基于BP方法的预测消息计算

根据BP准则，函数f_r(x_r)(其中x_r={x₁, x₂, …, x_i, …}是与该函数相关变量的集合) 到其相关变量x_i的消息是：f_r(x_r) 乘以除x_i外的其他相关变量到f_r(x_r) 的消息，再对这些变量求积分 (即对x_i边缘化)。状态转移函数f_i^k|k-1的相关变量为x_i^k和x_i^k-1，因此，f_i^k|k-1到x_i^k的消息为：

$ {m_{f_i^{k|k - 1} \to {\boldsymbol{x}}_i^k}}({\boldsymbol{x}}_i^k) = \int {p({\boldsymbol{x}}_i^k|{\boldsymbol{x}}_i^{k - 1}){m_{{\boldsymbol{x}}_i^{k - 1} \to f_i^{k|k - 1}}}({\boldsymbol{x}}_i^{k - 1})} d{\boldsymbol{x}}_i^{k - 1} $

(6)

其中：${m_{{\boldsymbol{x}}_i^{k{\rm{ - 1}}} \to f_i^{k|k - 1}}}({\boldsymbol{x}}_i^k)$是变量节点x_i^k-1到函数节点f_i^k|k-1的消息。由于不计算从当前时刻到过去时刻的消息，因此${m_{{\boldsymbol{x}}_i^{k{\rm{ - 1}}} \to f_i^{k|k - 1}}}({\boldsymbol{x}}_i^k)$等于x_i^k-1的置信b(x_i^k-1)。假设b(x_i^k-1) 的高斯近似记作$\hat b({\boldsymbol{x}}_i^{k - 1})$，其均值和协方差矩阵分别为${\boldsymbol{\hat x}}_i^{k - 1}$和${\boldsymbol{\hat V}}_i^{k - 1}$，则有：

$ {m_{{\boldsymbol{x}}_i^{k - 1} \to f_i^{k|k - 1}}}({\boldsymbol{x}}_i^{k - 1}) = {\rm{N}}({\boldsymbol{x}}_i^{k - 1};{\kern 1pt} {\kern 1pt} {\boldsymbol{\hat x}}_i^{k - 1},{\boldsymbol{\hat V}}_i^{k - 1}{\kern 1pt} ) $

(7)

将式 (2) 和式 (7) 和代入式 (6)，可得预测消息${m_{f_i^{k|k - 1} \to {\boldsymbol{x}}_i^k}}({\boldsymbol{x}}_i^k)$为：

$ \begin{array}{l} {m_{f_i^{k|k - 1} \to {\boldsymbol{x}}_i^k}}({\boldsymbol{x}}_i^k){\kern 1pt} {\kern 1pt} = \int {p({\boldsymbol{x}}_i^k|{\boldsymbol{x}}_i^{k - 1})} \cdot \\ \;\;\;\;{\rm{N}}({\boldsymbol{x}}_i^{k - 1};{\kern 1pt} {\kern 1pt} {\kern 1pt} {\boldsymbol{\hat x}}_i^{k - 1},{\kern 1pt} {\kern 1pt} {\boldsymbol{\hat V}}_i^{k - 1})d{\boldsymbol{x}}_i^{k - 1} \propto \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \exp \left( { - \frac{1}{2}{{({\boldsymbol{x}}_i^k - {\boldsymbol{\tilde x}}_i^k)}^{\rm{T}}}{{({\boldsymbol{\tilde V}}_i^k{\kern 1pt} {\kern 1pt} )}^{ - 1}}({\boldsymbol{x}}_i^k - {\boldsymbol{\tilde x}}_i^k)} \right) \end{array} $

(8)

其中：${\boldsymbol{\tilde x}}_i^k = {\kern 1pt} {\boldsymbol{\hat x}}_i^{k - 1} + {\boldsymbol{s}}_i^k \cdot \Delta T$，${\boldsymbol{\tilde V}}_i^k{\kern 1pt} = {\boldsymbol{\hat V}}_i^{k - 1} \cdot {\kern 1pt} {\kern 1pt} \sigma _{{\boldsymbol{n}}_i^k}^2$。

2.2 基于MF方法的协作消息计算

k时刻，x_i^k接收来自f_ij^k和f_ji^k的消息${m_{f_{ij}^k \to {\boldsymbol{x}}_i^k}}({\boldsymbol{x}}_i^k)$和${m_{f_{ji}^k \to {\boldsymbol{x}}_i^k}}({\boldsymbol{x}}_i^k)$，其中j∈N_i^k。按照VMP准则，某一函数f_r(x_r)(其中x_r={x₁, x₂, …, x_i, …}是与该函数相关变量的集合) 到其相关变量x_i的消息定义为：f_r(x_r) 取自然对数后，与除x_i外的其他相关变量到f_r(x_r) 的消息相乘，再对x_i边缘化，最后取exp。函数f_ij^k的相关变量为x_i^k和x_j^k，因此，第t次迭代中, f_ij^k到x_i^k的消息$m_{f_{ij}^k \to {\boldsymbol{x}}_i^k}^t({\boldsymbol{x}}_i^k)$为：

$ m_{f_{ij}^k \to {\boldsymbol{x}}_i^k}^t({\boldsymbol{x}}_i^k){\kern 1pt} = \exp \left( {\int {{{\hat b}^{t - 1}}({\boldsymbol{x}}_j^k) \cdot \ln p(d_{i \leftarrow j}^k|{\boldsymbol{x}}_i^k,{\boldsymbol{x}}_j^k){\rm{d}}{\boldsymbol{x}}_j^k} } \right) $

(9)

其中${\hat b^{t - 1}}({\boldsymbol{x}}_j^k)$是上次迭代时x_j^k的置信b^t-1(x_j^k) 的高斯近似。

将式 (4) 带入式 (9) 可得：

$ \begin{array}{l} m_{f_{ij}^k \to {\boldsymbol{x}}_i^k}^t({\boldsymbol{x}}_i^k) = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \exp \left( {\int {{{\hat b}^{t - 1}}({\boldsymbol{x}}_j^k) \cdot \frac{{{{(d_{i \leftarrow j}^k - \left\| {{\boldsymbol{x}}_j^k - {\boldsymbol{x}}_i^k} \right\|)}^2}}}{{2\sigma _{ij,k}^2}}{\rm{d}}{\boldsymbol{x}}_j^k} } \right) \end{array} $

(10)

同理可得：

$ \begin{array}{l} m_{f_{ji}^k \to {\boldsymbol{x}}_i^k}^t({\boldsymbol{x}}_i^k) = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \exp \left( {\int {{{\hat b}^{t - 1}}({\boldsymbol{x}}_j^k) \cdot \frac{{{{(d_{j \leftarrow i}^k - \left\| {{\boldsymbol{x}}_i^k - {\boldsymbol{x}}_j^k} \right\|)}^2}}}{{2\sigma _{ji,k}^2}}{\rm{d}}{\boldsymbol{x}}_j^k} } \right) \end{array} $

(11)

2.3 置信的计算和近似

k时刻，到变量节点x_i^k的消息有预测消息${m_{f_i^{k|k - 1} \to {\boldsymbol{x}}_i^k}}({\boldsymbol{x}}_i^k)$和协作消息$m_{f_{ij}^k \to {\boldsymbol{x}}_i^k}^t({\boldsymbol{x}}_i^k)$、$m_{f_{ji}^k \to {\boldsymbol{x}}_i^k}^t({\boldsymbol{x}}_i^k)$(其中j∈N_i^k)。因此，由式 (8)、式 (10) 和式 (11) 可得第t次迭代得到的x_i^k的置信b^t(x_i^k) 为：

$ \begin{array}{l} {b^t}({\boldsymbol{x}}_i^k){\kern 1pt} \propto \exp \left( { - \frac{1}{2}({\boldsymbol{x}}_i^k - {\boldsymbol{\tilde x}}_i^k){{({\boldsymbol{\tilde V}}_i^k{\kern 1pt} {\kern 1pt} )}^{ - 1}}({\boldsymbol{x}}_i^k - {\boldsymbol{\tilde x}}_i^k) + } \right.\\ {\kern 1pt} \;\;\;\;\;\;\;\;\sum\limits_{j \in N_i^k} {\int {{{\hat b}^{t - 1}}({\boldsymbol{x}}_j^k) \cdot \frac{{{{(d_{i \leftarrow j}^k - \left\| {{\boldsymbol{x}}_j^k - {\boldsymbol{x}}_i^k} \right\|)}^2}}}{{2\sigma _{ij,k}^2}}{\rm{d}}{\boldsymbol{x}}_j^k} } + \\ {\kern 1pt} \;\;\;\;\;\;\;\;\left. {\sum\limits_{j \in N_i^k} {\int {{{\hat b}^{t - 1}}({\boldsymbol{x}}_j^k) \cdot \frac{{{{(d_{j \leftarrow i}^k - \left\| {{\boldsymbol{x}}_i^k - {\boldsymbol{x}}_j^k} \right\|)}^2}}}{{2\sigma _{iji,k}^2}}{\rm{d}}{\boldsymbol{x}}_j^k} } } \right) \end{array} $

(12)

由式 (12) 计算得到的b^t(x_i^k) 是关于x_i^k的非高斯函数，为便于在节点间传输，本算法将所有消息约束为高斯函数。由于高斯函数的均值和方差分别与指数项的一次项和二次项有关，而二阶泰勒级数展开可以很方便地得到二次型函数。因此，将指数项中的非线性项f(x_i^k, x_j^k)$ \buildrel \Delta \over = $‖x_j^k-x_i^k‖=‖x_i^k-x_j^k‖在上一次迭代点 (${\boldsymbol{\tilde x}}_i^k,{\boldsymbol{\tilde x}}_j^k$) 处进行二阶泰勒级数展开，即：

$ \begin{array}{l} f({\boldsymbol{x}}_i^k,{\boldsymbol{x}}_j^k) \approx f({\boldsymbol{\tilde x}}_i^k,{\boldsymbol{\tilde x}}_j^k) + \frac{{{\partial ^{\rm{T}}}f({\boldsymbol{\tilde x}}_i^k,{\boldsymbol{\tilde x}}_j^k)}}{{\partial {\boldsymbol{x}}_i^k}}({\boldsymbol{x}}_i^k - {\boldsymbol{\tilde x}}_i^k) + \\ {\kern 1pt} \;\;\;\;{\kern 1pt} \frac{{{\partial ^{\rm{T}}}f({\boldsymbol{\tilde x}}_i^k,{\boldsymbol{\tilde x}}_j^k)}}{{\partial {\boldsymbol{x}}_j^k}}({\boldsymbol{x}}_j^k - {\boldsymbol{\tilde x}}_j^k) + \\ {\kern 1pt} \;\;\;\;\frac{1}{2}{\left( \begin{array}{l} {\boldsymbol{x}}_i^k - {\boldsymbol{\tilde x}}_i^k\\ {\boldsymbol{x}}_j^k - {\boldsymbol{\tilde x}}_j^k \end{array} \right)^{\rm{T}}}{\nabla ^2}f({\boldsymbol{\tilde x}}_i^k,{\boldsymbol{\tilde x}}_j^k)\left( \begin{array}{l} {\boldsymbol{x}}_i^k - {\boldsymbol{\tilde x}}_i^k\\ {\boldsymbol{x}}_j^k - {\boldsymbol{\tilde x}}_j^k \end{array} \right) \end{array} $

(13)

将式 (13) 代入式 (12) 并整理后可得b^t(x_i^k) 的高斯近似函数${\hat b^t}({\boldsymbol{x}}_i^k) = {\rm{N}}\left( {{\boldsymbol{x}}_i^k;{\kern 1pt} {\kern 1pt} {\boldsymbol{\hat x}}_i^k{\kern 1pt} ,{\kern 1pt} {\kern 1pt} {\boldsymbol{\hat V}}_i^k} \right)$，其中${\boldsymbol{\hat x}}_i^k$和${\boldsymbol{\hat V}}_i^k$为：

$ \begin{array}{l} {\boldsymbol{\hat V}}_i^k = \left[ {{{({\boldsymbol{\tilde V}}_i^k{\kern 1pt} {\kern 1pt} )}^{ - 1}} + \sum\limits_{j \in N_i^k} {(\frac{1}{{\sigma _{ij,k}^2}} + \frac{1}{{\sigma _{ji,k}^2}}){\boldsymbol{I}} - } } \right.\\ {\left. {\;\;\;\;\sum\limits_{j \in N_i^k} {(\frac{{d_{i \leftarrow j}^k}}{{\sigma _{ij,k}^2}} + \frac{{d_{j \leftarrow i}^k}}{{\sigma _{ji,k}^2}}){\nabla ^2}f({\boldsymbol{\tilde x}}_i^k,{\boldsymbol{\tilde x}}_j^k)} } \right]^{ - 1}} \end{array} $

(14)

$ \begin{array}{l} {\boldsymbol{\hat x}}_i^k = {\boldsymbol{\hat V}}_i^k\{ {({\boldsymbol{\tilde V}}_i^k{\kern 1pt} {\kern 1pt} )^{ - 1}} + \sum\limits_{j \in N_i^k} {[1/\left( {\sigma _{ij,k}^2} \right) + 1/\left( {\sigma _{ji,k}^2} \right)] \cdot {\boldsymbol{\tilde x}}_j^k} + \\ \;\;\;\;\sum\limits_{j \in N_i^k} {\frac{{d_{i \leftarrow j}^k}}{{\sigma _{ij,k}^2}}[\frac{{\partial f({\boldsymbol{\tilde x}}_i^k,{\boldsymbol{\tilde x}}_j^k)}}{{\partial {\boldsymbol{x}}_i^k}} - {\nabla ^2}f({\boldsymbol{\tilde x}}_i^k,{\boldsymbol{\tilde x}}_j^k) \cdot {\boldsymbol{\tilde x}}_j^k)]} + \\ {\kern 1pt} \;\;\;\;\sum\limits_{j \in N_i^k} {\frac{{d_{j \leftarrow i}^k}}{{\sigma _{ji,k}^2}}[\frac{{\partial f({\boldsymbol{\tilde x}}_i^k,{\boldsymbol{\tilde x}}_j^k)}}{{\partial {\boldsymbol{x}}_i^k}} - {\nabla ^2}f({\boldsymbol{\tilde x}}_i^k,{\boldsymbol{\tilde x}}_j^k) \cdot {\boldsymbol{\tilde x}}_j^k)]} \} \end{array} $

(15)

最后一次迭代得到的${\hat b^t}({\boldsymbol{x}}_i^k)$即为节点i的边缘后验概率分布函数p(x_i^k|D^0:K)，最后根据最大后验估计准则得到节点i的位置估计。

2.4 算法执行流程

基于BP方法的预测过程是无环图，因此预测消息不进行迭代更新，基于MF方法的协作消息以及每个节点位置变量的置信需要进行迭代计算。MF方法总是保证良好的收敛性能^[15]，本文算法在先验较好的情况下收敛较快，先验较差时收敛较慢。

本文算法的执行流程如下所示。

1) 初始化先验p(x_i⁰)。

2) 各节点根据式 (8) 计算预测消息${m_{f_i^{k|k - 1} \to {\boldsymbol{x}}_i^k}}({\boldsymbol{x}}_i^k)$。

3) 迭代开始，各节点并行执行：

a) 广播自己的位置和测量距离并接收邻居点的信息；

b) 由式 (14) 和 (15) 计算${\hat b^t}({\boldsymbol{x}}_i^k)$的均值和协方差矩阵。

4) 迭代终止。

5) 各节点根据最大后验估计准则确定估计位置。

6) 算法结束。

3 仿真分析

以下通过数值仿真评估本文算法的性能。考虑与文献[9]相同的场景，即包含13个参考节点和100个待定位节点的100 m×100 m的二维平面。各节点均可独立随机移动，其中参考节点的位置已知，但存在一定的不确定性。仿真时，设置k时刻节点i状态转移噪声标准差${\sigma _{{\boldsymbol{n}}_i^k}}$和测距噪声标准差σ_{ij, k}均为1 m，两次定位之间的时间间隔ΔT=1 s，运动速度s_i^k=[s_i1^k, s_i2^k]^T(其中s_i1^k, s_i2^k∈[-1 m/s, 1 m/s])，迭代次数设置为10。由于测量设备的误差，参考节点的测量位置可能存在一定的误差，即不确定性。本文算法考虑到这种不确定性，将参考节点的位置看作高斯型变量，其不确定性用σ_{a, k}表示，并将参考节点位置先验概率分布函数的协方差矩阵记为${\boldsymbol{\tilde V}}_a^k = \sigma _{a,k}^2 \cdot {{\boldsymbol{I}}_{2 \times 2}}$。仿真中，使用误差累积分布函数 (Cumulative Distribution Function, CDF)，即定位误差小于r(单位: m) 的概率，来衡量算法的性能。

当σ_{a, k}=1 m时，即参考节点位置的误差在0~1 m时，考虑其不确定性与不考虑其不确定性两种情况下的CDF曲线如图 2所示。可以看出，考虑参考节点位置可能存在的误差能够避免误差累积，从而获得更好的定位精度。

不同σ_{a, k}下定位误差累积分布函数如图 3所示。

由图 3可知，σ_{a, k}越小，即参考节点测量位置的误差越小，定位性能越好。此外，随着定位时刻k的增加，定位精度有所提高，这是因为随着k的增加待定位节点的先验信息 (即预测信息) 准确度提高。但当k大于一定值时 (实验中k > 10)，定位精度基本稳定，不会有大幅度提高。

为更好地评估算法性能，在此将本文算法 (记作BP-MF) 与目前性能最优的SPAWN算法进行比较。仿真时设置σ_a^k=0.25 m，SPAWN算法中每个消息的粒子数记为N_p，并考虑N_p=1 000和N_p=2 500两种情况，分别记作SPAWN1000和SPAWN2500。由图 4可知，k=1时，本文算法的定位精度与SPAWN1000和SPAWN2500有较大差别；但随着k的增加，差别变得很小，k=20时BP-MF的定位精度高于SPAWN1000，且与SPAWN25000接近。

当粒子数足够多时，虽然SPAWN算法的估计性能优于本文算法，但其通信开销和计算复杂度均非常高。具体地，假设k时刻节点i有N_{i, k}个邻居节点。在通信开销方面，在基于粒子的SPAWN中，节点i要向邻居节点发送自身的位置，则需要发送N_p个位置向量和N_p个权重标量，并接收N_{i, k}个邻居节点的位置信息。而在本文算法中，节点i只需向邻居节点发送1个位置向量和1个位置协方差矩阵，接收来自邻居节点的N_{i, k}个位置向量和N_{i, k}个位置协方差矩阵。对于计算复杂度，每次迭代时，SPAWN算法和本文算法中每个节点的计算复杂度分别为O(N_{i, k}·N_p²) 和O(N_{i, k})。由于N_{i, k}≪N_p，因此，与SPAWN算法相比，本文算法大幅降低了通信开销和计算复杂度。

4 结语

本文针对参考节点位置存在一定模糊性的无线传感器网络，提出一种基于BP和MF的联合消息传递协作定位算法。仿真结果和分析表明，其估计性能与基于粒子的SPAWN接近，并避免了通信开销和计算复杂度过大的问题。在本文基础上，下一步将研究基于消息传递算法的目标跟踪算法。