2-D Direction Finding of Cyclostationary Signals with Polarization Smoothing Algorithm in Multipath Environment

1.   引言

随着微机电系统(Micro-Electro-Mechanical Systems, MEMS)和无线通信技术的快速发展，使低成本，低功耗的传感器大规模部署成为可能^[1-3]。无线传感器网络(Wireless Sensor Network, WSN)由大量传感器构成，能够检测各种物理信息，广泛应用于医疗(病人检测)、环境(大鸭岛实验)、家庭(用水检测)等领域^[2,4,5]。在这些应用中，传感器节点感知到的信息与其位置结合起来才有意义。因此，传感器节点的定位问题越来越受到人们的重视。尽管可以在传感器上安装全球定位系统(Global Positioning System, GPS)模块来定位。但GPS模块成本高，功耗大，会增加大规模部署成本。为了解决这一问题，只对部分节点安装GPS模块，其余节点采用合作定位方法来估计。合作定位所需的距离信息可以通过接收信号强度(Received-Signal-Strength, RSS)、到达时间差(Time-Difference-Of-Arrival, TDOA)、到达角(Angle-Of-Arrival, AOA)等技术来获取^[3,4,6]。

使用合作定位方法时，在WSN中随机选取部分传感器安装GPS模块，以获取精确位置。将事先已知位置的传感器节点称为已知位置(Location-Aware, LA)节点或锚节点，对于其它未知位置的传感器，称为未知位置(Location-Unaware, LU)节点。使用测距技术获得节点之间的相互距离，之后采用合作定位方法对LU节点进行定位。如：三边定位算法、三角定位算法、极大似然估计算法以及优化算法等^[7]。其中，优化算法是比较典型的定位算法之一。优化算法将定位问题归结为非凸的目标函数优化问题，针对这一问题，许多做法是将非凸的优化问题转化为凸优化问题进行求解^[1,2,8,9]。优化算法可以分为集中式算法和分布式算法。集中式算法使用一个中央处理器收集节点之间的距离信息并进行定位计算。文献[2]提出了一种基于半正定规划(Semi-Definite Programming, SDP)的集中式算法。文中通过对目标函数添加正则项来降低解的秩，然后使用梯度下降法细化节点的位置，提高了定位精度。但随着网络规模的增大，网络通信量和中心节点的计算量变高，会超出中央处理器的处理能力^[3]。因此，集中式算法扩展性较低^[10]，只适合小规模的WSN。与集中式算法相比，分布式算法将处理器分配到每个节点上，每个处理器仅收集局部节点之间的距离信息并进行计算，有效降低了通信量和计算量，有良好的扩展性^[11]，可用于无中央处理器或大规模的WSN。但缺点是会导致误差的传播，降低定位精度^[1,3,6,9]。文献[1]中使用2阶锥规划(Second-Order Cone Programming, SOCP)松弛技术将非凸的优化问题松弛为凸的优化问题。使用LU节点和局部距离信息将WSN划分为若干个子区域，并提出了一种分布式算法，有良好的扩展性，能用于大规模网络。但此算法对网络边界的LU节点定位效果不理想。文献[9]中使用凸松弛方法来近似非凸的最大似然公式，利用优化特性提出了一种分布式算法，每个节点使用最优梯度法进行收敛，降低了通信量。但这种方法需要LU节点位于LA节点的凸包中，才能获得较精确的定位。这在实际情况中难以满足，导致网络边界定位效果不理想。文献[10]中提出了一种分布式定位算法，将锚节点和与其直接相连的传感器节点构成一个子区域，对每个子区域使用SDP松弛方法进行求解。然后将估计误差较低的LU节点作为伪锚节点，用来估计剩余的LU节点。但这种做法会导致误差的积累和传播。

本文的主要目的是解决大规模WSN中节点难以定位问题，提出了一种基于改进牛顿法的分布式定位算法，本算法包括网络划分和分布式算法。其创新点包括：(1)本文将WSN划分为部分重叠的子区域，从而建立了各个子区域的定位问题，可以独立优化求解。降低了通信量和计算量，提高了扩展性，使其能够用于大规模WSN；(2)分布式算法的每一步迭代包括子区域节点的定位求解和融合平均。首先，各个子区域定位问题采用基于Hesse矩阵修正的改进牛顿法求解，使其能够快速收敛；然后，对重叠子区域内的节点，进行局部融合求平均，能够修正部分子区域求解不够准确的问题，提高了区域边界LU节点的定位精度；(3)实验结果表明，本文所提分布式算法，能满足大规模WSN中节点定位需求，而且在不同的情形下，与文献[1,9]中提出的分布式算法相比，在大规模网络中定位精度更高。

2.   定位问题的描述

