强混响环境下基于声矢量传感器的多源DOA估计

孙思雨，张海剑，孙 洪

(武汉大学 电子信息学院，湖北 武汉 430072)

0 引言

声矢量传感器(Acoustic Vector Sensor,AVS)是一种由2～3个正交的方向传感器和一个全向麦克风组成的阵列,方向传感器用来测量沿正交轴的粒子速度,全向麦克风用来测量该点的声压[1-3]。AVS中的各传感器元件十分靠近,故可假定它们被放置在空间中的同一点上,从而使得AVS具有频率无关的特性。利用AVS所提供的丰富声学信息以及其频率无关的特性,AVS在阵列信号处理方面优于传统的传感器阵列。此外,由于具有体积小、质量轻等特点,AVS在声学研究领域也受到不断关注,并且在一些新兴行业中得到实际应用,如目标跟踪[4-6]、水声通信[7]和声源定位[8]等。

波达方向(Direction of Arrival,DOA)估计是一个备受关注的研究热点[9-13],随着AVS在阵列信号处理领域的兴起,一些基于AVS的DOA估计方法也纷沓而至。一种常用的方法是将基于传统压力传感器阵列的DOA估计方法扩展到AVS阵列,如波束形成[14-18]、多重信号分类(Multiple Signal Classification,MUSIC)[19]、旋转不变参数估计(Estimation of Signal Parameters via Rotational Invariance Technique,ESPRIT)[20]以及基于稀疏表示的方法[21]。此外,利用AVS特有的阵列结构,过去三十年中许多学者提出了适用于AVS的声源测向算法。Nehorai等[22]提出了一种基于强度矢量和速度协方差矩阵的方法,用单个AVS处理DOA估计问题,这是AVS应用在DOA估计中的开创性工作。此外,Levin等[23]还在其中推导了均方角误差的克拉美罗下界用来评估DOA估计的性能。文献[24]中提出了基于单个AVS的最大导向响应能量估计器,根据最大功率方向来估计DOA,该方法是基于强度矢量和速度协方差矩阵的扩展,能够利用合适的参数来接近克拉美劳下界。然而,上述工作主要是为单源测向而设计的,不适用于多源存在的欠定情况。

上述策略大多是在时域内处理观测信号,而不是充分利用时频特征,仅利用信号的时域统计信息,难以处理源数超过麦克风数目的欠定情况。随着时频分析在信号处理领域的兴起[25-26],时频稀疏性理论逐渐应用于基于AVS的多源DOA估计中。通过检测单个源占主导的时频点,并利用这些时频点进行DOA估计,多源DOA估计任务可以转换为多个单源DOA估计任务[8]。

准确估计DOA的关键在于正确检测出具有特定特性的有利于测向的时频点。文献[27]中提出了一种基于传感器间数据比(Inner-Sensor Data Ratio,ISDR)的方法,通过提取局域信噪比(Signal to Noise Ratio,SNR)较高的时频点来实现多源DOA估计[28]。基于ISDR的方法在自由场假设的前提下取得了较好DOA估计性能,但随着混响的增加,反射信号会破坏ISDR,从而导致该方法整体性能下降。为提高在混响环境下的DOA估计性能,一般选择受混响干扰较少的时频点从而实现精确DOA估计。文献[29]中提出了一种基于低混响单源点(Low-Reverberant-Single-Source Point,LRSSP)的方法,利用单个AVS在混响环境中实现多源DOA估计,其中LRSSP指代时频域中受混响干扰较小的单源时频点。针对AVS特有的频率无关特性,该方法利用每个LRSSP的实部和虚部具有相似方向向量的特点,用来在时频域中检测LRSSP。尽管基于LRSSP方法具有较好性能且易于实现,但随着混响增加,该方法不可避免地会检测到一些包含错误方向信息的非LRSSP从而导致性能下降。

