基于混合引力搜索算法的立体传声器阵列优化设计方法

摘要

本发明提出一种基于混合引力搜索算法的立体传声器阵列优化方法，包括：构建立体传声器阵列的主要参数在混合引力搜索算法框架下的适应度函数，最优化立体传声器阵列在一定旁瓣级下具有较窄主瓣宽度的阵列结构。该方法优化设计的立体传声器阵列，能弥补传统平面阵定位结果无法分辨出前后波达方向的不足，抑制后向来波，实现对复杂声环境下目标噪声源的被动定位。

说明书

【技术领域】

本发明属于传声器技术领域，特别涉及一种立体传声器阵列优化设计方法。

【背景技术】

实际情况下，变电站的噪声源较复杂，平面阵的定位结果无法分辨出前后波达方向，而体 积阵能够实现对全空间的扫描，抑制后向来波。体积阵性能的优劣很大程度上受阵型结构，阵 元数量以及阵列孔径的影响。在相同的阵元数量和阵列孔径的前提下，不同形式阵列结构的波 束方向图的主瓣宽度和最大旁瓣级不一样，产生的波束响应也不同，对噪声源识别定位的性能 影响很大。主瓣宽度和旁瓣级，是相互矛盾的性能指标，无法得到同时具有最窄主瓣宽度和最 低旁瓣级的最优阵形，因此需要综合考虑，取得二者的折中。如何结合声源特征，优化设计体 积阵列结构是多传感器阵列声成像定位的一个难点，具有重要的意义。

就智能优化算法来说，目前应用比较广泛的有遗传算法、蚁群算法、神经网络、粒子群算 法。这些算法大多启蒙于自然现象或自然规律，各具优缺点，适用面也不同。遗传算法对参数 的依赖性较强，和粒子群算法一样容易早熟，陷入局部最优；蚁群算法和模拟退火虽然不受初 始条件的约束，但计算量大，优化过程较长。

【发明内容】

本发明的目的在于提出一种基于混合引力搜索算法的立体传声器阵列优化方法，以克服现 有技术所存在的问题；本发明方法优化设计的立体传声器阵列，能弥补传统平面阵定位结果无 法分辨出前后波达方向的不足，抑制后向来波，实现对复杂声环境下目标噪声源的被动定位。

为了实现上述目的，本发明采用如下技术方案：

基于混合引力搜索算法的立体传声器阵列优化设计方法，包括以下步骤：

1)、构建立体传声器阵列优化参数模型：

待优化多臂立体传声器阵列有N个阵臂，类似伞骨架呈角度张开，每个臂上有M个麦 克风，共计N*M个麦克风；阵臂之间的夹角α＝360°/N，固定每个阵臂第一个阵元到阵中心 的距离l₁＝0.2m，最后一个阵元到阵中心的距离只考虑单个臂上阵元的距离 L＝[l₂,...,l_M-1]的变化，其余阵臂依次复制；与此同时，为防止阵元重叠，限制每条臂上相邻两 个阵元的最小间距为l₀，即

l_m+1-l_m≥l₀,m＝1,2,…,M-1

2)、设定引力系统粒子形式及其优化目标：

假设一个包含有N₀个粒子的引力系统，系统维度为D，这个引力系统的所有粒子用向量 表示为经过步骤1)处理之后，引力系统粒子的形式设置为

设频率为f₀的声信号所对应的常规波束旁瓣级指标为SLL₀，要求旁瓣级小于该指标，若 阵列结构为P，则适应度函数Fit(P,f₀)可以表示为

$Fit(P,f_0)=\left(\begin{array}{cc}MLW(P,f_0)& SLL(P,f_0)\le {\mathrm{SLL}}_0\\ +\infty & SLL(P,f_0)>{\mathrm{SLL}}_0\end{array}\right)$

当旁瓣级满足指标要求时，则以主瓣宽度为适应度函数返回值；若旁瓣级不满足指标要求， 则返回给适应度函数的值为无穷大；因此在一定旁瓣级的要求下，优化目标函数为

minFit(P,f₀),s.t.SLL(P,f₀)≤SLL₀

3)、具体的迭代优化过程：

