CN101350635A - 稀疏测量集上基于最短路径的传感器网络节点自定位方法 - Google Patents

Abstract

稀疏测量集上基于最短路径的传感器网络节点自定位方法属于无线传感器网络自定位技术领域，其特征在于，依次含有节点布撒、路由建立和数据传输、网络连接信息提取、节点相对坐标定位、绝对坐标变换、计算机输出结果等步骤，其中，在节点相对坐标定位时对于邻居节点，采用射频功率检测方法进行节点间距离测量，对于非邻居节点，用弗洛伊德最短路径距离进行距离估计，从而得到包括邻居和非邻居节点在内的最短路径距离矩阵，据此用多维标度分析方法得到进行相对定位的节点位置拓扑估计的初始值，接着求取未知距离的近似概率分布，再以完全距离矩阵的似然函数作为期望目标函数进行优化，以消除其随机性，得到相对定位结果。

Description

稀疏测量集上基于最短路径的传感器网络节点自定位方法

技术领域

本发明属于无线传感器网络自定位领域，具体包含：无线传感器网络定位技术、统计信号处理方法、数据推断技术、多维标度分析技术。

背景技术

在过去的80多年里，无线网络技术取得了突飞猛进的发展。从人工操作的无线电报网络到使用扩频技术的自动化无线局域、个域网络，无线网络的应用领域随着技术的进步不断地扩展。作为无线网络的一种典型应用和研究热点，无线传感器网络由大量震动、(地)磁、热量、视觉、红外、声音和雷达等多种不同类型传感器节点所构成，可以用于监控温度、湿度、压力、土壤构成、噪声、机械应力等多种环境条件。利用“无处不在”的传感技术，用户可以更加深入地了解和把握周围的世界。

在无线传感器网络的各种应用中，位置信息对传感器网络的监测活动至关重要。事件发生的位置，或者获取信息的节点位置都是传感器节点监测消息中所包含的重要信息。传感器网络的自定位就是根据少数已知位置的节点，按照某种定位机制确定网络中其他所有节点的位置。只有在传感器节点自身定位之后，传感器节点监测到的事件的具体位置才能够被确定，因此，在传感器网络中，传感器节点自身的正确定位是大量实际应用的前提。

目前已有的定位机制大都基于传感器节点间物理距离的测量列出多元方程组从而进行拓扑反解，例如TOA技术利用节点间的超声传播延时，RSSI利用节点间的射频信号传播强度。

利用无向图G_r，M(X，E，Z)表示无线传感器网络，其中顶点集合X代表网络中分布的传感器节点位置；E是无向边集合，代表可以互相通信的节点对之间的通信链路，在本发明中，我们假设节点间的通信能力同测距能力是等价的，即只有当节点对可以互相通信时，它们之间距离才可以被测量到；Z对应E中每条边的距离，r代表这个网络所处的空间维数，一般情况下r＝2或者3，M是网络中传感器节点的总数。在下文中，我们分别使用M阶对称矩阵W和D来表示E和Z，其中

而d_ij如果不为空项，其值就是n_i和n_j之间的距离测量值；如果为空项，则说明n_i和n_j无法互相通信。

为了得到较高精度的节点定位，需要有充足的节点对距离测量结果提供足够的信息。但是在实际传感器网络中，诸如电磁环境恶劣、节点分布稀疏、网络中存在大量障碍物等情况都会导致传感器节点间互相通信的困难，从而造成对应的节点对距离无法测量。以图1为例。在图1(a)中，50个节点分布在一个边长10米的正方形无遮挡网络中；在图1(b)中，网络中部多了射频信号无法穿透的障碍物(如阴影所示)。网络中的所有节点共享一个无线信道，节点的传输功率相同，传输距离都是2米，同时只允许视线传输(line-of-sight)。星号表示传感器节点的位置，它们之间的连接线代表对应的节点间通信链路。可以看到，此时两个网络中可以互相通信的节点对均非常少，从而使得各自的距离测量集合是都是稀疏集，如果用距离矩阵D表示，则图1所对应的D中均存在大量空项，是稀疏矩阵。

我们采用射频信号强度检测(RSS)的方法进行节点间距离测量，对于任一节点对(i，j)，假设节点i以固定发射功率P_i向节点j发射信号，通过无线信道传输模型的经验公式，可以将节点j的接收功率P_r转换成物理距离d_ij的经验数值。由于无线信道的随机性，根据[T.S.Rappaport.2001.“Wireless Communications：Principles and Practice 2nd ed.”，Pearson Education.]中的经典信道公式，该种方法得到的d_ij满足(2)式所示的概率分布，其期望值为(i，j)间的真实欧式距离δ_ij，其中方差σ_d代表测量中的不确定性。

$f\left(d_{\mathrm{ij}}\right|{\mathrm{\δ}}_{\mathrm{ij}})=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_d{d}_{\mathrm{ij}}}{e}^{-\frac{{(\log d_{\mathrm{ij}}-\log {\mathrm{\δ}}_{\mathrm{ij}})}^2}{2{\mathrm{\σ}}_d^2}}---\left(2\right)$

所以，如果获得所有的节点对距离测量值，根据(2)式的概率分布，就可以使用最大似然估计的方法，通过优化拓扑估计X使得如下似然函数最大

$L\left(X\right|D)=\underset{i,j}{\mathrm{\Σ}}\ln f\left(d_{\mathrm{ij}}\right|{\mathrm{\δ}}_{\mathrm{ij}})---\left(3\right)$