在一片区域内随机部署一定数量的传感器，假设WSN中的传感器有相同的发射功率，则它们的通信半径相同^[1]。如果不同传感器节点之间的距离小于$r$，则认为它们之间可以相互通信，能够获得传感器之间的相互距离。使用距离信息，可将WSN连接成一个网状图。因此，一个区域内的WSN可以构成一个无向图$G = ({\cal{V}}, {\cal{E}})$。其中，集合${\cal{V}}$表示传感器节点集合；集合${\cal{E}}$表示节点之间边的集合。

无线传感器网络的定位问题可以描述为：在一片区域${R^d}(d \ge 1)$内随机部署$N$个传感器，${R^d}$表示$d$维欧几里得空间^[8]。在区域内有$m$个LU节点，$n$个LA节点。采用合作定位方法^[12-14]，使用LA节点位置和节点之间的距离信息，来估计出$m$个LU节点的位置^[1]。本文中，传感器分布在${[ - 0.5, 0.5]^2}$的单位正方形区域内^{[1, 6]}。LU节点表示为${{\text{p}}_1}, {{\text{p}}_2}, ·\!·\!· , {{\text{p}}_m}$，其中，${{\text{p}}_i} = {[{x_i}, {y_i}]^{\rm{T}}} \in {R^2}$，表示LU节点$i$的坐标；LA节点表示为${{\text{a}}_1}, {{\text{a}}_2}, ·\!·\!· , {{\text{a}}_n}$，其中，${{\text{a}}_k} = {[{x_k}, {y_k}]^{\rm{T}}} \in {R^2}$，表示LA节点$k$的坐标。用${d_{ij}}$表示LU节点$i$和节点$j$之间带噪声的测距，${d_{ij}} = {d_{ji}}$^[9]; ${d_{ik}}$表示LU节点$i$和LA节点$k$之间带噪声的测距^{[2, 6]}; ${\omega _{ij}}$和${\omega _{ik}}$是根据距离的反比取的归一化权重。$(i, j) \!\in\! {\cal{A}}$，表示LU节点之间可测得距离的节点集合；$(i, k) \in {\cal{B}}$，表示LU与LA节点之间可测得距离的节点集合。因此，使用带噪声的测距信息将定位问题归结为一个无约束的高度非线性非凸的优化问题^[2]

$$ \begin{align} {\rm{}}&\mathop {{\rm{min}}}\limits_{{{\text{p}}_1}, {{\text{p}}_2}, \cdots , {{\text{p}}_m}} {\sum\limits_{(i, j) \in {\cal{A}}} {\omega _{ij}^2\left| {\left\| {{{\text{p}}_i} - {{\text{p}}_j}} \right\|_2^2 - d_{ij}^2} \right|} ^2} \\ {\rm{}}& \quad + \sum\limits_{(i, k) \in {\cal{B}}} {\omega _{ik}^2{{\left| {\left\| {{{\text{p}}_i} - {{\text{a}}_k}} \right\|_2^2 - d_{ik}^2} \right|}^2}} \end{align} $ $

(1)

3.   分布式定位算法的设计

3.1   划分子区域

通过观察式(1)可以看出，目标函数仅取决于通信半径$r$内与LU节点直接相连的邻居节点位置以及节点之间的距离信息。因此，以LU节点为中心，可将WSN构成的无向图$G = ({\cal{V}}, {\cal{E}})$分解为重叠的子图。

$$G = \bigcup\limits_{s = 1}^m {{G_s}} $ $

(2)

其中，${G_s} = ({{\cal{V}}_s}, {{\cal{E}}_s})$, $s = 1, 2, ·\!·\!· , m$表示第$s$个子图；${{\cal{V}}_s}$表示子图$s$中所有节点的集合；${{\cal{E}}_s}$表示子图$s$中节点之间边的集合，包含LU节点之间直接相连的边和LU节点与LA节点之间直接相连的边。

3.2   子区域优化

将WSN划分为重叠的子区域后，将目标函数式(1)分解为子区域优化问题。之后使用分布式算法对子区域进行优化求解和局部融合。本文采用基于单位步长的改进牛顿法对每个子区域进行优化求解，这种方法需要一个好的初始值^[15]，以便沿着目标函数下降的方向进行搜索，从而更快达到目标函数的局部最小值。因此，本文采用三边定位算法的粗略估计值作为改进牛顿法的初始值。

3.2.1   子区域问题的描述

使用子区域中节点位置和节点之间的距离信息，构造出子区域的目标函数，进而采用分布式算法进行计算。与集中式算法相比，提高了扩展性，降低了通信量和计算量，可以有效解决大规模WSN中的定位问题。

子区域的目标函数归结为

