高 玥 卢 铃 吴 鸣 杨 军 曹 寅

(1 中国科学院大学 北京 100049)

(2 中国科学院噪声与振动重点实验室(声学研究所)北京 100190)

(3 国网湖南省电力有限公司电力科学研究院 长沙 410007)

(4 国网电力设施噪声与振动实验室 长沙 410007)

(5 机器视觉 语音及信号处理研究中心 吉尔福德 GU2 7XH)

0 引言

传声器阵型的配置在解决声成像问题中起着重要的作用。传声器阵列的设计在很大程度上取决于所期望的应用。Meyer 等[1]在2002年提出基于球谐分解的刚性球阵波束形成器，可在不改变波束模式的情况下将观测方向指向三维任意位置。Rafaely[2]提出了一种基于球谐的球阵列设计和分析框架，还提出了球面上传声器定位的备选空间采样方案。Balmages等[3]提出了开放双球面阵，解决了由于贝塞尔函数零点引起的空间模态对应频率上出现的数值病态问题。Carious 等[4]和Lamotte等[5]利用了双球面阵列，其内部为刚性阵列，非常适合中高频范围，而外部是一个开放球面阵列，目的是提高低频性能，平衡了内部阵列方向性好但分辨率差和外部阵列的分辨率好但方向性差的问题，从而实现较好性能。上述研究在低频段均存在着所需阵列孔径较大不易制作和运输的问题。

基于传声器阵列的声成像方法以其简单和高效的优势有着广阔的应用背景，可分为阵元域和模态域两种。传统阵元域的算法中一类为间接法，其利用不同阵元接收信号的到达时间差来估计声源方位。其中，广义互相关相位变换方法使用广义互相关方法进行时间延迟估计[6]，并对幅度谱引入加权函数做相位变换[7]进行声源方位估计，由于其复杂度较低因此得到了广泛的发展和应用，但其受混响的影响较大。在此基础上，Michaud 等[8]研究了实时追踪移动声源的问题，Grondin 等[9]在噪声条件下通过时频掩膜技术来提高定位精度。另一类直接法是计算所有方位上的代价函数从而估计声源位置。其中可控响应功率法[10]的具体实现与波束形成器的设计方法相关，在智能音箱等领域有着较大应用，研究重点在于减少空间扫描的网络格点数；最大似然法[11]处理相干信号能力较好且统计特性优秀但是计算复杂度高；基于子空间的方法拥有超分辨率，其中非相干子空间法[12]的空间相干矩阵缺秩问题会导致算法性能下降，而相干子空间法则[13]需要一定的先验知识来获得聚焦矩阵。

传统的波束形成算法在低频段效果较差，且在相干声源、复杂声源存在的条件下算法性能会受限，因此本文采用广义逆波束形成(Generalized inverse beamforming, GIBF)算法来解决上述问题[14−15]。广义逆波束形成的改进主要以改进正则化方法为特征，求解反问题所涉及的优化参数对于准确的波束形成结果至关重要，Zavala 等[16]对此进行了讨论。为了验证算法并评估其性能，Bahr 等[17]将广义逆波束形成应用于与小型露天喷气式设施即NASA兰利静音流动设备相对应的实验基准数据集中，研究了由不同贡献者实施并应用于同一测试用例的几种阵列分析方法的可变性。与其他常见波束成形技术的比较可以显示广义逆波束形成在复杂应用中的潜力和优势。

在上述阵元域方法的基础上，模态域方法根据传声器阵列的拓扑结构，把接收到的信号从阵元域进一步展开到模态域，将环形阵列信号转换到环谐域，球形阵列信号转换到球谐域。由于模态域波束形成器的导向矢量具有频率无关的特性，可以直接实现频域平滑从而减少计算量，对于实际应用有很大的价值。Rafaely[18]在球谐域采用延时求和的方法设计波束形成器，Gao 等[19]将球谐域声源定位方法应用在水下环境中并分析其性能。

