小生境粒子群算法在多信号DOA估计中的应用*

朱 任，庄铭杰

(华侨大学 工学院，福建 泉州 362021)

0 引 言

空间谱估计是智能天线中阵列信号处理的一个重要研究方向，在现代通信、声呐、雷达和导航等领域中被广泛应用。波达方向(Direction of Arrival,DOA)是空间谱估计中的一项重要参数，快速求解高精度DOA对智能天线系统的应用具有重要意义。多重信号分类(Multiple Signal Classification,MUSIC)算法[1]具有高精度和高分辨的特点，是估计DOA参数的经典算法之一。但是MUSIC算法需要对空间谱函数进行网格搜索来寻找极大值点，即谱峰搜索。谱峰搜索通常是一个计算量大、耗时长的穷举最优化问题，在高精度搜索时会大幅降低算法的实时性，因此保持MUSIC算法的谱峰搜索高精度同时降低计算量是一项很有价值的研究。

近年来，以遗传算法(Genetic Algorithm,GA)、粒子群算法(Particle Swarm Optimization,PSO)等为代表的经典进化算法在许多领域的问题优化、求解方面表现出优异的性能[2]，有许多研究者将其改进后应用在谱峰搜索中：文献[3-6]分别将鸡群算法、GA、PSO和差分进化算法应用到DOA估计中，经过仿真证实了这类算法能有效降低谱峰搜索的计算量，大幅提升DOA估计的精度和实时性，在实际工程应用中具有潜在的应用价值。

但是，许多情况下空间谱函数并不只有一个“峰”。在多用户环境下，阵列同时接收到多个入射信号，此时空间谱函数的图形会出现多个与DOA相对应的极大值点，因此这是一个多峰函数优化问题。显然，经典进化算法求解此类多峰值优化问题却无能为力，它们只能搜索一个最优的极大值点，而无法找到其他的极大值点。对于这类问题，引入小生境机制是一种有效的解决办法，它模拟了自然环境中不同种群占据不同类型资源的生存方式，每个小生境都代表着解空间中包含一个局部解的区域以及该区域内的粒子。这些粒子只会在所属小生境内寻找最优解，这样有利于保持种群的多样性来防止收敛于全局最优解，非常适合用来解决多峰函数优化问题。文献[7-9]在几种进化算法的基础上，分别提出了基于并行免疫、排挤、隔离策略的小生境技术，并通过仿真验证了算法的有效性。文献[10]分析了PSO的拓扑结构并提出多种不同邻域结构的小生境PSO。适应值与欧氏距离比PSO(Fitness Euclidean Distance Ratio PSO,FER-PSO)[11]利用适应值与欧氏距离的比值来调整粒子位置，实验结果表明算法能可靠地定位多个局部最优解。邻域拥挤差分进化算法(Neighborhood based Crowding DE,NCDE)[12]是一种邻域变异策略的差分进化算法，变异是在每个基于欧氏距离的邻域内进行的，因而能保持进化过程中发现的最优解，实验证明能很好解决多峰搜索问题。文献[13]模拟了蒲公英种子的传播方式，提出了小生境蒲公英算法并应用于DOA估计中。虽然实验结果显示有效，但是该算法的运行前提是设置合适的小生境半径参数，而这个参数的选取则依赖于待求解问题的先验知识，其中包括解的分布间隔、疏密程度等。在没有相关先验知识的情况下，算法很难顺利求解问题，这也是许多进化计算算法普遍存在的缺点。在诸如雷达跟踪、声呐探测等许多工程应用中，尽管适用于多峰函数求解的小生境优化算法很多，但仍难以满足MUSIC算法要求的实时处理能力和多信号DOA估计能力，同时还要尽可能减少对先验知识的依赖等。总之，目前的小生境算法可以解决多峰函数搜索的优化问题，但难以完全满足应用中对DOA估计的全部要求。

1 基于MUSIC算法的DOA估计

1.1 均匀直线阵列和L形阵列