最大似然解X就可以作为拓扑估计的结果。

在传感器网络自定位中，如果距离测量集合D的稀疏程度过高，根据仅有的距离测量值利用式(3)进行拓扑反解，得到的定位结果会大大偏离节点的真实位置。因此，如果能够利用已有的网络观测信息，对未测量到的节点对距离进行推断，一定程度上补足不完全的距离矩阵D，可以为定位精度带来一定的增益。相比于传统方法，该种方法的优势在于：1.该方法是完全独立于测量模型的，即在任何测量模型假设下，都可以得到最优的定位结果；2.距离测量集合D非常稀疏的情况下，该方法可以保证较为稳定的定位性能。算法的具体细节将在下节详述。

发明内容

本发明中的算法称为基于最短路径的传感器网络自技术，英文缩写为IA-MLE，其目的在于设计一种在任何测量模型下，距离测量集合非常稀疏时均可以得到可靠定位结果的网络自定位技术。

本发明的特征在于，所述方法是在计算机中依次按以下步骤实现的：

步骤(1)，在需要进行环境监控和数据采集的区域中构建无线传感器网络：

在所述区域中布撒M个无线传感器和多个网关，所述无线传感器依次由物理传感器件、微处理器和射频模块依次串接而成，所述各无线传感器的传输功率相同，传输距离是2米，且只能视线传输，所述射频模块工作在300M～3GHz频段；

所述无线传感器网络用一个无向图G_r，M(X，E，Z)表示，r为该网络所处的空间维数，r＝2或者3，顶点集合X表示网络中分布的传感器节点位置，E为无向边集合，表示能互相通信的节点对之间的通信链路，其通信能力同测距能力等价，Z是无向边集合E中每条边的距离的集合；从而，以所述网关节点为根节点，以所述传感器节点为中间节点或叶子节点，形成多条树状路由，采用定向扩散协议建立所述无线传感器网络中的路由关系，把所述各节点的数据传回各网关节点；在所述树状路由建立阶段，在所述网关节点请求下，各传感器节点把自己下一跳传输的目的节点ID发往各网关节点，以在各网关节点内建立树状路由图，同时传输的还有该跳收发节点的距离测量值d_ij，i和j分别是两个相邻节点各自的编号，从而拼成一个距离矩阵D，是一个M阶对称矩阵；

步骤(2)，利用所述距离矩阵D依次按以下步骤对所述无线传感器网络进行相对定位：

步骤(2.1)，所述的每个节点均以固定发射功率向所述无线传感器网络广播数据包，若任一节点j收到来自节点i的数据包，表明两者之间存在一跳数据通路，两者互为邻居节点，用w_ij＝1表示，否则，为非邻居节点，用w_ij＝0表示，构成一个连接矩阵W，是一个M阶对称矩阵；

步骤(2.2)，依次按以下步骤进行相对定位：

步骤(2.2.1)，对于所述的邻居节点用射频信号强度检测测得互相之间的距离：所得的距离d_ij满足下式所示的概率分布：

$f\left(d_{\mathrm{ij}}\right|{\mathrm{\δ}}_{\mathrm{ij}})=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_d{d}_{\mathrm{ij}}}{e}^{-\frac{{(\log d_{\mathrm{ij}}-\log {\mathrm{\δ}}_{\mathrm{ij}})}^2}{2{\mathrm{\σ}}_d^2}}$

其中d_ij为节点i，j之间的距离测量值，

δ_ij是i，j间真实欧式距离的期望值，

σ_d为方差，代表测量中的不确定性，由方差σ_P根据下式所转换得到，

${\mathrm{\σ}}_d=\frac{{\mathrm{\σ}}_P}{10n}$

σ_P是节点测量接收射频信号的功率时所叠加的白噪声的方差，代表了射频功率测量的不确定性，n代表信道中的衰落指数，用于反映信号随着传输距离衰落的速度，在室内环境中，它的典型值为1.6-1.8；

上述公式表示节点间的距离测量值由于受到测量功率噪声的影响，概率密度满足以真实欧式距离为中心的对数高斯分布，

节点位置的最大似然解X在获得所有节点对距离测量值后，就是如下似然函数的最大值：

$L\left(X\right|D)=\underset{i,j}{\mathrm{\Σ}}\ln f\left(d_{\mathrm{ij}}\right|{\mathrm{\δ}}_{\mathrm{ij}})$

步骤(2.2.2)，当存在所述的非邻居节点时，用弗洛伊德最短路径算法求解它们之间的距离估计：通过遍历式搜索找出非邻居节点的所有多跳数据通路，并且选取其中单调路径距离总和最短的一条作为所述非邻居节点之间的最短路径，对应的距离就是该非邻居节点间的最短路径距离，用R(i，j)表示；

再通过对全网的遍历式搜索，得到任意两个相邻或非相邻节点间的最短路径距离，构成最短距离矩阵D_s；

步骤(2.2.3)，用多维标度分析定位方法，根据所述最短距离矩阵D_s按下式得到一个节点位置的初始拓扑估计X⁽⁰⁾作为迭代初值：

$X^{\left(0\right)}=\mathrm{diag}({\mathrm{\λ}}_1^{1/2},\·\·\·{\mathrm{\λ}}_{D_s}^{1/2})U^T$

定义e是M维的全1向量，，，，，定义H为中心化矩阵：1-ee^T/M，则可以得到一个对称矩阵B：

B＝-HD_sH＝HX^TXH

上述的

就是对矩阵B进行奇异值分解得到的特征值，U的列是对应的特征向量；