本文的核心思路为：提出了基于球谐分解的L2范数广义逆波束形成算法，并在分布式开放球阵的阵型基础上进行仿真及实验。得到如下结论：通过仿真及实验分析对比了不同算法在低频相干声源条件下声成像的准确性和对阵元位置误差的鲁棒性后，证明在此分布式开放球阵阵型下该算法具有较好的性能。

1 声成像算法原理

1.1 球谐系数估计法

Ueno 等[20]提出了适用于任意阵型和阵元指向性的全局谐波系数估计方法。通过贝叶斯准则，将截断后的后验均值作为区域中球谐波系数矢量的估计：

其中，区域中的任意一点用r=(r,θ,φ)表示，(r0,k)为在r0处展开的截断后的球谐波系数矩阵，k为波数，P为截断阶数；λ是与球谐波系数矢量的先验分布有关的常数，xe为接收的三维声场，Σ ∈CQ×Q是噪声的互功率谱密度矩阵。阵列由Q个阵元构成，分别分布在r1,···,rQ处，阵元指向性为c1,···,cQ，c是由球谐波分解系数cmn构成的矢量。ΞP(r0,k)是与不同展开中心展开球谐波系数之间的转换矩阵T及阵元指向性系数c(θ,φ)有关的矩阵，Ψ(k)可由ΞP(r0,k)表示[20]，即：

其中，(·)∗表示复共轭转置。

考虑幅度为S(w)的单频平面波从(θ0,φ0)方向入射，则在r处的频域观测信号的球谐波展开为[21]

其中，波矢k=(k,θ0,φ0)，jn(kr)为n阶球贝塞尔函数，Ymn(θ,φ)为球谐波函数。由于任意三维声场x(r,k)可以表示为球谐波展开的形式：

其中，αmn为球谐波系数，基函数定义为φmn(r,k)=则由公式(4)、公式(5)可得到若是D个平面波入射时，声场中的球谐波系数具有如下形式：

其中，Sd(ω)和(θd,φd)分别为dth 平面波的幅度和入射方向。

根据球谐波函数的正交性， 定义滤波器h(θ,φ)∈C∞，令滤波器系数为

则该滤波器输出为

可以看出，该滤波器具有理想的空间选择特性。因此通过改变滤波器的扫描角度则可得到方位谱的估计值，从而确定信号的入射角度。

1.2 广义逆波束形成

广义逆波束形成可以归类为逆方法，因为它需要对声源-传声器间的传播问题进行反演，可以建模为[22]

其中，A ∈CQ×N，q ∈CN×1，p ∈CQ×1，N为扫描网络格点数。由于一般情况下Q ≪N，则问题转化为求欠定方程的解。本文采用一范数及二范数两种方法来求解。

1.2.1 L1范数广义逆波束形成

L1范数的广义逆波束形成(L1-Generalized inverse beamforming, L1-GIBF)算法是用压缩感知(Compressive sensing, CS)来对稀疏的噪声源进行声成像的算法[23]。CS 为求解欠定问题Aq=p提供了一种解决方案，其中测量矩阵p、声源振幅q和线性传递矩阵A(即CS 中的感知矩阵)分别为Q×1 的向量、N ×1 的向量以及Q×N维矩阵，其中Q为组成传声器阵列的阵元数，N为格点数，且一般情况下Q ≪N。阵元处声信号可以由向量和矩阵形式表示为

其中，ε是由电噪声、实验误差、测量数据与模拟声压不匹配等引起的误差。A的列向量an等价于球面波束形成的导向矢量。由于Q远小于N，导致求解x的欠定问题没有唯一解。由于声源位于特定的局部区域内，因此解决不适定问题的一种可行方法是将稀疏度放在x上。则可将求解式(10)转化为约束问题：

但是L0范数的最小化问题是一个非凸问题，在计算上难度较大。当有足够稀疏的声源时，公式(11)等价于L1范数的最小化问题，则可以转化为凸优化问题从而有效地实现计算处理。又由于q中的一个元素乘以具有较大欧氏范数的A中的一列，在计算L1 范数最小化问题时趋向于一个非零解。为了防止这种偏差，将阵元处测量数据p和an分别按照欧几里得范数进行归一化。基于以上可得到