$$ \begin{align} {\rm{}}&\mathop {{\rm{min}}}\limits_{{{\text{p}}_{{s_1}}}, {{\text{p}}_{{s_2}}}, \cdots , {{\text{p}}_{{s_l}}}} {\sum\limits_{({s_i}, {s_j}) \in {{\cal{A}}_s}} {\omega _{{s_i}{s_j}}^2\left| {\left\| {{{\text{p}}_{{s_i}}} - {{\text{p}}_{{s_j}}}} \right\|_2^2 - d_{{s_i}{s_j}}^2} \right|} ^2} \\ {\rm{}}& \quad+ \sum\limits_{({s_i}, {s_k}) \in {{\cal{B}}_s}} {\omega _{{s_i}{s_k}}^2{{\left| {\left\| {{{\text{p}}_{{s_i}}} - {{\text{a}}_{{s_k}}}} \right\|_2^2 - d_{{s_i}{s_k}}^2} \right|}^2}} \end{align} $ $

(3)

其中，$s = 1, 2,·\!·\!· , m$，表示第$s$个子区域；${{\text{p}}_{{s_i}}} = {[{x_{{s_i}}}, {y_{{s_i}}}]^{\rm{T}}}{\rm{ }}(i = 1, 2, ·\!·\!· , l)$，表示子区域$s$中LU节点$i$的坐标；${{\text{a}}_{{s_k}}} = {[{x_{{s_k}}}, {y_{{s_k}}}]^{\rm{T}}}{\rm{ }}(k = 1, 2, ·\!·\!· , h)$，表示子区域$s$中LA节点$k$的坐标；$({s_i}, {s_j}) \in {{\cal{A}}_s}$，表示子区域$s$中LU节点之间可测得距离的节点集合；$({s_i}, {s_k}) \in {{\cal{B}}_s}$，表示子区域$s$中LU与LA节点之间可测得距离的节点集合；${d_{{s_i}{s_j}}}$和${d_{{s_i}{s_k}}}$表示带噪声的测距，噪声模型如下^{[2, 6]}。${\rm{nf}}$是噪声因子，用于控制噪声强度，取值$[0, 1]$; ${\varepsilon _{{s_i}{s_j}}}$和${\varepsilon _{{s_i}{s_k}}}$是产生的随机噪声，是$( - 1, 1)$之间均匀的随机数。

$${d_{{s_i}{s_j}}} = {\left\| {{{\text{p}}_{{s_i}}} - {{\text{p}}_{{s_j}}}} \right\|_2} \cdot \left| {1 + {\rm{nf}} \cdot {\varepsilon _{{s_i}{s_j}}}} \right|$ $

(4)

$${d_{{s_i}{s_k}}} = {\left\| {{{\text{p}}_{{s_i}}} - {{\text{a}}_{{s_k}}}} \right\|_2} \cdot \left| {1 + {\rm{nf}} \cdot {\varepsilon _{{s_i}{s_k}}}} \right|$ $

(5)

$$ \left| {{\varepsilon _{{s_i}{s_j}}}} \right|{\text{ 或 }}\left| {{\varepsilon _{{s_i}{s_k}}}} \right| < 1 \hspace{71pt} $ $

(6)

式(3)中，${\omega _{{s_i}{s_j}}}$和${\omega _{{s_i}{s_k}}}$是根据子区域$s$中距离的反比取的归一化权重，因为传感器节点之间随着距离的增加，测得的距离的可信度降低，误差会增加，其权重应该越小；反之，传感器节点之间距离越近，可信度越高，相应的权重应该越大，权重分别为

$$ {\omega _{{s_i}{s_j}}} = \frac{{d_{{s_i}{s_j}}^{ - 1}}}{{\displaystyle\sum\limits_{({s_i}, {s_j}) \in {{\cal{A}}_s}} {d_{{s_i}{s_j}}^{ - 1}} + \displaystyle\sum\limits_{({s_i}, {s_k}) \in {{\cal{B}}_s}} {d_{{s_i}{s_k}}^{ - 1}} }} $ $

(7)

$$ {\omega _{{s_i}{s_k}}} = \frac{{d_{{s_i}{s_k}}^{ - 1}}}{{\displaystyle\sum\limits_{({s_i}, {s_j}) \in {{\cal{A}}_s}} {d_{{s_i}{s_j}}^{ - 1}} + \displaystyle\sum\limits_{({s_i}, {s_k}) \in {{\cal{B}}_s}} {d_{{s_i}{s_k}}^{ - 1}} }} $ $

(8)

