CN101354435A - 基于距离大小顺序关系的传感器网络节点自定位方法 - Google Patents

Abstract

基于距离大小顺序关系的传感器网络节点自定位方法属于无线传感器网络自定位技术领域，其特征在于，依次含有节点布撒、路由建立和数据传输、网络连接信息提取、节点相对坐标定位、绝对坐标变换、计算机输出结果等步骤，在节点相对坐标定位时，利用最短路径距离中的距离大小顺序关系得到节点位置的初始估计拓扑，以此为基础建立优化节点距离矩阵和作为初始估计拓扑的反解欧式距离的对数似然函数，再用相邻模块比较方法交替迭代优化距离矩阵和所述欧式距离；本发明不依赖于距离的具体测量值，只利用节点间距离的相对顺序关系，具有自己独特的优势。

Description

基于距离大小顺序关系的传感器网络节点自定位方法

技术领域

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

背景技术

微型电子集成化(MEMS)技术的迅速发展奠定了设计和实现片上系统(SOC)的基础，使得将多种传感器集成为一体，制造小型化、低成本、多功能的传感器节点成为可能。由大量的MEMS传感器节点组成的无线传感器网络，已被广泛应用于军事作战、环境监测、医疗诊断、家庭娱乐、空间探索以及商业制造等多个领域。区别于传统的无线自组织网络，无线传感器网络以监控物理世界为主要目标，是一种测控网络，并且具有超大规模、无人值守、易受物理环境影响(动态性强)等特征。

在无线传感器网络的研究中，节点的定位问题是一个热点的研究领域，这是因为节点的准确定位是无线传感器网络应用的重要条件。例如在各种监控网络中，都需要知道传感器节点的位置信息，从而获知信息来源的准确位置。并且，利用节点的位置信息还可用于提高路由效率，向部署者报告网络的覆盖质量，实现网络的负载均衡和网络拓扑的自配置等。

目前已有的定位技术主要基于节点间物理距离的测量列出多元方程组从而进行拓扑反解，例如到达时间指示(TOA)技术利用节点间的超声传播延时，接收信号强度指示(RSSI)利用节点间的射频信号传播强度。这些方法应用的信号测量手段虽然各有差异，但是定位的基本思想都是相同的，就是对一个特定形式的代价函数的最小化过程。这些代价函数可以有一个通用的表示，如式(1)

$\mathrm{Stress}(X,\hat{D})=\underset{i}{\mathrm{\Σ}}\underset{j}{\mathrm{\Σ}}g({\mathrm{\δ}}_{\mathrm{ij}},{\hat{d}}_{\mathrm{ij}})---\left(1\right)$

其中，X是网络中传感器节点的位置坐标，也就是定位问题中需要求解的变量集合；{δ_ij}即是由X计算得到的节点间欧式距离； $\hat{D}=\left\{{\hat{d}}_{\mathrm{ij}}\right\}$ 是由距离测量值D＝{d_ij}构建的节点间距离矩阵，

和D通常有很紧密的联系，但是在很多情况下，两者并不需要严格相等。

基于式(1)，所有的定位问题都可以被抽象为以下三个步骤。

1)建立节点间距离矩阵

2)定义代价函数g和Stress

3)优化节点位置估计X，使得{δ_ij}同

的距离在Stress函数的意义下尽可能地接近，使得Stress函数最小化。

在传统的定位技术中，

和D中对应元素的值通常是通过一个线性关系直接映射的。例如在RSSI或者TOA技术中，距离测量值是连续值，所以可以直接应用，即 ${\hat{d}}_{\mathrm{ij}}=d_{\mathrm{ij}};$ 在基于跳数估计的方法如DV-Hop中，距离测量值d_ij是离散跳数，所以需要有一个线性关系将其转化为距离，即 ${\hat{d}}_{\mathrm{ij}}=b*d_{\mathrm{ij}},$ 其中b对应着网络中的平均每跳距离估计。