步骤(2.2.4)：用对非邻居节点间的距离估计补足距离矩阵D，得到完全距离矩阵

每个未知距离d_ij的近似概率分布为：

$f^{\left(t\right)}\left(d_{\mathrm{ij}}\right|D,C,X^{(t-1)})=\frac{\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_d{d}_{\mathrm{ij}}}{e}^{-\frac{{[\log d_{\mathrm{ij}}-\log {\mathrm{\δ}}_{\mathrm{ij}}\left(X^{(t-1)}\right)]}^2}{2{\mathrm{\σ}}_d^2}}}{{\∫}_{R(i,k)}^{R(i,j)}\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_d{d}_{\mathrm{ij}}}{e}^{-\frac{{[\log d_{\mathrm{ij}}-\log {\mathrm{\δ}}_{\mathrm{ij}}\left(X^{(t-1)}\right)]}^2}{2{\mathrm{\σ}}_d^2}}}$

C为所述最短路径约束集合，t为迭代次数，D为前述的只包含邻居节点间距离关系的距离矩阵；

所有的未知距离d_ij组成的集合D_mis在第t次迭代后的近似概率分布为：

$f^{\left(t\right)}\left(D_{\mathrm{mis}}\right|D,C,X^{(t-1)})=\underset{\mathrm{mis}}{\mathrm{\Π}}{f}^{\left(t\right)}\left(d_{\mathrm{ij}}\right|D,C,X^{(t-1)})$

步骤(2.2.5)，建立相对定位时需要优化的目标函数：

$Q\left(X\right|X^{(t-1)})=\underset{-\∞}{\overset{\∞}{\∫}}L\left(X\right|D,D_{\mathrm{mis}},C)f^{\left(t\right)}\left(D_{\mathrm{mis}}\right|D,C,X^{(t-1)}){\mathrm{\ΔD}}_{\mathrm{mis}}$

其中Δ代表微分记号，完全距离矩阵 $\hat{D}=D\∪D_{\mathrm{mis}},$ L(X|D，D_mis，C)＝ln f(D，D_mis，C|X)，代表完全距离矩阵在完全意义下的似然函数；

步骤(2.2.6)，在|Q^(t)-Q^(t-1)|＞ε的条件下，ε为设定值，使所述Q(X|X^(t-1))最大化，得到所述节点位置X的优化结果：

$X^{\left(t\right)}=\underset{X}{\max }Q\left(X\right|X^{(t-1)})$

作为最终相对定位结果X_re；

步骤(3)，采用多维标度分析方法把所述相对定位结果X_re转变为绝对定位结果X_ab：设Y为网关节点的真实坐标，Y_re代表网关节点的相对位置坐标，所述网关节点均为已知位置的节点，依次执行以下步骤：

步骤(3.1)，用多维标度分析中的普鲁克方法求取缩放系数s，旋转镜像变化矩阵K，K^TK＝I，I为单位矩阵，以及坐标平移向量t：

对Y_re进行平移和旋转变化使下式中Tr(s，t，K)最小，以使旋转、平移后的Y_re与已知位置Y的均方差最小，t^T为此时的坐标平移向量，1为全l的向量，

Tr(s，t，K)＝tr[Y-(sY_reK+1t^T)]^T[Y-(sY_reK+1t^T)]

步骤(3.2)，按下式计算X_ab

X_ab＝sX_reK+1t^T。

本发明的效果在于克服网络稀疏性对于传感器网络定位精度的影响。本发明中所提出的算法IA-MLE不仅可以根据信道模型自适应地调整优化策略，同时引入网络约束条件，保证了定位信息的充足，使得算法在稀疏网络中仍然可以有令人满意的表现。仿真结果表明，我们的方法相比于传统的未考虑基于最短路径距离的数据推断的方法(此处以经典的最大似然估计定位方法CMLE为例)，在测距信息充足的情况下可以保证不差于后者的表现，在测距信息不足的情况下则明显优于后者。为了观察基于最短距离的数据推断对于定位的实际效果，在数据仿真软件Matlab R2007上对于本发明中的算法IA-MLE进行对比仿真。为了深入分析说明IA-MLE的性能，我们同时观察经典的最大似然估计定位方法CMLE在相同拓扑上的性能，并选取了两种典型拓扑：(a)一个10米×10米的正方形无遮挡网络，如图1(a)所示；(b)一个C型走廊，每条边长10m，宽1m，如图1(b)所示。每个网络中均随机布撒了50个传感器节点。

因为IA-MLE的主要目的是克服网络稀疏性对于传感器定位的影响，所以我们的实验主要研究传感器节点的射频通信距离R，网络的连通度，以及全网平均定位误差三者之间的关系。

在正方形网络中，网络的连通度完全受R的控制，当R增大时，连通度也随之单调上升。取噪声方差 ${\mathrm{\σ}}_d^2=\mathrm{0.1,0.3},$ 已知位置的网关节点数目分别为4个和8个，此时网络连通度以及CMLE和IA-MLE的定位误差随着R的变化如图2所示。当R＞3.5m时，由于网络中的通信链路足够多，此时采用CMLE方法也可以得到较理想的定位结果，利用IA-MLE引入网络约束并不能得到很明显的增益。

但是当R＜3m时，节点间通信链路逐渐减少，网络趋向稀疏。此时，可观测到的距离测量值集合不足以为定位提供充分的信息，所以IA-MLE中引入的网络约束信息的重要性就显现出来的。可以明显地看到，随着R的减少，虽然CMLE和IA-MLE的性能同步恶化，但是后者的速度明显更慢，网络的稀疏程度越高，两者的性能差异越发明显。