使用${\text{p}} = {[{\text{p}}_{{s_1}}^{\rm{T}}{\rm{ }}{\text{p}}_{{s_2}}^{\rm{T}}{\rm{ }} ·\!·\!· {\rm{ }}{\text{p}}_{{s_l}}^{\rm{T}}]^{\rm{T}}}$，表示子区域$s$中所有LU节点的坐标构成的列向量，其中${\text{p}}_{{s_i}}^{\rm{T}} = [{x_{{s_i}}}, {y_{{s_i}}}]$; ${\text{a}} = {[{\text{a}}_{{s_1}}^{\rm{T}}{\rm{ }}{\text{a}}_{{s_2}}^{\rm{T}}{\rm{ }} ·\!·\!· {\rm{ }}{\text{a}}_{{s_h}}^{\rm{T}}]^{\rm{T}}}$，表示子区域$s$中所有LA节点的坐标构成的列向量，其中${\text{a}}_{{s_k}}^{\rm{T}} = [{x_{{s_k}}}, {y_{{s_k}}}]$。令${{\text{e}}_{i1}}$和${{\text{e}}_{i2}}$分别表示$2l \times 2l$单位矩阵的第$2i - 1$列和$2i$列；${{\text{e}}_{k1}}$和${{\text{e}}_{k2}}$分别表示$2h \times 2h$单位矩阵的第$2k - 1$列和$2k$列。则，${x_{{s_i}}} = {\text{e}}_{i1}^{\rm{T}}{\text{p}}$, ${y_{{s_i}}} = {\text{e}}_{i2}^{\rm{T}}{\text{p}}$, ${x_{{s_k}}} = {\text{e}}_{k1}^{\rm{T}}{\text{a}}$, ${y_{{s_k}}} = {\text{e}}_{k2}^{\rm{T}}{\text{a}}$，传感器节点之间的真实距离可写为

$$ \left\| {{{\text{p}}_{{s_i}}} - {{\text{p}}_{{s_j}}}} \right\|_2^2 = {{\text{p}}^{\rm{T}}}{\text{Ap}} \hspace{89pt} $ $

(9)

$$ \begin{align} \left\| {{{\text{p}}_{{s_i}}} - {{\text{a}}_{{s_k}}}} \right\|_2^2 =& {{\text{p}}^{\rm{T}}}{\text{Bp}} - {{\text{p}}^{\rm{T}}}{\text{Ca}} - {{\text{a}}^{\rm{T}}}{\text{Dp}} \\ & + {{\text{a}}^{\rm{T}}}{\text{Ea}} \end{align} $ $

(10)

其中，

$$ \begin{align} {\text{A}} =& ({{\text{e}}_{i1}} - {{\text{e}}_{j1}}){({{\text{e}}_{i1}} - {{\text{e}}_{j1}})^{\rm{T}}} \\ & + ({{\text{e}}_{i2}} - {{\text{e}}_{j2}}){({{\text{e}}_{i2}} - {{\text{e}}_{j2}})^{\rm{T}}} \end{align} $ $

(11)

$$ {\text{B}} = {{\text{e}}_{i1}}{\text{e}}_{i1}^{\rm{T}} + {{\text{e}}_{i2}}{\text{e}}_{i2}^{\rm{T}} \hspace{44pt} $ $

(12)

$$ {\text{C}} = {{\text{e}}_{i1}}{\text{e}}_{k1}^{\rm{T}} + {{\text{e}}_{i2}}{\text{e}}_{k2}^{\rm{T}} \hspace{43pt} $ $

(13)

$$ {\text{D}} = {{\text{e}}_{k1}}{\text{e}}_{i1}^{\rm{T}} + {{\text{e}}_{k2}}{\text{e}}_{i2}^{\rm{T}} \hspace{42pt} $ $

(14)

$$ {\text{E}} = {{\text{e}}_{k1}}{\text{e}}_{k1}^{\rm{T}} + {{\text{e}}_{k2}}{\text{e}}_{k2}^{\rm{T}} \hspace{41pt} $ $

(15)

目标函数式(3)可以写为

$$ \begin{aligned} f({\text{p}}) =& \mathop {{\rm{min}}}\limits_{{{\text{p}}_{{s_1}}}, {{\text{p}}_{{s_2}}}, \cdots , {{\text{p}}_{{s_l}}}} {\sum\limits_{({s_i}, {s_j}) \in {{\cal{A}}_s}} {\omega _{{s_i}{s_j}}^2({{\text{p}}^{\rm{T}}}{\text{Ap}} - d_{{s_i}{s_j}}^2)} ^2} \\ {\rm{}}&+\sum\limits_{({s_i}, {s_k}) \in {{\cal{B}}_s}} \omega _{{s_i}{s_k}}^2({{\text{p}}^{\rm{T}}}{\text{Bp}} - {{\text{p}}^{\rm{T}}}{\text{Ca}} - {{\text{a}}^{\rm{T}}}{\text{Dp}} \\ {\rm{}}& + {{{\text{a}}^{\rm{T}}}{\text{Ea}} - d_{{s_i}{s_k}}^2)}^2 \end{aligned} $ $

(16)

本文使用基于单位步长的改进牛顿法对子区域进行优化，在给出具体算法之前，先计算出目标函数$f({\text{p}})$的梯度向量$\nabla f({\text{p}})$和Hesse矩阵${\nabla ^2}f({\text{p}})$

