亲和传播聚类型电磁地图全覆盖重构算法

林栋明,王红军

(国防科技大学 电子对抗学院,合肥 230037)

1 引 言

电磁地图是一种可反映复杂地理环境中电磁信号的空域、时域和频域信息,进而为电磁频谱资源管理和电磁态势掌控提供有力支撑的图形或图像[1,2].其内涵覆盖无线电、雷达、卫星等电磁信号辐射源.

以往,该方面的研究基本以电磁地图中的无线电环境地图为主[3].通常在构建无线电环境地图之前,必须要先获得每个地理位置的无线电信息.传统方法主要以人力路测为主,但是,这种做法需要消耗大量人力物力,尤其是在无法到达的地区更是没有可实现性.随着分布式无线网络技术的发展,利用小型无人平台搭载感知设备实现目标区域内无线电信号探测的方式成为了数据获得的首选方式[4].然而,感知节点所得到的数据中,只有感知节点自身所处位置的接收信号功率值是准确的,而感知范围内其他位置的数据均是来源于传播方程的推算,如果对目标区域内环境的先验信息掌握较少,就无法准确选择传播模型以及对应的参数,那么这些数据便会存在较大的误差,导致无法使用.认为感知节点可精确获得探测范围内的数据,并将其用于无线电环境地图的构建这一做法,仅适合于先验信息充足的应用场景,但在实际应用时,先验信息往往难以获得,因此必须要研究先验信息缺少的情况下无线电环境地图的构建方法.

近年来,地质学上的插值算法如反距离加权插值法[5]、克里金插值法[6]等被用于无线电环境地图的构建.Denkovski等[5]利用经典的反距离加权插值法和改进的谢别德插值法(二者均属于反距离加权插值法)实现了相应无线电环境地图的构建,所用方法有效利用样本数据且具有较强的鲁棒性.但是,当感知节点与辐射源相距较近时,插值误差明显增加.Suchański等[6]着眼于无线电环境地图的构建问题,并兼顾了传感器的战术运用,用克里金插值法进行了无线电环境地图的构建,实验结果较为平滑,准确度较高;此外,该文还将最近邻法应用于地图的构建.总体来说,克里金插值法具有所得插值精度高,数据平滑性好的特点,但是对异常值敏感,少量虚假信息就会造成所得数据精度急剧下降[7],而最近邻法虽然对于无数值区域的填充效果较好,但是存在数值变化不连续的问题.另外,局部多项式法[8]、最小曲率法[9]也是常用的插值算法,也可以应用于无线电环境地图的构建.利用插值法构建无线电环境地图的方法又被称作直接法,此外,还有间接法以及混合法.间接法是通过构建电波传播模型实现无线电环境地图构建的方法[10].Gajewski[10]讨论了传播模型的作用以及构建的无线电环境地图的精度,但是这类方法往往需要目标地区的信道特点,除了自由空间的传播损耗外,还需考虑阴影衰落、电波衍射损耗等参数,以及选取合适的经验模型,参数的确定和模型的选取难度较大.另外,基于辐射源位置的无线电环境地图构建方法也可归为此类[11].混合法则是结合了传播模型或辐射源位置和插值法从而构建无线电环境地图的方法[12,13],由于该方法结合了直接和间接两种方法,所以精度高,但是方法复杂度也相应增加了.

凸集投影最初是用于图片超分辨率重建的,现在也被引入用于地震数据的重建.闫浩飞等[14]为了满足反演等工作需求,利用凸集投影算法通过二维傅里叶变换和反变换实现空缺数据的预测重建.该方法不需要和实际坐标保持网格一一对应,可直接进行数据迭代重建,能有效实现航磁数据的测量.该方法复杂度低,精度较高,但是针对不规则缺失数据的重建时,适用性不强.凸集投影在地震波数据的重建上,其基本思路都是将数据的存在域进行映射,在另一个领域对数据进行处理[15].