将粒子的位置用向量表示为i＝1,2,…,d,…,D；D为引力系统的 维度；根据牛顿万有引力定律：任意两个质点在连心线方向上有力的相互吸引，设第i个粒子 受到第j个粒子的引力则第d维方向上的引力分量

$F_{ij}^d\left(t\right)=G\left(t\right)·\frac{M_i\left(t\right)\times M_j\left(t\right)}{\left|\right|x_i\left(t\right),x_j\left(t\right)|{|}_2+ϵ}·(x_j^d\left(t\right)-x_i^d\left(t\right))$

其中，M_i(t)、M_j(t)是第i个粒子和第j个粒子被动引力质量；是第j个粒子 和第i个粒子在d维方向上的坐标，||x_i(t),x_j(t)||₂是这两个粒子的欧式距离；ε是一个很小的 常数，一般设为0；G(t)是随迭代时间递减的万有引力常系数，且

$G\left(t\right)=G_0·\exp (-k·\frac{t}{T})$

G₀是引力常系数的初始值，k是一个递减的系数，令G₀＝1，k＝20；t为当前迭代次数， T为设定的最大迭代次数；

给任意两个粒子在d维方向上的引力分量加一个在(0,1)范围内变化的随机权值rand()， 则第i个粒子在d维方向上受到的其他粒子的引力随机加权总和

$F_i^d\left(t\right)=\underset{j=1,j\ne i}{\overset{D}{\Sigma }}rand\left(\right)·F_{ij}^d\left(t\right)$

这一迭代时刻，粒子在d维方向的引力加速度

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

其中，M_ii(t)是第i个粒子的惯性质量；粒子被动引力质量和惯性质量的计算依赖于适应 度函数Fit返回的函数值；假设被动引力质量、惯性质量和个体质量三者相等，记为M_i，M_i根 据以下几个公式进行计算：

$M_i\left(t\right)=\frac{m_i\left(t\right)}{\underset{j=1}{\overset{N}{\Sigma }}{m}_j\left(t\right)}$

$m_i\left(t\right)=\frac{{\mathrm{Fit}}_i\left(t\right)-{\mathrm{Fit}}_{worst}\left(t\right)}{{\mathrm{Fit}}_{best}\left(t\right)-{\mathrm{Fit}}_{worst}\left(t\right)}$

其中Fit_best＝min(Fit_i(t)),i∈[1,N]，Fit_worst＝max(Fit_i(t)),i∈[1,N]；Fit_i(t)是第i个粒子在 t时刻的适应度函数值。

混合引力搜索算法的粒子速度和位置按以下两个公式进行更新：

$v_i^d(t+1)=w·v_i^d\left(t\right)+c_3·rand\left(\right)·a_i^d\left(t\right)+c_4·rand\left(\right)·(g_{best}^d\left(t\right)-x_i^d\left(t\right))$

$x_i^d(t+1)=x_i^d\left(t\right)+v_i^d(t+1)$

权重因子c₃、c₄取0.5和1.5；惯性权值w取(0,1)内的随机数；是指当前迭代时刻 粒子在d维方向上的最优解；当循环达到最大次数T时，终止迭代。

N＝7，M＝9；D＝1m。

阵列投影直径为2m，阵列中心高度为1.5m到2.5m可调，阵列中心配置高分辨率摄像头 1个，阵列工作时张开并锁定，不工作时能够折叠存放。

步骤2)中粒子个数为25个，系统维度D＝8。

基于混合引力搜索算法的立体传声器阵列优化设计方法，具体包括以下步骤：

1)、目标声源的频率为1000Hz，要求优化后阵列的旁瓣级不高于-8dB，优化目标为满足 旁瓣级的情况下获得最窄的主瓣宽度；阵列是投影直径为2米的7臂阵，固定每个阵臂第一个 阵元到阵中心的距离l₁＝0.2m，每个臂阵元数为9，最后一个阵元到阵中心的距离

2)、设定引力系统粒子形式及其优化目标：

将体积阵的参数设置为引力系统粒子，粒子个数为25个，系统维度D＝8；这个引力系统 的所有粒子可以用向量表示为X＝[x₁,x₂,…,x₂₅]；设频率为1000Hz的声信号所对应的常规波 束旁瓣级指标为-8dB，要求旁瓣级小于该指标；当旁瓣级满足指标要求时，则以主瓣宽度为适 应度函数返回值；若旁瓣级不满足指标要求，则返回给适应度函数的值为无穷大；因此在一定 旁瓣级的要求下，优化目标函数为

