基于粒子群优化的宽带同心圆环阵型优化方法

康雅聪，魏明洋，田巳睿，丁林宁

(1.南京理工大学 电子工程与光电技术学院 江苏 南京 210094；2.中国科学院 声学研究所，北京 100190)

0 引言

声学成像技术在航空气动声学和机械设备噪声检测[1-2]等领域具有广泛的应用。它结合了摄像头获取的实时画面和平面传声器阵列得到的声源能量辐射分布，能直观地确定目标声源的位置。阵列点扩散函数(PSF)[3-4]是指声源位于画面正中心时的阵列波束响应，它指示了阵列对空间声源的响应情况，且与阵元位置有关。此外，衡量声成像性能的一个重要指标是传声器阵列的抗干扰能力，即期望阵列处理后形成的波束图具有很低的旁瓣水平。自适应波束形成、反卷积波束形成等阵列算法能有效消除波束图中的旁瓣，但是这些算法的效果受限于原始波束图的质量，即点扩散函数的旁瓣水平。

在实际工程应用中，阵列的布放范围、阵元数量、阵元间距等都有一定的要求。为了在这些实际限定的因素下得到更低的旁瓣水平，需要对阵列的阵元分布进行优化设计。阵型优化是指以阵元空间位置为自变量，点扩散函数的旁瓣水平为因变量，获取最低旁瓣水平时阵元位置的优化过程。考虑到阵型优化是一个高维度、非线性、不连续、不可微的优化问题，当前的解决方案是利用全局优化的智能算法进行求解，如遗传算法[5]、粒子群算法[6]等。同心圆环阵列是一簇以坐标原点为圆心，不同圆半径的圆环阵组成的阵列。其圆周对称特性使阵列的方向图形状相对稳定，且在俯仰角方向具有理想的辐射特性[7]，因而受到越来越多研究者的关注。具体的相关研究包括但不限于：Haupt[8]利用混合遗传算法，以圆环半径和圆环上阵元数量作为优化对象对同心圆环阵列进行综合，取得了较好的效果。陈客松等[9]通过约束同一圆环上阵元间距相等,利用修正遗传算法优化圆环的半径,获得最小的峰值旁瓣电平。栾晓明等[7]提出了一种基于差分进化算法的稀布同心圆环阵列半径和阵元间距的联合优化方法，获得了阵列响应低峰值旁瓣电平的阵元分布。这些研究都是基于单一频率的研究，目前关于宽带的阵列优化方法的研究还很少。

本文提出了一种基于平面同心圆环阵列的阵型优化设计方法，构造了一种度量宽带点扩散函数旁瓣水平的目标函数，采用粒子群全局优化算法，有效降低了波束扫描区域的整体旁瓣水平。同时，目标函数的低维度流形表示以及粒子群的高效搜索使得优化算法快速收敛，有效地降低了优化设计耗时。

1 同心圆环阵列波束扫描模型

如图1所示，考虑一个由不同半径的圆环阵列构成的同心圆环传声器阵列，以同心圆环的圆心建立三维坐标系。阵列阵元位于xOy平面内，假设同心圆环阵列共有L个圆环，每个圆环的半径分别为r1、r2、…、rL、每个环上的阵元数量分别为H1、H2、…、HL，总的传声器个数为N。

图1 同心圆环阵列扫描模型

假设声源分布在距传声器阵列中心点为z0的扫描平面上，声源满足阵列远场条件，则第n个传声器接收到扫描点(θ,φ)处的声压信号表示为

p(xn)=ejk·xn

(1)

式中：xn=[xn,yn](n=1,2,…,N)是第n个阵元的位置向量；k为波导矢量，即

(2)

式中：c为声音在介质中的传播速度；f为入射声波频率；θ为入射波方位角；φ为入射波俯仰角。

假设扫描平面上有一单位幅度的声源位于(0,0,z0)，其波矢量为k0，则传声器阵列接收到的声压为

(3)

一个扫描点沿着方向k的导向矢量为

(4)

则基于常规波束形成得到的原始波束图为

(5)

由于传声器阵列位于z=0平面，则式(5)改写为