$$ \begin{aligned} \nabla f({\text{p}}) =& \sum\limits_{({s_i}, {s_j}) \in {{\cal{A}}_s}} {4\omega _{{s_i}{s_j}}^2({{\text{p}}^{\rm{T}}}{\text{Ap}} - d_{{s_i}{s_j}}^2)({\text{Ap}})} \\ {\rm{}}&+\sum\limits_{({s_i}, {s_k}) \in {{\cal{B}}_s}} 2\omega _{{s_i}{s_k}}^2({{\text{p}}^{\rm{T}}}{\text{Bp}} - {{\text{p}}^{\rm{T}}}{\text{Ca}} - {{\text{a}}^{\rm{T}}}{\text{Dp}} \\ {\rm{}} & + {{\text{a}}^{\rm{T}}}{\text{Ea}} - d_{{s_i}{s_k}}^2)(2{\text{Bp}} - {\text{Ca}} - {{\text{D}}^{\rm{T}}}{\text{a}}) \end{aligned} $ $

(17)

$$ \begin{aligned} {\nabla ^2}f({\text{p}}) = & \sum\limits_{({s_i}, {s_j}) \in {{\cal{A}}_s}} 4\omega _{{s_i}{s_j}}^2\left[ 2{\text{Ap}}{{\text{p}}^{\rm{T}}}{{\text{A}}^{\rm{T}}}\right. \\ {\rm{}} & \left.+ ({{\text{p}}^{\rm{T}}}{\text{Ap}} - d_{{s_i}{s_j}}^2){{\text{A}}^{\rm{T}}} \right] \\ {\rm{}} & + \sum\limits_{({s_i}, {s_k}) \in {{\cal{B}}_s}} 2\omega _{{s_i}{s_k}}^2\left[ (2{\text{Bp}} - {\text{Ca}} - {{\text{D}}^{\rm{T}}}{\text{a}})\right.\\ {\rm{}} & \cdot {{(2{\text{Bp}} - {\text{Ca}} - {{\text{D}}^{\rm{T}}}{\text{a}})}^{\rm{T}}} + 2({{\text{p}}^{\rm{T}}}{\text{Bp}} \\ {\rm{}}& \left.- {{\text{p}}^{\rm{T}}}{\text{Ca}} - {{\text{a}}^{\rm{T}}}{\text{Dp}} + {{\text{a}}^{\rm{T}}}{\text{Ea}} - d_{{s_i}{s_k}}^2){{\text{B}}^{\rm{T}}} \right]\\ {\rm{}} & \hspace{168pt}(18) \end{aligned} $ $

3.2.2   三边定位获取初始值

在已有的定位算法中，三边定位算法是比较简单的，可以用作LU节点初始位置的粗略估计。具体做法为：假设待定的LU节点D，坐标为$(x, y)$，周围有3个不共线的LA节点A, B, C，坐标分别为$({x_1}, {y_1})$, $({x_2}, {y_2})$, $({x_3}, {y_3})$。3个LA节点到LU节点的距离分别为${R_1}$, ${R_2}$, ${R_3}$，均小于通信半径$r$。分别以A, B, C为圆心，${R_1}$, ${R_2}$, ${R_3}$为半径，画3个相交的圆，相交点即为LU节点D的位置。

对应的方程组为

$$\left. {\begin{array}{*{20}{c}} {{{(x - {x_1})}^2} + {{(y - {y_1})}^2} = R_1^2} \\ {{{(x - {x_2})}^2} + {{(y - {y_2})}^2} = R_2^2} \\ {{{(x - {x_3})}^2} + {{(y - {y_3})}^2} = R_3^2} \end{array}} \right\}$ $

(19)

方程的解为LU节点D的坐标

$$ \begin{align} \left[ \begin{gathered} x \\ y \\ \end{gathered} \right] =& {\left[ \begin{gathered} 2{\rm{(}}{x_3} - {x_1}{\rm{) }}2({y_3} - {y_1}) \\ 2{\rm{(}}{x_3} - {x_2}{\rm{) 2}}({y_{\rm{3}}} - {y_2}) \\ \end{gathered} \right]^{ - 1}}\\ {\rm{}}& \cdot\left[ \begin{gathered} x_3^2 - x_1^2 + y_3^2 - y_1^2 + R_1^2 - R_3^2 \\ x_3^2 - x_2^2 + y_3^2 - y_2^2 + R_2^2 - R_3^2 \\ \end{gathered} \right] \end{align} $ $

(20)

然而，在实际情况中，传感器节点的部署是随机的，且受到$r$的限制，无法保证每个LU节点周围都有3个或3个以上的LA节点。因此，本文按照以下方法来获得LU节点的初始位置，作为基于单位步长改进牛顿法的初始值：当LU节点周围有3个或3个以上LA节点时，利用最近的3个LA节点，用三边定位方法粗略定出其初始位置；当LU节点周围有1～2个LA节点时，则将离LU节点最近的LA节点作为其初始位置；当LU节点周围没有LA节点时，则将传感器分布区域的中心作为其初始位置。

3.2.3   基于单位步长的改进牛顿法