已有的电磁地图构建方法中,克里金插值法复杂度较高且所得数据精度受异常值影响大,反距离加权插值法受孤立点影响较大,所得结果精度低,在大面积数据重建时,结果精度不高.最近邻法、局部多项式法等其他算法的适用场景局限性较大.并且区域内的信号传播环境常常不可知,所以无法将感知节点感知范围内所得的数据均视为准确数据.为实现先验信息缺少情况下电磁地图的高精度重构,并尽可能减少所需要的感知节点数量,考虑到克里金插值法的高准确度[6]、参考将重磁数据映射至另一存在域进行处理的做法[14,15],本文提出了一种基于亲和传播聚类的电磁地图全覆盖重构算法,该算法由3部分组成:亲和传播聚类算法,(Affinity Propagation Clustering Algorithm,APCA);克里金插值法(Kriging Algorithm,KGA);域映射(Domain Mapping Algorithm,DPA).由于电磁地图包含的地图内容十分丰富,限于篇幅,本文仅考虑其中的电磁信号接收强度地图的重构(以下简称电磁地图),并将参考信号接收功率(Reference Signal Receiving Power,RSRP)作为研究对象.

2 算法实现

该算法的主要框架为:依据实际情况预先对目标区域进行栅格划分,得到每个栅格的地理位置并确定电磁数据的总量,合理确定数据量既能有效避免数据的冗余,又便于实验.划分依据为:将目标区域的长宽按照一定间隔进行分段,该间隔长度取决于感知节点的感知半径,由此可得一定数目的栅格,并将每个栅格中心处的电磁数据作为该栅格电磁数据的代表.若目标区域长宽为4000米,间隔长度设为20米,则可得40000个栅格,即该区域的电磁数据总量为40000.然后在目标区域内部署一定数目的感知节点获取对应位置上的电磁数据,即样本点,节点可以是随机部署,也可以是按照一定的部署方案进行部署;之后选定一定数量的样本点对其地理位置进行亲和传播聚类得到聚类中心,但是需要注意,选择用于聚类的样本点总量不宜超过电磁数据总量的2.5%,点数过多会导致聚类中心过多从而增大计算量以及导致聚类中心冗余,不利于电磁数据的重构.如果电磁数据总量为40000,则选择用于聚类的样本点数目不应超过1000,一般情况下,用于聚类的样本点数目为400即可满足大部分场景的要求;然后利用克里金插值法重构出以聚类中心的位置为圆心、一定半径的小面积电磁地图;这里,小面积地图的重构中心必须由聚类算法来确定,如果聚类得到n个聚类中心,则表示利用克里金插值法得到的小面积电磁地图数量为n,小面积地图的重构中心不宜由人工确定,并且半径的选择不宜过大,一般为栅格长度的10～30倍,要保证克里金插值法所得数据量和采样点所得数据量之和小于电磁数据总量;之后将电磁数据进行二维傅里叶变换,将数据从时间—空间域映射至频率—波数域进行处理,最终将数据逆映射回原始存在域.运行算法时,3个过程得到的数据量均不能为零,通过将原始数据、克里金插值法所得数据、以及域映射所得数据进行融合并重构出目标区域最终的整体电磁地图,由此实现对无线通信网络的优化和战场电磁态势的掌控.

2.1 亲和传播聚类算法

利用分布式感知节点获取目标区域内无线电数据的做法是目前较为便捷有效的方法[16-18].在目标区域中随机抛洒或按照特定方式部署一定数量的感知节点,依托这些感知节点对电磁环境进行采样.然后选择一定数量的采样点作为聚类点,对其位置进行亲和传播聚类,进而获得对应聚类中心的地理位置.

亲和传播聚类算法[19]是一种通过传递信息实现聚类的算法.该算法把数据集视作节点网络,在各节点间传递代表度r(i,k)以及适应度a(i,k)两种信息,然后利用该信息不断对聚类中心作出调整,使得节点网络的相似度达到最大,并得到相应的聚类中心.

算法通过下列公式实现信息的更新:

(1)

(2)

