0 引言

信息技术、集成电路及传感器等技术的迅猛发展，使得价廉、功耗低、体积微小并具备一定计算、感知、存储及通信能力的设备得以实现，并且部署于各种物理环境之中，形成了无线传感器网络(Wireless Sensor Network, WSN)。它主要是为了对覆盖区域中的目标对象进行感知、收集和处理，并且把消息发送给监测者。监测者根据接收的数据了解监控区域情况，并根据数据信息对监控区域作相应的改变。由于传感器节点的能量非常有限，而且通常情况下能量不能及时得到补充，因此能量效率在无线传感器网络中非常重要，它影响了整个网络的寿命^[1]。为了提高能量效率，路由问题成为无线传感器网络最具有挑战性的课题。为了解决此难题，专家学者提出了多种路由协议。

低功耗自适应集簇分层型(Low Energy Adaptive Cluster Hierarchy, LEACH)路由算法^[2]是一种典型的分簇路由协议，以循环的方式随机选择簇头节点，有效降低了节点能耗，延长了网络的整体生存时间，但其簇头选举的随机性导致簇头分布不合理，且簇间单跳通信的方式导致簇头的能量消耗过快。文献[3]提出了能量有效的非均匀成簇(Energy-Efficient Uneven Clustering, EEUC)算法，通过使用非均匀的竞争半径来构造大小不等的簇，使得靠近基站的簇的规模小于远离基站的簇，均衡了簇头能耗，然而簇头的选择只由概率和门限值决定，无法保证所选簇头最优。分布式能量均衡非均匀成簇(Distributed Energy-Balanced Unequal Clustering, DEBUC)算法^[4]也是采用不同的竞争半径来实现非均匀分簇，簇头选举参考候选簇头的剩余能量和邻居节点的剩余能量，采用基于时间的簇头竞争算法，在簇间运用贪婪算法选择中继节点，均衡了能耗，但存在离基站较远的簇规模过大、节点能量消耗过多的问题。基于博弈理论的有效分簇算法(Game-theoretic Approach for Efficient Clustering, GAEC)^[5]采用博弈的方式选择簇头，在分簇的过程中还引入了分区的思想，但区域数量是随机设置的，容易导致分簇不合理，且该方法并没有采用循环方式更换簇头，会使得簇头能量消耗过快而过早死亡。翟春杰等^[6]设计了一种优化的分区算法，将节点基于区域的划分而形成簇，解决了簇的个数和分布的随机性问题。为了缓解热点效应，数据转发时隔层选择下一跳簇头，但该过程仅仅考虑距离因素，容易导致能量不均。李双双等^[7]通过非均匀分区对网络进行划分，而后在各个区内综合考虑节点能量、距离和密集程度来选举簇头，最后构造负载均衡路径树实现簇间数据的转发，有效地均衡了能量消耗，但其分区未考虑节点的分布特征，任何情况下都采用此分区模型，因此适应性较差。基于非均匀成簇的多跳路由算法(Multi-hop Routing Algorithm based on Uneven Clustering, MRAUC)^[8]建立了扇型场景下的无线传感器网络非均匀成簇模型，解决了复杂、不规则场景下的高效能组网问题，但其依然根据传统的概率阈值来进行临时簇头的选举，引入了随机性；算法以最小能量消耗原则竞选最佳中继簇头，然而忽略了簇头能耗的均衡问题。文献[9]提出一种基于随机虚拟骨干树和改进退避分布式聚类协议的能量感知路由算法，节省了网络能耗，但算法在全局范围内通过定时广播竞选簇头时，其延时时长仅由节点的剩余能量决定，这样可能导致簇头分布不合理。在簇间路由阶段，基于节点位置和密度的非均匀分簇路由算法(node Location and node Density based Unbalance Clustering routing algorithm, LDUC)^[10]利用通信簇头节点来分担簇头的簇间数据转发任务，减少了簇头的能耗，但通信簇头需要进行额外的选举过程，增加了控制能耗。文献[11]提出一种新的基于预测能量消耗效率的分簇路由算法，选择簇头时充分考虑节点度和节点间的平均距离，因此簇内通信代价很小；在考虑预测能量消耗、跳数和传播延迟的基础上，利用蜂群优化(Bee Colony Optimization, BCO)算法来优化簇间数据传输，提高整体的网络性能，降低和平衡整个网络的能量消耗，延长网络的生存时间。