本文采用基于单位步长改进牛顿法对子区域进行优化的原因是：首先，目标函数的Hesse矩阵不一定是正定的，有可能是奇异的，可能会导致位置的更新往错误的方向移动，使定位出现偏差。因此，本文将Hesse矩阵修正为正定的。其次，在将Hesse矩阵修正为正定后，基于单位步长的改进牛顿法可以快速收敛到局部最优解^[16]。

观察目标函数的Hesse矩阵，包含4项求和。其中，$2{\text{Ap}}{{\text{p}}^{\rm{T}}}{{\text{A}}^{\rm{T}}}$和$(2{\text{Bp}} - {\text{Ca}} - {{\text{D}}^{\rm{T}}}{\text{a}}) {(2{\text{Bp}} -}$${{\text{Ca}} - {{\text{D}}^{\rm{T}}}{\text{a}})^{\rm{T}}}$是正定的，但$({{\text{p}}^{\rm{T}}}{\text{Ap}} - d_{{s_i}{s_j}}^2){{\text{A}}^{\rm{T}}}$和$2({{\text{p}}^{\rm{T}}}{\text{Bp}} - {{\text{p}}^{\rm{T}}}{\text{Ca}} - {{\text{a}}^{\rm{T}}}{\text{Dp}} + {{\text{a}}^{\rm{T}}}{\text{Ea}} - d_{{s_i}{s_k}}^2){{\text{B}}^{\rm{T}}}$不一定是正定的。因此，将Hesse矩阵修正为正定的

$$ \begin{align} {\nabla ^2}f({\text{p}}) =& \sum\limits_{({s_i}, {s_j}) \in {{\cal{A}}_s}} {8\omega _{{s_i}{s_j}}^2{\text{Ap}}{{\text{p}}^{\rm{T}}}{{\text{A}}^{\rm{T}}}} \\ {\rm{}}& + \sum\limits_{({s_i}, {s_k}) \in {{\cal{B}}_s}} 2\omega _{{s_i}{s_k}}^2(2{\text{Bp}} - {\text{Ca}} - {{\text{D}}^{\rm{T}}}{\text{a}})\\ {\rm{}}& \cdot{{(2{\text{Bp}} - {\text{Ca}} - {{\text{D}}^{\rm{T}}}{\text{a}})}^{\rm{T}}} \end{align} $ $

(21)

基于上述分析，子区域定位算法步骤如表1所示。

表  1  子区域定位算法

输入：使用三边定位方法，粗略获得LU节点的初始值，提取子区域$s$中LU节点的位置${{\text{p}}_k}$，作为初始值。设$k = 0$，表示第$k$次迭代；

(1) 基于子区域目标函数$f({{\text{p}}_k})$，计算梯度向量$\nabla f({{\text{p}}_k})$和修正后的Hesse矩阵${\nabla ^2}f({{\text{p}}_k})$；

(2) 计算搜索方向${{\text{q}}_k}$：${{\text{q}}_k} = - {\nabla ^2}f{({{\text{p}}_k})^{ - 1}}\nabla f({{\text{p}}_k})$；

(3) 更新子区域$s$中LU节点位置，其中，步长设置为单位步长：${{\text{p}}_{k + 1}} = {{\text{p}}_k} + {{\text{q}}_k}$； 　(4) 判断迭代终止条件。如果$\frac{{\left| {f({{\text{p}}_{k + 1}}) - f({{\text{p}}_k})} \right|}}{{1 + \left| {f({{\text{p}}_k})} \right|}} < \delta$($\delta$是一个很小的正数，本文设置为$\delta = 1e{\rm{ - }}9$)或迭代次数$k > 100$，则终止迭代，${{\text{p}}_{k + 1}}$即为 子区域优化最终结果；否则令$k = k + 1$，返回至第(1)步；

输出：${{\text{p}}_{k + 1}}$。

下载: 导出CSV 

| 显示表格

3.3   融合求平均

在对网络进行划分时，以每个LU节点为中心，使LU节点和与其直接相连的邻居节点为一个子区域，将WSN划分为若干个重叠的子区域。如图1所示，a, b, c分别代表3个LU节点，它们的通信半径均为$r$，以a, b, c为圆心，$r$为半径，每个虚线圆包含一个子区域。可以看出，节点c存在于3个不同的子区域中。子区域优化后，会有3个估计值，分别为${c_1}$, ${c_2}$, ${c_3}$。为了降低定位误差，对于同一个LU节点，取出重叠区域中的估计值，进行融合求平均，作为本次迭代最终LU节点估计值。图1中${c_4}$表示最终估计值。融合求平均公式为

图  1  融合求平均示意图

下载: 全尺寸图片 幻灯片

表  3  本文算法和文献[1,2,9]算法在不同参数下实验结果对比