这种依赖于距离测量具体数值的机制存在着几个弱点：1.定位算法的精度对距离测量值的准确程度依赖性很高；2.现实网络中的恶劣电磁传播环境会导致节点间通信困难，使得某些节点对的距离无法测量。本发明提出了一种全新的不依赖于距离测量值的定位算法，节点的定位只需要利用节点对距离互相之间的顺序关系。这种方法对于距离测量值不足的问题，以及测量噪声都具有更好的鲁棒性。

发明内容

我们在本发明中提出了基于顺序关系的最大似然估计定位方法RMLE(Ranking-basedMLE)。本发明的目的在于设计一种在任何测量模型下，不依赖于射频距离测量值，得到可靠定位结果的网络自定位技术。

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

步骤(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间真实欧式距离的期望值，

σ_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)，把步骤(2.2.2)中的最短距离矩阵D_s中的所有元素按值得大小进行排序，得到所述元素的距离大小的距离顺序关系，并用一个长度为(M-1)(M-2)/2的三维数组N存放所述顺序关系；其中，第一维记录对应元素的值，第二维及第三维分别记录该元素在D_s中的行和列编号，

步骤(2.3)，按下式从所建最短距离矩阵D_s得到一个节点位置的初始拓扑估计X⁽⁰⁾作为迭代的初信，

$X^{\left(0\right)}=\underset{X}{\min }\mathrm{Stress}(D_s,X)=\sqrt{\underset{i,j}{\mathrm{\Σ}}{(d_{\mathrm{ij}}-d_{\mathrm{ij}}^{\mathrm{euc}}\left(X\right))}^2}$

d_ij ^euc(X)为反解的欧式距离

$d_{\mathrm{ij}}^{\mathrm{euc}}\left(X\right)=\sqrt{\underset{k}{\mathrm{\Σ}}{(x_{\mathrm{ik}}-x_{\mathrm{jk}})}^2}$

步骤(2.4)，建立作为优化标准的节点距离矩阵和节点位置矩阵X的对数似然函数：

$\mathrm{Stress}(X,\hat{D})=\underset{i}{\mathrm{\Σ}}\underset{j}{\mathrm{\Σ}}(\ln \frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_d{\hat{d}}_{\mathrm{ij}}}-\frac{{(\log {\hat{d}}_{\mathrm{ij}}-\log d_{\mathrm{ij}}^{\mathrm{euc}})}^2}{2{\mathrm{\σ}}_d^2})$

再使

最小，其步骤依次如下：

步骤(2.4.1)，对于由X⁽⁰⁾构建的欧式距离矩阵{d_ij ^euc(X⁽⁰⁾)}，找到最优的节点间距离矩阵

使代价函数S⁽⁰⁾最小：

$S^{\left(0\right)}=\mathrm{Stress}({\hat{D}}^{\left(0\right)},D^{\mathrm{euc}}\left(X^{\left(0\right)}\right)),$ 并开始迭代循环，其步骤如下：

步骤(2.4.1.1)，使用牛顿优化法，在第t次循环中，优化X^(t)，使该X^(t)与前一次迭代的距离矩阵

之间的代价函数最小，即

$X^{\left(t\right)}=\underset{X}{\min }\mathrm{Stress}({\hat{D}}^{(t-1)},D^{\mathrm{euc}}\left(X\right))$

步骤(2.4.1.2)，利用相邻模块比较方法来优化

使该

与D^euc(X^(t))之间的代价函数最小，即：

${\hat{D}}^{\left(t\right)}=\underset{\hat{D}}{\min }\mathrm{Stress}(\hat{D},D^{\mathrm{euc}}\left(X^{\left(t\right)}\right))$

步骤(2.4.1.3)，在第t次迭代中，得到

此时的代价函数为S^(t)，

$S^{\left(t\right)}=\mathrm{Stress}({\hat{D}}^{\left(t\right)},D^{\mathrm{euc}}\left(X^{\left(t\right)}\right))$

若：|S^(t)-S^(t-1)|＞ε，ε为设定值，则进入下一次迭代，优化X^(t)，

若：|S^(t)-S^(t-1)|＜ε，则收敛，所得到的

和X^(t)就是相对定位阶段的最终结果

