引言

声源定位是音频信号处理中重要的研究方向之一, 经常被应用于空间滤波、声源分离、去混响、语音增强、视频会议、机器人听觉等领域中^[1-5]。目前声源波达方位(DOA)主要利用声源定位方法从传声器阵列采集的声场信息中进行估计, 其中球面传声器阵列可以对三维声场进行等分辨率分析, 且便于将采集的声场信息转换到球谐波域中进行处理而受到人们的关注^[6-10]。在球谐波域中导向矢量与阵列分布无关, 并且可以将导向矢量与频率解耦^[7-8], 这为声源的定位带来了便利。

在球谐波域中常用的声源定位方法有: 可控波束响应方法 ^[10]、子空间类方法^[11-12]、最大似然估计方法^[13-14]、伪声强(PIV)类方法^[15-23]等。其中PIV^[15]声源定位方法具有声源定位的闭式解, 计算量较低, 适用于实时的声源定位场合, 但是PIV方法只采用了0阶和1阶的声场球谐波特征波束, 使得其空间分辨率较低, 在高混响和多声源环境下定位性能较差^[16-17]。为了弥补PIV方法的缺陷提出了一些改进方法, 如将声场高阶特征波束引入到PIV方法中提升其空间分辨率的增广伪声强(AIV)方法^[18-19]、局部扫描可控波束响应方法^[20]、子空间伪声强方法^[21]和直接路径占优(DPD)检测伪声强(DPD-PIV)方法^[22]等。

将高阶特征波束引入到PIV方法中提升了定位性能, 但是相应地增大了计算量, 且低频段的高阶特征波束易受噪声污染^[7-8], 限制了定位方法性能的提升。为了减小该问题的影响常用方法是将低频特征波束舍去^[11-22], 但是噪声的随机性使得保留下来的高阶特征波束中仍存在被噪声严重污染的特征波束, 并且若舍去的低频信息较多将会影响能量主要位于低频段的声源的定位。针对该问题, 文献[23]提出了阶数感知(OA)定位方法, 利用各阶特征波束因噪声产生的差异感知每阶特征波束受噪声污染的程度, 从中筛选出受噪声影响较小的特征波束用于伪声强声源定位中, 提升了定位方法的抗噪能力。并且该方法采用分阶处理思想, 能够将低频段受噪声影响较小的低阶特征波束筛选出来, 从而保留了声源的低频信息, 提升了伪声强定位方法的鲁棒性。

现有阶数感知因子只利用了单个时频点的信息, 限制了其对于混响和多声源干扰程度的鉴别能力, 以及在高频段各阶特征波束因噪声产生的差异减小引起的噪声感知能力下降。为解决现有阶数感知存在的问题, 本文将相邻多个时频点的信息应用到阶数感知中, 利用相邻时频点中源自同一声源的特征波束之间存在的相关性, 构造了局部时频相关的阶数感知因子, 将受噪声污染较小的单声源直达声信号占优的特征波束筛选出来用于伪声强声源定位方法中, 提升定位方法的抗噪和抗混响能力。通过理论分析、仿真实验和实测实验验证, 所提方法相比已有的伪声强声源定位方法在降低计算量的同时具有更高的定位精度和鲁棒性, 且随信噪比的降低、混响时间的增大、声源间距的减小优势提升得越明显。

1.   理论基础

1.1   信号模型

考虑由$Q$个传声器构成的半径为${r_a}$的球面阵列, 以阵列中心为球坐标系原点, 第$q$个传声器的坐标为$({r_a},{\varOmega _q}) = ({r_a},{\theta _q},{\varphi _q})$, 其中${\theta _q}$为俯仰角、${\varphi _q}$为方位角。$p(\tau ,k,{r_a},{\varOmega _q})$为第$q$个传声器采集到的声压信号, 其中$k$为波数, $\tau$为帧序号。对声压信号进行离散球傅里叶变换(DSFT)得到声场的球谐波系数^[8]:

其中, ${\alpha _q}$为离散采样权值, 若阵列为均匀采样分布^[8,24] ${\alpha _q} = {{4\pi } \mathord{\left/ {\vphantom {{4\pi } Q}} \right. } Q}$, ${( \cdot )^*}$表示复共轭, ${\rm Y} _n^m(\varOmega )$为$n$阶$m$维球谐波函数, $- n \leqslant m \leqslant n$, 其定义为^[8,24]

其中, $( \cdot )!$为阶乘, ${\text{P}}_n^m( \cdot )$为$n$阶$m$维连带勒让德函数。若需要的最大球谐波阶数为N, 为避免式(1)离散球傅里叶变换出现空间混叠, 要求N和传声器个数Q之间满足${(N + 1)^2} \leqslant Q$以及N和波数$k$之间满足$kr < N$ ^[7-8]。球谐波函数在离散采样点处满足正交性^[24]: $\sum\nolimits_{q = 1}^Q {{\alpha _q}{\rm Y} _n^m({\varOmega _q}){\rm Y} _{n'}^{m'*}({\varOmega _q})} = {\delta _{n - n'}}{\delta _{m - m'}}$, $n \leqslant N$, 其中${\delta _{n - n'}}$与${\delta _{m - m'}}$为克罗内克$\delta$函数。

对于由$L$个声源辐射出的平面波在球谐波域中可以表示为^[10-15]

其中, ${n_{nm}}(\tau ,k) = \sum\nolimits_{q = 1}^Q {{\alpha _q}{n_q}(\tau ,k){\rm Y} _n^{m*}({\varOmega _q})}$为加性噪声的球谐波系数, ${n_q}(\tau ,k)$为第$q$个传声器的噪声幅度; ${s_l}\left( {\tau ,k} \right)$和${\varTheta _l} = ({\theta _l},{\phi _l})$分别为第$l$个声源的幅度和DOA; ${b_n}(k{r_a})$为模式强度, 与波数$k$、球阵半径${r_a}$以及球阵的配置有关。对于本文采用的远场刚性球其表达式为^[24]

其中, ${\text{i}}$为虚数符号, ${{{\rm{j}}} _n}( \cdot )$为第一类球贝塞尔函数, ${{\rm{h}}} _n^{(2)}( \cdot )$为第二类球汉克尔函数, $( \cdot )'$表示一阶微分。

对式(3)进行模式强度补偿, 即式(3)除以${b_n}(k{r_a})$, 可以得到特征波束^[6-8]:

将$n$阶所有特征波束表示为向量的形式${{\boldsymbol{a}}_{\boldsymbol{n}}}(\tau ,k) = {[{a_{n( - n)}}(\tau ,k), \cdots, {a_{n0}}(\tau ,k), \cdots, {a_{nn}}(\tau ,k)]^{\text{T}}}$, ${[ \cdot ]^{\text{T}}}$为转置。

1.2   伪声强定位方法

1.2.1   PIV方法

设声场在时频点$(\tau ,k)$的伪声强向量为${\boldsymbol{U}}(\tau ,k) = {[{U_x}(\tau ,k),{U_y}(\tau ,k),{U_{\textit{z}}}(\tau ,k)]^{\text{T}}}$, 其中任一项的表达式为^[15]

其中, ${{\rm{Re}}} ( \cdot )$表示取实部, $\varsigma \in \{ x,y,{\textit{z}}\}$。${D_{ - \varsigma }}(\tau ,k)$的表达式为^[15]

式中, ${\varPsi _{ - x}} = ({\pi \mathord{\left/ {\vphantom {\pi 2}} \right. } 2},\pi )$, ${\varPsi _{ - y}} = ({\pi \mathord{\left/ {\vphantom {\pi 2}} \right. } 2},{{ - \pi } \mathord{\left/ {\vphantom {{ - \pi } 2}} \right. } 2})$, ${\varPsi _{ - {\textit{z}}}} = (\pi ,0)$。则估计的DOA由下式给出^[15]:

其中$|| \cdot ||$为二范数。从式(6)—式(8)中可以看出PIV方法估计声源DOA具有闭式解, 计算量较低, 但只利用0, 1阶特征波束使其空间分辨率较差。为改善该问题, 出现了一些应用高阶特征波束进行PIV定位的方法, 增广伪声强(AIV) ^[18-19]为其中典型的一种方法。

1.2.2   AIV方法

0阶模式强度${b_0}(kr)$的模值比较大, 进行模式补偿得到的0阶特征波束${a_{00}}(\tau ,k)$相对受噪声影响较小, 可以用${a_{00}}(\tau ,k)$得到声场幅度的近似值, $s(\tau ,k) \simeq \sqrt {4\pi } {a_{00}}(\tau ,k)$, 构造${\widetilde a_{nm}}(\tau ,k,\varTheta ) = \sqrt {4\pi } {a_{00}}(\tau ,k){\rm Y} _n^{m*}(\varTheta )$来逼近${a_{nm}}(\tau ,k)$, 定义代价函数^[18-19]:

通过寻找使代价函数${{\varPi }}\left( {\tau ,k,\varTheta } \right)$最小的$\varTheta$作为声源定位结果^[18-19]:

AIV方法通常采用局部扫描寻找式(10)最小值, 将PIV定位结果作为扫描初始值, 在初始值附近一定范围内进行搜索。用${a_{00}}(\tau ,k)$近似声场幅度在噪声、混响的影响下会产生偏差, 且高阶特征波束在低频段易受噪声影响, 使得AIV方法的定位性能在低信噪比和高混响环境下会有所下降。

2.   时频相关球谐波阶数感知鲁棒声源定位方法

针对低信噪比和高混响环境下PIV和AIV方法定位性能下降的问题, 本节提出局部时频相关的阶数感知(TFCOA)因子以及基于该因子的TFCOA-PIV和TFCOA-AIV定位方法。本节首先介绍时频相关球谐波阶数感知因子的构造, 并对其在噪声和混响环境下的表现进行理论分析, 最后介绍TFCOA-PIV和TFCOA-AIV定位方法流程。

2.1   时频相关球谐域阶数感知因子的构造

根据语音信号的稀疏性^[25], 若由相邻时频点构成的时频块中的信号主要来自于同一方位, 在噪声影响较小时它们的特征波束之间存在相关性, 可以利用这种相关性将受噪声、混响影响较小的单声源直达声信号占优的时频点块筛选出来用于定位, 提升定位方法的性能。为表达的简洁, 在下文中与时频点有关的变量除第一次出现外将省略$(\tau ,k)$。

考虑以时频点$(\tau ,k)$作为右下角参考点的时频块, 包含连续${J_\tau }$帧, 每一帧连续${J_k}$个频点, 共${J_\tau }{J_k}$个时频点。将其中任一时频点$(\tau - {j_\tau },k + {j_k})$的$n$阶特征波束向量${{\boldsymbol{a}}_{\boldsymbol{n}}}(\tau - {j_\tau },k + {j_k})$用${\boldsymbol{a}}_{\boldsymbol{n}}^{{j_\tau },{j_k}}(\tau ,k)$表示, 其中的特征波束和信号幅度分别用$a_{nm}^{{j_\tau },{j_k}}(\tau ,k)$, $s_l^{{j_\tau },{j_k}}(\tau ,k)$表示, ${j_\tau } \in [0, {J_\tau } - 1]$, ${j_k} \in [0, {J_k} - 1]$。定义:

其中, ${\varGamma _n} = {{4\pi } \mathord{\left/ {\vphantom {{4\pi } {(2n + 1)}}} \right. } {(2n + 1)}}$。为衡量参考时频点$(\tau ,k)$的n阶特征波束向量${\boldsymbol{a}}_{\boldsymbol{n}}^{0,0}$与${\boldsymbol{a}}_{\boldsymbol{n}}^{{j_\tau },{j_k}}$的相关程度, 定义:

其中, $\langle {\boldsymbol{a}},{\boldsymbol{b}}\rangle$为向量${\boldsymbol{a}}$与${\boldsymbol{b}}$的内积。根据柯西−施瓦茨不等式有$|H_n^{{j_\tau },{j_k}}| \leqslant {{{\varGamma}}_n}||{\boldsymbol{a}}_{\boldsymbol{n}}^{0,0}||||{\boldsymbol{a}}_{\boldsymbol{n}}^{{j_\tau },{j_k}}||$, 即$|H_n^{{j_\tau },{j_k}}| \leqslant \sqrt {P_n^{0,0}P_n^{{j_\tau },{j_k}}}$, 当${\boldsymbol{a}}_{\boldsymbol{n}}^{{j_\tau },{j_k}} = \beta _{}^{{j_\tau },{j_k}}{\boldsymbol{a}}_{\boldsymbol{n}}^{0,0}$时等号成立, $\beta _{}^{{j_\tau },{j_k}}$为常数, 即$|H_n^{{j_\tau },{j_k}}|$取值与${\boldsymbol{a}}_{\boldsymbol{n}}^{0,0}$和${\boldsymbol{a}}_{\boldsymbol{n}}^{{j_\tau },{j_k}}$的的相关程度有关。

直达声中含有声源的真实DOA信息, 是声源定位的重要信息, 现对只包含直达声信号的情况进行分析。当时频点$(\tau ,k)$和$(\tau - {j_\tau },k + {j_k})$中只含有声源$l$的直达声信号时, 根据式(5)有$a_{nm}^{0,0} = s_l^{0,0}{\rm Y} _n^{m*}({\varTheta _l})$, $a_{nm}^{{j_\tau },{j_k}} = s_l^{{j_\tau },{j_k}}{\rm Y} _n^{m*}({\varTheta _l})$, 此时满足条件${\boldsymbol{a}}_{\boldsymbol{n}}^{{j_\tau },{j_k}} = \beta _{}^{{j_\tau },{j_k}}{\boldsymbol{a}}_{\boldsymbol{n}}^{0,0}$, $\beta _{}^{{j_\tau },{j_k}} = s_l^{{j_\tau },{j_k}}/s_l^{0,0}$。根据式(11)和球谐波函数求和公式^[8]有$P_n^{0,0} = |s_l^{0,0}{|^2}$, $P_n^{{j_\tau },{j_k}} = |\beta _{}^{{j_\tau },{j_k}}{|^2}|s_l^{0,0}{|^2}$, $P_n^{0,0}$, $P_n^{{j_\tau },{j_k}}$与阶数n无关, 即各阶$P_n^{0,0}$, $P_n^{{j_\tau },{j_k}}$相等。根据式(12)有$H_n^{{j_\tau },{j_k}} = \beta {_{}^{{j_\tau },{j_k}*} }|s_l^{0,0}{|^2}$, 同样与阶数n无关, 因此$H_n^{{j_\tau },{j_k}}H_{n - 1}^{{j_\tau },{j_k}*} = |\beta _{}^{{j_\tau },{j_k}}{|^2}|s_l^{0,0}{|^4}$为与阶数n无关的实数, 且$H_n^{{j_\tau },{j_k}}H_{n - 1}^{{j_\tau },{j_k}*} = P_n^{0,0}P_n^{{j_\tau },{j_k}} = P_{n - 1}^{0,0}P_{n - 1}^{{j_\tau },{j_k}}$。因此在$(\tau ,k)$和$(\tau - {j_\tau },k + {j_k})$中只含有单声源直达声信号时, 各阶$P_n^{0,0}$, $P_n^{{j_\tau },{j_k}}$相等, 各阶$H_n^{{j_\tau },{j_k}}H_{n - 1}^{{j_\tau },{j_k}*}$相等且等于最大值$P_n^{0,0}P_n^{{j_\tau },{j_k}}$或$P_{n - 1}^{0,0}P_{n - 1}^{{j_\tau },{j_k}}$。