如图1示，yoz平面的y轴上放置M根天线，按均匀直线阵列(Uniform Linear Array,ULA)方式排列。

图1 阵列模型

令相邻两阵元的间距d，达波方向与z轴的夹角为α，S(t)表示t时刻入射到阵列上的带限远场信号，则该达波信号到达第m个阵元时相对第1个阵元所产生的时延大小如式(1)所示，其中c为光速。

(1)

由时延易得阵元间的相位差。因此对于ULA，达波信号的导向向量为

(2)

ULA的阵列流形由关于各个入射信号的导向向量组成，当有k个入射信号时，ULA的阵列流形为

A=[a(α1),a(α2),…,a(αk)]，

(3)

则t时刻经ULA接收处理后信号X(t)的矩阵表达式为

X(t)=AS(t)+N(t)。

(4)

式中：N(t)是均值为0、方差为σ2的高斯白噪声，与入射信号不相关。阵列接收数据的协方差矩阵Rx可改写为

Rx=E[XXH]=AE[SSH]AH+σ2I=ARSAH+σ2I。

(5)

式中：RS是信号的相关矩阵，I是单位矩阵。由于信号与噪声相互独立，Rx可经特征值分解为信号和噪声两部分：

(6)

其中，US和UN分别是由较大特征值和较小特征值所对应的特征向量组成的信号子空间和噪声子空间。由于这两个子空间是正交的，因此导向向量aH(α)和噪声子空间UN同样也是正交的，即

aH(α)UN=0。

(7)

由于在实际应用中采集到的数据长度是有限的，无法直接得到接收数据的协方差矩阵，所以通过L次有限快拍样本来估算：

(8)

通过样本估算的协方差矩阵无法满足信号子空间与噪声子空间完全正交，即式(7)为趋近于0的数值。应用其倒数为极大值的原理，MUSIC算法以最小优化搜索的方式建立式(9)的功率谱函数来估计达波方向：

(9)

通过搜索式(9)函数的极大值点，就可以估算信号的入射角度α。

ULA无法同时检测三维空间中的俯仰角和方位角，因此本文采用如图1所示的L形阵列结构。令x轴和y轴各有M个间距均为d的均匀分布阵元，达波信号S与z轴夹角为α，x轴与信号在xoy平面上投影夹角为β，则L形阵列的导向向量为

(10)

当有k个入射信号时L形阵列的流形为

(11)

式中：x轴上的阵列流形Ax和y轴上的阵列流形Ay分别为

Ax=[ax(α1,β1),ax(α2,β2),…,ax(αk,βk)]，

(12)

Ay=[ay(α1,β1),ay(α2,β2),…,ay(αk,βk)]。

(13)

所以L型阵列的功率谱函数为

(14)

同样地，通过搜索式(14)的极大值点，就可以估算出达波信号的俯仰角α和方位角β。

1.2 谱峰搜索及复杂度

假设有一图1中yoz平面上的ULA，阵元数M=6，阵元间距d=λ/2；有三个用户信号的DOA分别为α=-60°、α=20°和α=40°，SNR=20 dB，快拍数为500。将式(9)在α∈[-90°,90°]范围内绘制得到图2(a)的空间谱图形。另假设有一图1中xoy平面上的L形阵列，在x轴和y轴各有阵元数M=6，阵元间距d=λ/2;有两个用户信号分别从α=30°、β=60°和α=60°、β=120°入射，按照式(14)的谱函数在α∈[0°,90°]、β∈[0°,180°]范围内绘制得出图2(b)的空间谱图。其中，所有用户接收信号的SNR=20 dB，快拍数均为500。图2中，空间谱图形出现了与信号源数量相同、高度相近的极大值点(即谱峰)，其对应的坐标角度就是用户信号的DOA，因此MUSIC算法的最后一个步骤就是使用网格搜索的方式寻找空间谱函数的峰值。网格的大小将决定DOA估计的复杂度和精度：一元和二元网格搜索的时间复杂度分别为O(n)和O(mn)，其中m和n是α和β轴的网格数量。网格搜索的误差小于网格的步长，例如步长为1°的网格搜索得出的结果误差在1°以内。

