张国昌 吴 鸣 韩欣宇 武帅兵 殷 兰 杨 军

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

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

(3北京信息科技大学自动化学院 北京 100192)

0 引言

近年来，使用麦克风阵列的声源方位估计方法引起了学者的广泛关注，其在机器人、电话会议、视频监控、助听器等[1-4]领域有着诸多的应用。通常声源波达方向(Direction of arrival,DOA)估计方法分为间接法和直接法两类[5-7]。间接法，该方法需要首先计算声源到各阵元对的相对时间延迟(Time difference of arrival,TDOA)而后使用几何方法确定声源方位[8]；直接法，其计算一系列候选方位上的损失函数，并选择最有可能的方位，其根据使用的技术不同分为可控响应功率(Steered response power,SRP)法[9]、最大似然方法[10]和子空间类方法[11]。

其中TDOA方法虽然计算复杂度较低，但是由于其只利用了各通道信号互相关函数最大值所对应的时延信息，使得其极易受混响的影响。

由广义互相关函数(Generalized cross correlation,GCC)计算SRP的DOA估计方法被广泛的应用分析，其中基于相位变换(Phase transform,PHAT)广义互相关函数的SRP算法(Steered response power-phase transform,SRP-PHAT)[12-15]最为流行。但是这种方法的空间分辨率较差，另外还需要进行逐个频点的空间谱扫描，计算复杂度高。

最大似然方法以其良好的一致性(Consistency)、极佳的统计特性(Statistical efficiency)和处理相干信号的能力被广泛的研究[10]，但需要进行高维空间搜索来求取最优解，计算复杂度高，目前很难实际应用。

另外一种基于子空间的方法以其超分辨特性也受到了广泛的关注，这类方法在多声源或存在混响的情况下具有优于TDOA和SRP方法的性能[11]。针对宽带信号，两种基于窄带分解的子空间方法得到广泛的研究。非相干信号子空间方法(Incoherent signal subspace method,ISSM)[8]融合各个子带的DOA估计结果；相干子空间方法(Coherent signal subspace method,CSSM)[16-19]先将各窄带信号聚焦至某一特定频点，而后只在聚焦频点处使用窄带DOA估计。两种方法各有优缺点：前者由于相干早期反射声的存在使得空间相关矩阵(Spatial correlation matrix,SCM)的信号子空间缺秩，这导致估计性能的急剧下降；后者利用频率聚焦解决了前者的缺秩问题，但需要对声源方位进行初步估计以获得聚焦矩阵，该方法的性能对初始DOA估计精度非常敏感[18]。

为了避免初步DOA估计和聚焦矩阵的实时求解，针对环形阵列(Circular array,CA)，环谐波(Circular harmonic,CH)分解可以解耦频率与估计方位角之间的相关性，使得聚焦可以通过各频点下SCM的平均来实现。本文后续的声源定位方法均采用CH域的CSSM。

CH域的DOA估计方法以其优良的特性被广泛的分析。Tianaroig等[20]使用环谐波波束形成(Circular harmonic beamforming,CHB)改善了延时求和波束形成器(Delay and sum beamformer,DSB)的分辨力；Teutsch等[21]提出了特征波束域的基于旋转不变技术的信号参数估计方法(Estimating signal parameter via rotational invariance techniques,ESPRIT)，避免了空间谱扫描；Torres等[22]首次使用逐时频点处理的方法来提升CHB定位的稳健性。但是CA的CH域展开系数存在零点，这会导致严重的噪声放大问题。为了解决这一问题，以下三种方法被广泛应用：将CA固定至刚性圆柱散射体上，使用吉洪诺夫正则化[23]和使用同心多环阵(Concentric circular array,CCA)[24]。与前两种方法相比，多环阵无需增加刚性散射体也不存在正则化引入的波束失真。

为了解决CH展开零点处的噪声放大问题，本文由CH域的理想波束出发，推导出了CCA频域补偿滤波器的最小模求解方法。这种方法设计的补偿滤波器具有最优的噪声增益。另外，最小模CCA环谐波域的波束响应也不存在理想CA贝塞尔函数零点处的失真和正则化CA的低频波束展宽问题。为了对比CA和最小模CCA在CH域的声源DOA估计性能，还设计了一套麦克风阵列系统。仿真和实验均表明，与使用同阵元数、同孔径的CA相比，最小模CCA可以有效地降低声源方位估计的误差。

1 同心多环阵环谐波展开

1.1 连续同心多环阵

图1展示了由P个连续同心单环阵组成的连续同心多环阵列，其中第p个环阵的半径为rp，信源s入射的方位角和俯仰角分别为θi和φi。第p个连续环阵接收到的声压信号可在极坐标下表示为

其中，S(ω)为入射信源信号在角频率ω下的分量，k=2πf/c=ω/c为波数，c为声音在空气中的传播速度。