(6)

具体到同心圆环阵列，阵元的位置坐标为

[xlm,ylm]T=[rlcosφlm,rlsinφlm)]T

(7)

式中：rl为第l个圆环的半径；φlm表示第l个圆环上的第m个阵元与x轴的夹角。则同心圆环阵列的扫描声源的原始波束图为

(8)

定义波束旁瓣区域为Θ，则最大旁瓣级峰值(PPSLL)是旁瓣区域的最大值为

(9)

波束主瓣宽度(PMLW)为

(10)

式中D为圆环阵列最大直径。文献[10]定义了以下表示旁瓣能量水平的函数(即M函数)：

(11)

M函数反映了图1中环形扫描区域辐射声源能量水平，Ω包含所有感兴趣的扫描方向，k=|k|，φmin和φmax分别代表环形扫描区域的最小和最大俯仰角。为了最小化旁瓣，选择的区域应包含所有感兴趣的扫描方向但不包含主瓣，故式(11)改写为

(12)

Amn=J1(kdmnsinφmax)sinφmax-

J1(kdmnsinφmin)sinφmin

(13)

式中：dmn是第m个阵元与第n个阵元之间的距离；J1是第一类的一阶贝塞尔函数。

2 基于粒子群算法的阵型优化

2.1 构造适应度函数

由式(8)可知，fsinφ可看做原始波束图的综合自变量，而sinφ与扫描区域的范围相关。综合自变量不变时，原始波束图不变，由此可以推断，低频时扫描区域的原始波束图与高频时缩小扫描区域的原始波束图相似。因此，可将宽带范围最低旁瓣水平约束等效为最高频率点局部扫描区域最低旁瓣水平约束。以需求频率范围上限确定目标函数，构造如下适应度函数：

J(x1…,xN)=w1M1(x1…,xN)+

w2M2(x1…,xN)

(14)

式中：M1(x1，…，xN)，M2(x1，…，xN)分别为不同积分区域内的M函数；w1,w2分别为两个函数值的权重。

另外，由于第一旁瓣存在于整个设计频带范围内，我们的设计目标是在减小第一旁瓣的同时最小化整体旁瓣水平，因此需要确定合适的权重以及M1、M2函数合适的积分边界φmin,1、φmax,1和φmin,2、φmax,2。对于M1函数，φmin,1为主瓣第一次过零处的位置，φmax,1为第一旁瓣第一次过零处的位置；对于M2函数，φmin,2为M1函数所确定的φmax,1，φmax,2为实际扫描区域的最大俯仰角。

2.2 阵形设计与优化

实际工程应用中，阵列的布局区域有一定限制。为了降低优化变量维度，本文假设环状区域中内圈圆半径为rmin，外圈圆半径为rmax，在约束区域内设计L个圆环组成同心圆环阵列，各圆环的半径分别为r1、r2、…、rL，每个环上的阵元数目分别为H1、H2、…、HL，并且每个环上的阵元等角度分布，区域内阵元数量为N，约束所有阵元的最小间距为dmin。由于阵元的最小间距限制，则相邻圆环之间的间隔至少为dmin。故采取文献[11]中提出的余量思想：选取各个圆环的半径余量Δr1、Δr2、…、ΔrL-1作为半径的优化参量以减少无效空间的搜索。半径余量为相邻两个圆环之间除最小间距dmin外的其他可优化的间距，如图2所示。图中阴影部分表示半径余量，空白部分表示最小阵元间距,则有：

图2 半径余量示意图

(15)

式中：rl为第l个圆环的半径；Δrl为第l环和第l+1环之间的可优化距离。

定义每个圆环上第一个阵元与水平方向的偏转角φ1、φ2、…、φL，并将其作为待优化变量。第m个阵元是指坐标系中从x轴起始，沿逆时针方向第m个阵元，则第l个圆环上的第m个阵元的偏转角可表示为

(16)

综合以上约束，对于粒子群优化算法，定义待优化粒子为X=(Δr1,Δr2，…，ΔrL-1,φ1,φ2…φL)，构建以下粒子群优化模型：

minJ(X)