(a)一元MUSIC算法空间谱图形

为了展示网格搜索的性能，仿真实验选择图2的谱函数图形，仿真环境:CPU为Ryzen 5 3600，平台为Matlab2018。表1和表2给出了网格搜索的结果和用时，结果显示搜索用时随着网格步长缩小而大幅增长，导致MUSIC算法不适用于实时性要求高的场合中，而若选择大步长则使得结果得误差较大。因此，改进和优化谱峰搜索，使其在保证高精度的前提下降低计算量，是本文的主要研究内容。

表1 一元等步长网格搜索性能

表2 二元等步长网格搜索性能

1.3 经典进化算法在谱峰搜索中的问题

以遗传算法、粒子群算法为代表的经典进化算法在单信号DOA估计中表现出优异的性能，但却难以估计多信号的DOA。以图2(a)中的谱函数为例，使用经典PSO进行优化，算法参数选取：C1=C2=2.05，iwt=0.729 8，粒子数50。搜索结果如图3所示，其中黑色圆圈表示粒子，红色三角表示最优粒子。

图3 经典PSO优化图2(a)的效果图

从图3可以看出，经典PSO仅能找出α=40°处的极大值点，而却无法发现α=-60°和α=20°两处的极大值点。所有粒子都趋向于最优粒子，降低了粒子的多样性，从而失去发现其他极大值点的机会。为了降低最优粒子对全局粒子的吸引，本文引入小生境思想将解空间划分成多个区域，避免PSO在多峰的谱峰搜索中出现上述问题。

2 一种改进的粒子群算法

针对网格搜索计算量大而经典进化计算无法优化多峰谱函数的情况，本文提出了小生境粒子群算法。算法主要由中点顺序聚类和基于迭代数的搜索策略构成。本节将描述算法的原理并给出伪代码。

2.1 中点顺序聚类算法

顺序聚类(Sequential Algorithmic Scheme,SAS)是一种简单有效的聚类方式，可以将对象的集合划分成由相似对象组成的数个类。因此，顺序聚类结合PSO应用时，可以根据粒子间距和函数值将随机分布的粒子划分成数个小粒子群，形成多个小生境。本文将粒子的坐标中点信息用于改进顺序聚类，提出了一种中点顺序聚类(Midpoint SAS,MSAS)算法。

在MSAS中，f(x)表示粒子x的适应值，将所有粒子按照适应值大小进行降序排列，随后从大到小遍历每个粒子。如图4(a)所示，计算粒子Pi和距离它最近的小生境nicheN中最优秀粒子headN的欧氏距离中点midi，若midi的适应值小于Pi和headN，则认为两者之间存在一个“山谷”，Pi不属于nicheN。算法会建立一个新的小生境nichenew，并将Pi归入其中，如式(15)所示；否则，如图4(b)、(c)所示，Pi属于离它最近的小生境nicheN，算法连同中点midi一起归入nicheN，并比较这三个点的适应值，更新headN。

(a)中点适应值最小的情况

(15)

图5的黑色曲线是式(16)在0.1≤x<3.4取值范围内的图形，图中的“*”和“o”是运行MSAS将30个随机分布粒子聚类的结果，其中不同小生境的成员以不同颜色区分，“*”代表一个小生境中的最优粒子，“o”代表其他粒子。通过观察图5可知式(16)具有5个高度不同的“山峰”，且每个“山峰”的范围不同。经过MSAS对随机粒子聚类后，搜索空间被划分为5个不同大小的小生境，可见MSAS具有无需相关参数、自适应调整小生境大小的能力，对待这类大小不一、分布不均匀的多峰函数十分有效。

图5 中点顺序聚类优化式(16)的结果

(16)

中点顺序聚类算法的伪代码如下：

输入:粒子总数m,目标函数function

根据适应值大小将所有粒子降序排列。

P1∈niche1,head1←P1(第一个粒子属于第一个小生境，并记录最优值)

遍历：从i=2到m