图1 多环阵信号入射示意图Fig.1 The geometry of CCA with P rings

利用贝塞尔函数展开公式

式(2)中，Jl(α)为一类l阶贝塞尔函数，式(1)可以展开成

这 里ejlθ被 称 为l阶 环 谐 波CH(或 模 态)，Cl(krp,θi,φi)为由(θi,φi)方向入射信源信号的ω频率成分在l阶CH上的展开系数，

至此，第p个连续圆环阵列在各方位角下的声压信号被分解到无穷阶的CH上，xp(krp,θ)与Cl(krp,θi,φi)为一傅里叶变换对。

1.2 连续阵列均匀采样

实际应用时，需要对环形连续阵列进行离散采样，这里考虑均匀采样场景，这种离散的单环阵和同心多环阵列分别被称为均匀环形阵列(Uniform circular array,UCA)和均匀同心环形阵列(Uniform concentric circular array,UCCA)。采样操作引入了CH系数的误差[23,25]：

其中，为采样引起的第p环l阶CH系数误差，假设第p环的CH分解阶数为lp，那么当第p环麦克风数目Np＞|2lp|时，式(5)中的第一项占主要成分。另外CH分解阶数lp＞krpsin(φi)时的CH分解系数较小可以忽略。考虑sin(φi)取最大值1，p环的CH最高阶数假定为Lp≈krp，由于CH最高阶数为krp，所以当第p环麦克风数满足Np＞2krp时较小可以被忽略。

2 多环阵环谐波域波束形成

多环阵CH域波束形成旨在将各环各阶CH成分进行合适的补偿，以形成具有空间筛选特性的波束。其关键在于补偿系数的求解。下面将从理想波束出发，推导出UCCA的最优补偿滤波器。

理想情况下，针对入射声源s，指向θi方向的期望波束形成器的响应输出可表示为

这里δ(θ-θi)可以由无穷阶的傅里叶级数组合获得，则

CHB的输出可以表示成各环各阶环谐波分量进行合适补偿后的累加，

为了使得CHB的输出可以逼近于式(7)所示的理想波束形成器的输出，应当满足以下条件：对任意给定环数P，存在一组或者多组可使得式(10)成立：

这里ωL至ωH为设计补偿滤波器的频率范围。式(10)中的在环数P=1时有唯一确定解，在P＞1时有无穷个可行解。

考虑均匀离散采样，第p环CH分解最高阶数Lp，并假设L1=L2=···=Lp=L，理想波束近似为

将式(11)写成向量形式得

此时，BCHB(ω,θ)可以表示如下：

CH域波束BCHB(ω,θ)可以看作是CH系数矢量和CH域导向矢量的夹角余弦。

2.1 单环阵补偿滤波器求解及其正则化

当环数P=1时，UCCA退化为UCA，此时补偿滤波器Hl(ω,φi)需要满足

在频率ω对应贝塞尔函数零点时，该频点对应的Hl(ω,φi)趋于无穷大，这导致了严重的噪声放大问题。为了缓解这一问题，一种基于吉洪诺夫正则化求解补偿滤波器的方法被提出[23]

式(16)中，α为吉洪诺夫正则化系数，正则化的引入提升了CHB的稳健性，但牺牲了波束指向性。增大正则化系数可以获得更加稳定的波束形成器。

2.2 多环阵最小模补偿滤波器求解

多环的引入，松弛了式(10)对应的约束条件，使得补偿滤波器由UCA的唯一解变为UCCA的无穷个解。可以从这无穷个解中挑选出最鲁棒的解。与阵元域的加权系数的特性类似，最鲁棒解对应的补偿滤波器趋向于具有最小的模值[26]。最小模补偿滤波器可以通过求解如下最优化问题获得：

Hl(ω,φi)=为各环补偿滤波器的向量形式，

式(8)对应的最小模解为

与UCA利用标量Hl(ω,φi)进行CH系数补偿不同，UCCA使用矢量滤波器同时补偿多个环阵的同阶CH。为了补偿某一确定(l,ω,φi)条件下的CH系数至期望值，UCCA可以利用不同半径下的多个同阶CH来完成，从而具有更大的灵活性。通过求解最小模补偿滤波器，在接近贝塞尔函数零点处，UCCA会偏向于赋给CH分解系数较小的环阵以较小的补偿系数，这一点与式(16)对应的UCA补偿滤波器的解恰恰相反。通过合理的设置半径，可以避免各环CH分解系数零点处于相同频率条件下，UCCA可以彻底解决UCA中贝塞尔函数零点处的噪声放大问题。