针对以上一些研究中分簇灵活性和合理性较差的问题，且为均衡能耗、提高节点能量利用效率、延长网络寿命，本文提出一种基于最优簇数和改进引力搜索的WSN路由算法(WSN routing algorithm based on Optimal Number of Clusters and Improved Gravitational Search, ONCIGS)。首先，利用最优簇数对改进的凝聚嵌套(AGglomerative NESting, AGNES)算法进行初始化，实现网络的合理分簇；然后，将簇头剩余能量的标准差作为目标函数，采用改进的引力搜索算法搜索能量均衡的簇间数据转发路径。实验结果表明，ONCIGS可以减少并均衡节点能耗，从而延长网络生命周期。

1 系统模型 1.1 网络模型

本文对网络模型作如下假设：1)N个传感器节点随机分布在M×M监测区域内，基站位于监测区域之外。2)基站和所有节点静止不动，节点能量有限，基站能量无限。3)节点同构，即具有相同的属性，且具有唯一的ID。4)无线信道满足对称性，即正向传输和反向传输的能量消耗相同。5)节点可感知得到自身的位置坐标。6)节点可根据通信距离调整无线发射功率。

1.2 能量模型

本文能量模型采用经典的一阶无线电模型^[12]。将lbit数据传输d距离所消耗的能量如下：