其中,r(i,k)表示第k个样本作为第i个样本中心点时的代表程度,r(i′,k)表示第k个样本作为第i′个样本中心点时的代表程度,a(i,k)表示第i个样本将第k个样本选作中心点时的适合程度,a(i,k′)表示第i个样本将第k′个样本选作中心点时的适合程度,s(i,k)表示第k个样本和第i个样本的相似度,s(i,k′)表示第k′个样本和第i个样本的相似度,这里取欧式距离平方的负数,s(i,i)则表示将第i个样本确定为中心点的参考度,通常为相似度的中值.两个参数通过式(3)进行更新:

(3)

λ是阻尼因子,一般为0.9,用于算法收敛.

聚类开始时,a(i,k)和r(i,k)取值均为零,之后根据式(1)～式(3)对二者进行循环更新,当算法寻找到的聚类中心方案使得网络相似度达到最大时,迭代停止,输出此时的聚类中心.聚类是根据各个样本地理位置之间的相关性进行的,目的是得到使网络相关性达到最大的聚类中心,该过程除了所处的地理位置以外,不受其他因素影响.聚类处理实质上是对电磁数据的预处理,该处理避免了直接对电磁数据进行处理,从而避免了电磁数据之间相关性的损失,同时有效地将聚类样点之间的地理相关性提取出来,可得到相关性最大的样点类群.此外,由于亲和传播聚类得到的聚类中心是在所用的样点中进行选取的,这样避免了引入其他数据的干扰;另外,由于无需预先设定聚类的总类数,因此自动化程度高,可有效避免人工干预导致的不准确性,可为后期工程化应用降低难度.

2.2 克里金插值法

克里金插值法可得到较高精度的RSRP输出如下[20]:

yout(x)=wT(x)α+e(x)

(4)

其中,yout(x)表示RSRP值,x表示地理位置,w(x)=[w1(x),w2(x),…,wN(x)]T是组成回归模型的N个所选择的函数,α=[α1,α2,…,αN]T则表示回归模型的参数,e(x)作为一个随机过程,满足以下条件:

(5)

E[·]表示求均值运算,xi,xm表示两个输入数据,R(θ,xi,xm)表示相关模型,τ2表示e(x)这个随机过程的方差,e(x)可以看做所得预测值和实际值的误差分布,θ为相关模型的参数.

若有一组样本点为{s1,…,sz},其对应样值为Y=[y1,…,yz]T,此时,可以得到下式:

(6)

定义相关函数为:

(7)

式中,n为样点的维数.同时要使得下式成立:

(8)

|R|代表计算矩阵的行列式.这里,本文采用球状模型,即:

(9)

其中,dj=μj-xj.

(10)

然后,以亲和传播聚类算法所得到的聚类中心为圆心,利用式(6)～式(10)对其一定半径长度范围内未采样位置的RSRP值进行预测重构,便可得到小面积地图数据,即实现由点数据到小面数据的获取.如果感知节点的数量较少,可以考虑增大重构半径以及增加聚类样点占总采样点的比例来优化重构效果.由于亲和传播聚类充分保留了电磁数据的相关性,同时得到了一种强相关性的中心点分布方案,并利用插值精度高的克里金插值法从现有的数据中充分提取电磁数据之间的相关性,这有利于提高所得结果的精度;另外,算法可自动选择相关性强的地理位置作为小面积地图重构的中心,有效避免了人工选择以及随机选择的不准确性,提高了算法运行效率.

2.3 域映射处理

将上述克里金插值得到的RSRP预测值填补至由感知节点采样得到的残缺数据的对应位置中,并将此数据用做域映射处理的输入数据.

1)对数据进行扩边处理,然后再进行填零处理,并得到提取矩阵:

(11)

d0(x,y)表示的是在输入残缺数据中位置为(x,y)的待处理数据,NaN表示该处数据残缺.

2)利用公式(12)处理该数据:

dt(x,y)=F-1{Tk[F(d0(x,y)+S(x,y)dt-1(x,y))]}
t∈[1,N]

(12)

t表示迭代的次数,dt(x,y)表示第t次迭代后对应位置为(x,y)的电磁数据,F、F-1表示二维傅里叶正变换以及反变换,Tk表示圆形低通滤波器,k∈[1,…,κ],κ表示截止波数,N表示总迭代次数.