和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；

根据权利要求，1所建的基于距离大小的顺序关系的传感器网络节点自定位方法，其特征在于，所建的相邻模块比较方法依次含有以下步骤：

步骤(a)，把D^euc(X^(t))赋给

步骤(b)，把最短路径矩阵D_s中的所有项按照数值大小的顺序形式排列成一个一维数组A，即是前面所述三维数组N中的第一维，其长度为(M-1)(M-2)/2，M为所建传感器网络中的传感器数，按照相同的下标顺序，也把优化的距离矩阵

中的元素排列成一个一维数组Q；

步骤(c)，对

进行搜索更新：

在任何一个位置q，若下式成立，则不对

和d_ij ^(t)做调整，继续向前搜索：

A(q-1)＝d_kl小于等于A(q)＝d_ij

$Q(q-1)={\hat{d}}_{\mathrm{kl}}^{\left(t\right)}$ 小于等于 $Q\left(q\right)=d_{\mathrm{ij}}^{\left(t\right)}$

否则，便在保证位置q-1之前的顺序关系下减小

增加d_ij ^(t)，在必要的时候调整

之前的元素；

步骤(d)，根据步骤(b)和步骤(c)得到的结果，得到下述最佳调整公式：

${\hat{d}}_{\mathrm{ij}}^{\left(t\right)}=Q^{\mathrm{new}}\left(q\right)=...=Q^{\mathrm{new}}(q-k)={(Q^{\mathrm{old}}\left(q\right)...Q^{\mathrm{old}}(q-k))}^{\frac{1}{k+1}}$

其中，Q^new代表调整后的数值，Q^old代表调整前的数值；

所述的k值满足如下条件：Q^new(q-k)≥Q^old(q-k-1)；

步骤(e)，重复步骤(c)到步骤(d)。

本发明的效果在于仅利用节点间距离大小顺序关系就可以得到较好的定位结果。RMLE定位方法虽然精度不如直接利用测量距离，但是减少了网络内的距离测量成本，保证在稀疏和高噪声的网络中，仍然可以得到较为准确的定位结果。仿真结果表明，我们的方法相比于传统的未考虑基于最短路径距离的数据推断的方法(此处以经典的最大似然估计定位方法CMLE为例)，在测距信息充足的情况下可以保证不差于后者的表现，在测距信息不足的情况下则明显优于后者。为了观察基于最短距离的数据推断对于定位的实际效果，在数据仿真软件MatlabR2007上对于本发明中的算法RMLE进行对比仿真。为了深入分析说明RMLE的性能，我们同时观察经典的最大似然估计定位方法CMLE在相同拓扑上的性能，并选取了两种典型拓扑：(a)一个10米×10米的正方形无遮挡网络，如图1所示；(b)一个C型走廊，每条边长10m，宽1m，如图2所示。每个网络中均随机布撒了50个传感器节点。

取信道噪声的方差σ² _d＝0.2，首先观察网络连通性对于定位性能的影响。图3反映了在正方形拓扑中的情况。RMLE的性能优势主要表现在稀疏网络中，当节点的射频距离R＝2m时，CMLE的平均定位误差是RMLE的两倍。

而在C型拓扑中，RMLE的性能优势更加明显，如图4所示。这就说明在这种拓扑不规则性造成的稀疏网络当中，有向顺序关系以外的其他顺序约束对于定位也有很重要的意义，特别是那些不同边上的节点之间的距离信息。

接下来观察一下距离测量噪声对于定位性能的影响。如图5和图6所示，由于RMLE没有直接使用距离信息，所以它们比CMLE具有对测量噪声更强的鲁棒性。

在正方形拓扑下，设置网络平均连通度为8.9，如图5所示，随着噪声方差σ² _d从0.1增加到0.5，CMLE的平均定位误差上升了122％，而RMLE只上升了48％。在C型拓扑下，情况也类似，CMLE的定位误差随着测量噪声上升的速度要比RMLE快得多，如图6所示。

附图表说明

图1.50个传感器节点随机布撒在正方形拓扑网络中；