找到Pi最近的小生境nicheN，并计算中点midi。

如果：满足式(15)：

Pi∈nichenew,headnew←Pi(Pi属于新的小生境，并记录最优值)

否则：

Pi∈nicheN,比较Pi、headN和midi适应值大小，更新headN。

2.2 粗略搜索与精细搜索

经过MSAS的处理后NPSO可以获取到大部分小生境的模糊位置和大小，并将每个已知小生境的年龄设置为1。之后，由于初始的随机粒子没能落到个别隐蔽的“山峰”中，导致这些小生境未能被建立，所以算法开始执行粗略搜索来对解空间进行充分探索。粗略搜索是发生在小生境之间的搜索，目的是模糊定位所有局部最优解。如式(17)所示，粗略搜索会遍历每个小生境，计算其最优粒子head与相邻小生境的最优粒子headN位置之间的欧氏距离，以该最优粒子为中心、距离一半作为半径，在目标函数上随机产生一个粒子，最后将产生的所有粒子通过MSAS划分到小生境中。

(17)

为了增强粗略搜索发现较隐蔽“山峰”的能力，每个小生境会在整个解空间中随机产生一个粒子。每经过一次粗搜的小生境年龄会增加1。

经过数轮的粗略搜索，NPSO已经能发现所有局部最优解的模糊位置，建立了对应数量的小生境，但是还无法精确定位每个小生境的最优粒子位置。此时，达到预设年龄参数的小生境会转为精细搜索策略：精搜是发生在小生境内部的搜索，目的是找出最优解的精确位置。如式(18)所示，随机选取小生境内一点P，以P至本小生境内最优粒子head的距离一半为半径，以head为中心，在目标函数上随机生成一个粒子，并和head比较适应值，判断是否需要更新最优粒子。

(18)

在算法初始阶段，小生境还未被完全发现或者其中的最优解还不能模糊定位时，粒子和已知的小生境头部之间的关系十分复杂。例如，当粒子和已知最近的小生境之间隔着一个未被发现的小生境时，MSAS会错误判断。因此，NPSO在每轮迭代的结尾都会判断两个相邻的小生境的最优粒子是否满足MSAS的规则，对错误的小生境进行合并。小生境粒子群算法(NPSO)伪代码如下：

输入:最大粒子数m，搜索范围，预设年龄，最大迭代数，目标函数。

在搜索范围内随机生成m个粒子。

运行MSAS对粒子聚类，划分小生境。

循环：从迭代数=1到最大迭代数。

循环：从i=1 到小生境总数。

如果：年龄<=预设年龄

根据式(17)产生粒子。

在搜索范围内随机产生粒子。

否则：

根据式(18)产生粒子。

MSAS对产生的粒子聚类，划分小生境。

判断两个相邻的小生境的最优粒子是否满足MSAS中的规则并合并误判的小生境。

保持粒子总数不超过最大粒子数m。

返回每个小生境的最优粒子。

3 仿真与分析

3.1 仿真环境与参数

假设有一图1中xoy平面的L形阵列，其中x和y轴上各有6个间距半波长的均匀分布阵元。阵列接收信号受到零均值、方差为σ2的加性高斯噪声信号干扰，信噪比为20 dB，快拍数500。仿真所采用的空间谱参数如表3所示，三个样例分别有2、3和4个谱峰，且非均匀分布，搜索范围均为α∈[0°,90°],β∈[0°,180°]。选择数量不同、分布不均的样例是为了能更好地考察算法对小生境范围的自适应能力。由于干扰噪声的影响，谱峰坐标和DOA之间存在一定误差，因此谱峰坐标由0.001°边长的网格搜索得出。本节中的误差均是指算法输出结果与谱峰坐标之间的误差。

表3 空间谱样例的参数

为了比较算法在相同计算资源下的性能表现，实验限制了目标函数(谱函数)的最大求值次数(Max Function Evaluation,MaxFEs)为10 000次。同时使用峰识别率(Peak Ratio,PR)来统计算法输出结果的准确度。式(19)是PR的定义式，由算法发现的解数量除以总的解数量得出，若输出结果和实际坐标间的欧氏距离达到精度要求，则认为算法成功获取到该解。