参数设置				时间					${\rm{\bar e\bar r\bar r}}$
$N$	$r$	${\rm{nf}}$		本文算法	文献[1]	文献[9]	文献[2]		本文算法	文献[1]	文献[9]	文献[2]
100	0.40	0		12.6 s	15.1 s	9.8 s	2.6 s		5.7e-4	0.0583	0.0489	6.5e-10
100	0.40	0.20		32.0 s	180.6 s	27.8 s	3.9 s		0.0114	0.0942	0.0746	0.0107
500	0.15	0		71.1 s	132.1 s	393.2 s	92.7 s		0.0011	0.0124	0.0025	4.0e-4
500	0.15	0.20		176.9 s	127.7 s	986.3 s	91.8 s		0.0048	0.0325	0.0164	0.0072
1000	0.11	0		309.2 s	348.7 s	51.0 min	*		0.0015	0.0097	0.0026	*
1000	0.11	0.20		522.5 s	339.3 s	2.4 h	*		0.0043	0.0249	0.0121	*
2000	0.07	0		369.1 s	760.5 s	–	*		0.0022	0.0209	–	*
2000	0.07	0.20		576.2 s	649.5 s	–	*		0.0055	0.0241	–	*
4000	0.05	0		555.9 s	23.1 min	–	*		0.0029	0.0266	–	*
4000	0.05	0.20		840.3 s	27.5 min	–	*		0.0062	0.0238	–	*
10000	0.03	0		23.3 min	1.1 h	–	*		0.0034	0.0172	–	*
10000	0.03	0.20		31.6 min	1.2 h	–	*		0.0064	0.0147	–	*

下载: 导出CSV 

| 显示表格

$$ {\bar{\text{v}}_i} = \frac{1}{{{L_i}}}\sum\limits_{s \in {M_i}} {{{\text{p}}_{s, i}}} {\rm{ , }}i = 1, 2, ·\!·\!· , m;{{ s = 1}}, {\rm{2}}, ·\!·\!· , {L_i} $ $

(22)

其中，${L_i}$表示包含LU节点$i$的子区域个数；${M_i}$表示包含LU节点$i$的所有子区域的集合；${{\text{p}}_{s, i}}$表示子区域$s$中LU节点$i$的估计值；${\bar {\text{v}}_i}$表示局部融合后LU节点$i$的估计值。这样做的优点是：每个子区域对LU节点的估计不需要太准确，即使某些子区域中的LU节点估计位置较差，也可以通过融合求平均的方式，对其进行改善，从而降低定位误差。

基于上述分析，基于改进牛顿法的分布式定位算法总结如表2所示。

表  2  基于改进牛顿法的分布式定位算法

准备工作：将$N$个传感器随机部署在${[ - 0.5, 0.5]^2}$的单位正方形区域内，区域中有$m$个LU节点，$n$个LA节点，节点的通信半径均为$r$。根据 LU节点位置及与其直接相连的邻居节点，将WSN划分为$m$个重叠子区域，$m$个子区域可以进行分布式计算。设$t = 0$；

输入：利用三边定位方法获得WSN中LU节点的粗略估计值${\tilde{\text{v}}^{(t)}}$，作为基于单位步长改进牛顿法的初始值；

(1) 利用表1中子区域定位算法，对$m$个子区域独立优化求解。得到每个子区域$s(s = 1, 2, \cdots , m)$中所有LU节点的估计位置${\text{p}}_s^{(t + 1)}$；

(2) 收集$m$个子区域中LU节点估计位置${\text{p}}_s^{(t + 1)}(s = 1, 2, \cdots , m)$，对重叠子区域中的LU节点进行融合求平均，作为LU节点的估计位 　　　置${\tilde{\text{v}}^{(t + 1)}}$； 　　　　(3) 判断迭代终止条件。如果${\rm{max}}{\left\| {\tilde {\text{v}}_i^{(t + 1)} - \tilde {\text{v}}_i^{(t)}} \right\|_2} \le \eta$($\eta$是一个小的正数，本文设置为$\eta = 1{\rm e}{\rm{ - } }2$), $i = 1, 2, ·\!·\!· , m$，则终止迭代，${\tilde{\text{v}}^{(t + 1)}}$ 　　　即为最终估计出的LU节点位置；否则将${\tilde{\text{v}}^{(t + 1)}}$作为基于单位步长的改进牛顿法的初始值，设$t = t + 1$并返回至(1)继续迭代；

输出：${\tilde {\text{v}}^{(t + 1)}}$。

下载: 导出CSV 

| 显示表格

3.4   计算复杂度分析

本文采用基于单位步长的改进牛顿法对子区域进行优化求解，求解搜索方向${\nabla ^2}f{({{\text{p}}_k})^{ - 1}}\nabla f({{\text{p}}_k})$的计算复杂度为$O(N_i^3)$，其中${N_i}$表示子区域$i$中的节点个数。由于本文算法是迭代算法且每个子区域中节点个数可能不同，设$d = \mathop {\max }\limits_i \{ {N_i}\}$。因此，本文算法的计算复杂度为$O({d^{\ 3}}mK)$，其中$d \ll m$, $m$表示LU节点个数，$K$表示算法收敛的迭代次数；文献[1]算法的计算复杂度为$O(m{n_{\rm c}}K)$，其中${n_{\rm c}}$表示网络连接度^[1]；文献[2]算法将$2 \times m$的变量转化为$(m + 2) \times (m + 2)$的变量，很大程度上增加了计算复杂度，具体的计算复杂度与使用的求解算法有关。