minFit(P,1000),s.t.SLL(P,1000)≤-8

3-1)、初始化引力系统粒子；

3-2)、通过初始粒子计算阵列的波束指标，进而得到适应度函数值；

3-3)、根据适应度函数值计算粒子的引力质量和惯性质量；计算粒子引力，引力加速度， 按照粒子的速度和位置进行更新；

3-4)、更新后的粒子按照步骤3-2)得到适应度函数值，挑选适应度函数值最小的引力子 为当前最优参数，若满足优化目标函数，迭代停止；若不满足，则重复步骤3-3)、3-4)；若达 到最大迭代次数2000，目标函数仍不满足，则终止迭代，更换初始化粒子，重复步骤3-2)、 3-3)、3-4)。

相对于现有技术，本发明具有以下有益效果：

本发明结合了粒子群算法的混合引力搜索算法既改善了传统引力搜索算法在求解复杂多 峰问题时容易陷入局部最优的现象，又比原始的粒子群算法具有更快的迭代收敛速度，大大提 高了阵列优化的有效性；所设计的立体传声器阵列成功抑制了后向噪声的干扰，实现了对复杂 声环境下目标噪声源精确定位的要求。

【附图说明】

图1：待优化的7臂体积阵；

图2：优化后的7臂体积阵；

图3：优化前后的定位结果比较；图3(a)为优化前随机产生的阵列参数得到的定位结 果示意图；图3(b)为优化后的定位结果示意图。

图4为本发明方法的流程图。

【具体实施方式】

一、伞形阵模型

假设一个近场常规波束模型，将M元伞形阵垂直放置在空间中，建立三维直角坐标系， 各个阵元的三维空间坐标可以表示为，且p_i(x_i,y_i,z_i),i＝1,2,…M，且

其中，l_i是第i个阵元到基阵原点的距离，是第i个阵元所在阵臂与z轴的夹角，θ_i是该 阵臂在xOy平面的投影与x轴的夹角。阵列结构可以用矩阵表示为

P＝[p₁,p₂,p₃,…,p_M]^T

在距离基阵z＝H₀的平面上设立一个频率为f₀的点声源，其三维坐标为B(x_B,y_B,H₀)。为 方便研究，忽略多途效应的影响。若s_i(t)为第i个阵元接收到的声源B的声信号，那么基阵的 输入信号x(t)可用向量表示为

x(t)＝[s₁(t),s₂(t),…,s_M(t)]

则基阵信号输出y(t)的向量表达式为

$y\left(t\right)={\Sigma }_{i=1}^M{w}_i{s}_i\left(t\right)=w^{H}x\left(t\right)$

上式中，w为球面波相幅补偿量的权向量形式，该权值能够很好的抑制后向声源A的影响，(·)^H表示复共轭转置。

故该体积阵在来波方向的信号波束图为

其中w₀和分别是阵列的加权向量和流形向量。对波束图进行归一化处理，使期 望方向上的波束响应为0dB，则主瓣宽度MLW可定义为波束图-3dB点之间的距离Δφ_-3dB，旁 瓣级SLL即为除主瓣以外第一旁瓣所对应的幅度响应值，显然与频率和阵列结构有关。

设频率为f₀的声信号所对应的常规波束旁瓣级指标为SLL₀，即要求旁瓣级小于该指标， 阵列结构为P，则适应度函数Fit(P,f₀)可以表示为

$Fit(P,f_0)=\left(\begin{array}{cc}MLW(P,f_0)& SLL(P,f_0)\le {\mathrm{SLL}}_0\\ +\infty & SLL(P,f_0)>{\mathrm{SLL}}_0\end{array}\right)$

当阵列结构为P的条件下，旁瓣级满足指标要求，则以主瓣宽度为适应度函数；若旁瓣 级不满足指标要求，则返回给适应度函数的值为无穷大。因此在一定旁瓣级的要求下，优化目 标函数为

minFit(P,f₀),s.t.SLL(P,f₀)≤SLL₀