(19)

为了确定NPSO中的预设年龄参数，选取粒子数为200，预设年龄为1～20，并选择表3中的样例C，运行500次蒙特卡洛仿真并统计达到预设年龄时的PR，结果如表4所示。因为NPSO在达到预设年龄前的搜索策略是粗搜，更注重搜索的广度，旨在发现解的大致位置，故此处精度要求设置为5。

表4 NPSO不同预设年龄的仿真结果

表4中PR随着预设年龄增大而逐渐增大并趋近于100%。观察发现当预设年龄大于10之后PR的增量逐渐减少，呈现边际收益递减的规律。为了平衡搜索的广度和深度，预设年龄取10较为适宜。为全面体现NPSO算法是一种高效的多峰问题优化算法，表5给出了6种同类算法及参数，它们是一组性能比较优秀的算法及其参数，NPSO将在多信号DOA估计中与它们进行比较。

表5 算法参数

3.2 算法在低计算资源下的性能

仿真实验采用的谱峰搜索对象是表3中给出的样例A、B和C，粒子数均取200，每种算法各运行500次蒙特卡洛仿真。表6统计了算法的运行时间,而表7～9统计了MaxFEs为10 000时各算法在三个不同精度(0.5、0.1和0.01)要求下的PR。图6是各算法输出结果与谱峰坐标之间误差的累积分布函数(Cumulative Distribution Function,CDF)图。

表6 所有算法的运行时间

表7 精度级别ε=0.5时各算法的峰识别率

表8 精度级别ε=0.1时各算法的峰识别率

表9 精度级别ε=0.01时各算法的峰识别率

图6 各算法输出结果误差的CDF图

表6～9的结果表明，r2pso、r2psolhc、r3pso和r3psolhc四种算法因为策略相似，所以在耗时及各个精度要求下的PR都接近；NPSO以外的其他算法中，FER-PSO在0.5和0.1精度范围内下具有较高的PR，但所有算法在0.01精度条件下PR都非常低，只有NPSO算法在各个精度条件下均能保持高PR，这意味着其获取到了更精确的谱峰坐标值。图6中CDF曲线越往左移表明DOA的估计误差越小，其中NPSO的线条具有最大的斜率，表明绝大部分误差(大约95%)都分布在极小的范围内，估值结果的方差小、数值稳定；其余算法的曲线符合表7～9的数据，FER-PSO在线条起始阶段与其余算法重合，表明其低误差结果的数量与其他算法相近，之后斜率增大，意味着较低误差结果的数量会比其他算法多，而其他五种算法之间的误差辨识度不高。

通过上述实验可以看出，NPSO算法具有计算速度快的特点，相较于表2中0.01°网格算法所需的2 177 s,NPSO用时仅约其1/7 000。在保证高精度的前提下，新算法可大大提升谱峰搜索的速度。

3.3 算法在高计算资源下的性能

上节仿真实验均运行在低算力资源的环境中，且最大函数求值次数都选取10 000次，可以看出除NPSO外，在低计算资源条件下其余算法均表现出很低的峰识别率，这可能是因为算法受到算力的限制无法发挥出真实的性能。为此，本节实验中放宽算力资源的限制，使算法能对解空间充分地进行探索和挖掘，以此来检验算法的最大性能。选取表3样例C，算法参数不变，最大函数求值次数取值从1 000次开始，之后为2 500次，然后每隔2 500次取一个数值，最大为50 000次。在每个MaxFEs取值下都对7种算法运行500次蒙特卡洛仿真。图7统计了结果的均方根误差(Root Mean Square Error,RMSE)，它衡量了输出结果(观测值)与谱峰坐标(真值)之间的偏差，其解析式为

(20)