在C型网络中，情况又有所不同，如图3所示。此时，网络的连通度并不但受R的控制，拓扑的不规则性也会阻碍位于网络不同边上的节点之间互相通信。可以看到此时IA-MLE的定位误差只有CMLE的约50％，并且这个增益不随着节点通信能力R的增加而下降。

最后，我们比较两种定位方法结果的稳定性，这在传感器网络的定位中，甚至比平均结果的精确程度更加重要。将CMLE和IA-MLE在相同拓扑上分别独立实施30次，所得单次定位误差结果分布如图4所示。可以看到，不管在正方形还是C型拓扑中，IA-MLE定位所产生的结果一致性都更加良好，CMLE产生的结果分布很分散而没有规律。所以，IA-MLE不但比能得到比CMLE更加精确的结果，同时置信区间也更小，一致性更高。

附图说明

图1.(a).50个传感器节点随机布撒在正方形拓扑网络中；

图1.(b).50个传感器节点随机布撒在C形拓扑网络中；

图2.(a).σ² _d＝0.1时正方形网络中CMLE和IA-MLE的性能变化曲线：蓝线星形记号代表CMLE在4个网关节点时的结果；蓝线圆记号代表CMLE在8个网关节点时的结果；红线星形记号代表IA-MLE在4个网关节点时的结果；红线圆记号代表IA-MLE在8个网关节点时的结果；

图2.(b).σ² _d＝0.3时正方形网络中CMLE和IA-MLE的性能变化曲线：蓝线星形记号代表CMLE在4个网关节点时的结果；蓝线圆记号代表CMLE在8个网关节点时的结果；红线星形记号代表IA-MLE在4个网关节点时的结果；红线圆记号代表IA-MLE在8个网关节点时的结果；

图3.(a).σ² _d＝0.1时C型网络中CMLE和IA-MLE的性能变化曲线：蓝线星形记号代表CMLE在4个网关节点时的结果；蓝线圆记号代表CMLE在8个网关节点时的结果；红线星形记号代表IA-MLE在4个网关节点时的结果；红线圆记号代表IA-MLE在8个网关节点时的结果；

图3.(b).σ² _d＝0.3时C型网络中CMLE和IA-MLE的性能变化曲线：蓝线星形记号代表CMLE在4个网关节点时的结果；蓝线圆记号代表CMLE在8个网关节点时的结果；红线星形记号代表IA-MLE在4个网关节点时的结果；红线圆记号代表IA-MLE在8个网关节点时的结果；

图4.(a).正方形拓扑中有4个网关节点时，CMLE的定位误差分布：横轴代表平均定位误差的大小，纵轴代表在每个定位误差上出现的实验重复的次数；

图4.(b).正方形拓扑中有4个网关节点时，IA-MLE的定位误差分布：横轴代表平均定位误差的大小，纵轴代表在每个定位误差上出现的实验重复的次数；

图4.(c).C形拓扑中有3个网关节点时，CMLE的定位误差分布：横轴代表平均定位误差的大小，纵轴代表在每个定位误差上出现的实验重复的次数；

图4.(d).C形拓扑中有3个网关节点时，IA-MLE的定位误差分布：横轴代表平均定位误差的大小，纵轴代表在每个定位误差上出现的实验重复的次数；

图5.用于大范围数据采集的无线传感器网络实例图；

图6.基于最短路径的传感器网络节点自定位方法定位算法步骤框图。

具体实施方式

要对于未测量距离进行补足，除了已有的距离测量值外，我们还需要利用网络中的其他观测信息，在下文中，这些被统称为网络约束。

通过观察距离测量集合D中的空项可以发现，在无线传感器网络中，这些空项的缺失并不是随机的。例如节点i同节点j之间的距离d_ij没有被观测到，这个事实本身提供的信息虽然无法让我们推断出一个具体值，但我们仍然可以断定d_ij超过节点的射频信号传输距离，或者节点i、j之间存在遮挡，没有直接的视线通信路径。上述断定就可以被视为一个由观测得到的网络约束，用来限制节点i和j之间的相对位置关系；在定位中，我们除了使得拓扑反解的结果尽量满足已有的节点间距离观测值之外，也要保证节点i和j间的距离尽可能远，从而使得它们无法直接通信的可能性最大化。

在无线传感器网络中，类似的约束还有很多种，我们在算法中引入的约束越多，则定位精度就越高。在下文中，采用“最短路径约束”为例来详细说明算法。

首先简单介绍“最短路径约束”。

采用传统的弗洛伊德(Floyd)算法计算出整个网络中所有节点对之间的最短路径。对于没有直接通信链路的节点i和节点j，选取它们之间所有路径中权重和最短的一条，作为最短路径，用P_s(i，j)表示；对应的权重和就是最短路径距离R(i，j)。类似，定义P_s(i，j)上节点j的前跳节点k和i的最短路径及其对应距离分别为P_s(i，k)和R(i，k)，则我们假设节点i同节点j之间的距离d_ij的值落在R(i，k)同R(i，j)之间。在常规的传感器网络中，这样的假设是合理且以很大概率成立的。根据这一假设，d_ij的概率分布可以根据式(4)进行更新