二、混合引力搜索算法

假设一个包含有N个粒子的引力系统，系统维度为D，这个引力系统的所有粒子可以用 向量表示为X＝[x₁,x₂,…,x_N]，则第i个粒子的位置可以用向量表示为i＝1,2,…,d,…,D。根据牛顿万有引力定律：任意两个质点在连心线方向上有力的相互吸引， 设第i个粒子受到第j个粒子的引力则第d维方向上 的引力分量为

$F_{ij}^d\left(t\right)=G\left(t\right)·\frac{M_i\left(t\right)\times M_j\left(t\right)}{\left|\right|x_i\left(t\right),x_j\left(t\right)|{|}_2+ϵ}·(x_j^d\left(t\right)-x_i^d\left(t\right))$

其中，M_i(t)、M_j(t)是第i个粒子和第j个粒子被动引力质量；||x_i(t),x_j(t)||₂是这两个粒 子的欧氏距离；ε是一个很小的常数；G(t)是随迭代时间递减的万有引力常系数，

$G\left(t\right)=G_0·\exp (-k·\frac{t}{T})$

上式中，G₀为引力常系数的初始值，k是一个递减的系数，t为当前迭代次数，T为设定 的最大迭代次数。

为保证算法的随机性，我们给任意两个粒子在d维方向上的引力分量加一个在(0,1)范围 内变化的随机权值，则第i个粒子在d维方向上受到的其他粒子的引力随机加权总和

$F_i^d\left(t\right)=\underset{j=1,j\ne i}{\overset{D}{\Sigma }}rand\left(\right)·F_{ij}^d\left(t\right)$

因此，这一迭代时刻，该粒子在d维方向的引力加速度

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

其中，M_ii(t)是第i个粒子的惯性质量。粒子被动引力质量和惯性质量的计算依赖于适应 度函数F返回的函数值。假设被动引力质量、惯性质量和个体质量三者相等，记为M_i，M_i可 以根据以下几个公式进行计算：

$M_i\left(t\right)=\frac{m_i\left(t\right)}{\underset{j=1}{\overset{N}{\Sigma }}{m}_j\left(t\right)}$

$m_i\left(t\right)=\frac{{\mathrm{Fit}}_i\left(t\right)-{\mathrm{Fit}}_{worst}\left(t\right)}{{\mathrm{Fit}}_{best}\left(t\right)-{\mathrm{Fit}}_{worst}\left(t\right)}$

若是求解最小值，则Fit_best＝min(Fit_i(t)),i∈[1,N]，Fit_worst＝max(Fit_i(t)),i∈[1,N]；若求 解最大值，则Fit_best与Fit_worst内容互换。

在标准的GSA算法中，第i个粒子的速度v_i和位置x_i按照如下公式进行更迭：

${\nu }_i^d(t+1)=rand\left(\right)·{\nu }_i^d\left(t\right)+a_i^d\left(t\right)$

$x_i^d(t+1)=x_i^d\left(t\right)+v_i^d(t+1)$

标准PSO算法按照如下公式进行更迭：

v_i(t+1)＝w·v_i(t)+c₁·rand()·[p_best(t)-x_i(t)]+c₂·rand()·[g_best(t)-x_i(t)]

x_i(t+1)＝x_i(t)+v_i(t+1)

其中，c₁、c₂为权重因子，通常等于2；rand()是在(0,1)内取值的随机数；w是惯性权值， 起到平衡全局和局部搜索能力的作用，一般取0.9到0.2递减；p_best(t)和g_best(t)是通过比较粒 子自身和粒子群整体的适应度函数值后，选出的当前迭代时刻t的局部最优位置和全局最优位 置。

与GSA算法不同的是，PSO算法中粒子速度和位置的更新迭代是将粒子的每一维度合为 一个整体来参与迭代计算，而GSA是对粒子速度和位置每个维度方向的分量进行迭代计算， 所以，GSA的计算复杂度比PSO要高一些。

综合PSO和GSA的粒子速度位置更迭公式，得到基于粒子群算法的混合引力搜索算法(H PSOGSA)的粒子速度和位置迭代公式如下：