Tk的确定方法如下:

(13)

α,β分别表示数据点相对于中心数据点的横纵坐标值,δ表示设定的半径.步骤2)最后得到初步的结果ds(x,y).

3)对数据进行处理得到最终结果为:

(14)

broad表示扩边的长度,xori、yori分别表示域映射处理输入数据的单行横向、单列纵向的数据量,dall(x,y)表示所需要的完整的电磁地图数据.

图1为域映射处理数据的流程.

图1 域映射数据处理流程

在本算法中,电磁采样数据是被保存在一个数组中,因此,如果某个位置没有采集到电磁数据,那么该处对应的电磁数据即设为NaN(非数值).所谓扩边处理,即在该数组周围扩充一定量的数据,例如,电磁数据原先被保存在一个100×100的数组中,现在进行扩边处理,同时要求扩边长度为20个单位,那么所得的数组大小即为140×140,并且将新填充的值均设为NaN.所谓填零处理,则是将上述得到的经扩边处理后的数组中所有NaN值用零替换.

本文求取径向平均功率谱的步骤如下:

1.利用二维傅里叶变换将数据映射到波数域,并求得功率谱;

2.设置基频,以数据中心为圆心,并将基频的整数倍(也为波数)作为半径构建数个同心圆环,并计算各圆环范围内的功率谱均值;

3.将基频倍数作为横坐标,各圆环内的功率谱平均值取对数后的数值作为纵坐标即得径向平均功率谱.

这里,本文提出了1/5功率截止波数确定法,即将径向平均功率谱下降至小于或等于最大值1/5时的第一个波数作为截止波数.

滤波器处理是将数据通过一个圆形低通滤波器,通过滤波器提取电磁数据之间的相关性,从而更好地实现未知电磁数据的估计.

数据整合则是用域映射的输入数据中的非NaN数据替换掉达到截止波数时的初步输出的数据中对应位置的数据,也就是说,此时初步得到的数据并不是最终需要的数据.在实现数据整合后再对所得的整合数据进行去扩边处理,即去掉扩边处理时所增加的数据,输出所得结果.此结果就是最终的输出结果,也就是所求的完整的电磁地图数据.

3 仿真实现及结果分析

目前,许多研究无线电环境地图的学者都将接收信号强度(Received Signal Strength Indicator,RSSI)作为预测重建的对象,该指标表示无线网络中各终端所接收到的信号的强度.而参考信号接收功率RSRP表示无线通信网络终端中某个符号所承载的参考信号的功率平均值.由于参考信号接收功率衡量电磁信号强度的结果较接收信号强度更为精准,因此,本文将参考信号接收功率作为描述区域接收信号强度分布的相关参数.

不失一般性,仿真以LTE移动通信网作为目标网络,提取Brussels地图中4000米×4000米LTE移动通信网络的无线电覆盖情况作为实验数据,其原始电磁地图如图2所示.数据产生所用的损耗模型为奥村-哈达模型,且由于该区域为实际地理区域,内有山脉隆起和河流等地物的存在,因此可较大程度的反映真实的电磁传播情况,将该目标区域进行栅格化分,可得到40000个栅格,即该区域内共40000个电磁数据.奥村-哈达模型如式(15)所示:

图2 完整电磁地图

Lb=69.55+26.16lgf-13.82lgHb-α(hm)+(44.9-6.55lgHb)lgd

(15)

其中,Lb表示传播损耗,f=2110MHz,表示辐射源频率,α(hm)为修正因子,接收端的等效高度设为hm=1.5米,此时可忽略该因子,Hb=30米,代表辐射源的等效天线高度,d代表接收端和辐射源之间的距离.那么某点的接收功率Te可表示为:

Te=Tr-Lb

(16)

Tr为辐射源辐射功率.

随机在该区域抛洒4000个低成本的感知节点,即采样点(感知节点)占比为10%,该节点均可准确采集到数据.实验平台为core i7,所用软件为MATLAB R2018b,Surfer 14和Atoll 2.8.0.