$f^{\mathrm{new}}\left(d_{\mathrm{ij}}\right|{\mathrm{\δ}}_{\mathrm{ij}})=\frac{f^{\mathrm{old}}\left[d_{\mathrm{ij}}\right|{\mathrm{\δ}}_{\mathrm{ij}}]}{{\∫}_{R(i,k)}^{R(i,j)}{f}^{\mathrm{old}}\left[d_{\mathrm{ij}}\right|{\mathrm{\δ}}_{\mathrm{ij}}]}---\left(4\right)$

通过上述更新，采用两个最短路径值作为上下界，我们缩小了d_ij的可能分布范围，从而引入了“最短路径约束”信息。

本发明中的算法：基于最短路径的传感器网络节点自定位方法结合了最大似然估计和多维标度分析方法进行拓扑反解，通过迭代得到最优解。算法分为两阶段，首先是相对定位阶段，使用最大似然估计得到网络的相对拓扑；然后在绝对定位阶段中再用多维标度分析方法将相对拓扑转化为绝对拓扑。

相对定位阶段包含以下4个步骤。

步骤一：使用传统的定位方法，例如多维标度分析方法(MDS-MAP)，根据原始的不完全距离矩阵D得到一个初始拓扑估计X⁽⁰⁾作为迭代初值。

进入迭代。

步骤二：这一步的工作是补足距离矩阵，即利用网络约束对未观测到的距离进行推断。同传统意义上的推断不同，本发明中我们并没有给出一个具体的推断值d_ij，而是利用网络约束来推断d_ij的概率分布更新。例如，在第t次迭代中，我们可以得到上一次迭代所得到的拓扑估计X^(t-1)。由于网络中节点的真实位置X不可知，所以，我们只能用X^(t-1)作为X的近似估计，求得此时每个未知距离的近似概率分布。

$f^{\left(t\right)}\left(d_{\mathrm{ij}}\right|D,C,X^{(t-1)})=\frac{r^{(t-1)}\left(d_{\mathrm{ij}}\right)}{{\∫}_{d(i,k)}^{d(i,j)}{r}^{(t-1)}\left(d_{\mathrm{ij}}\right)}---\left(5\right)$

$r^{(t-1)}\left(d_{\mathrm{ij}}\right)=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_d{d}_{\mathrm{ij}}}{e}^{-\frac{{[\log d_{\mathrm{ij}}-\log {\mathrm{\δ}}_{\mathrm{ij}}\left(X^{(t-1)}\right)]}^2}{2{\mathrm{\σ}}_d^2}}$

其中C代表的是定位中所使用的网络约束集合，在此处，即为“最短路径约束”集合。

如果将所有未知距离组成的集合表示为D_mis，那么D_mis在该次迭代后的近似概率分布可以表示为

$f^{\left(t\right)}\left(D_{\mathrm{mis}}\right|D,C,X^{(t-1)})=\underset{\mathrm{mis}}{\mathrm{\Π}}{f}^{\left(t\right)}\left(d_{\mathrm{ij}}\right|D,C,X^{(t-1)})---\left(6\right)$

步骤三：这一步建立定位问题中需要优化的目标似然函数。由于未知距离只已知概率分布f^(t)(D_mis|D，C，X^(t-1))，所以“完全”的似然函数其实是一个期望似然函数，通过对D_mis部分的概率密度求期望以去掉D_mis的随机性。即需要优化的最终目标函数为

Q(X|X^(t-1))＝∫L(X|D，D_mis，C)f^(t)(D_mis|D，C，X^(t-1))dD_mis    (7)

其中

L(X|D，D_mis，C)＝ln f(D，D_mis，C|X)        (8)

步骤四：最后一步中，我们优化X，使得(7)中得到的目标函数Q(X|X^(t-1))最大化，这一步中所得到的优化结果

$X^{\left(t\right)}=\underset{X}{\max }Q\left(X\right|X^{(t-1)})$

即是该次迭代中所得到的全网拓扑估计的更新值。

如果|Q^(t)-Q^(t-1)|＞ε，则t自加1，回到步骤二开始下一次迭代；否则，迭代终止，X^(t)就是用基于最短路径的传感器网络节点自定位方法得到的最终相对定位结果X_re。

在基于最短路径的传感器网络节点自定位方法的每次迭代中，由于不知道真实的网络拓扑X，D_mis的概率密度以及对应的似然函数都是采用了X^(t-1)作为近似。但是通过理论推导可以证明，上述迭代过程最终可以收敛到一个使得似然函数L(X|D，C)最大化的点，即

$X^{\left(t\right)}\stackrel{t\→\∞}{\→}X=\underset{X}{\max }L\left(X\right|D,C)---\left(9\right)$

如果将(D，C)看成在网络中所有可以被我们观测到并用于定位的信息，那么显然基于最短路径的传感器网络节点自定位方法可以最充分地利用这些已有信息以得到尽可能高的定位精度。而对于这些已有信息的利用，正是通过对未观测到的节点间距离进行数据推断而实现的。

相对定位结果X_re只能反映节点间位置的相互关系，对X_re做任意正交变换，包括旋转、平移和镜像变换，所得结果仍然是满足式(9)的解。所以我们采用多维标度分析技术将X_re转换为绝对坐标X_ab。假设网络中有若干已知位置的网关节点，多维标度分析技术将X_re进行旋转、平移和镜像变换，使得X_ab中网关节点的绝对坐标同其已知的真实值的均方误差最小。

设Y_re代表锚节点的相对位置坐标，Y代表锚节点的真实坐标，寻找变换Y_ab＝sY_reK+1t^T，其中K为旋转镜像变化矩阵，要求K^TK＝I，s为缩放系数，1为全1的向量，t为坐标平移向量，使得Y_ab与已知位置Y的均方误差最小，即要求最小化Tr(s，t，K)