补偿滤波器模值与CH波束形成器的稳健性是负相关的[26]。下面对UCA和UCCA在俯仰角为90°下的前两阶补偿滤波器的模值进行了仿真，其中，正则化参数配置为0.00065，两种UCA和UCCA的阵列参数如表1所示。结果如图2所示，可知，UCA低阶CH下存在模值无穷大的补偿系数(图中展示的是模值的倒数)，UCCA通过组合两个无公共贝塞尔函数零点的UCA完全解决了这一问题。另外值得注意的是，与没有进行正则化处理的两UCA相比，UCCA在任一频点下都具有更小的补偿系数模值，也就是说，UCCA在任一频率下都具有更加稳健的性能。图3对比了正则化后UCA和UCCA的波束响应，其中，CH分解阶数为3阶，单环阵半径选择为60 mm。正则化虽然有效地提升了单环阵CHB的稳健性，但并不能有效地消除贝塞尔函数零点导致的波束图失真。正则化UCA的波束响应在2180 Hz、3440 Hz、4630 Hz等频点处存在明显的畸变，另外，正则化还引起了低频波束的展宽。而UCCA则不存在上述问题。由于CH分解阶数的限制，在高频区域，DSB方法具有最窄的主瓣。

表1 阵列参数Table1 Configurations of microphone arrays

图2 两种理想UCA、正则化UCA和UCCA补偿滤波器模值与频率的关系曲线Fig.2 The compensation filters’modulus values curves are shown as a function of frequency for UCAs with no regularization,UCA with regularization,and UCCA

图3 单环DSB、CHB和同心双环CHB的波束图Fig.3 Beampattern for DSB,CHB using UCA with regularization and CHB using UCCA

3 DOA估计

3.1 逐时频区域的构建CH域SCM

与传统UCA构建SCM的方法类似[27]，UCCA的CH域SCM可以由式(19)计算得到：

3.2 使用MUSIC方法进行DOA

在每个时频区域，DOA都会被使用多重信号子空间分类(Multiple signal classification,MUSIC)的方法估计出来。在假设每个TF区域只存在一个声源的前提下，声源方位角可以通过搜索空间谱最大值来获得：

这里A(θ)为方位角θ对应的导向矢量，UN∈C(2L+1)×2L为由SCM最小的2L个特征值对应的特征矢量组成的噪声子空间，Θ为待扫描的方位区间。

4 仿真及实验

4.1 仿真

本小节，为了验证UCCA具有更优的定位鲁棒性，该文评估了不同信噪比和混响条件下的DOA估计的统计性能。评估算法均采用CH域的CSSM方法。声源语料来自于TIMIT[28]语料库。UCA半径和阵元数分别设置为60 mm和16个，UCCA配置参数如表1所示，其他的仿真参数设置如表2所示。UCA的最高CH展开阶数为7阶，UCCA的CH展开阶数受限于阵元数最少的环，最高为3阶。在接下来的性能评估实验中，麦克风阵列被放置于距离地面1 m高的房间水平面中心位置，单个播放声源的扬声器被依次放置于方位角间隔30°的6个方位，扬声器高度同为1 m，声源到麦克风阵列中心位置距离为2 m。以下仿真结果会统计平均所有声源方位。虚源法[29]被用于生成带有混响的多通道信号，多通道不相关的高斯白噪声被添加至带混响信号，加入噪声后的信号被看作是阵列采样信号。另外，需要使用语音激活检测(Voice activity detection,VAD)模块[30]对阵列采样信号进行预处理，截取语音段，舍去非语音段。语音段数据每10帧(0.16 s)获得一个平均的声源方位估计结果。每一帧数据选取离散傅里叶变换后的88个子带(500～6000 Hz)，考虑频域步进间隔为5个子带，此时每一帧数据会进行16次声源方位估计。那么10帧一共获得160个定位结果，对该160个结果进行中值平均进而估计出一个较为稳健的声源方位。仿真分析中使用的软件为Matlab 2018b，处理器为Intel Core i7-6700HQ，主频2.6 GHz，处理器的核心数为4。

表2 仿真参数设置Table2 Parameters of simulations

首先仿真不同信噪比条件下的UCA和UCCA的DOA估计的统计性能，将房间的混响时间(T60)设置为0.5 s，信噪比依次设置为5 dB、10 dB、15 dB和20 dB。仿真结果如表3所示，其中最后一行的实时计算系数的计算公式为实时系数等于声源定位算法的执行时长除以VAD检测后的数据时长。

由表3可知，UCA在CH展开阶数为3阶或4阶的时候取得了较优的平均绝对离差性能。过低或者过高的分解阶数均导致了鲁棒性的退化，这是因为过低的分解阶数会使得过多的空间信息被舍弃，过高的分解阶数又会引起较差的白噪声增益。UCCA在阶数为3阶的时候具有较优的鲁棒性。对比UCA和UCCA，仅仅在高信噪比(20 dB)CH展开阶数为4阶这一个条件下，UCA具有优于UCCA的性能。在该条件下，因为UCCA构建的CH域空间相关矩阵维度较低，所以具有更低的计算复杂度。