为解决上述估计问题,进一步提高强混响情况下DOA估计的准确性和鲁棒性,本文提出了一种有效的方法用来检测直达声单源点(Direct-Path-Single-Source Point,DP-SSP),指代由单个实际源信号而非反射源所占主导的时频点。DP-SSP中包含了相关源的正确方向信息,能防止估计的方向与真实DOA产生较大偏差。本文受文献[30]直接路径优势(Direct Path Dominance,DPD)测试的启发,提出的DP-SSP检测技术首先使用频率平滑来减轻反射信号的影响,然后利用子空间投影技术来进一步提取DP-SSP并消除干扰项。检测出的DP-SSP随后被用来进行源数目估计和聚类。最后将MUSIC方法用于每一类DP-SSP点簇中用来计算该点簇对应源的DOA。下文将详细介绍本文所提出的基于频率平滑和子空间分解的DP-SSP检测方法。

1 信号模型

假设一个密闭混响室内中存在M个声源和一个由3个方向传感器和一个全向麦克风组成的AVS阵列,该AVS阵列接收到的混合信号经过短时傅里叶变换(Short-Time Fourier Transform,STFT),得到时频域信号可表示为:

(1)

式中:t和f表示分别表示时间轴和频率轴索引,X(t,f)=[Xp(t,f),Xvx(t,f),Xvy(t,f),Xvz(t,f)]T中包含了AVS四个通道信号的STFT系数,Xp(t,f)、Xvx(t,f)、Xvy(t,f)、Xvz(t,f)分别表示全向麦克风和3个方向麦克风的输出,Hq(f)表示第q个源到AVS之间的声传输函数,Sq(t,f)表示第q个源的STFT系数,N(t,f)表示噪声的STFT系数。

由于混响环境下的观测信号一般由直接路径声和反射声组成,因此第q个源到AVS之间的声传输函数可以写为:

(2)

图1 信号反射模型Fig.1 Signal reflection model

这里Iq表示强反射信号等效产生的虚拟镜像源的个数,ai(f)表示第i个源到AVS之间的声传递函数,Ψi=(θi,φi)表示第i个源的DOA,则基于AVS的阵列流形第i个源的导向矢量可以表示为[24]:

d(Ψi)=[1,cosφicosθi,sinφicosθi,sinθi]T,

(3)

以AVS所处位置为笛卡尔坐标系原点,φi∈[0°,360°)和θi∈(-90°,90°]分别表示第i个源到原点的方位角和仰角。

基于上述分析,本文的镜像模型可以等效为一个不含混响的远场模型:

(4)

式中:J表示实际源和虚拟源总数,其中包含Q个实际源和J-Q个由强反射产生的虚拟镜像源,E(t,f)表示建模误差,包含弱混响或者低能量的干扰源点。因此,在混响环境下的信号模型可以近似转化为自由空间中存在J个源的远场模型。在低混响情况下,J的值与Q接近。

2 基于DP-SSP检测的DOA估计方法

2.1 问题描述

根据式(4)中描述的信号模型,可以将时频点分为DP-SSP、强反射单源点(Strong Reflection Single Source Point,SR-SSP)以及多源时频点(Multi-Source Point,MSP)。DP-SSP是一种能量由一个实际源的直接路径信号占主导,受其他信号源和混响影响较少的时频点。SR-SSP表示由强早期反射信号等效的镜像源占主导的时频点,强早期反射信号通常指第一次反射的信号。MSP表示由多个能量值较大的源占主导的时频点,其中包含实际源和虚拟镜像源。由于SR-SSP和MSP会对最终的DOA估计结果产生不利影响,算法设计时需要考虑去除这类离群值。

从理论上来说,少量的DP-SSP已经足以获得准确的DOA估计。然而DP-SSP的检测通常会伴随着一系列的干扰时频点,这些时频点在检测过程中不可避免地会被错误识别成DP-SSP。因此,为了精确估计DOA,显然希望正确的DP-SSP点在所有检测出的时频点中占比越高越好。文献[29]提出了基于LRSSP的方法,致力于利用AVS接收到的信号存在的独特数学特性来设计算法提取DP-SSP,从而获得准确的源测向结果。然而随着混响增加,大量满足LRSSP检测规则的异常值被错误检测为DP-SSP,最终导致真实的DP-SSP在所有检测出的时频点中占比下降,进而对最终DOA估计造成影响。