Tr(s，t，K)＝tr[Y-(sY_reK+1t^T)]^T[Y-(sY_reK+1t^T)]        (10)

这个问题可以由多维标度分析中的普鲁克分析技术求解，将接得的(s，t，K)用于X_re，我们就得到了基于最短路径的传感器网络节点自定位方法的最终绝对定位结果

X_ab＝sX_reK+1t^T         (11)

无线传感器网络的应用场景各异，对定位的需求也各不相同。因此，在进行定位算法的设计前，必须选定应用场景进行有针对性的设计。

如图5所示，本文选用传感器网络中广泛应用的大范围数据采集场景(例如土壤温湿度监测、森林火险预警、智能大厦人员数据采集等)作为研究前提。这类应用主要是利用大量传感器节点，对某一较大区域中的数据进行周期性的采集，并将数据通过多跳路由传回到监测中心进行分析处理。这类网络的布撒形式多样，但网络一旦形成，网络拓扑将随之固定，不再改变。网络路由一般采用树状路由，数据收集节点是路由树的根结点。

由于传感器数量众多，基于成本、体积和功耗方面的考虑，在本发明的算法实施中，所有传感器节点不装配GPS、超声收发器、有向天线等额外的定位和测距设备，节点射频模块只具备射频信号强度检测能力(RSS)。

在微处理器模块上，我们使用美国德州仪器(TI)公司的16位超低功耗单片机MSP430F149。该款处理器的特点是功耗低、性能强，同时片上资源丰富，在1.8-3.6V下都可以工作，具有超低功耗模式、待机模式、关闭模式和活动模式，并且所有低功耗模式唤醒到运行模式的唤醒时间小于6us，非常适合传感器网络的超低功耗和对功耗控制的需求。

在射频端，由于300M-3GHz UHF频段的无线电波具有很强的直射、绕射、反射以及抗干扰能力强等特性，所以我们在该频段中选取无线传感器网络工作的信道。最终，我们选择了挪威Nordic Semiconductor新推出的单片射频收发芯片nRF905。与其他类似芯片相比，nRF905功耗较低，工作电压低至1.9～3.6V，以-10dBm的输出功率发射时电流只有11mA，工作于接收模式时的电流为12.5mA，支持空闲模式与关机模式，易于实现节能，并且具有发送和接收状态的多级功率控制，可以方便传感器网络射频功率控制相关技术的研究。nRF905采用优化的GMSK调制技术，工作于433/868/915MHz三个ISM频道，共有79个频道，频道之间的转换时间小于650us，这一特性对研究传感器网络的多信道协议十分有利。

在上述的应用场景和硬件平台上，基于最短路径的传感器网络节点自定位方法可以按照图6中所示流程由以下6个步骤实现。

下面将详细介绍每个步骤的实现过程。

1.节点布撒

将传感器节点布撒在需要进行环境监控和数据采集的区域中；同时放置网关节点，与电脑终端或者外界网络相连。

节点可以布撒在网络中事先指定的位置，也可以随机布撒。图5是一个典型的布撒后传感器网络。

2.路由建立和数据传输

在大范围数据采集传感器网络中，由于传感器节点分布范围广，所以在本方法中，我们采用树状路由令节点将数据传回网关节点。网络中的所有节点到任意一个网关节点的路由就形成一颗以该网关节点为根节点的树。图5就给出了网络中大量传感器节点到上方网关节点的一个树状路由图，图中线段表示实际的路由链路。

从上图可以发现，从所有传感器节点到一个网关节点的树状路由图中，已经可以获得网络中大量节点之间连接的信息。同时由于传感器节点均具有射频信号检测能力，对于存在通信链路的节点对，就可以获得这些链路的距离测量值。在本发明中，我们采用定向扩散路由协议(Directed Diffusion)来建立网络中的路由关系。

在路由的建立阶段，网关节点向所有传感器节点发送对于任务描述的“请求”(interest)，“兴趣”会逐渐在全网中扩散，最终达到所有传感器节点，与此同时也建立起了从网关节点到传感器节点的“梯度”。每个传感器节点都有自己对网关节点的最大“梯度”方向，即下一跳传输的目的节点编号(ID)。只需每个节点都将此下一跳节点ID打入传感器数据包，按照选择的路由发往网关节点，即可以在网关节点重建树状路由图。同时，在每一跳数据的传输中，射频信号强度测距技术就可以得到该跳收发节点间的的距离测量值。

3.网络连接信息的提取

这一步的目的是提取多个网关收到数据中关于节点测距的信息。在上一步骤中，每个节点都将自己的下一跳节点ID发往网关节点，在每个网关节点都可以得到全网节点到其自身的完整树状路由。然后所有的网关节点都将这些信息发往一个网关节点或者专用的定位服务器，将信息进行综合。此时，我们已经掌握了相当多可以用于定位的节点间连接或测距信息。

4.相对坐标定位

利用上一步中的节点间测距信息，进行全网相对定位。定位方法首先利用最短路径算法补足距离矩阵D，然后采用期望最大化方法进行循环优化。详细说明四个步骤。

步骤一：每个节点均以固定发射功率向全网广播数据包，如果任一节点i收到来自节点j的数据包，则表明节点i和j之间存在一跳数据通路，两者互为邻居节点，可以利用RSS功率检测法测得互相之间的距离。对于非邻居节点i和j，我们采用弗洛伊德最短路径算法求解它们之间的距离估计，通过遍历式的搜索找到i和j之间的所有多跳数据通路，并且选取其中单跳路径总和最短的一条作为i和j之间的最短路径，对应的距离即是非邻居节点间的最短路径距离。通过对全网的遍历式搜索，我们可以得到任意两个节点(相邻或非相邻)的最短路径距离，即最短路径距离矩阵D_s。