$ {E_{{\rm{Tx}}}}\left( {l,d} \right) = \left\{ \begin{array}{l} l{E_{{\rm{elec}}}} + l{\varepsilon _{{\rm{fs}}}}{d^2},\;\;\;\;\;\;d < {d_0}\\ l{E_{{\rm{elec}}}} + l{\varepsilon _{{\rm{amp}}}}{d^4},\;\;\;\;d \ge {d_0} \end{array} \right. $

(1)

接收lbit数据所消耗的能量为：

$ {E _{{\rm{Rx}}}}{\rm{(}}l{\rm{)}} = l{E_{{\rm{elec}}}} $

(2)

其中：E_Tx表示发射节点能耗；E_Rx表示接收节点能耗；E_elec表示发送或接收每比特数据电路消耗的能量；ε_fs和ε_amp分别表示节点中放大器在自由空间以及多径衰落模型下的能耗系数。若通信距离d < d₀，采用电波传输的自由空间模型；若通信距离d≥d₀，则采用多径衰落模型。阈值${d_0} = \sqrt {{\varepsilon _{{\rm{fs}}}}-{\varepsilon _{{\rm{amp}}}}} $。

1.3 数据融合模型

本文通过数据融合技术来减少网络传输的数据量，进而降低网络能量消耗。由于不同簇中采集到的数据具有较大差异性，本文仿真不考虑簇间数据的融合，只考虑簇内的数据融合。假设每个簇内成员向簇头发送的数据包大小为lbit，那么无论簇内有多少个成员，簇头均将接收到的数据压缩为lbit。

2 ONCIGS

网络初始化时，每个节点将自身ID以及感知到的位置坐标发送给基站，以便基站掌握整个网络中节点的分布情况，为之后的运算提供条件。

2.1 非均匀分簇 2.1.1 ONCIGS最优簇数

为了合理分簇，需要先确定最优簇数。目前大多数涉及最优簇数的研究利用在均匀分簇的前提下推导出的最优簇数来指导非均匀分簇，这样有失合理性。因此，本文在非均匀分簇的条件下确定最优簇数。

设将网络区域划分为k+1个簇，簇半径在R_c/2到R_c之间等步长选取，那么第i+1个簇的半径为：

$ {R_{i + 1}} = (k + i){R_c}/\left( {2k} \right) $

(3)

第i+1个簇的面积为：

$ {S_{i + 1}} = \pi {R^2} = {(k + i)^2}{R_c}^2\pi /\left( {4{k^2}} \right) $

(4)

为了保证网络覆盖率，需要满足条件$\sum\limits_{i = 0}^k {{S_{i + 1}} \geqslant {M^2}} $，即：

$ {R_c} \geqslant \sqrt {\frac{{24k{M^2}}}{{\pi (14k + 1)(k + 1)}}} $

(5)

网络节点密度λ=N/M²，则第i+1个簇内的节点数为：

$ {N_{i + 1}} = \lambda {S_{i + 1}} = \frac{N}{{{M^2}}} \cdot \frac{{{{(k + i)}^2}{R_c}^2\pi }}{{4{k^2}}} = \frac{{{{(k + i)}^2}{R_c}^2\pi N}}{{4{k^2}{M^2}}} $

(6)

设第i+1个簇内普通节点到簇头的距离为d_toCH，则有：

$ E{\rm{[}}d_{{\rm{toCH}}}^2{\rm{]}} = \iint\limits_S {{\rm{(}}{x^2} + {y^2}{\rm{)}}\rho {\rm{(}}x, y{\rm{)d}}x{\rm{d}}y} = \frac{{{{{\rm{(}}k + i{\rm{)}}}^2}{R_c}^2}}{{8{k^2}}} $

(7)

其中ρ(x, y)为簇内普通节点的分布密度，计算式如下：

$ \rho {\rm{(}}x, y{\rm{)}} = \frac{1}{{{S_{i + 1}}}} = \frac{{4{k^2}}}{{{{{\rm{(}}k + i{\rm{)}}}^2}{R_c}^2\pi }} $

(8)

那么第i+1个簇内的成员节点平均消耗的能量为：

$ {E_{{\rm{member}}}} = l{E_{{\rm{elec}}}} + l{\varepsilon _{{\rm{fs}}}}{d_{{\rm{toCH}}}}^2 $

(9)

其簇头消耗的能量为：

$ {E_{{\rm{CH}}}} = ({N_{i + 1}}-1)l{E_{{\rm{elec}}}} + {N_{i + 1}}l{E_{{\rm{DA}}}} + l{E_{{\rm{elec}}}} + l{\varepsilon _{{\rm{amp}}}}{d_{{\rm{toBS}}}}^4 $

(10)

则第i+1个簇的总能耗为：

$ \begin{gathered} {E_{i + 1}} = {E_{{\rm{CH}}}} + ({N_{i + 1}}-1){E_{{\rm{member}}}} = l{E_{{\rm{elec}}}}(2{N_{i + 1}}-1) + \hfill \\ \;\;\;\;\;\;\;\;\;{N_{i + 1}}l{E_{{\rm{DA}}}} + l{\varepsilon _{{\rm{amp}}}}{d_{{\rm{toBS}}}}^4 + ({N_{i + 1}}-1)l{\varepsilon _{{\rm{fs}}}}{d_{{\rm{toCH}}}}^2 \hfill \\ \end{gathered} $

(11)

因此，整个网络每一轮的总能耗为：

$ \begin{gathered} {E_{{\rm{total}}}} = \sum\limits_{i = 0}^k {{E_{i + 1}}} = l{E_{{\rm{elec}}}}(2N- k- 1) + Nl{E_{{\rm{DA}}}} + \hfill \\ {\rm{ }}(k + 1)l{\varepsilon _{{\rm{amp}}}}{d_{{\rm{toBS}}}}^4 + l{\varepsilon _{{\rm{fs}}}} \times \left[{\frac{{R_c^4\pi N}}{{32{M^2}}}\left( {\frac{{55}}{{3k}} + 8-\frac{1}{{30{k^3}}}} \right)-} \right. \hfill \\ \;\;\;\left. {\frac{{R_c^2}}{8}\left( {\frac{{7k}}{3} + \frac{5}{2} + \frac{1}{{6k}}} \right)} \right] \hfill \\ \end{gathered} $

(12)

为求得最优k值，令$\frac{{{\rm{d}}{E_{{\rm{total}}}}}}{{{\rm{d}}k}} = 0 $，得：

$ {k_{{\rm{opt}}}} = \sqrt {\frac{{-b + \sqrt {{b^2}-4ac} }}{{2a}}} $

(13)

其中：

$ a =-3840{M^2}{E_{{\rm{elec}}}} + 3840{M^2}{\varepsilon _{{\rm{amp}}}}{d_{{\rm{toBS}}}}^4-1120{M^2}{\varepsilon _{{\rm{fs}}}}R_c^2 $

(14)

$ b =-280{\varepsilon _{{\rm{fs}}}}R_c^4\pi N + 80{M^2}{\varepsilon _{{\rm{fs}}}}R_c^2 $

(15)

$ c = 12{\varepsilon _{{\rm{fs}}}}R_c^4\pi N $

(16)

2.1.2 基于改进AGNES聚类的非均匀分簇

确定了最优簇数，本文采用凝聚嵌套(AGNES)算法^[13]进行网络节点的聚类。AGNES是一种采用自底向上聚合策略的层次聚类算法，它先将数据集中的每个样本看作一个初始聚类簇，然后在算法运行的每一步中找出距离最近的两个聚类簇进行合并，该过程不断重复，直至达到预设的聚类簇个数^[13]。为了缓解能量空洞，需要对AGNES进行必要的改进，实现分簇的非均匀。

定义1  相对邻近度RC(C_i, C_j)为：

$ {RC} {\rm{(}}{C_i}, {C_j}{\rm{)}} = \frac{{{{{\rm{[}}{ED} {\rm{(}}{C_i}{\rm{)]}}}^2} + {{{\rm{[}}{ED} {\rm{(}}{C_j}{\rm{)]}}}^2}}}{{2E{D_{{\rm{max}}}}^2}} $

(17)

其中：C_i与C_j分别表示簇i和簇j；ED(C_i)为C_i到基站的距离，具体表示如式(18)；ED_max为当前所有簇到基站距离的最大值。

$ \begin{gathered} {ED} {\rm{(}}{C_i}{\rm{)}} = \hfill \\ {\rm{ }}\sqrt {{{\left( {\frac{1}{{\left| {{C_i}} \right|}}\sum\limits_{k \in {C_i}} {{x_k}}-{x_{{\rm{BS}}}}} \right)}^2} + {{\left( {\frac{1}{{\left| {{C_i}} \right|}}\sum\limits_{k \in {C_i}} {{y_k}}-{y_{{\rm{BS}}}}} \right)}^2}} \hfill \\ \end{gathered} $

(18)

其中：|C_i|表示C_i内的传感器节点个数；x_k、y_k分别表示C_i内节点k的横、纵坐标；x_BS、y_BS分别表示基站的横、纵坐标。相对邻近度衡量了簇对到基站的距离，簇对到基站距离越大，相对邻近度越大。

定义2  合并开关因子TH(C_i, C_j)为：

$ {TH} {\rm{(}}{C_i}, {C_j}{\rm{)}} = \left\{ {\begin{array}{*{20}{c}} {1, }&{{{d} _{{\rm{max}}}}{\rm{(}}{C_i}, {C_j}{\rm{)}} \leqslant 2\varphi {{R} _c}{\rm{(}}{C_i}, {C_j}{\rm{)}}} \\ {0, }&{{{d} _{{\rm{max}}}}{\rm{(}}{C_i}, {C_j}{\rm{)}} > 2\varphi {{R} _c}{\rm{(}}{C_i}, {C_j}{\rm{)}}} \end{array}} \right. $

(19)

其中：d_max(C_i, C_j)表示C_i与C_j的最大距离，具体表示如式(20)；φ为限制因子；R_c(C_i, C_j)为限制簇半径，表达式如式(21)所示。

$ {{d} _{{\rm{max}}}}{\rm{(}}{C_i}, {C_j}{\rm{)}} = \mathop {{\rm{max}}}\limits_{x \in {C_i}, z \in {C_j}} {dist} {\rm{(}}x, z{\rm{)}} $

(20)

其中：dist(x, z)表示节点x与节点z的欧氏距离。

$ {{R} _c}{\rm{(}}{C_i}, {C_j}{\rm{)}} = \left( {1-c\frac{{{D_{{\rm{max}}}}-{ED} {\rm{(}}{C_i} \cup {C_j}{\rm{)}}}}{{{D_{{\rm{max}}}}-{D_{{\rm{min}}}}}}} \right)R_c^0 $

(21)

其中：R_c⁰表示网络最大限制簇半径；c用于控制R_c(C_i, C_j)的取值范围；D_max、D_min分别为网络中节点与基站的最大和最小距离。设计TH(C_i, C_j)的目的是对簇半径加以限制，防止由于极大簇的产生而造成簇头负载过重。

关于两个簇之间的距离，本文采用平均距离：

$ {{d} _{{\rm{avg}}}}{\rm{(}}{C_i}, {C_j}{\rm{)}} = \frac{1}{{|{C_i}||{C_j}|}}\sum\limits_{x \in {C_i}} {\sum\limits_{z \in {C_j}} {{dist} {\rm{(}}x, z{\rm{)}}} } $

(22)

综合上述定义，设计簇对C_i与C_j对应的权值S如下：

$ S = \frac{{{RC} {{{\rm{(}}{C_i}, {C_j}{\rm{)}}}^\alpha }}}{{{{d} _{{\rm{avg}}}}{\rm{(}}{C_i}, {C_j}{\rm{)}}}} \times {TH} {\rm{(}}{C_i}, {C_j}{\rm{)}} $

(23)

其中：α由用户指定，用于调节相对邻近度的权重。权值S综合考虑了簇对距离、簇对到基站距离以及簇对合集的几何尺寸。在聚合的过程中，簇对是否合并不再像原始算法那样仅仅由簇对之间的距离决定，而是将最大的S值所对应的两个簇合并。那么算法运行时，在簇对合集的规模不超过限定值即TH(C_i, C_j)=1的前提下，距离基站较远的相邻簇对具有较大的相对邻近度，而距离基站较近的簇对相对邻近度较小，那么如果此时的d_avg(C_i, C_j)值大小相当，则距离基站较远的相邻簇对权值S较大，即具有较大的可能性合并在一起，而距离基站较近的簇对合并的可能性较小；从而基站附近的簇几何规模较小，远离基站的簇几何规模较大，实现了非均匀分簇，基站附近的簇头将预留出更多能量来完成数据转发任务。如果某簇对合集的簇几何尺寸超出限定值，则TH(C_i, C_j)=0，该簇对的权值S被置零，不具备合并可能性。

非均匀分簇过程的流程如图 1所示。

考虑到降低能耗，本文避免在每一轮都进行重新分簇，而只是在网络初始化阶段和死亡节点比例每上升10%时进行分簇操作，簇头选举和路由选择则是每轮都需要进行。分簇完成之后，基站将分簇结果发送给网络节点。

2.2 簇头选举

为了降低算法复杂度，本文簇头的选举采用文献[14]的定时广播方法，等待时间的长短取决于节点剩余能量的大小，拥有较多剩余能量的节点有更大的机会当选为簇头。普通节点接收到其所属簇的簇头消息后，向簇头发送加入请求消息，簇头回复请求并分配时隙，而后开始数据的采集与上传。

考虑到节省能耗，如果当选簇头没有改变，依然为上一轮的簇头，则簇内节点无需发送加入请求，簇头也无需重新分配时隙，可直接按照上一轮的拓扑结构和时间表来运行。

2.3 簇的路由

为了均衡簇间能耗，簇间路由结合改进的引力搜索算法(Gravitational Search Algorithm, GSA)，寻找各簇头与基站之间的最优路径，簇头沿着最优路径将数据转发到基站。

2.3.1 改进的引力搜索算法

引力搜索算法(GSA)^[15]是一种基于牛顿万有引力定律和运动定律的启发式搜索算法。假设在一个D维的搜索空间中有N个搜索粒子，其中第i个粒子在该搜索空间中的位置为：

$ {{\boldsymbol{X}}_i} = (x_i^1, x_i^2, \cdots, x_i^k, \cdots, x_i^D);\;\;i = 1, 2, \cdots, N $

(24)

以第t次迭代的执行过程为例，对算法的实现描述如下。

GSA首先初始化各粒子的位置，并计算粒子的适应度值。适应度值决定了粒子惯性质量的大小，使用以下计算式更新粒子的惯性质量：

$ {{m} _i}{\rm{(}}t{\rm{)}} = \left[{{{f} _i}{\rm{(}}t{\rm{)}}-{w} {\rm{(}}t{\rm{)}}} \right]/\left[{{b} {\rm{(}}t{\rm{)}}-{w} {\rm{(}}t{\rm{)}}} \right] $

(25)

$ {{M} _i}{\rm{(}}t{\rm{)}} = {{m} _i}{\rm{(}}t{\rm{)/}}\sum\limits_{j = 1}^N {{{m} _j}{\rm{(}}t{\rm{)}}} $

(26)

其中：f_i(t)表示粒子i的适应度值；b(t)、w(t)分别表示种群中的最优适应度值和最差适应度值。

粒子i受到粒子j在第k维上的作用力大小为：

$ {F} _{ij}^k{\rm{(}}t{\rm{)}} = {G} {\rm{(}}t{\rm{)}}\frac{{{{M} _{pi}}{\rm{(}}t{\rm{)}} \times {{M} _{aj}}{\rm{(}}t{\rm{)}}}}{{{{R} _{ij}}{\rm{(}}t{\rm{)}} + \varepsilon }}{\rm{[}}{x} _j^k{\rm{(}}t{\rm{)}}-{x} _i^k{\rm{(}}t{\rm{)]}} $

(27)

其中：G(t)是引力常量；R_ij(t)表示粒子i和粒子j的欧氏距离；ε是一个很小的常量。

粒子i在第k维上受到的合力为：

$ F_i^k\left( t \right) = \sum\limits_{j = 1, j \ne i}^N {ran{k_j}F_{ij}^k\left( t \right)} $

(28)

其中：rank_j是0到1之间的随机数。

依据牛顿第二定律，粒子i在第k维上的加速度为：

$ a_i^k\left( t \right) = F_i^k\left( t \right)/{M_{ii}}\left( t \right) $

(29)

粒子i在第k维上的速度和位置使用式(30)更新：

$ \left\{ {\begin{array}{*{20}{l}} {v_i^k\left( {t + 1} \right) = ran{k_i} \times v_i^k\left( t \right) + a_i^k\left( t \right)}\\ {x_i^k\left( {t + 1} \right) = x_i^k\left( t \right) + v_i^k\left( {t + 1} \right)} \end{array}} \right. $

(30)

针对基本GSA易陷入局部最优的问题，本文首先引入反向学习机制^[16]，以增强种群的多样性；随后结合精英策略思想，来提高种群的质量。

设一个D维搜索空间中有N个搜索粒子，构成种群PP。每个粒子的位置表示如式(24)所示，其中，x_i^k∈[x_min^k, x_max^k]，x_min^k和x_max^k分别表示种群在第k维上位置的最小值和最大值。通过引入反向学习机制，可产生PP中粒子X_i的反向粒子X′_i，进而构成PP的反向种群RP，X′_i的位置如下：

$ {\mathit{\boldsymbol{X'}}_i} = \left( {{{x'}_{i1}}, {{x'}_{i2}}, \cdots, {{x'}_{ik}}, \cdots, {{x'}_{iD}}} \right);\;\;\;i = 1, 2, \cdots, N $

(31)

其中，x′_ik的计算式如下：

$ {x'_{ik}} = x_{\min }^k + x_{\max }^k-x_i^k $

(32)

将种群PP与反向种群RP进行合并，组成合集，并将合集中的2N个粒子按照适应值由优到劣的顺序进行排列，接下来结合精英策略思想来更新种群。设合集中前20%较优的粒子所处的位置为X_{i_best}=(x_{i_best}¹, x_{i_best}², …, x_{i_best}^k, …, x_{i_best}^D)(i=1, 2, …, 2N/5)，通过以下规则来产生新粒子：

$ F = {R_m} \times rand\left( {-0.5, 0.5} \right)/D $

(33)

$ {\mathit{\boldsymbol{X}}_{i\_new}} = {\mathit{\boldsymbol{X}}_{i\_best}} \times F $

(34)

其中：因子F描述了新粒子的产生规则；R_m表示最优粒子与离其最近的粒子之间的欧氏距离；rand(-0.5, 0.5)表示[-0.5, 0.5]范围内的随机数；D为搜索空间的维数；X_{i_new}为产生的新粒子的位置。利用X_{i_new}(i=1, 2, …, 2N/5)替代合集中后20%较差的粒子，并将合集中的粒子重新排序，最后保留前N个粒子，其余舍弃，至此，种群更新完成。

GSA迭代后期，粒子在最优解附近聚集，种群密度较大，因此粒子间相互作用力会变得非常大，导致粒子在最优解附近高速振荡而减慢收敛速度。针对这一问题，本文首先定义种群密度因子。

定义3  种群密度因子δ为：

$ \delta {\rm{ = }}\frac{1}{N}\sum\limits_{i = 1}^N {{\delta _i}} $

(35)

其中：N为粒子个数；δ_i表示粒子i和其他所有粒子的平均距离，表达式如式(36)所示。

$ {\delta _i} = \frac{1}{{N-1}}\sum\limits_{j = 1, j \ne i}^N {\sqrt {\sum\limits_{k = 1}^D {{{\left( {x_i^k-x_j^k} \right)}^2}} } } $

(36)

根据式(35)对式(27)进行改进：

$ F_{ij}^k\left( t \right) = G\left( t \right)\frac{{{M_{pi}}\left( t \right) \times {M_{aj}}\left( t \right)}}{{{R_{ij}}{{\left( t \right)}^{Ep\left( \delta \right)}} + \varepsilon }}\left[{x_j^k\left( t \right)-x_i^k\left( t \right)} \right] $

(37)

其中，Ep(δ)采用如下设计：

$ Ep\left( \delta \right) = \left\{ \begin{array}{l} E{p_{{\rm{min}}}} + \left( {E{p_{{\rm{max}}}}-E{p_{{\rm{min}}}}} \right){{\rm{e}}^{1-1/\delta }}, \;\;\;\delta < 1\\ E{p_{{\rm{max}}}}, \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\delta \ge 1 \end{array} \right. $

(38)

其中：Ep_min和Ep_max分别表示Ep(δ)的最小值和最大值。如此设计，当种群密度过大时，R_ij(t)的指数变小，那么粒子间的距离对相互作用力的影响也变小，相距很近的粒子间作用力不至于过大，算法的收敛能力得以提升。

2.3.2 基于改进引力搜索的路由

设网络中有7个簇头节点，簇头集合为{CH₁, CH₂, …, CH₇}。在[0, 1]范围内随机产生N个粒子，粒子的维数为簇头节点的个数，每个粒子代表一种多跳路径。

统计每个簇头通信范围内的前向簇头节点，构建候选下一跳簇头列表，假设如表 1所示。

簇头CH_k在其候选下一跳簇头集合中选择第n_k个簇头作为最终的下一跳，n_k的取值如下：

$ {n_k} = \left\lceil {x_i^k \times NextNum\left( k \right)} \right\rceil $

(39)

其中：x_i^k为粒子i的第k维分量，其表示簇头k在选择最终下一跳时的概率；NextNum(k)表示簇头k的候选下一跳簇头个数。假设搜索结束后最优粒子的位置为(0.72, 0.44, 0.17, 0.25, 0.86, 0.21, 0.93)，那么结合表 1和式(39)可知${n_1} = \left\lceil {0.72 \times 2} \right\rceil = 2 $，即CH₁选择CH₆作为下一跳，以此类推，CH₂~CH₇的下一跳分别为CH₅、BS、CH₃、BS、CH₃、CH₆，由此得出多跳路径。

为均衡簇头能耗、延长网络寿命，本文将簇头剩余能量的标准差作为改进GSA的适应值函数，具体表示如下：

$ \sigma = \sqrt {\frac{1}{m}\sum\limits_{i = 1}^m {{{\left[{E\left( i \right)-\bar E} \right]}^2}} } $

(40)

其中：m表示簇头节点的个数；E(i)为簇头i在每轮通信过后的剩余能量；E表示簇头剩余能量的均值。簇间通信中簇头的能耗包括接收其他簇头数据包的能耗、发送自身数据包和其他簇头数据包的能耗，所以计算E(i)时需要扣除上述能耗。E(i)的计算方式如下：

$ \begin{array}{l} E\left( i \right) = {E_{\rm{b}}}\left( i \right)-{E_{{\rm{Rx}}}}\left( {{l_{{\rm{data}}}} \times PackNum\left( i \right)} \right)-\\ \;\;\;\;\;\;\;\;\;\;{E_{{\rm{Tx}}}}\left( {{l_{{\rm{data}}}} \times \left( {PackNum\left( i \right) + 1} \right), {d_{{\rm{next}}}}\left( i \right)} \right) \end{array} $

(41)

其中：E_b表示簇间通信开始前簇头i的剩余能量；l_data为数据包大小；PackNum(i)表示簇头i接收到的数据包个数；d_next(i)表示簇头i与下一跳簇头的距离。

优化过程的详细步骤如下：

步骤1  在[0, 1]区间内随机初始化种群。

步骤2  依据粒子位置和候选下一跳簇头列表得到各粒子所代表的多跳路径，计算各粒子的适应度值，并更新b(t)和w(t)。

步骤3  计算粒子的惯性质量。

步骤4  计算粒子所受合力。

步骤5  计算粒子的加速度。

步骤6  更新粒子的速度和位置。

步骤7  如果rand(0, 1) < T，则利用反向学习机制和精英策略思想来更新种群，其中T∈[0, 1]，由用户指定。

步骤8  跳转至步骤2继续执行，直到达到最大迭代次数。

3 实验仿真与分析

为了充分验证ONCIGS的性能，本文利用Matlab仿真软件在不同的环境下对ONCIGS、LEACH和DEBUC算法进行仿真分析，网络环境分为两种：网络覆盖面积100m×100m，节点数量100个，基站位置(50, 0)m；网络覆盖面积200m×200m，节点数量400个，基站位置(100, 0)m。节点初始能量设置为0.5J，E_Tx=E_Rx=E_elec=50nJ/bit，ε_fs=10pJ/(bit·m²)，ε_amp=0.001 3pJ/(bit·m⁴)，数据融合能耗设定为E_DA=5nJ/bit，阈值d₀=87m，数据包大小4 000bit，控制包大小200bit。对于ONCIGS算法，取c=0.5、φ=1、α=0.5、Ep_min=0.5、Ep_max=1。

3.1 最优簇数

两种网络规模中簇数与每轮总能耗的关系曲线如图 2所示。从图 2中看出，随着簇数的增加，网络一轮总能耗呈现出先降低后增长的趋势，且两条曲线分别在簇数为4和9时总能耗达到最低。这是因为簇数较少时，簇的几何规模较大，则簇内节点的通信距离较大，导致了网络总能耗较大；而簇数较多会增加簇间通信的任务量，也将导致总能耗较大。因此，由图 2可知，网络规模100m×100m时的最优簇数为4个，网络规模200m×200m时的最优簇数为9个。本文采用此数值对AGNES进行初始化。

3.2 网络生命周期

在两种网络规模下三种算法的网络生命周期对比如图 3所示。本文将从网络开始运行到20%节点死亡的轮数定义为网络的生命周期。由图 3可知，网络规模为100m×100m时，LEACH、DEBUC和ONCIGS的生命周期分别为465轮、624轮和660轮，相比于LEACH和DEBUC，ONCIGS分别将生命周期延长41.94%和5.77%；网络规模为200m×200m时，三种算法的生命周期分别为359轮、588轮和634轮，分别将生命周期延长76.60%和7.82%。

以上数据表明，相对于LEACH和DEBUC算法，ONCIGS可有效延长网络生命周期。这是因为ONCIGS首先求出非均匀分簇下的最佳簇数，而后采用改进的AGNES算法进行聚类，并将上述最佳簇数作为迭代的截止条件，分簇更加合理，簇内能耗更低。

3.3 网络总能耗

随着网络的运行，三种算法网络总能耗的变化情况如图 4所示。从图 4中可以看出，从网络开始运行到网络能量全部耗尽之前，ONCIGS的总能耗始终低于LEACH和DEBUC算法，这说明ONCIGS能够有效降低网络的能量消耗。

首先，ONCIGS在基于最优簇数进行合理分簇的前提下，通过改进的GSA搜索出理想的多跳路径，依此进行数据包的转发，簇间通信能耗得以减少；其次，ONCIGS避免了其余两种算法频繁分簇的问题，本文方法只在网络初始化和死亡节点比例每上升一定百分点时才进行分簇操作，这进一步降低了能耗。

3.4 网络节点平均剩余能量

三种算法的节点平均剩余能量对比如图 5所示。从图 5中可以看出，在两种网络规模下，ONCIGS的节点平均剩余能量均明显高于其余两种算法，而较高的节点平均剩余能量能够反映出较均衡的能量消耗，因此ONCIGS有效均衡了网络能耗。

在路由阶段，ONCIGS采用改进的GSA进行多跳路径的搜索，且将簇头剩余能量的标准差作为搜索的目标函数，均衡了簇头能耗。

4 结语

针对现有一些协议分簇机制或路由机制不完善而导致网络节点能量利用效率较低的问题，本文提出一种基于最优簇数和改进引力搜索的WSN路由算法ONCIGS。首先，推导出非均匀分簇条件下的最优簇数，并通过改进的AGNES算法实现最优数目的非均匀分簇；然后，采用改进GSA搜索能量均衡度较高的簇间多跳路径。仿真结果表明，ONCIGS能够在有效降低网络能耗的同时均衡节点能耗，提高节点能量的利用效率，进而延长网络生命周期。本文虽取得一定成果，但研究是在节点位置固定的假设下进行，且对于大规模网络的路由算法研究仍显不足。因此，进一步的工作将着重于节点移动的大规模无线传感器网络路由算法的研究。