图2.50个传感器节点随机布撒在C形拓扑网络中；

图3.CMLE，RMLE在正方形无遮挡网络中的性能比较，σ² _d＝0.2：红线为CMLE的结果，蓝线为RMLE的结果；

图4.CMLE，RMLE在C型网络中的性能比较，σ² _d＝0.2：红线为CMLE的结果，蓝线为RMLE的结果；

图5.CMLE，RMLE在正方形无遮挡网络中的性能比较，连通度＝8.9：红线为CMLE的结果，蓝线为RMLE的结果；

图6.CMLE，RMLE在正方形无遮挡网络中的性能比较，连通度＝5.7：红线为CMLE的结果，蓝线为RMLE的结果；

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

图8.基于距离大小顺序关系的传感器网络自定位算法RMLE步骤框图

表1.UDP算法典型工作工程示例

具体实施方式

在RMLE中，由于使用的顺序关系信息只是直接测量值的函数，所以可以在一定程度上减小σ² _d的影响。举个例子，对于任意两个节点间距离δ_ij和δ_kl，只要噪声导致的测量误差小于/δ_ij-δ_kl/，它们之间的顺序关系就可以同实际保持一致。/δ_ij-δ_kl/这个误差容限就是保证RMLE对于测量噪声鲁棒性的真正原因。

下面进行详细说明。

RMLE定位分为两个阶段：相对定位阶段和绝对定位阶段。相对定位利用最大似然估计得到统计意义上最优的网络相对拓扑，绝对定位阶段再利用多维标度分析方法(MDS)将相对拓扑转化为绝对拓扑。

相对定位阶段由以下8个步骤所组成：

步骤一：全网中的节点之间进行互相通信，测量可以直接通信的节点对距离。

在本发明中，我们采用射频信号强度检测(RSS)的方法来测量节点间距离。这种测量方法的核心思想是在一个特定的无线信道中，射频信号的功率衰减同信号的传输距离有一个一一映射的函数关系。如果发射节点按照恒定的功率发射射频信号，接收节点就可以根据接收信号的强度通过函数关系求得信号传输的物理距离。在实际网络中，对于信道建模的最准确方法是根据实际测量，利用经验数据来确定功率衰减同传输距离的关系。由于本发明中的算法是完全独立于信道模型的，即可以工作在任意信道模型上，所以为了叙述方便，在下文中我们采用Rapport在[T.S.Rappaport.2001.“Wireless Communications：Principles andPractice 2nd ed.”，Pearson Education.]中提出的经典对数信道模型进行描述。在这种信道模型的假设下，信号功率衰减量的dB数同传输距离的对数成正比，即

$P\left[\mathrm{dB}\right]=\stackrel{\‾}{P}\left(d\right)\left[\mathrm{dB}\right]+X_{{\mathrm{\σ}}_P}\left[\mathrm{dB}\right]$

$\stackrel{\‾}{P}\left(d\right)\left[\mathrm{dB}\right]{=P}_0\left[\mathrm{dB}\right]-10\mathrm{nlg}\left(\frac{d}{d_0}\right)---\left(2\right)$

其中d代表信号传播距离，d₀是一个固定值的参考距离，P₀是在固定发射功率下，d₀处接收功率的典型值。P代表在距离d处接收功率的典型值，可以看到它的衰减量同d的对数成正比。P代表在d接收功率的实际值。在实际信道中，由于遮蔽和多径等效应的影响，P不总是和典型值P相等，而表现为在P附近有一些波动，所以在Rapport的模型中，该种波动被表示成一个正态分布的随机变量

其均值为0，方差为σ_P。

在这种信道模型假设下，通过简单推导可以得出，距离的测量值满足以真实欧式距离为中心的对数正态分布，即

$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}}-{\mathrm{long\δ}}_{\mathrm{ij}})}^2}{{2\mathrm{\σ}}_d^2}}---\left(3\right)$

其中δ_ij为发射节点和接收节点之间的真实欧式距离，d_ij为距离测量值。而σ_d由σ_P折算得到，代表了接收信号功率值的随机性。