借助数学工具中凸优化算法建模软件包对式(12)进行求解。由于它的稀疏性，即使在单频情况下，此算法也能提供高分辨率的声成像结果。然而，当测量数据受到噪声污染时，稀疏解的非零分量会出现在与实际源位置不同的位置。

1.2.2 L2范数广义逆波束形成

L2 范数的广义逆波束形成算法(L2-Generalized inverse beamforming, L2-GIBF)是用Tikhonov 正则化[24]修正问题的方法来求解公式(9)欠定方程的解。通过添加残差范数项和正则化参数λ2修正等式后，得到正则化的解为

其中，A=USV∗为奇异值分解结果，U ∈CQ×Q是酉矩阵，S ∈RQ×Q是一个对角矩阵，其中包含的奇异值按降序s1≥s2≥···≥sQ≥0 排列，V ∈CN×Q满足V∗V=IQ，正则化系数λ2定义了滤波器的强度。

1.3 基于球谐分解的L2范数广义逆波束形成

基于球谐分解的L2 范数广义逆波束形成(Spherical harmonics decomposition-L2-Generalized inverse beamforming, SHD-L2-GIBF)算法的原理是先通过球谐分解(SHD)的预处理来稳定后续的求逆数值问题，尤其在中低频范围内效果更好，而后通过广义逆波束形成(GIBF)处理截断后的阵元接收数据进行声成像。

1.3.1 算法原理

pnm为复数球谐波系数，是由球傅里叶变换(Spherical Fourier transform, SFT)计算得到的球谐函数的谱，表示为

其中，p(θq,φq)为第q个阵元的声压分布，αnmq是由积分转换到累加和形式的近似，对第q个阵元的m阶n次项的加权系数。由于球形阵列传声器个数Q是有限的，因此p(θq,φq)可以通过截断P阶后的逆球傅里叶变换(Inverse spherical Fourier transform,ISFT)近似表示为[25]

则可将式(15)转换成矩阵形式表示为

其中，p为各个阵元采集声压数据经过短时傅里叶变换的频谱，pnm为球谐域谱，YP矩阵的阶数为M ×(P+1)2，包含在每个阵元上计算得到P阶的球谐函数Ymn(θq,φq)。

对式(16)中的线性系统进行反演，得到了SFT的数值解为式(17)：

其中，TP为YP的逆矩阵。

短时傅里叶变换后得到的信号频谱先经过ISFT 转换到模态域，然后通过SFT 转换回来，则实现了球谐分解步骤[25]。但是SFT和ISFT只能用于单球阵阵列条件下，为了克服这一限制，本文提出将球谐系数估计法引入此步球谐分解操作中。因此公式(16)、公式(17)可以用公式(1)、公式(8)的形式表示，则可在任意阵型下实现球谐分解预处理步骤。将预处理后的数据再进行本文1.2.2 节中叙述的二范数广义逆波束形成算法，实现基于球谐分解的L2范数广义逆波束形成算法。

1.3.2 正则化参数选择

对于公式(13)中正则化参数的选取方法，本文采用的为L-曲线准则[26]，在二维上，以lg-lg 为尺度，(‖p −Aq‖,‖q‖)的值构成了如图1所示形状类似于L 的曲线，故被称为L 曲线。为了平衡欠正则化和过正则化，在L曲线曲率最大的点(即由竖直到水平的转角)处选择正则化参数。以λ2为参数的L曲线的曲率函数可以表示为

图1 L 曲线图Fig.1 L–curve

其中，u(λ2)= lg‖p −Aq‖，v(λ2)= lg‖q‖，且“′”，“′′”表示对λ2求解一阶导数和二阶导数。

2 声成像算法性能仿真分析

2.1 仿真阵型及条件