为了分析基于LRSSP[29]的方法在高混响条件下识别DP-SSP的性能,给出了每个检测到的LRSSP点的强度向量分布,并计算其中真实DP-SSP的占比。检测到的LRSSP的强度向量计算方法如下[29]:

(5)

(6)

式中:{·}*表示共轭操作,uq表示第q个源的单位方向向量,当强度向量η(t,f)近似与一个源的单位方向向量平行时,则其对应的时频点(t,f)被认为是该源的一个DP-SSP点,其判别规则可以表示为:

(7)

(8)

式中:Card(·)表示一个集合中元素的数量,集合Ωd中包含了所有符合检测规则的时频点,Ωd中所有满足式(7)中的时频点包含在集合ΩrD中。在本文中,实际DP-SSP占比P被用来作为评估DP-SSP检测性能的一个标准。

本文实现一个简单的仿真来辅助说明上述分析,对比结果展示在图2中,不同颜色的圆圈表示对应于不同源的真实DP-SSP。仿真房间大小为8 m×6 m×4 m,AVS被放置在房间中心,3个声源来自于TIMIT数据集[32],均被放置于离声矢量传感器1.5 m处,3个源的DOA分别为(50°,0°)、(110°,15°)和(170°,30°)。

(a)基于LRSSP的方法(RT60=0.3 s)

(b)基于LRSSP的方法(RT60=0.6 s)图2 检测到的时频点与真实DP-SSP点对应的强度矢量分布图(SNR=20 dB)。Fig.2 The intensity vector distribution of the detected time-frequency points and real DP-SSP (SNR=20 dB)

从图2中可以明显看到,检测到的LRSSP中仅有一小部分是对DOA估计有利的真实DP-SSP点。通过对比图2(a)和图2(b)可以发现,随着混响的增加,强度向量分布的方向特征变得不明显,且实际DP-SSP的比例逐步下降。尽管基于LRSSP[29]的方法会对检测到的时频点做进一步的离群值去除使强度矢量分布具有更加明显的方向性,但是这仅适用于低混响的情况。随着混响的增加,大量出现的离群值会使得每个源对应的强度矢量簇中心与真实的方向向量发生偏移,进而在进行离群值去除时导致大量有利于DOA估计的时频点被去除,从而影响最终DOA估计的精度。为解决上述实际DP-SSP占比较低的问题,本文提出了一种有效的DP-SSP检测方法,通过准确检测DP-SSP并消除离群值,来提升有利DP-SSP的整体占比。

2.2 DP-SSP 检测

基于上述分析可以得到,在高混响环境下,DP-SSP的检测对DOA估计来说是至关重要的,因为DP-SSP中包含了对应源的关键方向信息。随着混响的增强,DP-SSP的检测变得越来越有挑战。本文提出了一种有效的两阶段策略来检测有用的DP-SSP并去除干扰时频点,在第一阶段应用频率平滑技术来尽可能去除潜在的早期反射强源点SR-SSP,在第二阶段应用子空间投影进一步有效去除多源点MSP并提取DP-SSP。

2.2.1 采用频率平滑消除SR-SSP

早期反射信号被认为是对直达声信号做延迟和衰减。随着混响增强,早期反射信号的能量逐渐增加变得不可忽略,因此第一次反射的信号通常被建模为一个与实际源信号高度相关的虚拟镜像源发出的信号。SR-SSP指代那些由虚拟源占主导的时频点。为了减少SR-SSP对最终DOA估计的影响,首先对所有的时频点做一次频率平滑以消除潜在的SR-SSP点,并粗略提取不受混响干扰的DP-SSP[33]。

受DPD测试[30]的启发,本文对分别对语谱图在时间和频率轴上平滑取平均,构建局部时频相关矩阵:

(9)

式中:Jt和Jf分别表示滑动窗口的尺寸大小。鉴于AVS的频率无关特性,导向矢量在频率平滑的过程中不会受到影响,因此上述时频相关矩阵可表示为:

(10)

(11)

(12)

式中:σi(t,f),i=1,2,…,J表示矩阵σi(t,f),i=1,2,…,J的一个特征值。如果第i个源在该时频区域内不活跃,则其对应的特征值将接近于0。将式(12)带入式(10)得到:

σJ(t,f)d(ΨJ)Td(ΨJ)。

(13)

(14)

因此,潜在DP-SSP的检测规则可以写为:

(15)

式中:σmax和σsubmax分别表示最大和第二大的特征值。

通过频率平滑和特征值分解操作去除了潜在的SR-SSP,减轻了混响的影响。由于上述操作是一个“区域级”的操作,会不可避免地检测到许多非DP-SSP点,比如多源点MSP,因此本文将第一阶段检测得到的点重新定义为低混响自项源点(Low Reverberant Auto Source Points,LR-ASP)。所有满足式(15)中检测规则的时频点被包含在集合ΩA中。为了进一步精确提取DP-SSP,将在下一小节对第一阶段检测到的LR-ASP进行“点级”的处理。

2.2.2 采用子控件分解消除MSP

由于集合ΘA中的LR-ASP主要包含DP-SSP和MSP,因此需要对LR-ASP进行提纯,在“点”级别上提取DP-SSP点。首先利用信号时频稀疏性,对每个LR-ASP做如下假设:在每个LR-ASP出现的强能量源点不超过2个。这个假设在大多数情况下是可行的,因为涉及3个或更多高能量源的时频点数量太少可忽略不计。

因此,每个LR-ASP可以写作:

(16)

式中:q1,q2∈{1,2,…,Q}表示出现在该时频点上的2个源的序号。当2个源之间能量值比值大于一个远大于1的阈值时,能量值高的源被认为主导该时频点,可以被看作是对应源的一个DP-SSP点,该阈值一般为经验选取。计算并归一化每个时频点对应的空间向量可得:

Cq1(t,f)d(Ψq1)+Cq2(t,f)d(Ψq2),

(17)

(18)

然而,当LR-ASP是一个由2个源占主导的MSP,上述条件将不会被满足。接着对所有的归一化空间向量进行k-means聚类,得到N簇向量,每个簇心对应的向量构成了一个过完备导向矢量字典。由于归一化空间向量是一个方向向量,对应的是与单个源相关的DP-SSP或者与多个源相关的MSP,因此将N设置为一个大于Q的值。每个簇心对应的方向向量计算为:

(19)

(20)

根据前述假设,每个LR-ASP上活跃的源信号不超过2个,则在时频点(t,f)上最可能出现的2个源的STFT系数可以计算为:

(21)

式中:[·]†表示Moore-Penrose伪逆操作符,aq1、aq2表示2个源对应的导向矢量。在字典ΘA中搜索,找到最优的2个导向矢量aq1、aq2,满足下式:

(22)

(23)

A2=[ar1ar2],

(24)

式中:P表示矩阵A2噪声子空间上的正交投影矩阵,ar1和ar2分别表示从字典ΘA中随机挑选的2个列向量,I表示单位矩阵,(·)-1表示求逆操作。

2个源的STFT系数可由式(21)计算得到。通过计算2个源之间的能量比,并设定阈值,即可挑选出一组DP-SSP点:

(25)

式中:λ2为一个经验阈值。所有满足上述不等式的LR-ASP点都被认为是DP-SSP点,包含在集合ΩD中。通过对集合ΩA中的所有LR-ASP点进行子空间分解,可以有效去除2个强能量源占主导的多源点MSP。