步骤二：采用最短路径算法，得到全网中所有节点对的最短路径和最短路径距离，建立最短路径距离矩阵D_s。在本发明中，我们采用弗洛伊德(Floyd)算法来计算最短路径。该算法的基本思想是通过穷尽搜索，找到任意两个节点之间的最佳路径，使得这条路径上的权重(对应为这两个节点之间的距离)最小。

有了网络中任意两个节点之间的最短距离，我们可以将这些距离进行排序，得到它们的顺序关系。

此处我们以最短路径距离为基准来确定节点间距离的顺序关系，这是基于这样的假设，即最短路径距离同节点间的实际距离有同样的单调顺序关系，两个节点间的最短路径距离越长，它们的真实距离也就越远。在实际的传感器网络中，虽然不是完全准确的，这样的假设在很大程度上是成立的。

步骤三：使用传统的定位方法，例如MDS-MAP，根据D_s得到一个初始拓扑估计X⁽⁰⁾作为迭代初值。该种方法优化节点坐标估计X，使得反解的欧式距离

$d_{\mathrm{ij}}^{\mathrm{euc}}\left(X\right)=\sqrt{\underset{k}{\mathrm{\Σ}}{(x_{\mathrm{ik}}-x_{\mathrm{jk}})}^2}---\left(4\right)$

同对应最短路径距离的均方差最小。即

$X^{\left(0\right)}=\underset{X}{\min }\mathrm{Stress}(D_s,X)=\sqrt{\underset{i,j}{\mathrm{\Σ}}{(d_{\mathrm{ij}}-d_{\mathrm{ij}}^{\mathrm{euc}}\left(X\right))}^2}---\left(5\right)$

步骤四：建立代价函数f和Stress。此处我们采用最大似然估计的方法，根据式(6)寻找X和

使得对数似然函数

$\mathrm{Stress}(X,\hat{D})=\underset{i}{\mathrm{\Σ}}\underset{j}{\mathrm{\Σ}}g(d_{\mathrm{ij}}^{\mathrm{euc}},{\hat{d}}_{\mathrm{ij}})=\underset{i}{\mathrm{\Σ}}\underset{j}{\mathrm{\Σ}}\ln f\left({\hat{d}}_{\mathrm{ij}}\right|d_{\mathrm{ij}}^{\mathrm{euc}})$

$=\underset{i}{\mathrm{\Σ}}\underset{j}{\mathrm{\Σ}}(\ln \frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_d{\hat{d}}_{\mathrm{ij}}}-\frac{{(\log {\hat{d}}_{\mathrm{ij}}-{\log d}_{\mathrm{ij}}^{\mathrm{euc}})}^2}{{2\mathrm{\σ}}_d^2})---\left(6\right)$

最小。此处由于同时优化X和

的复杂度太高，在本发明中，我们采用循环优化的做法，轮流优化X和

使得它们趋向最优值。整个循环分为以下四步。

步骤五：确定初值。对于由X⁽⁰⁾构建的欧式距离矩阵 $D^{\mathrm{euc}}\left(X^{\left(0\right)}\right)=\left\{d_{\mathrm{ij}}^{\mathrm{euc}}\left(X^{\left(0\right)}\right)\right\},$ 找到最优的节点间距离矩阵

使得代价函数

$S^{\left(0\right)}=\mathrm{Stress}({\hat{D}}^{\left(0\right)},D^{\mathrm{euc}}\left(X^{\left(0\right)}\right))---\left(7\right)$

最小。

迭代循环开始，交替优化X和

步骤六：在第t次循环中，优化X^(t)，使其与前一次迭代的距离矩阵之间的代价函数最小，即

$X^{\left(t\right)}=\underset{X}{\min }\mathrm{Stress}({\hat{D}}^{(t-1)},D^{\mathrm{euc}}\left(X\right))---\left(8\right)$

在这一步中，由于对X^(t)是无约束优化，所以可以采用常用的牛顿优化算法，令Stress函数对于X^(t)中的每一项求偏导均为零，从而列出联立多元方程组求解。

步骤七：优化