式中:N表示进行的蒙特卡洛仿真的轮数，S是入射信号总数，αn(i)表示第i轮仿真中算法输出的第n个入射信号的俯仰角α与真值之间的误差，βn(i)表示第i轮仿真中算法输出的第n个信号的方位角β与真值之间的误差。

观察图7可以看出，在较低算力资源下，NPSO较其余算法能取得更高的精度。随着MaxFEs的增大，算法获得了足够的计算资源，RMSE变化逐渐稳定，此时NPSO、FER-PSO和NCDE的RMSE相较其他算法能领先一个数量级。实验证明，无论算力资源是否充足，NPSO算法较其他算法均能保持精度的优势。

图7 RMSE随函数求值次数变化曲线

3.4 SNR对算法性能的影响

为了检验算法在不同信噪比(Signal-to-Noise Ratio,SNR)下的性能，SNR分别设置为-20 dB、-10 dB、0 dB、10 dB和20 dB，其余参数不变，运行500次蒙特卡洛仿真。图8给出了8种方法的RMSE随SNR变化曲线。当SNR较低时空间谱会出现许多假峰、干扰峰，且峰值的坐标与DOA之间会存在较大误差，因此为了准确比较结果，图中给出了网格搜索的曲线，与该线重合度越高则说明结果越准确。

图8 RMSE随SNR变化曲线

从图8中可以看出，NPSO的曲线与网格搜索的曲线高度重合，表明无论是在信噪比高还是低的情况下，本文算法不受“假”峰、“干扰”峰等的影响，均能准确地识别出谱函数的真正谱峰，有助于保持空MUSIC算法在低信噪比条件下性能维持稳定。

3.5 不同高度谱峰对算法性能的影响

在实际应用中，各信号的强度存在强弱之分，由此产生的谱峰高度就会有较大反差。为了检验算法在信号强度反差较大情况下的性能，本节仿真在样例C的基础上，PMUSIC采用功率表示，其余参数不变，得到的空间谱图形如图9所示。图中四个信号的功率强度反差较大，表现为谱峰高度大小不一。

图9 不同谱峰高度的空间谱图形

实验使用表5的算法参数，MaxFEs设置为10 000次，运行500次蒙特卡洛仿真。统计算法在三个精度级别下的PR，结果如表10所示。

表10 不同精度级别时各算法的峰识别率

仿真结果显示，在信号强度反差较大情况下，本文提出的NPSO不受谱峰之间落差大的影响，依旧能够保持极高的精度。

4 结束语

谱峰搜索耗时长、实时性低是制约MUSIC算法应用因素之一。在单谱峰中，研究者将进化计算算法应用于该领域已取得不错的效果，但是面对多谱峰尚未有十分满意的解决方案。本文提出了一种改进的小生境粒子群算法(NPSO)，并将其应用于MUSIC算法的多信号DOA估计中，通过与同类小生境算法和网格搜索进行比较，得到如下研究结论：

(1)精度高。NPSO能选择不同的搜索策略，在提高搜索广度的同时，也兼顾了搜索精度。实验结果显示，在低算力资源情况下可以取得较高精度，在算力资源丰富的情况下精度能达到10-3，满足了高精度应用场景的需求，并且不易受低信噪比以及信号强度反差大的影响。

(2)速度快。小生境的并行性使得在10-2精度级别下，NPSO仅需网格搜索用时的1/7 000，大大提高了DOA估计速度，对提升通信、雷达等的性能具有十分重要的意义。

(3)本文所提出的小生境粒子群算法NPSO适用于求解多峰函数优化问题，不仅仅是应用于多信号的MUSIC算法谱峰搜索当中，也可推广到其他具有相似问题的工程中，具有潜在的应用价值。

(4)本文仅仅提升了谱峰搜索性能，并未优化MUSIC算法的其他环节，而DOA估计的整体误差由所有环节的误差组成。因此，下一步的研究内容是将PSO用于提升MUSIC算法在接收数据处理等环节的性能。

综上所述，本文提出的方法不仅可以很好地应用于MUSIC算法的多信号DOA估计中，而且也可推广到其他需要求解多峰优化问题的工程应用中。