在所投放的分布式节点中随机选取400个节点,并对其地理位置进行亲和传播聚类,得到地理相关性最强的一组聚类结果如图3所示,可以看出一共得到16个聚类中心,各类中所有的非中心点均同中心点相连.

图3 聚类结果

获得聚类中心后,以各个聚类中心为圆心,半径为200米,重构出小面积的电磁地图,即利用所有分布式节点采集的电磁数据重构出16个小面积电磁地图,完成由点及面的数据获取,所得结果如图4中所示的圆形区域.

图4 小面积地图构建

将克里金插值所得的数据填补至原始采样数据的对应空缺位置后,再对该数据进行域映射处理,即:先对电磁数据进行扩边处理,扩边长度为20个单位长度,然后将其中的非数值用零来填补,之后将数据由时间—空间域映射至频率—波数域进行循环处理,直至循环结束,对数据进行整合和去扩边,便得到最终的结果.

域映射完成由点和小面到大面的数据获取,如图5除小面积圆和分散点以外的部分所示.

图5 域映射处理结果

在验证所提算法的作用效果时,本文从构建的电磁地图质量、均方根误差(Root Mean Square Error,RMSE)、决定系数以及鲁棒性等方面对算法的作用效果进行评价.

(17)

(18)

本文将克里金插值法、反距离加权插值法(Inverse Distance Weighting,IDW)、改进谢别德插值法(Modified Shepard′s Method,MSM)、最近邻法(Nearest Neighbor,NN)和局部多项式法(Local Polynomial,LP)应用于电磁地图的重构,并与本文算法进行了比较,得到图6的仿真结果.其中图6(a)为克里金插值法重构结果、图6(b)为反距离加权插值法重构结果、图6(c)为改进谢别德插值法重构结果、图6(d)为最近邻法重构结果、图6(e)为局部多项式法重构结果,图6(f)为本文算法重构结果.由于本文研究的是缺少先验信息情况下电磁地图的构建问题,而基于传播模型的构建方法作用效果较差,故这里不再赘述比较.

图6 重构结果对比

克里金插值法所得的结果在辐射源所在区域的构建效果较差,等磁线的形态十分杂乱,异常值较多,无法清楚的显示该区域的接收功率分布情况,地图质量低;从反距离加权插值法所得的结果来看,整个地图的等磁线均非常杂乱,异常值很多,牛眼现象严重,地图质量很低;改进谢别德插值法得到的结果较为清晰、等磁线形态较好,但是该结果存在牛眼现象,且存在部分异常值,地图质量较好;最近邻法得到的地图中等磁线较为杂乱,且等磁线存在波动的问题,数值变化也存在不连续的问题,该地图的质量很低;局部多项式法得到的结果虽然等磁线形态良好且无牛眼现象,但是异常值成片存在,辐射源所处区域的电磁分布情况无法准确展示出来,地图的质量很低.而本文方法所得结果中等磁线形态较好,尤其是辐射源所在区域,可较好地描述该区域接收功率的变化,且异常值较少,不存在牛眼现象,构建的地图质量较好.

除了地图的构建质量外,所得结果的精度同样是衡量算法作用效果的关键.表1为6种算法在10%采样点占比情况下的均方根误差以及决定系数.

表1 数据精度衡量1

从表1中可以看出,本文所提算法具有最小的均方根误差和最大的决定系数,因此该方法所得结果和真实值的平均偏离程度最小,且预测重构结果的数据分布同原始数据的分布最接近,是6种方法中精度最高的方法.

从理论上来说,可用于采样的感知节点越多,即感知节点的数量占电磁数据总量的比例越高,地图构建的精度也会越高.但是,如果因为成本等原因导致可用节点数量比预定数量少,所用算法必须要能保证地图的构建精度较高且和预定情况相比,其变化幅度要足够小,即算法要有较好的鲁棒性.

图7和图8分别显示当采样点占比从50%降低至10%以及采样点占比从5%降低至1%时各算法所得结果的均方根误差值.