$v_i^d(t+1)=w·v_i^d\left(t\right)+c_3·rand\left(\right)·a_i^d\left(t\right)+c_4·rand\left(\right)·(g_{best}^d\left(t\right)-x_i^d\left(t\right))$

$x_i^d(t+1)=x_i^d\left(t\right)+v_i^d(t+1)$

公式中的权重因子c₃、c₄通常设为0.5和1.5；惯性权值w取(0,1)内的随机数；是 指当前迭代时刻粒子在d维方向上的最优解。

一种基于混合引力搜索算法的立体传声器阵列优化方法，包括以下步骤：

1)、构建多臂立体传声器阵列优化参数模型：

待优化多臂立体传声器阵列有N个阵臂，类似伞骨架呈俯仰角度张开每个臂上有M个麦克风，共计N*M个。阵臂之间的夹角α＝360°/N，固定每个阵臂第一个阵 元到阵中心的距离l₁≈0.2m，最后一个阵元到阵中心的距离R＝1m，只考虑 单个臂上阵元的距离L＝[l₂,l₃,...,l_M-1]的变化，其余阵臂依次复制。与此同时，为防止阵元重 叠，限制每条臂上相邻两个阵元的最小间距为l₀，即

l_m+1-l_m≥l₀,m＝1,2,…,M-1

其中l₀是一个很小的正数，令l₀＝0.05m。

2)、设定引力系统粒子形式及其优化目标：

将体积阵的参数设置为引力系统粒子，系统维度D＝8。假设一个包含有 N个粒子的引力系统，系统维度为D，这个引力系统的所有粒子可以用向量表示为 X＝[x₁,x₂,…,x_N]。设频率为f₀的声信号所对应的常规波束旁瓣级指标为SLL₀，要求旁瓣级 小于该指标，若阵列结构为P，则适应度函数Fit(P,f₀)可以表示为

$Fit(P,f_0)=\left(\begin{array}{cc}MLW(P,f_0)& SLL(P,f_0)\le {\mathrm{SLL}}_0\\ +\infty & SLL(P,f_0)>{\mathrm{SLL}}_0\end{array}\right)$

当旁瓣级满足指标要求时，则以主瓣宽度为适应度函数返回值；若旁瓣级不满足指标要求，则 返回给适应度函数的值为无穷大。因此在一定旁瓣级的要求下，优化目标函数为

minFit(P,f₀),s.t.SLL(P,f₀)≤SLL₀

3)、具体的迭代优化过程：

将粒子的位置用向量表示为i＝1,2,…,d,…,D。根据牛顿万有 引力定律：任意两个质点在连心线方向上有力的相互吸引，设第i个粒子受到第j个粒子的引 力$F_{ij}\left(t\right)=[F_{ij}^1\left(t\right),F_{ij}^2\left(t\right),...,F_{ij}^d\left(t\right),...,F_{ij}^D\left(t\right),],$则第d维方向上的引力分量

$F_{ij}^d\left(t\right)=G\left(t\right)·\frac{M_i\left(t\right)\times M_j\left(t\right)}{\left|\right|x_i\left(t\right),x_j\left(t\right)|{|}_2+ϵ}·(x_j^d\left(t\right)-x_i^d\left(t\right))$

其中，M_i(t)、M_j(t)是第i个粒子和第j个粒子被动引力质量；是第j个粒子和第 i个粒子在d维方向上的坐标，||x_i(t),x_j(t)||₂是这两个粒子的欧式距离；ε是一个很小的常数， 一般取值为0；G(t)是随迭代时间递减的万有引力常系数，且

$G\left(t\right)=G_0·\exp (-k·\frac{t}{T})$

G₀是引力常系数的初始值，k是一个递减的系数，令G₀＝1，k＝20；t为当前迭代次数，T为 设定的最大迭代次数。为保证算法的随机性，我们给任意两个粒子在d维方向上的引力分量加 一个在(0,1)范围内变化的随机权值rand()，则第i个粒子在d维方向上受到的其他粒子的引力 随机加权总和

$F_i^d\left(t\right)=\underset{j=1,j\ne i}{\overset{D}{\Sigma }}rand\left(\right)·F_{ij}^d\left(t\right)$

这一迭代时刻，粒子在d维方向的引力加速度

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

