分布式传声器阵列的低频宽带信号方位估计

霍茹静，魏明洋，许志勇

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

0 引 言

随着国内民用车辆数量的增加，车辆私自改装、大型卡车超载、车辆老化不报废等问题日益突出，其带来的噪声污染、妨碍安全行车等问题亟待解决。这些问题车辆在道路上行驶时都会产生异常声音，所以现阶段的解决方案是通过异常声检测[1]与方位估计(Direction of Arrival, DOA)算法[2]确定问题车辆的位置和违法情况。由于异常声源往往处于移动状态，对位置信息的确定要求具有实时性；同时通过对这类异常声进行分析统计发现，异常声的声源级高、传播距离远、声信号的能量集中于低频段(500～1 000 Hz)且带宽较宽[3]，但是低频信号的高精度DOA估计需要使用大尺寸传声器阵列对声音进行采集。在实际应用中，需要结合上述特点选择合适的DOA算法与传声器阵型，确保DOA估计的高精度和低运算复杂度。

20世纪70年代末，多重信号分类(Multiple Signal Classification, MUSIC)[4]算法的出现开启了超分辨DOA估计算法的篇章。该算法通过谱峰搜索得到目标方位。当对精度要求较高时，谱峰搜索步长很小，这会导致算法的运算量很大。旋转不变子空间技术(Estimation of Signal Parameters via Rotational Invariance Techniques, ESPRIT)[5-6]算法相比MUSIC算法无需进行谱峰搜索，计算速度得到较大提高且DOA估计误差小，因而得到了广泛应用。

分布式阵列[7-8]是由空间散布的多个子阵列构成的阵列系统，其灵活的空间布置特性使得分布式阵列广泛应用于复杂安装环境的场景，子阵的分开布放使得分布式阵列具有很大的阵列尺寸，从而具有很好的低频信号DOA估计性能。分布式阵列的阵元总数多，阵列具有很强的鲁棒性以及很高的阵增益。这些优势使得基于分布式阵列的DOA估计算法研究具有重要意义。但是大尺寸阵列的阵元间距大于信号半波长，不满足空间采样定理，在对信号进行DOA估计时会同时出现多个估计值，即测角模糊。为了得到无模糊且高精度的DOA估计结果，需要针对阵型采取解模糊[9-11]策略。

本文提出了一种用于公路异常声监测系统的分布式阵列DOA估计算法，利用ESPRIT算法无需进行谱峰搜索的特点降低了计算复杂度。选取三个均匀矩形阵(Uniform Rectangular Array, URA)呈三角形布放在公路监控杆上，阵列整体孔径大且具有空间平移不变性。该算法利用宽带聚焦[12-13]方式将经典ESPRIT算法的适用范围从窄带信号拓展至宽带信号，并利用阵列两种尺度下的方位信息旋转不变特性估计声源方位，再根据最小误差准则解模糊法实现公路异常声源的高精度DOA估计。本文分析了分布式阵列较单个阵列对低频信号DOA估计精度的提升，研究了子阵间基线距离对DOA估计精度的影响。仿真结果验证了本文算法的有效性。

1 分布式阵列数学模型

如图1所示，子阵呈三角形分布的传声器阵列摆放在监控杆上，假设阵列分布在yOz平面上，坐标原点O位于子阵S1的第一个传声器处，子阵S1，S2，S3为孔径相等的URA，均包含M2个传声器，单个子阵的传声器间距d≤λ/2，λ为声源信号分析频段中最大频率对应的波长。在y轴方向上两个子阵的基线距离为Dy1，且满足Dy1≫Md，在z轴方向上的子阵S3相对于子阵S1的基线距离在y轴和z轴方向上的投影分别为Dy2和Dz，且满足Dy2≫Md,Dz≫Md。

图1 位于交通监控杆上的分布式传声器阵列示意图Fig.1 Schematic diagram of the distributed microphone array on traffic monitoring pole