s.t. J(X)=w1M1(X)+w2M2(X)

0≤Δrl≤rmax-rmin-(L-1)dmin

(17)

粒子群优化算法(PSO)过程中每个优化问题的解都是搜索空间中的一只鸟，称为“粒子”。所有的粒子都有一个由被优化的适应度函数决定的适应度值(Fitness Value)，每个粒子根据飞行速度决定其飞行的方向和距离，群体中的所有粒子追随当前的最优粒子在解空间中搜索。对于本文的阵型优化，粒子的解就是能代表阵元位置的一组参数向量Xopt,适应度函数就是2.1节所提出的函数。本文需要根据选取的(2L-1)个待优化参数布置阵元，用一个(2L-1)维的向量表示一个粒子，随机生成I个粒子，其中第t次迭代第i个粒子的位置Xi和速度vi分别表示为

(i=1,2,…,I)

(18)

(i=1,2,…,I)

(19)

第i个粒子当前t次迭代为止搜索到的最优布阵方式即为个体极值，记为

(20)

整个粒子群当前t次迭代为止搜索到的最优布阵方式为全局极值，记为

(21)

根据下式更新粒子的速度和位置：

(22)

(23)

式中：w为惯性权重；c1,c2为加速常数；s1,s2为两个在[0,1]范围内变化的随机数。本文w采用随迭代次数不断减小的惯性权重：

(24)

式中：wmax为最大的惯性权重；wmin为最小的惯性权重；T为总迭代次数；t为当前迭代次数。当w较大时，粒子搜索倾向于全局搜索；当w较小时，粒子可以在局部区域获得更精确的极值。因此，采用随迭代次数不断减少的惯性权重使得搜索过程首先进行全局搜索，然后在逐步减小的搜索区域进行细化搜索以获得更精确的解。整个算法流程如图3所示。

图3 粒子群算法流程图

3 仿真实验

3.1 宽带转最高频点可行性分析实验

为验证2.1节提出的将宽带范围最低旁瓣水平约束等效为最高频率点局部扫描区域最低旁瓣水平约束方法的可行性，本文设计的同心圆环阵参数如表1所示。每个环的半径相差0.007 2 m，阵列平面距离扫描平面2 m，扫描平面为2.4 m×2.4 m的方形网格点，声源频率范围为8～24 kHz，得到的原始波束图沿x轴的频率切片图如图4所示。

表1 同心圆环参数

图4 原始波束图沿x轴的频率切片图

由图4可以看出，随着频率的增大，主瓣宽度逐渐减小，旁瓣向主瓣靠拢，频率越高，则出现在视野范围内的旁瓣越多。因此，低频时扫描区域的原始波束图与高频时缩小扫描区域的原始波束图相似的推论是正确的。

3.2 两种适应度函数对比实验

利用粒子群工具包，分别以本文所提出的适应度函数和最高旁瓣级(PSLL)为优化目标，对同心圆环阵列进行优化设计，同心圆环阵的相关参数如表2所示。粒子群算法的参数如表3所示。两种适应度函数都在相同条件下进行仿真，阵列平面距离扫描平面2 m，扫描平面为2.4 m×2.4 m的方形网格点，网格分辨率为0.01 m。根据环形的尺寸可以求得主瓣边界对应的俯仰角φmin,M1选择最大俯仰角为8°，M2选择最小俯仰角为8°，最大俯仰角为33°。设w1=5，w2=0.1，声源频率为48 kHz，声源强度为30 dB。

表2 同心圆环参数

表3 粒子群算法参数

针对两种适应度函数分别进行10次优化求解实验，取仿真中最好的结果，算法最优值随迭代次数的变化曲线如图5所示。图中黑色曲线表示全局最优值随迭代次数的变化情况，蓝色点线代表种群的平均适应度值随迭代次数的变化情况，该曲线很好地反映了粒子群算法根据全局最优值和局部最优值不断迭代，并最终收敛到最优解的过程。

图5 两种适应度函数迭代对比图