为了验证所提的DP-SSP检测方法的有效性,首先在对基于LRSSP[29]的方法应用频率平滑,去除潜在的SR-SSP点。结合LRSSP[29]检测和SR-SSP去除后得到的时频点强度矢量分布如图3所示。虽然最后检测到的DP-SSP分布的方向性仍不明显,但与原始基于LRSSP[29]检测的方法相比,真实DP-SSP所占比例有明显提高。本文所提的两阶段DP-SSP检测方法检测到的DP-SSP强度矢量分布如图3(b)所示,其分布具有明显的方向性。虽仍然存在一些方向向量偏离真实DOA的较为分散的离群点,但总体强度矢量的分布非常接近真实源方向。此外,采用所提出的DP-SSP检测方法,相比基于LRSSP[29]的方法,真实DP-SSP占比大大提高,有力地说明了该方法在检测DP-SSP上的优越性。

(a)LRSSP检测+SR-SSP去除 (RT60=0.6 s)

(b)本文所提DP-SSP检测算法 (RT60=0.6 s)图3 2种检测规则检测到的时频点与真实DP-SSP点对应的强度矢量分布图(SNR=20 dB)Fig.3 The intensity vector distribution of the detected time-frequency points and real DP-SSP using different detection rules (SNR=20 dB)

随后利用平滑直方图[29]方法对检测到到DP-SSP进行源数估计,然后根据估计出的源数目将所有检测到的DP-SSP聚类为多组。最后对每组DP-SSP利用MUSIC方法进行DOA估计,整合每组DOA估计结果得到最终多源DOA估计结果。

3 仿真实验分析

给出了一系列实验结果用来评估所提算法在仿真环境中的DP-SSP检测、源数估计以及DOA估计3个方面的性能,并与基于LRSSP[29]的方法进行了比较。为了突出所提两阶段DP-SSP方法每一步的有效性,本文还将SR-SSP去除与LRSSP[29]检测结合作为另外一种对比方法,用来评估DP-SSP的检测性能。

仿真混响室内为一个理想的封闭长方体,尺寸为8 m×6 m×4 m,AVS被放置在房间中心。采用镜像原理[31]生成混响室内的房间脉冲响应,声源被放置在以AVS为圆心,半径1.5 m的圆上。声源从TIMIT语料库[32]中挑选,采样率为16 kHz。麦克风接收到的信号进行STFT变换,FFT点数为1 024,帧间重叠为50%。加性噪声为高斯白噪声,其他实验参数设置如表1所示,表1中的阈值均为经验选取。

表1 参数设置Tab.1 Parameter settings

与文献[29]类似,本文采用均方角误差(Root-Mean-Square Angular Error,RMSAE)作为评估指标来量化所提方法在所有实验中的整体性能,RMSAE定义为:

(26)

(27)

3.1 DP-SSP 检测

不同混响环境下利用3种不同的检测规则检测到的时频点中真实DP-SSP占比如图4所示。

图4 采用3种方法得到的时频点中真实DP-SSP所占百分比随混响时间变化Fig.4 Variation of real DP-SSP percentage with reverberation using three detection rules

从图3中可以看到,随着混响的增强,基于3种检测规则的真实DP-SSP占比呈现下降趋势。但随着混响增加,所提的两阶段DP-SSP检测方法始终保持着明显优势。结合SR-SSP去除和LRSSP[29]检测的方法在DP-SSP检测上优于原始的LRSSP[29]方法,这说明对去除SR-SSP可以有效辅助DP-SSP的检测。通过对比也可发现,除了第一阶段的SR-SSP去除,第二阶段的子空间投影算法也在DP-SSP检测中展现极大优势。真实DP-SSP占比的提高使得所提算法能在高混响情况下实现更高精度的源计数和DOA估计。

3.2 源数目估计

表2给出了所提DP-SSP检测方法与基于LRSSP[29]的方法相比,不同混响环境下正确估计源信号数目的概率。从表中可以发现,2种检测规则在RT60=0.3 s的情况下均能准确估计源数,信噪比对估计结果的影响较小。当RT60<0.45 s,SNR=20 dB时,基于LRSSP[29]的方法基本可以正确估计源数,准确率大于90%,且当RT60<0.45 s时,基于LRSSP[29]的方法估计精度断崖式下降。此外,从表中可以看到,基于LRSSP[29]的方法性能受SNR影响很大,这说明该方法对噪声更加敏感。本文所提方法随着混响增强和信噪比降低,仍然能保持稳定有效的源数估计性能。