使其与D^euc(X^(t))之间的代价函数最小，即

${\hat{D}}^{\left(t\right)}=\underset{\hat{D}}{\min }\mathrm{Stress}(\hat{D},D^{\mathrm{euc}}\left(X^{\left(t\right)}\right))---\left(9\right)$

同步骤六中的优化不同，此处

是有约束优化，

中的元素必须满足对应最小路径距离之间所呈现的顺序关系。牛顿优化方法方法并不适用于这类问题。所以，此处我们采用相邻模块比较算法(UDP：up-and-down-blocks algorithm)，利用单调回归的思想来求解。Kruscal已经证明了该算法可以在满足顺序关系的前提下，找到使得(9)中的代价函数最小的

简单描述UDP算法，这是一个搜索更新式的方法。

首先，将D^euc(X^(t))赋给此时(9)中的代价函数达到最小值为0，但是D^euc(X^(t))中元素的顺序关系同D_s中并不完全一致，所以需要继续调整。

然后，将最短路径矩阵D_s中的所有项按照数值大小以升序排列成一个一维数组A，如果网络中共有M个传感器节点，那么A的长度就是

按照相同的下标顺序也将

中的元素排列成一个一维数组Q。

接下来，开始对于

搜索更新过程。

在任何一个位置q，如果顺序关系成立，即

A(q-1)＝d_kl≤A(q)＝d_ij

$Q(q-1)={\hat{d}}_{\mathrm{kl}}^{\left(t\right)}\≤Q\left(q\right)=d_{\mathrm{ij}}^{\left(t\right)}---\left(10\right)$

则不对

和d_ij ^(t)做调整，继续向前搜索。

如果在位置q，顺序关系被破坏了，即

A(q-1)＝d_kl≤A(q)＝d_ij

$Q(q-1)={\hat{d}}_{\mathrm{kl}}^{\left(t\right)}>Q\left(q\right)=d_{\mathrm{ij}}^{\left(t\right)}---\left(11\right)$

此时就需要对

和d_ij ^(t)做一定的调整更新，

要减小同时d_ij ^(t)要增大，使得更新后的数值满足相应的顺序关系。此处需要注意到，如果我们单纯调整

和d_ij ^(t)的大小，那么

更新后的数值可能小于它之前的元素，即Q数组上位置在q-1之前的节点对距离值。所以，为了保证Q在位置q-1之前的顺序关系仍然得以保持，一些

之前的元素也需要进行调整。

最后，根据以上要求，我们得出的最佳调整公式为

${\hat{d}}_{\mathrm{ij}}^{\left(t\right)}=Q^{\mathrm{new}}\left(q\right)=...=Q^{\mathrm{new}}(q-k)={(Q^{\mathrm{old}}\left(q\right)...Q^{\mathrm{old}}(q-k))}^{\frac{1}{k+1}}---\left(12\right)$

其中Q^new代表调整后的数值，Q^old代表调整前的数值。而k是满足式(13)的最小值

Q^new(q-k)≥Q^old(q-k-1)                        (13)

此时，更新过的数组Q在位置q之前的元素已经完全满足了对应最短路径距离之间的顺序关系，算法继续往前搜索。

当UDP算法搜索达到了数组A和Q的末尾，对于

的更新就完成了。表1中的例子说明了在一个4个节点(6条距离)的网络中，UDP如何更新

得到最优解，每一步中被调整了的节点间距离用斜体标出。

步骤八：更新完得到

计算此时的代价函数

$S^{\left(t\right)}=\mathrm{Stress}({\hat{D}}^{\left(t\right)},D^{\mathrm{euc}}\left(X^{\left(t\right)}\right))---\left(14\right)$

如果|S^(t)-S^(t-1)|＞ε，回到步骤六进入下一次迭代；否则，迭代收敛，

和X^(t)就是相对定位阶段的最终结果

和X_re。我们可以保证是满足最短路径距离D_s中的所有顺序关系的，在这个前提下，我们找到了使得代价函数最小化的相对定位结果X_re。

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

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

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

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

X_ab＝sX_reK+1t^T                     (16)

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

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