假设一个平面波宽带声源入射到该分布式阵列上，记来波方向在xOy平面上的投影与x轴正半轴的夹角θ为方位角，来波方向与z轴正半轴的夹角φ为俯仰角，y轴与z轴上的方向余弦分别为u=sinθsinφ，v=cosφ，方位角与俯仰角的取值范围分别为θ∈[0，π2]，φ∈[0，π2)。

由于阵列接收的公路声源信号通常是非平稳的，因此需要先对接收信号进行分帧，然后逐帧变换到频率域。定义信号能量主要分布频段的带宽为B，将该频段划分为J个频点，可以得到第i帧的宽带信号模型为

式中：Xi(fj)、Si(fj)、Ni(fj)分别为阵列接收数据、声源信号、噪声在fj频点处的频谱向量，j=1，2，…，J。Ai(fj)为阵列流形矩阵，本文假定声源个数为1，则矩阵退化为列向量a(u，v，fj)。令子阵S1的第一行与第一列的导向矢量分别为

则子阵S1，S2，S3的导向矢量分别为

可得：

式中：⊗表示克罗内克(Kronecker)积。

在第i帧的频点fj下，X(fj)的协方差矩阵为

其中，信号能量σS2(fj)=Si(fj)Si*(fj)，Rn，i(fj)是噪声协方差矩阵。假设声源信号在连续2K+1帧时间内是广义平稳的，则第i帧协方差矩阵可用前后K帧信号协方差矩阵的统计平均来估计，即：

2 分布式二维DOA估计算法

2.1 宽带聚焦算法

公路车辆异常声信号为宽带信号，在利用ESPRIT算法进行DOA估计时，需要先将频带内各频点的信号子空间聚焦到参考频点，得到基于参考频点的样本协方差矩阵，进而根据子空间估计波达方向。本文采用双边相关变换(Two-sided Correlation Transform, TCT)聚焦算法[14]，该算法在理想条件下满足：

而实际情况中，样本协方差矩阵中包含噪声，聚焦变换矩阵要满足归一化约束，即要满足：

其中：T(fj)=U(f0)UH(fj)，U(f0)和U(fj)分别是和的最大奇异值对应的奇异向量。

2.2 二维经典ESPRIT算法

应用ESPRIT算法的核心思想是构造两个相同的子阵，一子阵的导向向量乘以关于方位信息的旋转因子可以得到另一子阵的导向向量(即空间平移不变)，二维ESPRIT算法要求在两个方向上构造两组具有空间平移不变性的子阵。

分布式均匀矩形阵布放示意图如图2所示。在图2中，在y轴方向选取子阵S1左边M-1列传声器构成子阵SL，右边M-1列传声器构成子阵SR，SL与SR之间的位移量为d。

图2 分布式均匀矩形阵布放示意图Fig.2 Layout diagram of the distributed URA

则阵列SL与SR的导向矢量满足关系

其中，JSL与JSR分别表示阵列SL与SR的选择矩阵，JSL=I2×3⊗(IM⊗[IM-10(M-1)×1])，I为单位矩阵，JSR=I2×3⊗(IM⊗[0(M-1)×1IM-1])，Φ=ej2πduλ表示阵列SL与SR的平移关系，其中波长λ=c/f0，c为声速，f0为聚焦频点。

聚焦后的信号协方差矩阵R(f0)进行特征值分解：

其中，ΣS和ΣN分别表示R(f0)的最大特征值和小特征值组成的对角阵，最大特征值对应的特征向量为信号子空间US。由于导向矢量a(u，v)与子空间JSLUS以及子空间JSRUS张成的空间相同，即span{a(u，v)}=span{JSLUS}=span{JSRUS}，span{}·表示生成向量张成的空间，则信号子空间的旋转不变特性可以表示为

根据最小二乘法得到旋转矩阵：

由于信号源个数为1，Ψ矩阵退化为标量Ψ。定义空间角频率β=angle(Ψ)，angle(·)表示求相位角。则y轴的方向余弦估计结果为

同理在z轴方向，选取子阵S1和S3下方M-1行传声器构成子阵SD，上方M-1行构成子阵SU，SD与SU之间的位移量为d，可以求得z轴的方向余弦估计结果v。