其中，M_ii(t)是第i个粒子的惯性质量。粒子被动引力质量和惯性质量的计算依赖于适应度函 数Fit返回的函数值。假设被动引力质量、惯性质量和个体质量三者相等，记为M_i，M_i可以 根据以下几个公式进行计算：

$M_i\left(t\right)=\frac{m_i\left(t\right)}{\underset{j=1}{\overset{N}{\Sigma }}{m}_j\left(t\right)}$

$m_i\left(t\right)=\frac{{\mathrm{Fit}}_i\left(t\right)-{\mathrm{Fit}}_{worst}\left(t\right)}{{\mathrm{Fit}}_{best}\left(t\right)-{\mathrm{Fit}}_{worst}\left(t\right)}$

其中Fit_best＝min(Fit_i(t)),i∈[1,N]，Fit_worst＝max(Fit_i(t)),i∈[1,N]，Fit_i(t)是第i个粒子在 t时刻的适应度函数值。

混合引力搜索算法的粒子速度和位置按以下两个公式进行更新：

$v_i^d(t+1)=w·v_i^d\left(t\right)+c_3·rand\left(\right)·a_i^d\left(t\right)+c_4·rand\left(\right)·(g_{best}^d\left(t\right)-x_i^d\left(t\right))$

$x_i^d(t+1)=x_i^d\left(t\right)+v_i^d(t+1)$

权重因子c₃、c₄通常设为0.5和1.5；惯性权值w取(0,1)内的随机数；是指当前迭 代时刻粒子在d维方向上的最优解。当循环达到最大次数T时，终止迭代。

下面结合具体实例，对本发明做进一步的描述。

请参阅图4所示，一种基于混合引力搜索算法的立体传声器阵列优化方法，包括以下步骤：

1)、根据阵列的使用目的构造阵列优化目标函数。假设目标声源的频率为1000Hz，要求 优化后阵列的旁瓣级不高于-8dB，优化目标为满足旁瓣级的情况下获得最窄的主瓣宽度。阵列 是投影直径为2米的7臂阵，固定每个阵臂第一个阵元到阵中心的距离l₁≈0.2m，每个臂阵元 数为9，最后一个阵元到阵中心的距离如图1所示。

2)、设定引力系统粒子形式及其优化目标：

将体积阵的参数设置为引力系统粒子，粒子个数为25个，系统维度D＝8。这个引力系统 的所有粒子可以用向量表示为X＝[x₁,x₂,…,x₂₅]。设频率为1000Hz的声信号所对应的常规波 束旁瓣级指标为-8dB，要求旁瓣级小于该指标。当旁瓣级满足指标要求时，则以主瓣宽度为适 应度函数返回值；若旁瓣级不满足指标要求，则返回给适应度函数的值为无穷大。因此在一定 旁瓣级的要求下，优化目标函数为

minFit(P,1000),s.t.SLL(P,1000)≤-8

3)、初始化引力系统粒子；

4)、通过初始粒子计算阵列的波束指标，进而得到适应度函数值；

5)、根据适应度函数值计算粒子的引力质量和惯性质量。计算粒子引力，引力加速度，按 照粒子的速度和位置进行更新。

6)、更新后的粒子按照步骤4)得到适应度函数值，挑选适应度函数值最小的引力子为当 前最优参数，若满足优化目标函数，迭代停止；若不满足，则重复步骤(5)(6)。若达到最大 迭代次数2000，目标函数仍不满足，则终止迭代，更换初始化粒子，重复(4)(5)(6)。

得到如表1所示的最优引力子后，计算当前参数情况下，阵列在其他频率的性能指标，如 表2所示。图2是优化后的体积阵阵型。图3是优化前后阵列的声源定位效果：图3(a)表 示优化前，随机产生的阵列参数得到的定位结果，图3(b)表示优化后的定位结果。对比两 个基阵的定位结果可以发现，优化后的伞形阵对后向噪声的抑制作用更强；而且优化之后 1000Hz下的-3dB主瓣宽度由3.4m缩小到2m，旁瓣级虽然略有升高，但仍然满足性能指标。 可见经过优化后，阵列的空间分辨率有较大提升。

表1优化后的伞形阵的详细参数

表2部分频率下的性能指标