表2 在不同信噪比和混响程度下正确估计源数目的概率Tab.2 Percentage of the correct estimation of the number of sources with different SNRs and reverberation

3.3 DOA估计

图5中展示了SNR为20 dB,在RT60=0.6 s和RT60= 0.75 s的仿真混响环境下,不同算法进行DOA估计的可视化结果。3个源信号所在位置分别为(70°,0°)、(110°,15°)、(170°,30°)。然而,当混响时间增加到0.75 s时,基于LRSSP[29]的方法无法准确定位出每个源,如图5(b)所示,一些源的DOA估计与对应的真实DOA有很大偏差。与基于LRSSP[29]的方法不同,随着混响增加,所提算法仍能实现高精度的DOA估计,且误差较小。此外,图6展示了在多个源间隔较近的情况下DOA的估计结果。3个源的DOA分别为(70°,0°)、(110°,15°)、(150°,25°)。结果表明,基于LRSSP[29]的方法无法处理距离较近的多源DOA估计,而所提方法能实现准确估计,具有明显优势。

(a)基于LRSSP[29]的方法(RT60=0.60 s) (b)所提DP-SSP检测的方法(RT60=0.60 s)

(c)基于LRSSP[29]的方法(RT60=0.75 s) (d)所提DP-SSP检测的方法(RT60=0.75 s)图5 仿真环境下2种方法估计DOA的可视化结果图(SNR=20 dB)Fig.5 DOA estimation results in the simulated environment using different DOA estimation algorithms (SNR=20 dB)

(a)基于LRSSP[29]的方法(RT60=0.60 s) (b)所提DP-SSP检测的方法(RT60=0.60 s)图6 源相隔较近时2种方法估计DOA的可视化结果图(SNR=20 dB)Fig.6 DOA estimation results using different DOA estimation algorithms in the case of closely spaced sources (SNR=20 dB)

表3展示了利用2种不同算法计算的RMSAE随混响时间变化的情况。由表3可以看出,在RT60<0.6 s时,基于LRSSP[29]的方法在3个源DOA估计上表现较好,但当混响时间超过0.6 s时,该方法的DOA估计精度明显下降。值得注意的是,尽管所提方法的性能不可避免地随着混响的增加而下降,但当RT60增加到0.9 s时,在SNR=20 dB的条件下,对比算法已经完全无法进行多源DOA估计,然而所提算法仍然具有一定性能。当信噪比降到15 dB,且RT60=0.9 s时,所提方法性能略有下降。

表3 RMSAE随混响和SNR变化Tab.3 Variation of RMSAE with the reverberation and SNR

图7展示了对不同数量的源进行DOA估计后计算得到的RMSAE随混响变化的趋势。结果表明,随着源数量的增加,本文所提方法性能不可避免地下降,特别是存在4个源的情况下,性能下降明显,这是因为DP-SSP检测的准确性随着源数量的增加而下降。所提算法在2个源的情况下,在不同的环境中都能展现出较好的性能,即使当RT60增加到0.9 s时,仍然能获得较高精度。

图7 不同源数情况下使用所提方法计算的RMSAE随混响时间变换图(SNR=20 dB)Fig.7 Variation of RMSAE with reverberation time in case of different number of sources by the proposed algorithm (SNR=20 dB)

4 结束语

本文针对传统算法利用单个AVS在高混响情况下多源DOA估计性能差的问题,提出了一种新的多源DOA估计算法。该算法的关键在于检测DP-SSP,这些DP-SSP为DOA的估计提供了有利信息。本文所提出的DP-SSP检测算法有两方面的贡献。首先,针对单个AVS设计了一种基于频率平滑的SR-SSP去除方法,有效地消除了混响带来的不利影响,为准确检测DP-SSP提供帮助。其次,采用子空间投影的算法进一步细化DP-SSP的检测,消除大部分MSP,提取正确的DP-SSP。实验结果表明,该算法对DP-SSP的高精度检测使得本文算法在源数目估计和DOA估计方面取得较好性能。