2.3 二维双尺度ESPRIT算法

采用图2中的分布式传声器阵列，利用子阵内相邻阵元的间距作为阵列平移不变尺度可以得到精度较低但无模糊的方向余弦估计值，称为粗估计；利用两两子阵间的间距作为阵列平移不变尺度可以得到精度高但有模糊的估计值，称为精估计。双尺度ESPRIT算法先对阵列接收信号进行粗估计和精估计，再依据最小误差准则并结合粗估计结果对精估计结果解模糊，得到最终DOA估计结果。

2.3.1 粗估计与精估计

根据定义，由2.2节所得的方向余弦估计值u和v可记为方向余弦粗估计uc和vc。类似地，精估计的空间平移不变性可表示为

信号子空间的旋转不变特性可表示为

其中，JSf1与JSf2表示精估计的选择矩阵，Ψf为精估计的平移关系。由于信号源个数为1，Ψf退化为标量Ψf，精估计相位为βf=angle(Ψf)。

定义子阵S1，S2，S3的选择矩阵分别为JS1= [IMM,02MM×MM]，JS2= [0MM,IMM,0MM]，JS3= [0MM,0MM,IMM]，0表示元素全为0的矩阵，则利用两两子阵进行精估计的选择矩阵为(JSf1，JSf2)=(JS1，JS2)，(JSf1，JSf2)=(JS1，JS3)，(JSf1，JSf2)=(JS2，JS3)，对 应 的Φf分 别 为Φ12=ej2πDy1uλ，Φ13=ej2π(Dy2u+Dzv)λ，Φ23=ej2π((Dy1+Dy2)u+Dzv)λ，将精估计选择矩阵依次代入式(21)求得角频率βf12、βf13、βf23。

2.3.2 最小误差准则解模糊

由于Dy1、Dy2、Dz均远大于d，在估计相位βf时会产生以2π为周期的测角模糊，则两两子阵间角频率βf12、βf13、βf23满足关系：

其中：n1，n2，n3∈Z，uf与vf为存在模糊的方向余弦精估计。联立式(22)中前两式即可求解无模糊的方向余弦精估计udf与vdf。

由图3所示在平面(uf，vf)上，式(22(a))表示平行于vf轴的均匀直线簇，(22(b))表示斜率相等而截距不同的均匀直线簇，方向余弦精估计为单位圆内两直线簇的交点，粗估计结果在该平面上表示参考点(uc，vc)。对方向余弦精估计解模糊，即根据最小误差准则，找到距参考点最近的直线交点。为实现正确解模糊，参考点应满足：

图3 解模糊法示意图Fig.3 Schematic diagram of defuzzification method

若直接按式(23)的约束遍历所有直线交点，计算复杂度较高。为此，本文提出一种结合几何意义的解模糊方法，首先确定周期模糊数n1，即找到距参考点最近的斜率为0的直线L1：

此时，参考点(uc，vc)在直线L1上的映射为(udf，vc)。

然后，确定周期模糊数n2，即找到距映射点最近的斜率不为0的直线L2：

最后得到声源DOA估计的方位角φDOA和俯仰角θDOA：

2.4 算法流程

综上所述，基于双尺度ESPRIT的宽带声源DOA估计算法流程总结如下：

(1) 将信号频带内各频点处的信号子空间聚焦到参考频点f0，并对聚焦后的协方差矩阵R(f0)进行特征值分解，得到信号子空间。

(2) 根据式(16)分别构造沿y轴和z轴粗估计的旋转不变性等式，求出方向余弦粗估计uc和vc。

(3) 对于S1，S2，S3中的两两子阵，根据式(20)构造精估计的旋转不变性等式，得到存在周期模糊的相位βf12、βf13、βf23。

(4) 由式(22)结合几何意义，根据最小误差准则求得无模糊的方向余弦精估计udf和vdf。

(5) 由式(28)求出异常声源DOA估计的方位角和俯仰角φDOA和θDOA。

3 仿真分析