图7 鲁棒性分析1

图8 鲁棒性分析2

从图7中可以看出,随着采样点占比下降,论文所提算法的均方根误差总增加幅度最小,为0.3274,并且,除了采样点占比为30%的情况以外,均方根误差均为最小,而当采样点占比为30%时,其均方根误差值也和最小的均方根误差值相差无几;6种算法中,变化幅度最大的是反距离加权插值法,为0.7145.而从图8中可以看出,本文所提算法的均方根误差变化幅度为0.4851,是6种算法中变化幅度第2小的算法,但是该值和最小值0.4244相差不大,且本文所提算法的均方根误差小于其他任何算法.因此,本文所提的算法具有较好的鲁棒性,作用效果稳定.

此外,从图7和图8中可以看出,要使得所得结果均方根误差小于1,最近邻法至少需要20%的采样点占比,而改进谢别德插值法需要30%,局部多项式法需要40%,克里金插值法需要10%,而反距离加权插值法在50%的采样点占比时也无法使得结果均方根误差小于1,而本文算法仅需要5%的采样点占比,因此,实际应用中,在保证同等作用效果的情况下,本文算法所需要的节点数目少.

当采样点数量较少时,算法考虑扩大重构半径以及增加聚类样点占采样点的比例来优化作用效果.经过大量实验结果的验证可以发现,当采样点数目占电磁数据总量的5%及以上时,重构半径约为划分栅格长度的10倍时可有较好的作用效果;当采样点占电磁数据总量的5%以下时,重构半径约为栅格长度的25倍时可有较好的作用效果.半径选择过大或过小,会造成所得地图质量恶化,数据精度下降.聚类样点数目一般不超过电磁数据总量的2.5%.

下面,本文针对采样点占比为1%的情况从数据精度方面进行研究.

表2显示了在采样点占比为1%时的均方根误差和决定系数.

表2 数据精度衡量2

从表2可知,本文所提出的算法有最小的均方根误差值以及最大的决定系数,因此该方法所得结果和原始值的平均偏离程度最小,且所得结果的数据分布更接近真实情况,即本算法所得结果要优于其他算法所得结果.此外,在此种情况下,本文算法仍可确保均方根误差小于1.5,而要达到这种精度,局部多项式法和反距离加权插值法等算法均需要更多的采样节点.

综上所述,本文提出的基于亲和传播聚类的电磁地图全覆盖重构算法的作用效果优于现有的克里金插值法、反距离加权插值法、改进谢别德插值法、最近邻法以及局部多项式法.

4 结 论

1)在对目标区域内信号传播环境情况未知的情况下,依托分布式节点网络获得采样数据,利用本文所提出的基于亲和传播聚类的电磁地图全覆盖重构算法可有效实现区域电磁地图的构建,且所需感知节点数目少.

2)本文所提的算法具有高精度、强鲁棒性的作用效果,作用效果优于现有的克里金插值法、反距离加权插值法等多种算法.

3)后期可以进一步减少所需的样本数据占总数据的比例,并改进本文所提算法以优化作用效果.

本文所提3种算法无法用同类型的算法进行替代.利用感知节点进行电磁信号采样后,需要对少部分数据进行处理以便下一步的数据利用.此时选择亲和传播聚类,并得到相关性高的一种地理位置聚类结果,这样既避免人工设置参数导致的结果不确定性,又避免了直接对数据进行操作,且该方法无需预先规定聚类类数,中心点为真实点,适合于先验信息较少的应用场景.立足于这些聚类中心,通过精度较高的克里金插值法重构出以聚类中心为圆心的一个个小面积电磁地图,即实现由点及面的数据获取,有效提高了可用的数据量并保证了一定的数据精度.由于傅里叶变换蕴含着信号的时频关系,因此,本文通过域映射利用数据间的相关性重构出剩余的未知数据,即实现由点、面至大面的数据获取.尽管3个阶段都有同类的方法,但是,任何一个阶段方法的改变都会造成负面影响,如准确率降低、数据的重建需要充足的先验知识等.