将D_s作为节点间真实欧式距离的估计，采用经典的多维标度分析MDS-MAP算法来进行相对定位，得到拓扑估计的初值X⁽⁰⁾。任意两个节点i和j(它们的坐标分别用x_i和x_j表示)之间的欧式距离d_ij可以被表示为

$d_{\mathrm{ij}}=\sqrt{{(x_i-x_j)}^T(x_i-x_j)}---\left(12\right)$

定义

则平方距离矩阵 $D_s^2=[{d_{\mathrm{ij}}^2]}_{i=1,j=1}^{i=M,j=M}$ 可以被表示为：

其中e是n维的全1向量。定义H为中心化矩阵：1-ee^T/N，则可以推出：

B＝-HD_sH＝HX^TXH         (14)

因此，由平方距离矩阵可以得到矩阵B，然后就可以对矩阵B进行分解得到坐标矩阵X_o因为B是实对称矩阵，则可以将B写成：

B＝Udiag(λ₁，…λ_p)U^T        (15)

即为矩阵B的奇异值(SVD)分解，其中p为坐标点的维数，U的列是B的特征向量，则：

$X=\mathrm{diag}({\mathrm{\λ}}_1^{1/2},\·\·\·{\mathrm{\λ}}_D^{1/2})U^T---\left(16\right)$

即为MDS-MAP算法所得到的初始定位结果X⁽⁰⁾。

步骤二：这一步的工作是补足距离矩阵，即对非邻居节点间的距离进行估计。步骤一采用最短路径距离作为非邻居节点间的距离估计，在通常的网络中，是非常不精确的。所以我们在这一步中采用网络约束，对距离进行更精确的推断。采用所有的未知距离d_ij作为变量集合，可以得到它们的联合概率密度如式(6)所示。

步骤三：建立目标似然函数如式(7)所示。所得结果是一个以坐标位置X为变量的表达式。此处只需要记录关键变量σ_d以及已知值X^(t-1)，D和C。

步骤四：采用牛顿法，优化X_o在每个优化点寻找Q函数上升最快的方向，最终达到式(7)的最大值。

5.绝对坐标变换

采用多维标度分析技术将相对定位结果转变为绝对定位结果。

用数组分别保存相对定位结果X_re，锚节点的相对定位结果Y_re，锚节点的真实坐标Y。首先根据式(10)求得旋转镜像变化矩阵K和缩放系数s，最后利用式(11)求得最终的绝对定位结果，保存在数组X_ab中。

6.结果输出

定位服务器将全网定位结果发回各个网关节点，并由网关节点沿着各条路由的逆方向将每个节点自身位置传回对应的节点。

Claims (1)

1.稀疏测量集上基于最短路径的传感器网络节点自定位方法，其特征在于，所述方法是在计算机中依次按以下步骤实现的：

步骤(1)，在需要进行环境监控和数据采集的区域中构建无线传感器网络：

在所述区域中布撒M个无线传感器和多个网关，所述无线传感器依次由物理传感器件、微处理器和射频模块依次串接而成，所述各无线传感器的传输功率相同，传输距离是2米，且只能视线传输，所述射频模块工作在300M～3GHz频段；

所述无线传感器网络用一个无向图G_r，M(X，E，Z)表示，r为该网络所处的空间维数，r＝2或者3，顶点集合X表示网络中分布的传感器节点位置，E为无向边集合，表示能互相通信的节点对之间的通信链路，其通信能力同测距能力等价，Z是无向边集合E中每条边的距离的集合；从而，以所述网关节点为根节点，以所述传感器节点为中间节点或叶子节点，形成多条树状路由，采用定向扩散协议建立所述无线传感器网络中的路由关系，把所述各节点的数据传回各网关节点；在所述树状路由建立阶段，在所述网关节点请求下，各传感器节点把自己下一跳传输的目的节点ID发往各网关节点，以在各网关节点内建立树状路由图，同时传输的还有该跳收发节点的距离测量值d_ij，i和j分别是两个相邻节点各自的编号，从而拼成一个距离矩阵D，是一个M阶对称矩阵；

步骤(2)，利用所述距离矩阵D依次按以下步骤对所述无线传感器网络进行相对定位：

步骤(2.1)，所述的每个节点均以固定发射功率向所述无线传感器网络广播数据包，若任一节点j收到来自节点i的数据包，表明两者之间存在一跳数据通路，两者互为邻居节点，用w_ij＝1表示，否则，为非邻居节点，用w_ij＝0表示，构成一个连接矩阵W，是一个M阶对称矩阵；

步骤(2.2)，依次按以下步骤进行相对定位：

步骤(2.2.1)，对于所述的邻居节点用射频信号强度检测测得互相之间的距离：所得的距离d_ij满足下式所示的概率分布：

$f\left(d_{\mathrm{ij}}\right|{\mathrm{\δ}}_{\mathrm{ij}})=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_d{d}_{\mathrm{ij}}}{e}^{-\frac{{(\log d_{\mathrm{ij}}-\log {\mathrm{\δ}}_{\mathrm{ij}})}^2}{2{\mathrm{\σ}}_d^2}}$

其中d_ij为节点i，j之间的距离测量值，