由于传感器数量众多，基于成本、体积和功耗方面的考虑，在本专利的算法实施中，所有传感器节点不装配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，这一特性对研究传感器网络的多信道协议十分有利。

在上述的应用场景和硬件平台上，RMLE可以按照图8中所示流程由以下6个步骤实现。

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

1.节点布撒

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

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

2.路由建立和数据传输

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

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

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

3.网络连接信息的提取

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

4.相对坐标定位

利用上一步中的节点间测距信息，进行全网相对定位。定位方法如上一章所示。

步骤一：每个节点均以固定发射功率向全网广播数据包，如果任一节点i收到来自节点j的数据包，则表明节点i和j之间存在一跳数据通路，两者互为邻居节点，可以利用RSS功率检测法测得互相之间的距离。对于非邻居节点i和j，我们采用网络约束，对它们之间的距离进行精确推断，可以得到任一未知距离的概率密度函数如式(13)所示。由于所有分布的形式均相同，所以只需记录每一个概率密度分布的参数σ_d，供后面的优化使用。

步骤二：建立最短路径矩阵。对于非邻居节点i和j，我们采用弗洛伊德最短路径算法求解它们之间的距离估计，通过遍历式的搜索找到i和j之间的所有多跳数据通路，并且选取其中路径总和最短的一条作为i和j之间的最短路径，对应的距离即是非邻居节点间的最短路径距离。通过对全网的遍历式搜索，我们可以得到任意两个节点(相邻或非相邻)的最短路径距离，即最短路径距离矩阵D_s。

将最短路径矩阵中的所有元素按照大小进行完全的排序，就得到了该定位方法中所需要的距离顺序关系，可以用一个长度为(M-1)(M-2)/2的三维数组N存放该顺序关系，其中N的第一维记录对应元素的值，第二维记录对应元素在D_s中的行编号，第三维记录列编号。

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

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

定义

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

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

B＝-HD_sH＝HX^TXH            (19)

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

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

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

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

即为多维标度分析算法所得到的初始定位结果X⁽⁰⁾。

步骤四：建立代价函数如式(6)所示。由于代价函数形式已知，此处只需要记录关键变量σ_d即可

步骤五：根据前述的UDP算法建立数组A和Q，通过完整的搜索和位置交换，找到最优的节点间距离矩阵

使得其         在满足N中的顺序大小关系的前提下，令式(7)最小。

步骤六：采用牛顿优化法，优化节点位置坐标X^(t)。在每个优化点寻找Q函数上升最快的方向，最终达到式(8)的最大值。

步骤七：同步骤五，找到最优的节点间距离矩阵

使得其在满足N中的顺序大小关系的前提下，令式(9)最小。

步骤八：根据计算得到的和X^(t)代入式(14)，判断迭代是否收敛，以决定是否继续进行相对定位循环。

5.绝对坐标变换

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

6.结果输出

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

表1

Claims (2)

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 }_{\mathrm{ij}}-{\log \mathrm{\δ}}_{\mathrm{ij}})}^2}{{2\mathrm{\σ}}_d^2}}$

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

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

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

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

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

$L(X/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)，把步骤(2.2.2)中的最短距离矩阵D_s中的所有元素按值得大小进行排序，得到所述元素的距离大小的距离顺序关系，并用一个长度为(M-1)(M-2)/2的三维数组N存放所述顺序关系；其中，第一维记录对应元素的值，第二维及第三维分别记录该元素在D_s中的行和列编号，

步骤(2.3)，按下式从所建最短距离矩阵D_s得到一个节点位置的初始拓扑估计X⁽⁰⁾作为迭代的初值，

$X^{\left(0\right)}=\underset{X}{\min }\mathrm{Stress}(D_s,X)=\sqrt{\underset{i,j}{\mathrm{\Σ}}{(d_{ij}-d_{\mathrm{ij}}^{\mathrm{euc}}\left(X\right))}^2}$

d_ij ^euc(X)为反解的欧式距离