接下来仿真不同混响时间下的DOA估计性能，设定信噪比为15 dB，混响时间设置范围为0.2～1.4 s，绝大部分房间的混响时间也均处在这一范围。表4展示了不同混响时间下的声源定位的统计性能。由表4可知，当T60＞0.2 s时，3阶CH展开下的UCCA具有最优的鲁棒性，当T60为0.2 s时，3阶CH展开的UCCA与4阶CH展开的UCA具有接近的性能。另外，与UCA相比，相同CH展开阶数的UCCA的鲁棒性具有显著的优势。这也就是说，在计算复杂度相同的情况下，UCCA性能优势明显。

表3 不同信噪比和阵列参数下的DOA估计平均绝对离差Table3 Mean absolute deviation performance for various SNRs and different array configurations(单位：°)

表4 不同混响时间和阵列参数下的DOA估计平均绝对离差Table4 Mean absolute deviation performance for various T60s and different array configurations(单位：°)

值得注意的是，目前UCCA各环的半径设置只遵循以下准则：没有同阶CH展开系数的零点处在相同频率点下。这一设计准则不是最优的，有理由认为，存在某种最优准则下的环半径设计方法可以使得UCCA具有更优的声源定位性能。

4.2 实验

仿真未考虑实际使用场景下的麦克风一致性、阵元位置误差等因素，因此设计一套UCCA用于评估其在真实环境下的定位鲁棒性具有重要意义。采用楼氏电子微机电系统(MEMS)模拟麦克风设计一套参数如表1所示的UCCA麦克风阵列硬件，如图4所示。麦克风型号为SPH1642，其幅度灵敏度误差为±1 dB，信噪比为65 dB。使用该阵列，在一个经过装修的房间中进行实验。房间尺寸为5.8 m×4.6 m×2.4 m。实验时扬声器和麦克风阵列布置方法如图5所示，每次实验只有一个扬声器发声音。其中扬声器和麦克风阵列的放置高度均为距离地面1.2 m，扬声器距离麦克风阵列中心1.7 m。

图4 含有16个麦克风的均匀同心双环阵实物图，半径分别为60 mm、40 mmFig.4 Photograph of 16-element uniform concentric circular microphone array with radius of[60 mm,40 mm]

图5 实验配置Fig.5 Experimental setups

在实验开始前，首先将扬声器放置于图5所示的0°方向测量了房间的脉冲响应，并使用脉冲反向积分方法估计了房间的T60，房间实物图和房间脉冲响应如图6所示，脉冲反向积分法[31]估计出T60约为0.377 s。

图6 实验房间实物图及其对应的房间脉冲响应Fig.6 Photograph of the experiment room and room impulse response of this room

接下来选用与仿真相同的参数和语料数据进行了实验。实验时依次将扬声器放置在图5所示的7个位置，每个方位播放30 min的语料。麦克风阵列采样信号首先进行VAD处理，处理后的数据每0.16 s获得一次中值平均的定位结果。表5展示了不同CH展开阶数和不同声源方位下DOA估计的统计性能。表格中的UCA对应为UCCA的外环。

表5 不同声源位置和阵列参数下的DOA估计平均绝对离差Table5 Mean absolute deviation performance for various source localizations and different array configurations(单位：°)

与仿真结果类似，在CH展开阶数相同时，UCCA的定位鲁棒性在所有实验布放的声源位置下均具有明显优势。分析CH展开阶数，在大部分的声源方位下，UCA在CH展开阶数为4阶时定位性能较优，但与3阶CH展开的UCCA相比，仍有不小的差距。

综合上述的仿真和实验结果来看，UCCA在以下条件下均具有显著优于 UCA的定位性能：(1)相同CH展开阶数，上述所有的T60和SNR条件下；(2)强混响或者低信噪比时，任意的CH展开阶数下。在高噪声和弱混响的部分CH展开阶数下，UCA与UCCA具有接近的DOA估计性能。

5 结论

本文阐述了一种针对多环阵列的环谐波展开技术。环谐波展开系数经补偿滤波器补偿后被用于声源方位估计。为了提升CH域方位估计的稳健性，针对UCA贝塞尔函数零点处噪声放大问题，提出了一种最小模准则设计UCCA补偿滤波器的方法。针对混响环境，利用环谐波展开解耦频率与方位角相关性的特点，CSSM被用于UCA和UCCA。仿真结果表明，合理设置孔径的最小模UCCA可以解决UCA零点问题，同时也避免了正则化引入的低频波束展宽。统计性能的仿真和实验显示，在相同麦克风数目和阵列孔径的前提下，与UCA相比，使用UCCA可以显著提升混响环境下的DOA稳健性。