若时频点$(\tau ,k)$和$(\tau - {j_\tau },k + {j_k})$中含有多个声源的直达声信号时, 不妨假设包含声源$l$和$l'$的直达声信号, 则$a_{nm}^{0,0} = s_l^{0,0}{\rm Y} _n^{m*}({\varTheta _l}) + s_{l'}^{0,0}{\rm Y} _n^{m*}({\varTheta _{l'}})$, $a_{nm}^{{j_\tau },{j_k}} = s_l^{{j_\tau },{j_k}}{\rm Y} _n^{m*}({\varTheta _l}) + s_{l'}^{{j_\tau },{j_k}}{\rm Y} _n^{m*}({\varTheta _{l'}})$, 若${\varTheta _l} \ne {\varTheta _{l'}}$且两个声源辐射的语音信号不同时, 则由各${\rm Y} _n^{m*}({\varTheta _l})$与${\rm Y} _n^{m*}({\varTheta _{l'}})$之间的差异, 以及语音信号随机性导致的$s_l^{{j_\tau },{j_k}}/s_l^{0,0} \ne s_{l'}^{{j_\tau },{j_k}}/s_{l'}^{0,0}$, 使得${\boldsymbol{a}}_{\boldsymbol{n}}^{{j_\tau },{j_k}} \ne \beta _{}^{{j_\tau },{j_k}}{\boldsymbol{a}}_{\boldsymbol{n}}^{0,0}$。此时$P_n^{{j_\tau },{j_k}}$和$H_n^{{j_\tau },{j_k}}$可分别表示为

其中, ${{{\rm{P}}} _n}( \cdot )$为n阶勒让德多项式, ${\varTheta _{ll'}}$为声源$l$和$l'$的角度间距。从式(13)和式(14)中可以看出, 此时各阶$P_n^{{j_\tau },{j_k}}$, $H_n^{{j_\tau },{j_k}}$的值与${{{\rm{P}}} _n}\left( {\cos \left( {{\varTheta _{ll'}}} \right)} \right)$有关, 而当${\varTheta _l} \ne {\varTheta _{l'}}$时各阶${{{\rm{P}}} _n}\left( {\cos \left( {{\varTheta _{ll'}}} \right)} \right)$的值存在差异, 使得各阶$P_n^{{j_\tau },{j_k}}$, $H_n^{{j_\tau },{j_k}}$不再相等, 同理各阶$P_n^{0,0}$也不再相等。可以推广到声源个数大于2的情况。

另外在实际室内环境中存在噪声和混响, 在噪声的随机性、各阶${b_n}(k{r_a})$的差异、不同混响信号传输路径的不同、语音信号的随机性等因素的影响下, 综合导致${\boldsymbol{a}}_{\boldsymbol{n}}^{{j_\tau },{j_k}} \ne \beta _{}^{{j_\tau },{j_k}}{\boldsymbol{a}}_{\boldsymbol{n}}^{0,0}$, 各阶$P_n^{0,0}$, $P_n^{{j_\tau },{j_k}}$, $H_n^{{j_\tau },{j_k}}$, $H_{n - 1}^{{j_\tau },{j_k}}$不再相等。因此在多声源混响噪声环境下, $H_n^{{j_\tau },{j_k}}H_{n - 1}^{{j_\tau },{j_k}*}$实部与最大值之间的差距, $P_n^{0,0}$, $P_n^{{j_\tau },{j_k}}$与$P_{n - 1}^{0,0}$, $P_{n - 1}^{{j_\tau },{j_k}}$之间的差异, 反映了${\boldsymbol{a}}_{\boldsymbol{n}}^{0,0}$, ${\boldsymbol{a}}_{{\boldsymbol{n - }}1}^{0,0}$, ${\boldsymbol{a}}_{\boldsymbol{n}}^{{j_\tau },{j_k}}$, ${\boldsymbol{a}}_{{\boldsymbol{n - }}{\text{1}}}^{{j_\tau },{j_k}}$受噪声、混响以及多个声源间相互干扰的影响程度。

若时频点$(\tau ,k)$和$(\tau - {j_\tau },k + {j_k})$中包含的信号以经同一方位反射的混响信号为主, 则可能会出现误选情况, 为减小这种可能性, 考虑时频块中每个时频点的特征波束向量与参考时频点特征波束向量之间的相关程度。将除$({j_\tau },{j_k}) = (0,0)$外的$H_n^{{j_\tau },{j_k}}$和$P_n^{{j_\tau },{j_k}}$表示为向量形式, ${{\boldsymbol{H}}_{\boldsymbol{n}}}(\tau ,k) = {[H_n^{0,1}, \cdots, H_n^{{j_\tau },{j_k}}, \cdots, H_n^{{J_\tau } - 1,{J_k} - 1}]^{\text{T}}}$, ${P}_{n}(\tau ,k)={[{P}_{n}^{0, 1},\cdots, {P}_{n}^{{j}_{\tau },{j}_{k}},\cdots, {P}_{n}^{{J}_{\tau }-1,{J}_{k}-1}]}^{\text{T}}$, 同理可以得到${{\boldsymbol{H}}_{{\boldsymbol{n - }}{\text{1}}}}(\tau ,k)$和${{\boldsymbol{P}}_{{\boldsymbol{n - }}{\text{1}}}}(\tau ,k)$。定义时频相关阶数感知因子(TFCOA)为

将式(15)中分子的内积展开, 由前面的分析可知, 分子中的任一项都小于分母中对应的项, 即$2{{\rm{Re}}} [H_n^{{j_\tau },{j_k}}H_{n - 1}^{{j_\tau },{j_k}*}] \leqslant 2\sqrt {P_n^{0,0}P_n^{{j_\tau },{j_k}}} \sqrt {P_{n - 1}^{0,0}P_{n - 1}^{{j_\tau },{j_k}}}$, 而$2\sqrt {P_n^{0,0}P_n^{{j_\tau },{j_k}}} \sqrt {P_{n - 1}^{0,0}P_{n - 1}^{{j_\tau },{j_k}}} \leqslant P_n^{0,0}P_n^{{j_\tau },{j_k}} + P_{n - 1}^{0,0}P_{n - 1}^{{j_\tau },{j_k}}$, 因此$|K_n^{}| \leqslant 1$, $K_n^{}$反映了时频块中其他时频点的特征波束向量与参考时频点$(\tau ,k)$的特征波束向量的相关程度。只有当时频块中只包含源自同一方位的信号时, 每一时频点都有$P_n^{0,0}P_n^{{j_\tau },{j_k}} = P_{n - 1}^{0,0}P_{n - 1}^{{j_\tau },{j_k}}$, ${\boldsymbol{a}}_{\boldsymbol{n}}^{{j_\tau },{j_k}} = \beta _{}^{{j_\tau },{j_k}}{\boldsymbol{a}}_{\boldsymbol{n}}^{0,0}$, 此时$K_n^{} = 1$。

2.2   时频相关阶数感知因子抗噪抗混响能力分析

2.2.1   抗混响能力分析

为简化分析和突出时频相关阶数感知因子在多声源和混响条件下的表现, 本节只考虑多声源和混响, 忽略噪声的影响。

先考虑$K_n^{}$分子中的任一项${\rm Re} [H_n^{{j_\tau },{j_k}}H_{n - 1}^{{j_\tau },{j_k}*}]$。在多声源和混响环境下, 令$a_{nm}^{0,0} = \sum\nolimits_l {s_l^{0,0}{\rm Y} _n^{m*}({\varTheta _l})}$, $a_{nm}^{{j_\tau },{j_k}} = \sum\nolimits_d^{} {s_d^{{j_\tau },{j_k}}{\rm Y} _n^{m*}({\varTheta _d})}$。根据式(12)和球谐波函数求和公式, 可得

其中, ${\rm P} _n^{l,d} = {{\rm P} _n}(\cos ({\varTheta _{ld}}))$, ${\rm P} _{n - 1}^{l',d'} = {{\rm P} _{n - 1}}(\cos ({\varTheta _{l'd'}}))$, ${\varTheta _{ld}}$为声源$l$和声源$d$的DOA角度间距。式(16)中除$l = l'$, $d = d'$以外的每一项都存在一个对称项${{\rm{Re}}} (s_{l'}^{0,0}s_{d'}^{{j_\tau },{j_k}*} s_l^{0,0*}s_d^{{j_\tau },{j_k}}){\rm P} _{n - 1}^{l,d}{\rm P} _n^{l',d'}$, 将两项合并在一起, ${{\rm{Re}}} (s_l^{0,0}s_d^{{j_\tau },{j_k}*}s_{l'}^{0,0*}s_{d'}^{{j_\tau },{j_k}}) ({\rm P} _n^{l,d}{\rm P} _{n - 1}^{l',d'} + {\rm P} _{n - 1}^{l,d}{\rm P} _n^{l',d'})$。即${{\rm{Re}}} [H_n^{{j_\tau },{j_k}}H_{n - 1}^{{j_\tau },{j_k}*}]$由$(\tau ,k)$中声源和$(\tau - {j_\tau },k + {j_k})$中声源的DOA角度间距的函数${\rm P} _n^{l,d}{\rm P} _{n - 1}^{l',d'} + {\rm P} _{n - 1}^{l,d}{\rm P} _n^{l',d'}$和每个声源的幅度共同决定, 而${\rm P} _n^{l,d}{\rm P} _{n - 1}^{l',d'} + {\rm P} _{n - 1}^{l,d}{\rm P} _n^{l',d'}$随声源DOA角度间距${\varTheta _{ld}}$, ${\varTheta _{l'd'}}$的变化情况如图1所示。

图  1  各阶${\rm P} _n^{l,d}{\rm P} _{n - 1}^{l',d'} + {\rm P} _{n - 1}^{l,d}{\rm P} _n^{l',d'}$随声源DOA角度间距变化情况 (a) n = 1; (b) n = 2; (c) n = 3

下载: 全尺寸图片 幻灯片

从图1中可以看出, 只有当${\varTheta _{ld}}$, ${\varTheta _{l'd'}}$为$0$或$2\pi$时, 即声源$l$和声源$d$、声源$l'$和声源$d'$的DOA相同时, ${\rm P} _n^{l,d}{\rm P} _{n - 1}^{l',d'} + {\rm P} _{n - 1}^{l,d}{\rm P} _n^{l',d'}$取最大值2, 之后随着声源角度间距${\varTheta _{ld}}$, ${\varTheta _{l'd'}}$的增大${\rm P} _n^{l,d}{\rm P} _{n - 1}^{l',d'} + {\rm P} _{n - 1}^{l,d}{\rm P} _n^{l',d'}$逐渐减小, 阶数n越大减小的速度越快, 虽然n > 1时${\rm P} _n^{l,d}{\rm P} _{n - 1}^{l',d'} + {\rm P} _{n - 1}^{l,d}{\rm P} _n^{l',d'}$下降到0附近开始波动, 但是波动的幅度较小, 在${\varTheta _{ld}}$, ${\varTheta _{l'd'}}$增大到$\pi$时达到最小值−2。对于式(16)中的$l = l'$, $d = d'$项, 此时${\varTheta _{ld}}$, ${\varTheta _{l'd'}}$, ${\rm P}_n^{l,d}{\rm P}_{n - 1}^{l,d}$按图1中对角线的趋势变化, 幅度为其一半, 波动幅度更小。从整体上看${\rm P} _n^{l,d}{\rm P} _{n - 1}^{l',d'} + {\rm P} _{n - 1}^{l,d}{\rm P} _n^{l',d'}$随声源角度间距的增大是减小的趋势。即当$(\tau ,k)$和$(\tau - {j_\tau },k + {j_k})$中声源DOA角度间距增大时,${\rm Re} [H_n^{{j_\tau },{j_k}}H_{n - 1}^{{j_\tau },{j_k}*}]$的取值呈减小趋势, 反映了$(\tau ,k)$和$(\tau - {j_\tau },k + {j_k})$中声源信息的相似程度。

考虑一种最坏情况, 时频块中任一时频点都包含相同方位的直达声和混响信号, 且任一时频点$(\tau - {j_\tau },k + {j_k})$满足条件${\boldsymbol{a}}_{\boldsymbol{n}}^{{j_\tau },{j_k}} = {\beta ^{{j_\tau },{j_k}}}{\boldsymbol{a}}_{\boldsymbol{n}}^{0,0}$。根据式(11)和式(12)有$P_n^{{j_\tau },{j_k}} = |{\beta ^{{j_\tau },{j_k}}}{|^2}P_n^{0,0}$, $H_n^{{j_\tau },{j_k}} = {\beta ^{{j_\tau },{j_k}*}}P_n^{0,0}$, $P_{n - 1}^{{j_\tau },{j_k}} = |{\beta ^{{j_\tau },{j_k}}}{|^2}P_{n - 1}^{0,0}$, $H_{n - 1}^{{j_\tau },{j_k}} = {\beta ^{{j_\tau },{j_k}*}}P_{n - 1}^{0,0}$, 代入式(15), 有

从式(17)可以看出, $K_n^{}$在最坏情况下退化为单时频点的阶数感知(OA)检测因子^[23]$O_n^{n - 1}$, 此时$O_n^{n - 1}$的值由$P_n^{0,0}$和$P_{n - 1}^{0,0}$的差距决定。由于语音信号具有随机性, 要使每个时频点都满足条件${\boldsymbol{a}}_{\boldsymbol{n}}^{{j_\tau },{j_k}} = {\beta ^{{j_\tau },{j_k}}}{\boldsymbol{a}}_{\boldsymbol{n}}^{0,0}$几乎是不可能的, 除非时频块中每个时频点只含有来自同一声源的直达声信号, 此时$P_n^{0,0} = P_{n - 1}^{0,0}$, $K_n^{} = O_{n,n - 1}^{} = 1$。因此多声源混响环境下$K_n^{} < O_{n,n - 1}^{}$, 即TFCOA检测因子对多个声源和混响的区分能力优于OA检测因子, 当$K_n^{}$值较大时说明时频块中的信号主要为来自同一声源的直达声信号。

2.2.2   时频相关阶数感知因子抗噪声能力分析

本节分析时频相关阶数感知因子在噪声环境下的表现, 为简化分析只考虑噪声, 忽略语音和混响信号。设噪声的均值为$\mu$, 方差为$\sigma _{}^2$, 此时${K_n}$为

其中, ${\boldsymbol{n}}_{\boldsymbol{n}}^{{\text{0,0}}}$和${\boldsymbol{n}}_{\boldsymbol{n}}^{{j_\tau },{j_k}}$分别为时频点$(\tau ,k)$和$(\tau - {j_\tau },k + {j_k})$的第$n$阶噪声球谐波系数构成的向量。从上式可以看出, 当时频点块中只包含噪声时, $K_n$的值由两方面决定, 一方面由模式强度补偿${1 \mathord{\left/ {\vphantom {1 {(|{b_n}(k{r_a}){|^2}|{b_{n - 1}}(k{r_a}){|^2})}}} \right. } {(|{b_n}(k{r_a}){|^2}|{b_{n - 1}}(k{r_a}){|^2})}}$与${1 \mathord{\left/ {\vphantom {1 {|{b_n}(k{r_a}){|^4}}}} \right. } {|{b_n}(k{r_a}){|^4}}}$, ${1 \mathord{\left/ {\vphantom {1 {|{b_{n - 1}}(k{r_a}){|^4}}}} \right. } {|{b_{n - 1}}(k{r_a}){|^4}}}$的差距决定, 这与OA方法是相同的, 在本文中不再讨论; 另一方面由参考时频点$(\tau ,k)$的${\boldsymbol{n}}_{\boldsymbol{n}}^{{\text{0,0}}}\text{,}{\text{ }}{\boldsymbol{n}}_{{\boldsymbol{n - }}1}^{{\text{0,0}}}$与时频点$(\tau - {j_\tau },k + {j_k})$的${\boldsymbol{n}}_{\boldsymbol{n}}^{{j_\tau },{j_k}}\text{,}{\text{ }}{\boldsymbol{n}}_{{\boldsymbol{n - }}1}^{{j_\tau },{j_k}}$的内积决定。由于噪声是随机的, 因此采用统计平均观察$K_n^{}$在纯噪声环境下的取值情况, 为分析方便, 将$K_n^{}$的分子分母分开讨论。

先考虑$K_n^{}$的分子中的任一时频点$(\tau - {j_\tau },k + {j_k})$对应的项, 由于${b_n}(k{r_a})$为确定函数, 省略${b_n}(k{r_a})$不影响对均值的讨论。当球阵的传声器分布为均匀分布、高斯分布和等角度分布时, 有

关于式(19)的证明见附录A。即式(18)分子中任一时频点对应项的均值均为0, 将式(18)中$K_n^{}$的分子和分母同除以总项数$({J_\tau }{J_k} - 1)$, 则分子是$({J_\tau }{J_k} - 1)$项随机变量的平均, 接近于随机变量平均值0, 即纯噪声时$K_n^{}$的分子取值在0附近。

对于式(18)中的分母, 其均值不为0, 证明见附录B。因此当只有噪声时, $K_n^{}$取值在0附近, 且与频率无关, 因此$K_n^{}$在整个频段都具有衡量噪声大小的能力, 弥补了OA检测因子在高频段衡量噪声大小能力下降的问题。在信号与噪声共存时, 若$K_n^{}$值较大则说明该时频块中语音信号占优, 若$K_n^{}$值较小则说明噪声占优, 可以舍去。

2.3   球谐域时频相关阶数感知伪声强声源定位方法

将时频相关阶数感知因子(TFCOA)与PIV和AIV方法相结合, 形成鲁棒的伪声强声源定位方法, 分别记为TFCOA-PIV方法、TFCOA-AIV方法, 相应的声源定位流程图如图2所示, 具体步骤如下:

图  2  所提定位方法流程图

下载: 全尺寸图片 幻灯片

(1) 将阵列采集的声压信号转换到球谐域。把采集的声压信号${{\boldsymbol{p}}_{0,0}}(\tau ,k)$, ···, ${{\boldsymbol{p}}_{{J_\tau },{J_k}}}(\tau ,k)$代入到式(1)中进行离散球傅里叶变换, 再进行模式强度补偿得到特征波束向量${\boldsymbol{a}}_{{\boldsymbol{nm}}}^{0,0}(\tau ,k)$, ···, ${\boldsymbol{a}}_{{\boldsymbol{nm}}}^{{J_\tau },{J_k}}(\tau ,k)$。其中, ${{\boldsymbol{p}}_{{j_\tau },{j_k}}}(\tau ,k) = {[p_1^{{j_\tau },{j_k}},p_2^{{j_\tau },{j_k}}, \cdots ,p_Q^{{j_\tau },{j_k}}]^{\text{T}}}$, 为Q个传声器采集声压信号组成的向量, ${{p}}_q^{{j_\tau },{j_k}} = p_q^{}(\tau - {j_\tau },k + {j_k})$为第q个传声器采集的声压信号; ${\boldsymbol{a}}_{{\boldsymbol{nm}}}^{{j_\tau },{j_k}}$为时频点$(\tau - {j_\tau }, k + {j_k})$所有特征波束构成的向量, ${\boldsymbol{a}}_{{\boldsymbol{nm}}}^{{j_\tau },{j_k}} = {[a_{00}^{{j_\tau },{j_k}},a_{1( - 1)}^{{j_\tau },{j_k}}, \cdots ,a_{NN}^{{j_\tau },{j_k}}]^{\text{T}}}$。

(2) 计算$K_n^{}(\tau ,k)$并进行时频点和可靠特征波束筛选。将时频块各时频点的$N$阶和$N - 1$阶的特征波束向量代入到式(15)中计算$K_N^{}(\tau ,k)$, 并与阈值${T_N}$比较, 超过${T_N}$则可靠特征波束阶数$\kappa = N$; 否则计算次高阶的$K_{N - 1}^{}(\tau ,k)$, 判断是否超过阈值${T_{N - 1}}$, 超过则$\kappa = N - 1$; 否则处理下一阶, 以此类推, 直到$K_2^{}(\tau ,k)$, 若$K_2^{}(\tau ,k)$超过阈值${T_2}$则$\kappa = 2$, 否则丢弃该时频点。

(3) 计算局部时频相关矩阵并进行特征值分解。由筛选出的时频块的$\kappa$阶可靠特征波束构成局部时频相关矩阵^[11,21,22]:

其中, ${\boldsymbol{a}}_{{\boldsymbol{\kappa m}}}^{{j_\tau },{j_k}} = {[a_{00}^{{j_\tau },{j_k}},a_{1( - 1)}^{{j_\tau },{j_k}}, \cdots ,a_{\kappa \kappa }^{{j_\tau },{j_k}}]^{\text{T}}}$。对${{\boldsymbol{R}}_a}(\tau ,k)$进行特征值分解, 可以表示为^[21-22]

其中, ${{\boldsymbol{U}}_{\boldsymbol{S}}}$为信号特征向量构成的信号子空间矩阵, ${{\boldsymbol{U}}_{\boldsymbol{V}}}$为噪声特征向量构成的噪声子空间矩阵, ${{\boldsymbol{\varSigma}} _{\boldsymbol{S}}}$和${{\boldsymbol{\varSigma }}_{\boldsymbol{V}}}$为由信号和噪声特征值构成的对角矩阵。

(4) PIV定位和AIV定位。选取${{\boldsymbol{U}}_{\boldsymbol{S}}}$中最大特征值对应的特征向量${{{\widehat{{\boldsymbol{a}}}}}}_{{{{\boldsymbol{\kappa m}}}}}^{}(\tau ,k) = [\widehat a_{00}^{}(\tau ,k),\widehat a_{1( - 1)}^{}(\tau ,k), \cdots , \widehat a_{\kappa \kappa }^{}(\tau ,k)]^{\text{T}}$, 将其中的0阶和1阶代入到式(6)—式(8)中得到TFCOA-PIV定位结果${\varTheta _{\text{P}}}(\tau ,k)$; 然后把${{\widehat {\boldsymbol{a}}}}_{{\boldsymbol{\kappa m}}}^{}(\tau ,k)$代入到式(9)中, 以${\varTheta _{\text{P}}}(\tau ,k)$作为空间扫描初始值, 根据式(10)得到TFCOA-AIV定位结果${\varTheta _{\text{A}}}(\tau ,k)$。

(5) 由平滑处理过的直方图给出声源定位结果。将筛选出的每个时频块的TFCOA-PIV或TFCOA-AIV的定位结果构成直方图, 并进行平滑处理^[18,22], 将其最高$L$个峰对应的角度坐标作为$L$个声源的定位结果。

2.4   计算量分析

本节对所提方法的计算量进行分析, 并与PIV^[15]、AIV^[18-19]、OA-AIV^[23]、DPD-AIV^[11,18]方法的计算量进行对比。其中DPD-AIV方法是DPD与AIV方法的结合, 将通过DPD检测的最大特征值对应的特征向量代入到AIV方法中进行定位。计算量由定位每个时频点需要的实数乘法次数进行衡量, 其中复数乘法按照4次实数乘法计算, $| \cdot {|^2}$按照两次实数乘法计算。

各对比方法定位一个时频点所需要的实数乘法次数如表1所示。其中PIV、AIV、OA-AIV、DPD-AIV方法的计算量参考了文献[18, 23]。TFCOA检测因子和OA检测因子的实数乘法次数按照计算全部(N − 1)个因子所需的计算量统计, 其中${P_{{\text{OA}}}}$, ${P_{{\text{DPD}}}}$, ${P_{{\text{TFC}}}}$分别表示OA、DPD和TFCOA检测因子选中该时频点的概率; ${G_{\text{L}}}$表示局部扫描的平均栅格点数, $\overline \kappa$和${\overline \kappa _{{\text{OA}}}}$分别为TFCOA和OA检测因子筛选出的可靠特征波束的平均阶数, 其中$\overline \kappa = \sum\nolimits_{n = 2}^N {nP_{{\text{TFC}}}^n}$, ${\overline \kappa _{{\text{OA}}}} = \sum\nolimits_{n = 2}^N {nP_{{\text{OA}}}^n}$, $P_{{\text{TFC}}}^n$和$P_{{\text{OA}}}^n$分别为TFCOA和OA检测因子筛选出的最大阶数为n的可靠特征波束百分比。

表  1  各对比方法定位一个时频点所需的实数乘法次数

对比方法	分阶筛选	计算相关矩阵	特征值分解	预定位	空间扫描
PIV	0	0	0	48	0
TFCOA-PIV	$[6{(N + 1)^2} + 8N]{J_\tau }{J_k}$	$4{(\overline \kappa + 1)^2}{J_\tau }{J_k}{P_{{\text{TFC}}}}$	$3{(\overline \kappa + 1)^6}{P_{{\text{TFC}}}}$	$48{P_{{\text{TFC}}}}$	0
AIV	0	0	0	48	$8{\left( {N + 1} \right)^2}{G_{\text{L}}}$
OA-AIV	$4{(N + 1)^2} + 6$	0	0	$48{P_{{\text{OA}}}}$	$8{\left( {{{\overline \kappa }_{OA}} + 1} \right)^2}{G_{\text{L}}}{P_{{\text{OA}}}}$
DPD-AIV	0	$4{(N + 1)^2}{J_\tau }{J_k}$	$3{(N + 1)^6}$	$48{P_{{\text{DPD}}}}$	$8{\left( {N + 1} \right)^2}{G_{\text{L}}}{P_{{\text{DPD}}}}$
TFCOA-AIV	$[6{(N + 1)^2} + 8N]{J_\tau }{J_k}$	$4{(\overline \kappa + 1)^2}{J_\tau }{J_k}{P_{{\text{TFC}}}}$	$3{(\overline \kappa + 1)^6}{P_{{\text{TFC}}}}$	$48{P_{{\text{TFC}}}}$	$8{(\overline \kappa + 1)^2}{G_{\text{L}}}{P_{{\text{TFC}}}}$

下载: 导出CSV 

| 显示表格

接下来测试各对比方法定位${\text{1 s}}$时长信号所需的运行时间，如表2所示。 运行时间为各对比定位方法运行100次的平均值, 实验参数设置与仿真实验一致。实验所用计算机的处理器为英特尔i7-10875H, 内存为24 G。

表  2  各对比方法定位1s信号所需的平均运行时间

对比方法	PIV	TFCOA-PIV	AIV	OA-AIV	DPD-AIV	TFCOA-AIV
运行时间 (s)	0.2375	0.5779	1.3149	0.8265	0.7675	0.6097

下载: 导出CSV 

| 显示表格

从表1和表2中可以看出, PIV方法具有闭式解, 其计算量是最低的, 运行时间最短; 而TFCOA-PIV需要对通过TFCOA检测的时频块进行子空间分解, 因此计算量和运行时间高于PIV方法。对于采用高阶特征波束定位的方法, 其中AIV方法需要对每个时频点进行局部空间扫描搜索, 而OA-AIV方法则只需要对通过OA检测的时频点进行AIV定位, 其计算量约是AIV方法的${P_{{\text{OA}}}}$; TFCOA-AIV方法虽然需要对通过TFCOA检测的时频块进行相关矩阵计算和特征值分解, 但是其筛选时频点的概率${P_{{\text{TFC}}}}$比${P_{{\text{OA}}}}$低许多, 使得其计算量低于OA-AIV方法, 定位1 s时长信号所需时间比OA-AIV方法减小了约25%; 而DPD-AIV方法需要对每个时频块进行子空间分解, 因此在${P_{{\text{DPD}}}} = {P_{{\text{TFC}}}}$时其计算量高于TFCOA-AIV方法, 如表1和表2所示。

3.   仿真实验

本节通过仿真实验测试所提声源定位方法的性能, 测试多声源环境下定位性能随信噪比、混响时间以及声源角度间距变化的情况, 并与PIV、AIV、OA-AIV、DPD-AIV方法进行对比, 主要测试低信噪比、高混响条件下的性能表现, 以及阵列存在增益、相位误差时对定位性能的影响。

3.1   仿真参数设置

房间的尺寸设置为6 m × 5 m × 3 m, 阵列中心位于4.4 m × 2.5 m × 1.5 m, 阵列中心到每个声源的距离均为1 m。阵列为刚性球面阵, 半径4.2 cm, 表面均匀分布放置32个传声器。声源个数为3, 放置在48个不同的位置进行测试, 其中声源1的方位角以30°递增分别放置于${\text{\{ }}{0^{\circ}},{30^{\circ}}, \cdots ,{300^{\circ}},{330^{\circ}}{\text{\} }}$, 俯仰角位于$\{{60}^{\circ}, {80}^{\circ},{100}^{\circ}, {120}^{\circ}\}$, 共48种组合; 声源2和声源3的俯仰角与声源1相同, 相互之间方位角间隔30°。每个位置测试10组语音信号, 共进行$48 \times 10{\text{ = }}480$次蒙特卡罗实验。这10组语音信号从TIMIT语音库中随机挑选, 每组信号由3个不同人的语音组成, 从中截取1 s作为测试信号, 信号采样频率为16 kHz。短时傅里叶变换的帧长为256, 帧移为50%, 采用汉宁窗函数进行分帧。最大球谐波阶数N = 3, 为减小球傅里叶变换的空间混叠和模式强度补偿对低频噪声过度放大的影响, 将参与定位的信号频率范围设置为400~5000 Hz。房间混响由基于像−源法的刚性球室内声学单位冲激响应^[26]产生。噪声为高斯白噪声, 不同阵元之间噪声互不相关。

TFCOA因子的筛选阈值按照筛选出总时频点的10%设置, 其中三阶阈值为固定值${T_3} = 0.3$, 若三阶筛选出的时频点不足总时频点的10%时, 则按照${K_2}$从大到小的顺序从剩余时频点中筛选补足, 因此二阶阈值${T_2}$为动态的。DPD-AIV方法的阈值也是按照筛选出总时频点的10%设置^[11,22]。OA-AIV^[23]方法的三阶阈值为0.75, 二阶阈值为0.8, 与文献[23]一致。所有对比方法的定位结果由每个时频点定位结果形成的直方图经高斯平滑处理^[18,22]后的最高$3$个峰的角度坐标给出, 其中高斯函数的方差$\sigma = 4$。AIV、OA-AIV和TFCOA-AIV的区域扫描半径为10°, 扫描栅格的俯仰角和方位角间距均为1°。在参考文献[11, 18, 21, 22]中关于${J_\tau }$, ${J_k}$设置和大量仿真测试的基础上, 将${J_\tau }$设置为3, ${J_k}$设置为7。

声源定位误差用估计角度单位矢量与实际角度单位矢量之间的夹角计算, 若估计的声源DOA为$\{ \widehat \theta ,\widehat \phi \}$, 其单位矢量为$\widehat {\boldsymbol{u}} = {[\sin\widehat \theta \cos\widehat \phi ,\sin\widehat \theta \sin \widehat \phi ,\cos\widehat \theta ]^{\text{T}}}$; 声源真实DOA为$\{ \theta ,\phi \}$, 其单位矢量${\boldsymbol{u}} = [\sin\theta \cos\phi , \sin\theta \sin \phi ,\cos\theta ]^{\text{T}}$; ${\boldsymbol{u}}$和$\widehat {\boldsymbol{u}}$之间的夹角为^[8]

3.2   性能指标

为衡量定位方法的性能采用了两种指标: 平均声源检测个数^[11,21,22]用来衡量定位方法的鲁棒性, 均方根误差用来衡量定位方法的精度。

(1) 平均声源检测个数: 将每次仿真实验的声源定位结果代入到式(22)中计算每个声源的估计误差, 若估计误差小于15°, 则认为检测出声源, 否则认为没有检测出声源。平均声源检测个数由每次实验检测出的声源个数的均值给出, 定义为

其中, ${J_{{\text{mon}}}}$为总的蒙特卡罗实验次数, $L_{\text{D}}^j$为第$j$次实验检测出的声源个数。

(2) 均方根误差: 将成功检测出声源的定位结果代入到式(22)中计算估计误差, 根据估计误差可以得到均方根误差, 有

其中, ${F_{\text{D}}} = {\overline I_{\text{D}}}{J_{{\text{mon}}}}$为成功检测出的声源总数, ${\widehat {\boldsymbol{u}}_{j,l}}$为第$j$次实验检测出的第l个声源的估计DOA单位矢量, ${{\boldsymbol{u}}_{j,l}}$为其对应的真实DOA单位矢量。

3.3   信噪比对所提方法定位性能的影响

本节测试所提方法在不同信噪比条件下的性能表现, 信噪比设置为${\text{SNR}} = \{ 10,15,20,25,30\} {\text{ dB}}$, 混响时间设置为${T_{60}} = 0.7{\text{ s}}$保持不变。各对比声源定位方法的均方根误差和平均声源检测个数随信噪比变化情况如图3所示。

图  3  各对比方法在不同信噪比条件下的均方根误差和平均声源检测个数 (a) 均方根误差; (b)平均声源检测个数

下载: 全尺寸图片 幻灯片

从图3中可以看出在不同信噪比条件下, TFCOA-AIV方法的均方根误差和平均声源检测个数在所有对比方法中的表现是最好的, 且信噪比越低优势越明显, 反映出TFCOA-AIV方法相比其他对比方法具有更强的抗噪能力。其中TFCOA-AIV方法与AIV和OA-AIV方法相比在均方根误差上分别有1.3°~2.1°和1.1°~1.6°的改善, 定位精度提升比较明显。在平均声源检测个数上, TFCOA-AIV在10 dB时平均能检测出2.99个声源, 相比AIV、OA-AIV方法分别提升了8%和4.5%。TFCOA-AIV在信噪比增大到$20{\text{ dB}}$以后能够检测出全部声源, 而AIV和OA-AIV方法在信噪比为$30{\text{ dB}}$时也无法检测出全部声源, TFCOA-AIV具有更高的鲁棒性。TFCOA-AIV与DPD-AIV方法相比, 在均方根误差上约有0.3°~0.5°降低, 在信噪比为$10{\text{ dB}}$时平均声源检测个数提升了约2%。

TFCOA-PIV方法与PIV方法相比, 在均方根误差上有$2.7\text{°} \sim 2.8\text{°}$的降低, 甚至在信噪比为$10{\text{ dB}}$时均方根误差低于采用高阶特征波束定位的AIV方法; 在平均声源检测个数上, 仅次于TFCOA-AIV方法, 优于其他对比方法, 相比PIV有26%~21%的提升。综合可以看出TFCOA因子具有较强的抗噪能力, 能够将受噪声影响较小的可靠特征波束筛选出来。

TFCOA因子的时频点筛选和定位效果如图4所示。图4为信噪比为$10{\text{ dB}}$时一次仿真实验的OA-AIV、DPD-AIV、TFCOA-AIV方法筛选出的各声源时频点分布及纯净语音语谱图。图中按照筛选出的时频点的定位结果归属到每个声源, 若时频点的定位结果与某声源真实DOA的误差在15°以内则认为属于该声源。图中用不同颜色表示不同声源的时频点, 其中蓝色为声源1、紫色为声源2、红色为声源3, 绿色为不属于任何声源的误选时频点, 其中声源2能量主要位于低频部分。从图4(a)中可以看出, OA-AIV筛选出的时频点中误选时频点较多, 主要位于中高频部分, 从中反映出OA方法主要针对低频噪声, 对于高频噪声、混响和多声源相互干扰的鉴别能力有限。从图4(b)和图4(c)中可以看出, TFCOA-AIV和DPD-AIV方法的误选时频点很少, 反映出TFCOA检测因子和DPD检测因子能够有效地将受噪声影响较小和单声源直达声占优的时频点筛选出来用于定位。但是从图4(b)中可以看出, DPD检测因子筛选出的时频点中属于声源2的时频点较少, 原因是声源2主要位于低频段, 其高阶特征波束易受噪声影响, 很难通过DPD检测。而TFCOA检测因子采用分阶处理, 在低频段无法筛选出高阶特征波束的情况下, 将受噪声污染程度较小的低阶特征波束筛选出来用于声源定位, 尽可能地保留了声源2的信息, 如图4(c)中紫色时频点所示, 增加了声源2被定位出的可能性, 这也是低信噪比环境下TFCOA-PIV和TFCOA-AIV方法的平均声源检测个数优于DPD-AIV方法的主要原因, 即TFCOA检测因子具有更强的鲁棒性。

图  4  OA-AIV、DPD-AIV、TFCOA-AIV方法筛选出的各声源时频点分布及纯净语音语谱图 (a) OA-AIV; (b) DPD-AIV; (c) TFCOA-AIV; (d) 纯净语音语谱图

下载: 全尺寸图片 幻灯片

3.4   混响对所提方法定位性能的影响

本节测试所提方法定位性能在不同混响条件下的表现, 其中信噪比为${\text{SNR}} = 10{\text{ dB}}$保持不变, 混响时间为${T_{60}} = \{ 0.4,0.5,0.6,0.7,0.8,0.9,1.0\} {\text{ s}}$。各对比声源定位方法的均方根误差和平均声源检测个数随混响时间的变化情况如图5所示。

图  5  各对比方法在不同混响时间条件下的均方根误差和平均声源检测个数 (a) 均方根误差; (b) 平均声源检测个数

下载: 全尺寸图片 幻灯片

从图5可以看出各对比方法的定位性能均随混响时间的增大而有所降低, 其中性能最优的TFCOA-AIV方法的变化趋势最小, 反映出TFCOA-AIV方法与其他对比方法相比具有更强的抗混响能力。随着混响时间的增大, TFCOA-AIV方法的均方根误差从1.69°增大到1.88°, 只增大了0.2°, 平均声源检测个数基本没变; 而DPD-AIV方法的均方根误差从1.99°增大到了2.26°, 约增大了0.27°, 平均声源检测个数也基本没变; 而AIV和OA-AIV方法的均方根误差分别从3.5°增大到4.2°和从3.03°增大到3.78°, 均增大了$0.7^\circ$左右, 平均声源检测个数分别下降了约4%和2.5%。

对于只利用0, 1阶特征波束定位的方法, TFCOA-PIV方法的均方根误差从3.29°增大到3.93°, 增大了约0.64°, 增大趋势低于OA-AIV和AIV方法, 平均声源检测个数基本没变; PIV方法的均方根误差和平均声源检测个数随混响时间的增大变化是最大的, 说明其抗混响能力是最差的。

通过以上分析可以看出TFCOA-AIV和TFCOA-PIV方法的平均声源检测个数基本不随混响时间的增大而变化, 均方根误差变化也小于或与其他方法持平, 反映出TFCOA因子具有较强的抗混响能力, 这一点可以从图4中进行分析。对比图4(c)中TFCOA-AIV方法筛选出的时频点在图4(d)语谱图中的位置可以看出, 筛选出的时频点主要位于语音信号的起始位置, 根据优先效应可知这部分信号主要是经直接路径传输的直达声信号, 受混响的影响较小, 因此TFCOA因子具有较强的抗混响能力。DPD检测因子筛选出的时频点中有一部分位于语音信号的末尾, 这部分信号以混响为主, 即DPD检测因子会将部分混响占优的信号误选出来, 产生错误的定位结果, 影响了其抗混响能力, 如图4(d)语谱图中红色方框的位置。虽然TFCOA检测因子也会出现误选情况, 但是数量要少于DPD检测因子, 且比较分散, 对定位结果的影响较小, 因此具有更强的抗混响能力。

3.5   声源角度间距对所提方法定位性能的影响

接下来测试不同声源角度间距对各对比方法的定位性能影响, 主要测试空间分辨能力。将信噪比设置为$10{\text{ dB}}$, 混响时间为$0.7{\text{ s}}$, 保持仿真设置中各声源的俯仰角不变, 声源之间的方位角间距进行变化, 方位角间距取值为$\left\{ {30^\circ ,45^\circ ,60^\circ ,75^\circ ,90^\circ } \right\}$。不同声源角度间距条件下各对比方法的均方根误差和平均声源检测个数的变化情况如图6所示。

图  6  各对比方法在不同声源角度间距条件下的均方根误差和平均声源检测个数 (a) 均方根误差; (b) 平均声源检测个数

下载: 全尺寸图片 幻灯片

从图6中可以看出各对比方法的定位性能随声源角度间距的增大而有所提升, 其中TFCOA-AIV、TFCOA-PIV和DPD-AIV方法的均方根误差和平均声源检测个数随声源角度间距的增大变化得比较缓慢, 说明这三种方法受声源角度间距的影响较小, 具有较高的空间分辨率。原因是TFCOA检测因子和DPD检测因子筛选出的时频点绝大部分可以清晰地归属到每个声源, 误选时频点较少, 在形成直方图时每个声源对应的峰比较尖锐, 即使在声源角度间距较小时相互之间的干扰也比较小, 因此随角度间距的增大定位性能变化不大。对于PIV、AIV和OA-AIV方法, 由于采用全部时频点或误选时频点较多导致形成直方图时相邻声源的峰之间会产生相互干扰、淹没等现象, 这种现象会随着声源角度间距的增大而改善, 因此随声源角度间距的增大定位性能有较大程度提高。

3.6   阵列增益相位误差对所提方法定位性能的影响

接下来测试传声器阵列存在增益相位误差时各对比方法的定位性能表现。为突出增益相位误差对定位性能的影响, 将信噪比设置为20 dB, 混响时间为0.3 s, 相位误差和增益误差按均匀分布随机产生, 其中相位误差的取值范围为$\pm 5^\circ$, 增益误差设置了5种取值范围$\left\{ { \pm 0.02,{\text{ }} \pm 0.04,{\text{ }} \pm 0.06,{\text{ }} \pm 0.08,{\text{ }} \pm 0.10} \right\}$。不同增益相位误差条件下各对比方法的均方根误差和平均声源检测个数的变化情况如图7所示。

图  7  各对比方法在不同增益误差条件下的均方根误差和平均声源检测个数 (a) 均方根误差; (b) 平均声源检测个数

下载: 全尺寸图片 幻灯片

从图7中可以看出, 各对比方法的定位性能均随传声器阵列增益误差的增大出现不同程度的下降。其中TFCOA-AIV方法的定位性能仍然是最优的, 其均方根误差随着阵列增益误差的增大约增加了0.2°, 平均声源检测个数几乎没变; 与定位性能次优的OA-AIV方法相比在均方根误差上降低了0.5°~0.7°, 在平均声源检测个数上两者几乎相同; 而TFCOA-PIV方法在阵列增益相位误差的影响下, 其均方根误差和平均声源检测个数均略优于AIV方法。对于DPD-AIV方法, 在没有增益相位误差时其均方根误差和平均声源检测个数分别为1.47°和3, 仅次于TFCOA-AIV方法的1.22°和3, 而当存在增益相位误差时其定位性能迅速恶化, 在增益误差取值范围增大到$\pm 0.1$时, 其均方根误差增大到2.8°, 平均声源检测个数减小到2.91, 定位性能仅优于PIV方法。从中可以看出, TFCOA检测因子与DPD检测因子相比具有更强的抗传声器阵列增益相位误差的能力。

4.   实测实验

本节通过实测实验对所提方法的定位性能进行测试。实测实验环境及球面传声器阵列如图8所示。实验的房间尺寸为9 m × 7 m × 3 m, 房间的混响时间约为0.35 s。阵列采用MH声学公司的em32 Eigenmike刚性球面传声器阵列, 半径为4.2 cm, 阵列表面按均匀分布放置了32个全指向性传声器。阵列的采样频率为44.1 kHz, 为了与仿真实验保持一致, 将采样频率降到16 kHz。阵列中心位于房间的4.5 m × 3.5 m × 1.5 m, 阵列与3个声源的距离均为1.5 m, 共测试了3个位置, 每个位置的声源角度坐标如表3所示(俯仰角在前)。每个位置测试了10组从TIMIT语音库中随机选取的语音信号, 从阵列采集的信号中截取1 s信号进行性能测试。各对比方法的参数设置同仿真实验。各对比定位方法的均方根误差和平均声源检测个数如表4所示。

图  8  实测实验环境和球面传声器阵列

下载: 全尺寸图片 幻灯片

表  3  实测实验声源角度坐标

声源	测试位置1	测试位置2	测试位置3
声源1	(108°, −60°)	(108°, −20°)	(108°, −100°)
声源2	(78°, 40°)	(78°, 40°)	(78°, −20°)
声源3	(98°, 140°)	(98°, 100°)	(98°, 60°)

下载: 导出CSV 

| 显示表格

表  4  实测实验各对比方法定位结果的均方根误差和平均声源检测个数

	PIV	DPD-AIV	TFCOA-PIV	AIV	OA-AIV	TFCOA-AIV
均方根误 差 (°)	4.71	2.35	3.36	2.54	2.23	2.02
平均声源 检测个数	2.30	2.70	2.97	2.77	2.87	3.00

下载: 导出CSV 

| 显示表格

从表4中可以看出, 定位性能表现最好的是TFCOA-AIV方法, 均方根误差和平均声源检测个数均优于其他对比方法; 而TFCOA-PIV方法由于只采用了0阶和1阶特征波束, 使得其均方根误差高于OA-AIV、DPD-AIV和AIV方法, 但是其平均声源检测个数均优于这三种方法; 对于DPD-AIV方法, 其平均声源检测个数仅高于PIV方法, 均方根误差高于OA-AIV方法, 与仿真实验中阵列增益相位误差影响下的性能表现相似, 说明DPD检测因子在实际应用中由于阵列增益相位误差的存在使其筛选直接路径占优时频点的能力下降。PIV方法的定位性能最差。通过实测实验可以看出, TFCOA检测因子在实际环境中的定位表现与仿真实验中一致, 具有较高的定位精度和鲁棒性, 由于实测环境的混响较低、信噪比较高, 使得TFCOA-AIV和TFCOA-PIV方法与其他方法相比优势并不像仿真实验中那么明显。

图9为一次实测实验各对比声源定位方法的归一化直方图。该次实验中声源的真实方位为表3中的测试位置1, 其中声源2的音量较低, 且有语音的时长约为0.4 s, 可视为一个弱声源。在图9中符号“+”为声源的真实方位, 符号“o”为各对比方法估计出的声源方位。从图9中可以看出TFCOA-AIV、TFCOA-PIV和OA-AIV方法成功定位出了3个声源的方位, 而PIV、AIV和DPD-AIV在定位声源2时出现了较大的误差。其中PIV、AIV方法由于没有经过时频点筛选, 其归一化直方图中存在许多由混响噪声等因素产生的虚假峰, 这些虚假峰会影响对声源DOA的估计。PIV方法在定位声源2时就错误定位在最高虚假峰的方位(93°, 139°); AIV方法虽然在声源2的真实方位处形成明显的峰, 但是其归一化高度0.090低于最高虚假峰的归一化高度0.105, 声源2被错误定位在(70°, 154°)。而OA-AIV方法由于采用OA检测因子进行可靠特征波束和时频点筛选, 对噪声和混响进行了一定程度的抑制, 使得最高虚假峰的归一化高度降低到0.071, 低于声源2的峰的高度0.078, 因此将声源2成功检测出来。

图  9  一次实测实验各对比声源定位方法归一化平滑直方图 (a) PIV; (b) TFCOA-PIV; (c) AIV; (d) OA-AIV; (e) DPD-AIV; (f) TFCOA-AIV

下载: 全尺寸图片 幻灯片

对于DPD-AIV方法, 其归一化直方图相比PIV、AIV、OA-AIV方法要干净许多, 但是在定位声源2时也出现了错误, 定位在(92°, 154°), 原因是DPD检测因子在阵列增益相位误差的影响下筛选出许多位于低频段的误选时频点, 影响了DPD-AIV方法的鲁棒性, 如图10(a)所示。图10为该次实测实验DPD-AIV、TFCOA-AIV方法筛选出的各声源时频点分布图。从图10(a)中可以看出DPD检测因子误选的绿色时频点数量远远大于声源2的紫色时频点, 并且分布比较集中, 使得在图9中形成的伪峰高于声源2的峰, 产生了错误的定位结果。在图9中TFCOA-AIV、TFCOA-PIV方法的归一化直方图与其他方法相比非常清晰干净, 原因是TFCOA因子筛选出的时频点可以清晰地归类到每个声源, 如图10(b)所示, 因此形成的3个声源的峰非常尖锐, 具有较高的空间分辨率。另外误选时频点很少, 虽然数量与声源2紫色时频点数量相近, 但是分布较为分散, 形成的伪峰低于声源2的峰, 因此能够将声源2定位出来, 具有较高的鲁棒性。

图  10  一次实测实验DPD-AIV、TFCOA-AIV方法筛选出的各声源时频点分布图 (a) DPD-AIV; (b) TFCOA-AIV

下载: 全尺寸图片 幻灯片

TFCOA检测因子在实际环境中同样能够将受噪声和混响影响较小的单声源直达声占优的时频点和特征波束筛选出来, 具有较强的抗噪和抗混响能力, 使得TFCOA-PIV和TFCOA-AIV方法在实际环境中也具有较高的定位精度和鲁棒性。

5.   结论

本文提出一种利用时频相关的球谐域特征波束阶数感知因子以及基于该因子的伪声强声源定位方法, 理论分析了该阶数感知因子相比已有的阶数感知具有更强的抗噪和抗混响能力, 并通过仿真和实测实验测试了所提声源定位方法在混响噪声环境下的定位性能。针对球谐域高阶特征波束易受噪声污染的问题, 利用相邻多个时频点特征波束之间的相关性进行阶数感知的自适应分阶处理, 在降低计算量的同时提升了定位方法的抗噪能力和抗混响能力。仿真实验和实测实验数据处理结果显示, 利用时频相关性的球谐域阶数感知因子能够将受噪声影响较小的和主要位于语音信号起始位置的特征波束筛选出来, 减小了噪声和混响的影响, 同时对阵列的增益和相位误差也具有较强的抵抗能力, 将其应用于伪声强声源定位中提升了定位的性能, 并且随着信噪比降低、混响时间增大、声源角度间距减小, 性能提升得越明显。

附录A

将式(19)中的内积展开, 并根据不同时频点噪声之间互不相关, 可得

其中, $n' = n - 1$。

考查其中的一项${\text{E}}(n_{nm}^{{j_\tau },{j_k}*}n_{n'm'}^{{j_\tau },{j_k}})$, 根据式(3)将$n_{nm}^{{j_\tau },{j_k}}$, $n_{n'm'}^{{j_\tau },{j_k}}$展开, 可得

其中, $n_q^{{j_\tau },{j_k}}$为第$q$个阵元的$(\tau - {j_\tau },k + {j_k})$时频点的噪声幅度。将上式的求和分成两部分: $q = q'$和$q \ne q'$, 根据不同阵元噪声互不相关, 有

在式(A3)等号右侧第1项中添加$\mu _{}^2\sum\nolimits_{q = 1}^Q {\alpha _q}^2\cdot {{\rm{Y}} _n^m({\varOmega _q}){{\rm{Y}}} {{_{n'}^m}^{'*}}({\varOmega _q})}$, 并在第2项中减去该添加项, 可得

由于${{\rm{Y}}} _0^0(\varTheta ) = {1 / {\sqrt {4\pi } }}$, 则式(A4)等号右侧第1项可以改写为

根据离散采样球谐波函数的正交性可知, 当$n = 0,{\text{ }}m = 0$, $n' = 0,{\text{ }}m' = 0$时式(A5)才不为0, 而$n \ne n'$, 因此式(A5)等于0。

对于式(A4)第2项$\sigma _{}^2\sum\nolimits_{q = 1}^Q {{\alpha _q}^2{{\rm{Y}}} _n^m\left( {{\varOmega _q}} \right){{\rm{Y}}} {{_{n'}^m}^{'*}}\left( {{\varOmega _q}} \right)}$, 若采样点为均匀分布时有${\alpha _q} = {{4\pi } \mathord{\left/ {\vphantom {{4\pi } I}} \right. } Q}$, 同样由于$n \ne n'$有${\alpha _q}\sigma _{}^2\sum\nolimits_{q = 1}^Q {{\alpha _q}{{\rm{Y}}} _n^m\left( {{\varOmega _q}} \right){{\rm{Y}}} {{_{n'}^m}^{'*}}\left( {{\varOmega _q}} \right)} = 0$。若采样点为高斯分布或等角度分布时, 根据其分布特点: 同俯仰角的采样点权值相同、同俯仰角采样点的方位角是等角度间距分布, 由式(2)可知${{\rm{Y}}} _n^m\left( {{\varOmega _q}} \right){{\rm{Y}}} {{_{n'}^m}^{'*}}\left( {{\varOmega _q}} \right)$中包含虚指数信号${{\rm{e}}^{{\rm j}(m - m'){\varphi _q}}}$, 因此同俯仰角的采样点的${{\rm{Y}}} _n^m({\varOmega _q}){{\rm{Y}}} {{_{n'}^m}^{'*}}({\varOmega _q})$相加时, 只有当$m = m'$时不为0, 其他情况为0; 下面分析$m = m'$时, 由于等角度和高斯分布的采样点关于阵列赤道对称且对称采样点权值相同, 对称采样点${\varOmega _q}$和${\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\varOmega } _q}$的俯仰角存在关系${\theta _q} = \pi - {\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\theta } _q}$, 有$\mathrm{Y}_n^m\left(\breve{\varOmega}_q\right) \mathrm{Y}_{n^{\prime}}^{m^*}\left(\breve{\varOmega}_q\right)=(-1)^{n+2 m+n^{\prime}} \mathrm{Y}_n^m\left(\varOmega_q\right) \mathrm{Y}_{n^{\prime}}^{m^*}\left(\varOmega_q\right)$, 由于$n + n' = 2n - 1$, 因此${( - 1)^{n + 2m + n'}} = - 1$, 使得$\mathrm{Y}_n^m\left(\breve \varOmega_q\right) \mathrm{Y}_{n^{\prime}}^{m^*}\left(\breve\varOmega_q\right)+ \mathrm{Y}_n^m\left(\breve{\varOmega}_q\right) \mathrm{Y}_{n^{\prime}}^{m^*}\left(\breve{\varOmega}_q\right)=0$, 因此等角度和高斯分布时$\sigma _{}^2\sum\nolimits_{q = 1}^Q {{\alpha _q}^2{\rm Y} _n^m\left( {{\varOmega _q}} \right){\rm Y} {{_{n'}^m}^{'*}}\left( {{\varOmega _q}} \right)} = 0$。综合以上分析当采样点为均匀分布、高斯分布和等角度分布时有${\text{E}}(n_{nm}^{{j_\tau },{j_k}*}n_{n'm'}^{{j_\tau },{j_k}}) = 0$, 结论同样适用于${\rm{E}}(n_{nm}^{0,0*}n_{n'm'}^{0,0})$, 因此:

附录B

考察式(18)${K}_n^{}(\tau ,k)$的分母中的任一时频点中阶数为n的项, 省略${1 \mathord{\left/ {\vphantom {1 {|{b_n}(k{r_a}){|^4}}}} \right. } {|{b_n}(k{r_a}){|^4}}}$, 有

其中, $n_{nm}^{0,0},{\text{ }}n_{n'm'}^{0,0*}$与$n_{nm}^{{j_\tau },{j_k}*},{\text{ }}n_{n'm'}^{{j_\tau },{j_k}}$属于不同时频点, 是统计独立的, 因此上式可改写为

考查其中一项:

将上式的求和分成两部分: $q = q'$和$q \ne q'$, 根据不同阵元噪声互不相关, 有

由于$n > 0$, 根据式(A4)可知上式第1项为0, 因此${{{\varGamma }}_n}\sum\nolimits_{m = - n}^n {\text{E}}(n_{nm}^{{j_\tau },{j_k}}n_{nm}^{{j_\tau },{j_k}*}) = \sigma^2\sum\nolimits_{q = 1}^Q {{\alpha _q}^2}$; 同理${{{\varGamma }}_n}\sum\nolimits_{m = - n}^n {\text{E}}(n_{nm}^{0,0} n_{nm}^{0,0*}) = \sigma _{}^2\sum\nolimits_{q = 1}^Q {{\alpha _q}^2}$。因此:

对于${\text{E}}({{\varGamma }}_{n - 1}^2||{\boldsymbol{n}}_{{\boldsymbol{n - }}{\text{1}}}^{{\text{0,0}}}|{|^2}||{\boldsymbol{n}}_{{\boldsymbol{n - }}1}^{{j_\tau },{j_k}}|{|^2})$, 由于$n - 1$可以取0, 则n = 1时式(B4)中的第1项不为0。因此:

应用到整个分母, 因此分母的均值为