4.   仿真结果及分析

仿真实验中传感器节点随机分布在${[ - 0.5, 0.5]^2}$的单位正方形区域内。为了表明算法的鲁棒性，对于每次改变实验参数，都随机运行5次实验。为了观察算法整体的定位精度，本文采用与文献[1]相同的评价指标，平均定位误差(${\rm{\overline {err}}}$)。

$$ {\rm{{\overline{err}}}} = \frac{{\sum\limits_{i = 1}^m {{{\left\| {{{\text{v}}_i} - {{\tilde{\text{v}}}_i}} \right\|}_2}} }}{m} $ $

(23)

其中，$m$表示LU节点的数量，${{\text{v}}_i}$表示LU节点$i$的真实位置，${\tilde {\text{v}}_i}$表示LU节点$i$的估计位置。本文仿真在3.6 GHz的Intel i7-7700处理器和8GB RAM的PC上运行，使用软件MATLAB R2016a编写代码。

实验1　本次实验展示了锚节点数目对定位精度的影响。将500个传感器随机部署在${[ - 0.5, 0.5]^2}$的单位正方形区域内，$r = 0.15$, ${\rm{nf}} = 0.20$。假设LA节点占总节点总数的比例为$p$, $p$取值范围是0.10～0.20。 图2展示了本文算法和文献[1]、文献[9]两种分布式算法对比结果。可以看出，随着锚节点数目增加，各个算法的定位精度都会得到提高。但本文算法的性能${\rm{\overline {err}}}$在不同的锚节点比例下，明显优于另外两种算法。因此，本文算法可以用更少的锚节点数目，达到相似的定位效果，从而降低部署成本。

图  2  不同$p$情况下的${\rm{\overline {err}}}$

下载: 全尺寸图片 幻灯片

实验2　本次实验展示了通信半径$r$对定位精度的影响。由于$r$与定位精度和传感器能耗有关。因此，选择合理的$r$很重要。将500个传感器随机部署在${[ - 0.5, 0.5]^2}$的单位正方形区域内，$p \!=\! 0.15$, ${\rm{nf}} \!=\! 0.20$。本次实验$r$取值范围是0.13～0.17。由于$r$不同，子区域中的节点数目会发生变化，则定位精度也会发生变化。图3展示了本文算法和文献[1]、文献[9]两种分布式算法对比结果。可以看出，本文算法随着$r$的增加，定位精度得到提高，而且性能明显优于文献[1]和文献[9]的算法。因此，本文所提算法可以使用更短的通信半径，达到相似的定位效果，从而降低能耗，增加传感器的使用寿命。

图  3  不同$r$情况下的${\rm{\overline{err}}}$

下载: 全尺寸图片 幻灯片

实验3　本次实验与文献[1,2,9]在不同的实验参数下，进行定位时间和定位精度对比，$p = 0.15$。其中，文献[2]是集中式算法，文献[1]和文献[9]为分布式定位算法。由表3可以看出，在小规模网络($N \le 500$)中，当${\rm{nf}} = 0$时，文献[2]提出的集中式算法有很短的运行时间和很高的定位精度。但在大规模网络($N \ge 1000$)中，集中式算法便无法定位。而分布式算法能够用于大规模网络，具有良好的扩展性。在大规模网络中($N \ge 1000$)，文献[9]提出的算法运行时间很长，而本文算法和文献[1]算法运行时间明显低于文献[9]。当$N \ge 4000$时，本文算法的运行时间明显低于文献[1]所需时间。而且在不同的仿真参数下，本文算法的定位精度明显优于另外两种分布式算法。表3中，每个数据都是随机运行5次取平均所得。“—”表示由于所需时间过长，没有记录数据；“*”表示内存溢出，无法定位。综上所述，本文提出的算法能用于大规模网络，且能够在不同的参数设置下，使用较短的时间，更精确的定位。体现了本文算法定位精度，运行时间和扩展性的优势。

5.   结束语

本文针对大规模传感器网络中节点难以定位问题，提出了一种基于改进牛顿法的分布式定位算法。首先，将传感器区域划分为若干个重叠子区域，将每个子区域定位问题归结为一个无约束的优化问题。然后使用分布式算法估计子区域中LU节点的位置并进行局部融合。通过仿真实验表明，本文算法可以用于大规模传感器网络，有良好的扩展性。与现有分布式算法相比，在大规模网络中定位精度更高。但本文算法在局部融合之前需要每个子区域都优化完成。在实际情况中，子区域优化所需时间不同，会导致时间的浪费。后续工作将考虑异步分布式算法，在保证定位精度的同时，降低运行时间。