$d_{\mathrm{ij}}^{\mathrm{euc}}\left(X\right)=\sqrt{\underset{k}{\mathrm{\Σ}}{(x_{\mathrm{ik}}-x_{\mathrm{jk}})}^2}$

步骤(2.4)，建立作为优化标准的节点距离矩阵和节点位置矩阵X的对数似然函数：

$\mathrm{Stress}(X,\hat{D})=\underset{i}{\mathrm{\Σ}}\underset{j}{\mathrm{\Σ}}(\ln \frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_d{\hat{d}}_{\mathrm{ij}}}-\frac{{(\log {\hat{d}}_{\mathrm{ij}}-\log d_{\mathrm{ij}}^{\mathrm{euc}})}^2}{{2\mathrm{\σ}}_d^2})$

再使

最小，其步骤依次如下：

步骤(2.4.1)，对于由X⁽⁰⁾构建的欧式距离矩阵{d_ij ^euc(X⁽⁰⁾)}，找到最优的节点间距离矩阵

使代价函数S⁽⁰⁾最小：

$S^{\left(0\right)}=\mathrm{Stress}({\hat{D}}^{\left(0\right)},D^{\mathrm{euc}}\left(X^{\left(0\right)}\right)),$ 并开始迭代循环，其步骤如下：

步骤(2.4.1.1)，使用牛顿优化法，在第t次循环中，优化X^(t)，使该X^(t)与前一次迭代的距离矩阵

之间的代价函数最小，即

$X^{\left(t\right)}=\underset{X}{\min }\mathrm{Stress}({\hat{D}}^{(t-1)},D^{\mathrm{euc}}\left(X\right))$

步骤(2.4.1.2)，利用相邻模块比较方法来优化

使该

与D^euc(X^(t))之间的代价函数最小，即：

${\hat{D}}^{\left(t\right)}=\underset{\hat{D}}{\min }\mathrm{Stress}(\hat{D},D^{\mathrm{euc}}\left(X^{\left(t\right)}\right))$

步骤(2.4.1.3)，在第t次迭代中，得到

此时的代价函数为S^(t)，

$S^{\left(t\right)}=\mathrm{Stress}({\hat{D}}^{\left(t\right)},D^{\mathrm{euc}}\left(X^{\left(t\right)}\right))$

若：|S^(t)-S^(t-1)|＞ε，ε为设定值，则进入下一次迭代，优化X^(t)，

若：|S^(t)-S^(t-1)|＜ε，则收敛，所得到的

和X^(t)就是相对定位阶段的最终结果

和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；。

2.根据权利要求1所述的基于距离大小的顺序关系的传感器网络节点自定位方法，其特征在于，所建的相邻模块比较方法依次含有以下步骤：

步骤(a)，把D^euc(X^(t))赋给

步骤(b)，把最短路径矩阵D_s中的所有项按照数值大小的顺序形式排列成一个一维数组A，即是前面所述三维数组N中的第一维，其长度为(M-1)(M-2)/2，M为所建传感器网络中的传感器数，按照相同的下标顺序，也把优化的距离矩阵

中的元素排列成一个一维数组Q；

步骤(c)，对

进行搜索更新：

在任何一个位置q，若下式成立，则不对

和d_ij ^(t)做调整，继续向前搜索：

A(q-1)＝d_kl小于等于A(q)＝d_ij

否则，便在保证位置q-1之前的顺序关系下减小增加d_ij ^(t)，在必要的时候调整

之前的元素；

步骤(d)，根据步骤(b)和步骤(c)得到的结果，得到下述最佳调整公式：

${\hat{d}}_{\mathrm{ij}}^{\left(t\right)}=Q^{\mathrm{new}}\left(q\right)=...=Q^{\mathrm{new}}(q-k)={(Q^{\mathrm{old}}\left(q\right)...Q^{\mathrm{old}}(q-k))}^{\frac{1}{k+1}}$

其中，Q^new代表调整后的数值，Q^old代表调整前的数值；

所述的k值满足如下条件：Q^new(q-k)≥Q^old(q-k-1)；

步骤(e)，重复步骤(c)到步骤(d)。

Images