将本文提出的适应度函数称为第一种适应度函数，最高旁瓣级函数(PSLL)称为第二种适应度函数。图5(a)是第一种适应度函数的迭代曲线，由图5(a)可看出迭代约200次后开始收敛。图5(b)是第二种适应度的迭代曲线，迭代300次开始收敛。对于第一种适应度函数，平均每次迭代耗时1.436 1 s；而第二种适应度函数平均每次迭代耗时17.876 s。这是由于第一种适应度函数的计算量只与传声器的数量有关，第二种适应度的计算量与传声器的数量、扫描的网格点数有关，而实际应用中往往需要高精度的波束形成图，扫描的网格点数远大于传声器数量，故第一种适应度函数平均每次迭代的运算量远小于第二种适应度函数。由此可见，本文所提出的适应度函数对宽带信号可以有效降低旁瓣，说明优化方案可行。选取本文提出适应度函数进行寻优可以减少运算量，提高计算速度。

这两种适应度函数得到的优化阵型和原始波束形成图如图6所示。表4为两种阵型的参数对比。

图6 两种适应度函数得到的优化阵型比较

表4 两种阵型不同频率的最高旁瓣级对比

通过图6(b)、(d)和表4可以看出，以第一种适应度函数优化得到的阵列一和以第二种适应度函数优化得到的阵列二拥有相似的主瓣宽度，48 kHz时均为0.32 m。在48 kHz时，虽然阵列二的最高旁瓣级为-19.55 dB，低于阵列一最高旁瓣级，但是该最高旁瓣出现在离主瓣较近的第一旁瓣处，而阵列一的最高旁瓣出现在远离主瓣视野的四角处。另外，阵列一的第一旁瓣级比阵列二低10.71 dB，这表明第一种适应度函数能有效降低主瓣附近的第一旁瓣高度。由表4可见，阵列一低频部分的最高旁瓣级总是低于阵列二，这说明对于宽频波束图，阵列一的整体旁瓣水平低于阵列二。以上对比分析证实了本文适应度函数的有效性。

3.3 权重系数对第一种适应度函数优化结果的影响

设置不同的权重系数会对实验结果产生不同的影响。为了说明实验第3.2节设置的权重的必要性，本文设计了另外两组不同权重的仿真实验，其他条件同仿真实验第3.2节，得到的阵型原始波束图如图7所示。对比图7(a)和图6(b)可以看出，相比于两个权重系数相等的情况，w1大于w2时可使靠近主瓣的第一旁瓣级更低，而其他范围的旁瓣水平相似；对比图7(b)和图6(a)可以看出，相比于两个权重系数相等的情况下，w1

图7 两组权重得到的原始波束图

3.4 两种优化方式对比实验

为了说明本文所提出的加入偏转角的优化方式能有效降低阵列在宽频信号的旁瓣水平，本节设计了仅优化半径的仿真实验。采用本文提出的适应度函数，其他条件同仿真实验第3.2节，得到仅优化半径的优化阵型和原始波束图如图8所示，所得阵型参数与前文阵型一对比如表5所示。

图8 仅优化半径得到的阵型

表5 两种阵型参数对比

由图8和表5可以看出，对于仅优化半径得到的阵列波束形成图，其主瓣宽度和整体旁瓣能量水平都远不及优化半径联合偏转角度的优化方式，这表明本文提出优化方式的有效性。

4 结束语

针对实际工程中约束阵列阵元总数、最大孔径、最小阵元间距等条件下的阵形优化设计问题，提出了一种基于粒子群优化算法且适用于宽带信号的同心圆环阵列阵型优化方法。该方法以约束的最高频率构造了反映宽带信号在扫描区域内旁瓣水平的适应度函数，同时以圆环半径和阵元偏转角作为联合优化变量，基于粒子群优化算法对阵型优化问题进行求解。数值结果仿真表明，本文所提出的适应度函数在优化效率上优于传统的以最高旁瓣级作为适应度函数的方法，优化后的阵型在扫描区域内旁瓣能量水平更低，更适合于宽带信号场景。下一步的工作重点为声像仪信号采集系统的实现，将本文设计的同心圆环传声器阵型应用到该系统中，对优化后阵型的性能进行实验验证。