δ_ij是i，j间真实欧式距离的期望值，

σ_d为方差，代表测量中的不确定性，由方差σ_P根据下式所转换得到，

${\mathrm{\σ}}_d=\frac{{\mathrm{\σ}}_P}{10n}$

σ_P是节点测量接收射频信号的功率时所叠加的白噪声的方差，代表了射频功率测量的不确定性，n代表信道中的衰落指数，用于反映信号随着传输距离衰落的速度，在室内环境中，它的典型值为1.6-1.8；

上述公式表示节点间的距离测量值由于受到测量功率噪声的影响，概率密度满足以真实欧式距离为中心的对数高斯分布，

节点位置的最大似然解X在获得所有节点对距离测量值后，就是如下似然函数的最大值：

$L\left(X\right|D)=\underset{i,j}{\mathrm{\Σ}}\ln f\left(d_{\mathrm{ij}}\right|{\mathrm{\δ}}_{\mathrm{ij}})$

步骤(2.2.2)，当存在所述的非邻居节点时，用弗洛伊德最短路径算法求解它们之间的距离估计：通过遍历式搜索找出非邻居节点的所有多跳数据通路，并且选取其中单调路径距离总和最短的一条作为所述非邻居节点之间的最短路径，对应的距离就是该非邻居节点间的最短路径距离，用R(i，j)表示；

再通过对全网的遍历式搜索，得到任意两个相邻或非相邻节点间的最短路径距离，构成最短距离矩阵D_s；

步骤(2.2.3)，用多维标度分析定位方法，根据所述最短距离矩阵D_s按下式得到一个节点位置的初始拓扑估计X⁽⁰⁾作为迭代初值；

$X^{\left(0\right)}=\mathrm{diag}({\mathrm{\λ}}_1^{1/2},\·\·\·{\mathrm{\λ}}_{D_s}^{1/2})U^T$

定义e是M维的全1向量，，，，，定义H为中心化矩阵：1-ee^T/M，则可以得到一个对称矩阵B：

B＝-HD_sH＝HX^TXH

上述的λ₁，…

就是对矩阵B进行奇异值分解得到的特征值，U的列是对应的特征向量；

步骤(2.2.4)：用对非邻居节点间的距离估计补足距离矩阵D，得到完全距离矩阵

每个未知距离d_ij的近似概率分布为：

$f^{\left(t\right)}\left(d_{\mathrm{ij}}\right|D,C,X^{(t-1)})=\frac{\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_d{d}_{\mathrm{ij}}}{e}^{-\frac{{[\log d_{\mathrm{ij}}-\log {\mathrm{\δ}}_{\mathrm{ij}}\left(X^{(t-1)}\right)]}^2}{2{\mathrm{\σ}}_d^2}}}{{\∫}_{R(i,k)}^{R(i,j)}\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_d{d}_{\mathrm{ij}}}{e}^{-\frac{{[\log d_{\mathrm{ij}}-\log {\mathrm{\δ}}_{\mathrm{ij}}\left(X^{(t-1)}\right)]}^2}{2{\mathrm{\σ}}_d^2}}}$

C为所述最短路径约束集合，t为迭代次数，D为前述的只包含邻居节点间距离关系的距离矩阵；

所有的未知距离d_ij组成的集合D_mis在第t次迭代后的近似概率分布为：

$f^{\left(t\right)}\left(D_{\mathrm{mis}}\right|D,C,X^{(t-1)})=\underset{\mathrm{mis}}{\mathrm{\Π}}{f}^{\left(t\right)}\left(d_{\mathrm{ij}}\right|D,C,X^{t-1})$

步骤(2.2.5)，建立相对定位时需要优化的目标函数：

$Q\left(X\right|X^{(t-1)})=\underset{-\∞}{\overset{\∞}{\∫}}{L\left(X\right|D,D_{\mathrm{mis}},C)f}^{\left(t\right)}\left(D_{\mathrm{mis}}\right|D,C,X^{(t-1)})\mathrm{\Δ}{D}_{\mathrm{mis}}$

其中Δ代表微分记号，完全距离矩阵 $\hat{D}=D\∪D_{\mathrm{mis}},$ L(X|D，D_mis，C)＝ln f(D，D_mis，C|X)，代表完全距离矩阵

在完全意义下的似然函数；

步骤(2.2.6)，在|Q^(t)-Q^(t-1)|＞ε的条件下，ε为设定值，使所述Q(X|X^(t-1))最大化，得到所述节点位置X的优化结果：

$X^{(t-1)}=\underset{X}{\max }Q\left(X\right|X^{(t-1)})$

作为最终相对定位结果X_re；

步骤(3)，采用多维标度分析方法把所述相对定位结果X_re转变为绝对定位结果X_ab：设Y为网关节点的真实坐标，Y_re代表网关节点的相对位置坐标，所述网关节点均为已知位置的节点，依次执行以下步骤：

步骤(3.1)，用多维标度分析中的普鲁克方法求取缩放系数s，旋转镜像变化矩阵K，K^TK＝I，I为单位矩阵，以及坐标平移向量t：

对Y_re进行平移和旋转变化使下式中Tr(s，t，K)最小，以使旋转、平移后的Y_re与已知位置Y的均方差最小，t^T为此时的坐标平移向量，1为全1的向量，

Tr(s，t，K)＝tr[Y-(sY_reK+1t^T)]^T[Y-(sY_reK+1t^T)]

步骤(3.2)，按下式计算X_ab

X_ab＝sX_reK+1t^T。

Images