传声器阵列阵型的配置在解决声成像问题中起着重要的作用。本文的研究重点在于低频相干声源的声成像性能。对于固定尺寸的球形阵列，其定位极限也随之确定，这种限制尤其在低频段影响更为明显。增加阵列尺寸能够提高阵列的空间分辨率，但是在实际操作中随着频率降低不断增大阵列孔径是不现实的，而且大孔径球形阵列的系统过于庞大，在制作、运输和使用等方面都有很大困难。本文采用了一种分布式的开放球形阵列，子阵孔径较小适用于中高频范围，外部整体阵列孔径较大可以扩大频率范围到低频。具体的布放方案为：为实现多声源声成像的目的，选取截断阶数P为4。传声器的分布方案要求符合奈奎斯特采样定理的空间球谐形式Q≥(P+1)2，从而保证空间采样不混叠，因此选取传声器数量Q= 32。32 阵元整体分布式球形阵列示意图见图2(a)，它由8 个子阵组成，子阵的球心位于正六面体的8 个顶点上，正六面体的外接球半径为R= 0.6 m。子阵如图2(b)所示，阵元位于正四面体的4 个顶点上，正四面体的外接球半径为r=0.13 m。

图2 阵列结构示意图Fig.2 Geometry model of the array

仿真阵型为上述分布式开放球形阵列，仿真条件为空间中两个等能量白噪声声源分别位于(θ1=−135°,φ1= 75°)和(θ2= 60°,φ2= 120°)，其中θ为方位角，φ为俯仰角；空气中声速c= 343 m/s，采样率为16 kHz，帧长为640，在时域上添加一定信噪比(Signal-noise ratio, SNR)的高斯白噪声为仿真噪声。格点采集于方位角−180° ～180°、俯仰角0° ～180°的区间内，步长均为3°。仿真对比研究了6 种声成像算法，包括常规频域波束形成(Conventional frequency domain beamforming, CFDBF)[27]、函数波束形成(Functional beamforming, FUNBF)[28]、增强超分辨率空间声源相干CLEAN(Enhanced high resolution-CLEAN-source coherence, EHR-CLEAN-SC)[29]算法、L1-GIBF、球谐系数估计法以及SHD-L2-GIBF。

2.2 相干声源声成像

两个等能量相干白噪声声源， 在信噪比15 dB 条件下，300 Hz 处的6 种算法谱图结果见图3。可见在两相干声源存在的条件下，低频段处理时CFDBF 及以CFDBF 为基础的两种算法即FUNBF 和EHR-CLEAN-SC 的谱图具有一定的相似性。CFDBF 由于其波束主瓣宽导致声成像结果有虚峰且谱峰位置不准确。而FUNBF 虽然经过参数幂次处理后旁瓣降低，但仍然存留着CFDBF 的谱峰不准确问题。EHR-CLEAN-SC 算法是基于旁瓣与主瓣相干的假设从而降低旁瓣，因此在低频段相干声源存在时性能大大降低。

图3 在SNR=15 dB、300 Hz 仿真条件下相干声源声成像归一化谱图Fig.3 Normalized spectrogram of coherent sound source acoustic imaging under SNR = 15 dB,300 Hz simulation conditions

在SNR = 15 dB、300 Hz 条件下，通过50 次Monte-Carlo实验，6 种算法对两相干声源定位误差的均值结果见表1。

表1 SNR=15 dB、300 Hz 下，6 种算法对两相干声源50 次仿真定位误差的均值Table 1 The mean of the localization error of the six algorithms for 50 times simulation of two coherent sound sources at f =300 Hz, SNR = 15 dB

对于低频相干可成像的3 种算法进行进一步计算误差方差可得到表2。

表2 SNR=15 dB、300 Hz 下，3 种算法对两相干声源50 次仿真定位误差的方差Table 2 The variance of the localization error of the three algorithms for 50 times simulation of two coherent sound sources at f =300 Hz, SNR=15 dB

2.3 阵元位置误差

由于在实际应用中，阵元位置误差可能会引起一定程度上的算法性能损失。因此仿真引入每个阵元存在均值0、方差为1 cm的随机位置误差，再进行上述3 种低频相干声源声成像算法50 次仿真计算可以得到定位误差均值和方差，如表3和表4所示。