设M=4，三个子阵内传声器间距d=0.01 m，分布式阵列传声器总数为3M2=48。子阵间基线距离满足0.04 m＜Dy1+Dy2＜4 m，0.04 m＜Dz＜3.5 m。模拟车辆异常声宽带信号源的入射方向(θ，φ)=(30°，60°)，信噪比RSN=10 dB。采用联合角度估计与方向余弦的均方根误差(Root Mean Square Error,RMSE)作为衡量算法的估计精度性能，定义为

其中，K表示蒙特卡洛(Monte Carlo)实验次数本文中K=5 000，与为第k次实验信号的估计值。

3.1 分布式DOA估计的性能比较与验证

本次实验设置分布式阵列的基线距离分别为Dy2=0.3 m，Dy1=1.5 m，Dz=1.5 m，此外，设置阵元数和阵元间距与分布式阵列相同的6行8列的URA作为对比阵列。图4给出了两个阵列基于ESPRIT算法的联合DOA估计误差在不同频率范围(1/3倍频程带宽)的变化。图4中的RMSE根据式(29)计算，横坐标为1/3倍频程带宽对应的中心频率。

图4 各频带下分布式阵列的联合DOA估计均方根误差Fig.4 RMSE of joint DOA estimation of distributed array in each frequency band

由图4可以看出，DOA估计误差随信号频率的增大而减小，分布式阵列采用双尺度ESPRIT算法的估计精度比URA提高了约10倍，而仅采用粗估计尺度经典ESPRIT算法的分布式阵列估计精度比URA差，这是由于粗估计只利用了(24，24)个阵元估计u和v，而URA使用了(42，40)个阵元。该实验验证了本文提出算法以及通过增大阵列孔径的方式提高DOA估计精度的可行性。

3.2 估计误差随基线距离变化情况

下面探究分布式阵列ESPRIT算法的估计精度与基线距离的关系。仿真使用的低频宽带声源信号的带宽为562～708 Hz。

(1) 仿真假设Dy2=0.3 m，图5给出了方向余弦u的估计误差随Dy1/d的变化情况。图5中的RMSE根据式(30)计算。由图5可以看出，随着Dy1/d的增大，在一定范围内，双尺度ESPRIT算法的方向余弦u的估计精度逐渐提升，但当d＞1.9 m、Dy1/d＞190时，估计精度变差，此时的基线距离称为基线模糊门限。经典ESPRIT算法只利用了相邻阵元间距，基线距离Dy1对方向余弦u的估计精度影响不大。

图5 方向余弦u估计精度与Dy1/d的关系Fig.5 Relationship between accuracy of direction cosine estimation u and baseline distance Dy1

(2) 令Dy1=1.9 m，图6给出了在不同基线距离Dy2下，方向余弦v的估计精度随基线距离Dz的变化情况。图6中的RMSE根据式(31)计算。由图6可以看出，在基线模糊门限范围内，Dy2=0时的估计精度高于Dy2=0.3 m时的估计精度。

图6 方向余弦v估计精度与基线距离Dz/d的关系Fig.6 Relationship between accuracy of direction cosine estimation v and baseline distance Dz

由于解模糊以粗估计值为参考，当基线距离Dy1与Dz增大到一定程度，图3中的直线簇变密集，当粗估计值无法满足式(23)中的条件时，会导致解模糊错误率升高，DOA估计误差增大。当Dy2=0，图3中的直线簇相互垂直，式(27)中v不受u的误差影响，此时的DOA估计精度最高。在实际应用中，可根据公路监控杆的条件设计分布式阵型。

4 结 论

为解决公路异常声源的高精度DOA估计问题，并满足低计算复杂度、布阵灵活的需求，本文提出了一种适用于均匀矩形子阵呈三角形分布的分布式传声器阵列二维DOA估计算法。本文对信号采用宽频段聚焦和ESPRIT算法估计DOA，避免了逐频点估计和二维谱峰搜索造成的复杂计算，算法的双尺度策略与最小误差准则解模糊法能利用扩展孔径的优势，提高低频声源的DOA估计精度，且适用于任意三个相同均匀矩形子阵组成的分布式阵列，具有较好的工程应用前景。