表3 SNR=15 dB、300 Hz 下，在随机阵元位置误差下，对两相干声源50 次仿真定位误差的均值Table 3 The mean of the localization error of the three algorithms for 50 times simulation of two coherent sound sources at f =300 Hz,SNR=15 dB with random microphone position error

表4 SNR = 15 dB、300 Hz 下，在随机阵元位置误差下，对两相干声源50 次仿真定位误差的方差Table 4 The variance of the localization error of the three algorithms for 50 times simulation of two coherent sound sources at f = 300 Hz, SNR = 15 dB with random microphone position error

通过上述分析得到结论，在此分布式开放球形阵列阵型配置下，低频相干声源时CFDBF、FUNBF、EHR-CLEAN-SC算法由于在低频导向矢量的相干性较强，波束的主瓣较宽，导致声源位置估计误差较大，无法达到声成像要求。SHD-L2-GIBF、L1-GIBF 和球谐系数估计法有较好的性能。其中L1-GIBF 算法由于要进行最优化迭代，拥有较高分辨率，可用于区分距离较近的声源，但是计算复杂度较高，且算法稳定性较差。球谐系数估计法由于模态展开的频率无关性，可以直接实现频域聚焦，拥有较好的定位性能。SHD-L2-GIBF 的声成像由于球谐分解预处理的截断作用，高频的误差及噪声对算法的影响效果减弱，为算法提供了数值稳定性。通过多次实验得到误差均值及方差的表格对比，可以定量得到SHD-L2-GIBF 在此条件下比其他算法定位误差高，即算法的准确度较高。在阵元位置误差存在下，SHD-L2-GIBF 的误差均值及方差最小，即算法的鲁棒性较好；并且相较最优化迭代一范数方法，此算法步骤计算复杂度较低，物理意义明确，可以量化地衡量声场，为后续实现实时处理研究提供较大可能性。

2.4 半消声室实验结果

四阵元单球阵的实物图见图4(a)，由于实际实验中硬件设备只有一个四阵元阵，无法实现直接测量32 阵元数据，因此实验时分别对8 个位置测8 次四阵元数据，并在每次测试时，于建立的坐标系原点处放置参考传声器，通过互相关确定时延从而进行32 阵元采集数据的时间同步，得到完整实验数据。

在中国科学院声学研究所半消声室进行实验，实验搭建场景在图4(b)中展示。用两个索尼SRS-XB01 音箱作为声源，同时发出等能量的相同白噪声并用四阵元子阵进行记录，每次子阵阵列的位置示意图如图2(a)中蓝色星号所示。两个音箱分别位于(θ1=−135°, φ1= 75°)和(θ2= 60°,φ2= 120°)处。在300 Hz 处各个算法声成像结果见图5。

图4 子阵及实验场景实物图Fig.4 The actual picture of sub-array and experimental scene

图5 300 Hz 处相干声源声成像半消声室实验结果归一化谱图Fig.5 Normalized spectra of experimental results in a semi-anechoic chamber of coherent sound source acoustic imaging at 300 Hz

根据声成像谱图计算其中3 种可以定位算法在半消声室实验中的误差结果，见表5。

表5 在半消声室，300 Hz 处两相干白噪声声源3 种算法的定位误差Table 5 Localization errors of three algorithms for coherent white noise sources at 300 Hz in the semi-0033 cxanechoic room

与2.1 节仿真结果对比，可见实验结果与仿真结果类似，说明阵型和算法具有较高的可行性，且SHD-L2-GIBF 算法在此阵列基础上准确性和鲁棒性较好。

3 结论

本文提出了一种基于球谐分解的二范数广义逆波束形成(SHD-L2-GIBF)算法，通过在分布式开放球形阵列阵元设置下仿真多种声成像算法，对比分析可得到此算法可以很好地适配此阵型，在低频相干声源的声成像中提供较好的空间定位能力和算法稳定性。在半消声室中搭建的实验进一步证明了此算法的可行性。