信号源的定位是传感器阵列信号处理领域重要部分, 曾应用于雷达、声呐等大型通信设备中, 到20世纪80年代, 随着语音通信发展, 基于麦克风阵列的语音信号处理也逐步受到学者们关注, 在会议发言人定位、机器人听觉系统、安防等场景都得到了广泛的研究与应用^[1]. 目前基于小型麦克风阵列的声源定位大多是二维层面, 要实现三维定位则需要增加麦克风组成大规模麦克风阵列采集声源信息, 增加了设备成本和计算成本. 因此, 如何利用小型麦克风阵列进行声源的三维实时定位成了重要的研究课题.

在目前的研究中, 近场声源定位算法主要分为3类: (1)基于高分辨谱估计的定位算法. 该算法主要是针对窄带信号, 若要运用到语音信号中, 需要将宽带的语音信号先分解为窄带信号, 因转换损失会造成偏差^[2,3]. (2)基于可控响应功率的定位算法. 该算法对宽带信号和窄带信号都适用, 但需要对整个空间的网格点进行搜索和功率计算, 计算量大导致响应速度慢, 不适用于实时的声源定位系统^[4]. (3)基于到达时间差的定位算法. 该算法分为两个部分, 第1步是时延估计, 第2步是定位估计, 计算量小, 实时性高, 在时延估计阶段有较高要求, 否则会在定位估计中会放大时延估计的误差^[5].

针对SRP-PHAT算法响应速度慢的缺陷, 很多学者提出了改进方法. Cho等^[6]提出了针对小型麦克风阵列的SSC加速算法, 利用相同到达时间差聚类取质心的方法减少搜索网格, 但仅对俯仰角和方向角进行了计算, 缺少测距, 且聚类后的质心数目仍然; Yook等^[7]在SSC优化算法的基础上, 提出了针对大型麦克风阵列的两级搜索空间聚类(Two-Level Search Space Clustering, TL-SSC)算法, 将声源的候选位置划分成组, 找到少数可能包含最大功率的组, 进一步减少SSC算法的搜索网格, 但对于小型麦克风阵列的计算量减小效果不明显; Salvati等^[8]利用搜索网格上TDOA信息的非均匀积累提出了基于灵敏度的区域选择SRP (Region selection SRP, R-SRP)算法以减少声源定位误差, 但计算量较多, 实时性差.

本文介绍了阵列信号模型和传统定位算法, 改进了TDOA定位算法以提高定位准确性, 改进了SRP-PHAT-SSC定位算法以简化聚类逻辑, 同时针对小型麦克风阵列三维定位算法实时性不足的问题, 提出了TDOA定位算法与SRP-PHAT-SSC算法相结合的搜索空间收缩聚类算法, 并在Matlab平台对3种SRP-PHAT搜索方式进行了仿真比较. 算法拟针对五元小型麦克风阵列, 在SRP-PHAT-SSC加速算法的基础上进一步减少搜索网格数, 提高SRP-PHAT算法的实时性和三维定位精度.

1 阵列信号模型的建立

阵列模型的建立是声学研究的基础. 本文针对室内的声源定位, 声场是近场模型, 声波以球面波的形式扩散, 图1为近场中一个声源到两个麦克风的传播示意图^[9], 其中, Source表示声源; Mic₁和Mic₂表示两个麦克风; τ为声源信号到达两个麦克风的时延, 单位为s; c为声速, 单位为m/s; 声源到两个麦克风的距离差d为τ×c.

根据图1, 时延为τ的两个麦克风接收到的信号模型x_i(t)(i=1, 2)可表示为^[10]:

${x_i}\left( t \right) = {h_i}*s\left( {t - {\tau _i}} \right) + {n_i}\left( t \right)$

(1)

其中, s(t)为声源信号; h_i为声源到麦克风i之间的脉冲响应函数, 由式(2)表示; τ_i为声源到达麦克风i的时间, τ=τ₁−τ₂; n_i(t)为第i个麦克风的环境噪声信号, s(t)、n₁(t)、n₂(t)之间互不相关.

${h_i}\left( {r,{r_s},t} \right) \!\!=\!\!\! \mathop \sum \limits_{p \in P} \!\mathop \sum \limits_{\lambda \in \Lambda} \beta _{{x_1}}^{\left| {{\lambda_x} \!-\! q} \right|}\beta _{{x_2}}^{\left| {{\lambda_x}} \right|}\beta _{{y_1}}^{\left| {{\lambda_y} \!-\! j} \right|}\beta _{{y_2}}^{\left| {{\lambda_y}} \right|}\beta _{{{\textit{z}}_1}}^{\left| {{\lambda_{\textit{z}}} \!-\! k} \right|}\beta _{\textit{z}}^{\left| {{\lambda_{\textit{z}}}} \right|}\frac{{\delta \left( {t - R/c} \right)}}{{4\pi R}}$

(2)

其中, r表示麦克风i的位置坐标; r_s表示声源位置坐标; t表示时间; c为声速; β表示一个房间6面墙的反射系数; R表示虚源的位置坐标; δ{∙}表示虚源R到麦克风i的直达声压, 可替代为使用Hanning窗的理想低通滤波器, 窗长为4 ms, 截止频率为奈奎斯特频率; P={(q, j, k) | q, j, kϵ{0, 1}}; Λ={(λ_x, λ_y, λ_z) | λ_x, λ_y, λ_z∈Z: −N≤ λ_x, λ_y, λ_z≤N}, N表示虚源的级数; R=((1−2q)x_s+2λ_xL_x−x, (1−2j)y_s+2λ_yL_y−y, (1−2k)z_s+2λ_zL_z−z).

在实际使用中, 一般通过设定混响时间RT₆₀来改变声源和麦克风i之间的脉冲响应, 根据Sabin经验公式可以由混响时间求出式(2)中的墙壁反射系数β. 混响时间和墙壁反射系数的关系为:

$R{T_{60}} = \frac{{24\ln\left( {10} \right)V}}{{c\displaystyle\mathop \sum \nolimits_{i = 1}^6 {S_i}\left( {1 - \beta _i^2} \right)}}$

(3)

其中, V表示房间的体积, S_i表示第i面墙的面积.

本文使用一个五元全向型麦克风阵列进行仿真, 以一个中型大小的会议室为场景, 房间大小为9 m×8 m×3 m; 5个麦克风的位置为M₁(2, 2, 1.6)、M₂(2, 1.9, 1.5)、M₃(2, 2.1, 1.5)、M₄(1.9, 2, 1.5)、M₅(2.1, 2, 1.5), 阵列结构如图2所示. 声速c设定为340 m/s; 声源使用Noizeus中一段3 s的纯净男声, 采样频率为8 kHz; 声源坐标设为(1.6, 2.3, 1.7), 在Matlab平台根据式(1)模拟麦克风接收到的信号, 以加性高斯噪声作为环境噪声. 对麦克风模拟信号进行带通滤波和分帧加窗预处理, 窗函数选择典型的Hamming窗, 帧长取256点, 帧移为50%.

2 定位算法的分析与改进 2.1 TDOA定位算法 2.1.1 广义互相关法时延估计

TDOA定位算法的第一步是时延估计, 广义互相关法(Generalized Cross-Correlation, GCC)是一种最常用的时延估计算法^[11,12]. 通过最大化一对麦克风信号在频域上的互相关函数可估计出声音信号到达这对麦克风的时间差^[5].

信号x₁与x₂之间的时延 ${\tau _{{x_1}{x_2}}}$ 定义为:

${\tau _{{x_1}{x_2}}} = \mathop {{\rm{argmax}}}\limits_\omega {R_{{x_1}{x_2}}}\left( \tau \right)$

(4)

信号x₁与x₂之间的广义互相关函数 ${R_{{x_1}{x_2}}}\left( \tau \right)$ 为:

${R_{{x_1}{x_2}}}\left( \tau \right) = \mathop \int \limits_{ - \infty }^\infty {G_{{x_1}{x_2}}}\left( \omega \right){e^{j\omega \tau }}{\rm{d}}\omega $

(5)

${G_{{x_1}{x_2}}}\left( \omega \right) = {X_1}\left( \omega \right)X_2^*\left( \omega \right)$

(6)

其中, ω为频率; j=−2πi; ${G_{{x_1}{x_2}}}\left( \omega \right)$ 为x₁与x₂的互功率谱; X₁(ω)与X₂(ω)分别是x₁与x₂经过快速傅里叶变换(FFT)后的频域信号; (*)表示复共轭运算.

由于混响和噪声容易对广义互相关函数产生影响, 为了减少混响和噪声的影响, 本文引入PHAT(PHAse Transform, 相位变换)加权函数对互功率谱函数作白化处理, PHAT加权对混响有一定的抑制效果^[13,14]. 采用一对麦克风进行时延估计仿真, GCC-PHAT函数与GCC函数的对比如图3所示, GCC-PHAT函数只保留了相位信息, 从而锐化互相关函数的峰值. PHAT加权函数的公式为:

${\phi _{{\rm{PHAT}}}} = \frac{1}{{\left| {{G_{{x_1}{x_2}}}\left( \omega \right)} \right|}}$

(7)

则式(5)引入PHAT加权可写为:

${R_{{x_1}{x_2}}}\left( \tau \right) = \mathop \int \nolimits_{ - \infty }^\infty \frac{{{G_{{x_1}{x_2}}}\left( \omega \right)}}{{\left| {{G_{{x_1}{x_2}}}\left( \omega \right)} \right|}}{e^{j\omega \tau }}{\rm{d}}\omega $

(8)

2.1.2 球形插值法定位估计

定位估计是TDOA定位算法的第2步, 为了保证TDOA定位算法的实时性, 本文选用球形插值法(Spherical Interpolation, SI)用于定位估计. 球形插值法属于最小二乘(Least Square, LS)算法, 复杂度低, 基本思想是以麦克风阵列中的一个麦克风为参考麦克风, 通过最小化声源到参考麦克风与其他麦克风之间距离差的估计值与实际值的误差平方和, 求出声源位置^[15,16].

假设有N个麦克风的坐标(M₁, M₂,…, M_N), 若以M₁为参考麦克风, 并作为空间坐标原点, 此时, M_i的坐标为 ${{\boldsymbol{q}}_i}$ , 声源s到M_i与M₁的估计距离差为 $\hat {{d_i}}$ , $\hat {{d_i}}$ 是根据第一步的估计时延τ得到的值, M_i到坐标原点的距离为R_i, 声源位置坐标为 ${{\boldsymbol{q}}_s}$ , 则M_i到原点距离的平方可表示为: $d_i^2 = {({R_s} + \hat {{d_i}})^2} = ||{{\boldsymbol{q}}_i} - {{\boldsymbol{q}}_s}|{|^2} = R_i^2 + R_s^2 - 2{\boldsymbol{q}}_i^{\rm T}{{\boldsymbol{q}}_s}$ .

由于估计值 $\hat {{d_i}}$ 存在误差, 所以估计距离与实际距离的最小二乘误差函数可写为:

${\boldsymbol{\varepsilon }} = 2{{\boldsymbol{R}}_s}{\hat{\boldsymbol D}} + 2{\boldsymbol{M}}{{\boldsymbol{q}}_s} - {\boldsymbol{\delta }}$

(9)

其中, ${\boldsymbol{M}} = \left[\!\! {\begin{array}{*{20}{c}} {{x_2}}&{{y_2}}&{{{\textit{z}}_2}} \\ {{x_3}}&{{y_3}}&{{{\textit{z}}_3}} \\ \vdots& \vdots & \vdots \\ {{x_N}}&{{y_N}}&{{{\textit{z}}_N}} \end{array}} \!\!\right]$ , $\hat{{\boldsymbol{D}}}=\left[\!\!\begin{array}{c}\hat{{d}_{2}}\\ \hat{{d}_{3}}\\ \vdots \\ \hat{{d}_{N}}\end{array}\!\!\right]$ , ${\boldsymbol{\delta}} =\left[\!\!\begin{array}{c}{R}_{2}^{2}-{\hat{{d}_{2}}}^{2}\\ {R}_{3}^{2}-{\hat{{d}_{3}}}^{2}\\ \vdots \\ {R}_{N}^{2}-{\hat{{d}_{N}}}^{2}\end{array}\!\!\right]$ .

式(9)可简化为:

$ {\boldsymbol{\varepsilon }} = {\boldsymbol{A}\mu } - {\boldsymbol{b}} $

(10)

其中, ${\boldsymbol{A}} = [{\boldsymbol{M}}|{\hat{\boldsymbol D}}]$ , ${{\mu }} = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{q}}_s}} \\ {{{\boldsymbol{R}}_s}} \end{array}} \right]$ , ${\boldsymbol{b}} = \dfrac{1}{2}{\boldsymbol{\delta }}$ .

则式(10)的最小二乘解可表示为:

${\hat{\boldsymbol \mu }} = {({{\boldsymbol{A}}^{\rm T}}{\boldsymbol{A}})^{ - 1}}{{\boldsymbol{A}}^{\rm T}}{\boldsymbol{b}}$

(11)

定义投影矩阵 ${{\boldsymbol{P}}_{{d}}} = { I} - \dfrac{{{\hat{\boldsymbol D}}{{{\hat{\boldsymbol D}}}^{\rm T}}}}{{{{{\hat{\boldsymbol D}}}^{\rm T}}{\hat{\boldsymbol D}}}}$ , 代入式(11)可得估计位置:

${{\hat{\boldsymbol q}}_s} = {({{\boldsymbol{M}}^{\rm T}}{{\boldsymbol{P}}_{{d}}}{\boldsymbol{M}})^{ - 1}}{{\boldsymbol{M}}^{\rm T}}{{\boldsymbol{P}}_{{d}}}{\boldsymbol{b}}$

(12)

若估计声源的位置与实际位置中的方向角φ相差∆φ度、俯仰角θ相差∆θ度、径向距离r相差∆r米, 则记为异常值, 否则记为正确值. 本文使用准确率作为评判准则: 准确率=正确值的帧数和/总帧数.

由于GCC-PHAT算法的时延估计误差会在球形插值法中放大, 定位错误率较高, 所以本文引入离群值校正, 对基于TDOA算法的定位结果进行校正处理, 将离群值校正为中位数, 其中, 离群值为中位数绝对偏差(Median Absolute Deviation, MAD)与MAD中位数相差超过3倍的值.

设置方向角、俯仰角和径向距离的容错范围分别是∆φ=15、∆θ=15、∆r=1. 在5元锥形麦克风阵列的情况下, TDOA定位算法无法准确计算出俯仰角, 经过仿真实验得知, 俯仰角的准确率均在20%以下. 当无混响影响时, 如图4所示, 方向角和径向距离准确率随信噪比的减少而下降, 定位结果容易受到噪声影响, 做了离群值校正后, TDOA定位算法的鲁棒性增强, 准确率都在85%以上; 当只有混响影响时, 由于PHAT加权函数对混响的抑制作用, 所以定位结果都在容错范围内, TDOA在方向角和径向距离两个方面的准确率均为100%, 不受混响时间的影响.

2.2 SRP-PHAT定位算法

SRP-PHAT算法属于一步定位的算法, 其PHAT滤波被用于提高SRP算法鲁棒性^[17]. 通过对定位空间划分格子, 每一个格子进行所有麦克风对之间互相关函数的求和来精确功率的计算, 从而确定最大功率的空间点作为估计声源的位置^[18].

定义q为空间点位置坐标, ${R_{lm}}\left[ {{\tau _{lm}}\left( q \right)} \right]$ 为麦克风l和麦克风m所接收信号的 GCC-PHAT 函数, 则可控响应功率P(q)为:

$P\left( q \right) = \mathop \sum \limits_{l = 1}^N \mathop \sum \limits_{m = l + 1}^N {R_{lm}}\left[ {{\tau _{lm}}\left( q \right)} \right]$

(13)

那么, 估计声源的位置就在可控响应功率值最大的空间网格:

${{\hat{\boldsymbol q}}_s} = \mathop {{\rm{argmax}}}\limits_q P\left( q \right)$

(14)

可控响应功率取决于阵列中每个麦克风对的接收信号在频域上的相位差, 而频谱的相位差在时域上是离散的, 所以如果有两个空间点的位置很接近, 则这两个空间点到各个麦克风对的时间差相同, 那么就无法依靠式(14)求出声源位置.

SSC算法是一种针对基于SRP-PHAT的小型麦克风阵列声源定位而提出的加速算法^[6]. SSC算法将具有相同TDOA的空间点集合在一起形成集群, 然后求出每个集群的质心, 将每个集群的质心作为代表坐标存储在一个查找表中, 从而避免了上述问题的出现, 同时进行声源定位时大大减少了功率计算的时间, 因为对于小型麦克风阵列, 最终组成的集群的数目远小于初始的空间网格总数.

文献[6]中使用的是自顶向下的聚类方法构造查找表, 本文使用逻辑与实现更为简单的自底向上聚类方法, 其算法流程如算法1.

算法1. 自底向上聚类算法

1)以球面坐标(φ, θ, r)划分空间网格, 形成网格点坐标集合{b};

2)遍历{b}, 计算每个网格点中每组麦克风对的TDOA, 存储到TDOA查找表;

3)将相同TDOA值的网格点划分成簇;

4)计算每个簇的质心, 形成质心集合{c};

5)遍历{c}, 计算每个集群质心的TDOA, 存储到TDOA_SSC查找表;

6)计算每个质心的可控响应功率;

7)功率最大的质心位置为估计声源位置.

2.3 搜索空间收缩聚类算法

本文改进后的TDOA算法提高了定位计算的鲁棒性, 在五元锥形麦克风阵列且信噪比大于10 dB的情况下, 可在一定容错范围内得到正确的方向角和径向距离, 但无法计算得到俯仰角.

SRP-PHAT-SSC算法需要对室内空间的每个集群质心进行可控功率计算, 数量相比SRP-PHAT的全网格搜索虽然已经减少, 但算法的实时性仍然不高, 可进一步缩减搜索网格数目以减少可控响应功率的计算量, 提高算法实时性.

因此, 针对小型麦克风阵列的三维定位, 本文提出一种基于TDOA和SRP-PHAT-SSC的组合算法−搜索空间收缩聚类算法(Search Space Shrinking Clustering, SSSC), 通过引入改进的TDOA算法以进一步减少SRP-PHAT-SSC算法中的候选声源位置, 增加算法的鲁棒性并提高实时性. 该算法主要分为两个阶段, 粗定位阶段和细定位阶段, 整体流程如图5所示.

首先使用自底向上聚类的SSC算法得到整个空间的集群质心位置集合B并保存, 用于提供网格查询的功能, 采取5 cm的细粒度网格以提高SRP-PHAT-SSC算法的定位精度. 首先是粗定位阶段, 对音频信号帧进行基于TDOA的定位计算, 先利用式(8)求得估计时延, 再利用式(12)求得估计声源坐标, 得到估计声源位置集合Q_s, 对Q_s进行离群值校正以提高TDOA定位的鲁棒性, 由于基于五元锥形麦克风阵列的TDOA定位算法无法准确判断俯仰角, 所以取其方向角 ${\hat \varphi _{\rm TDOA}}$ 和径向距离 ${\hat r_{\rm TDOA}}$ 作为空间收缩的参考. TDOA定位的容错范围设为∆φ=15和∆r=1, 根据设置的TDOA定位容错范围进行搜索空间收缩以减少可控功率的计算量, 收缩后的空间B_s为 $\{ (\varphi ,\theta ,r)|\varphi \in [{\hat \varphi _{{\rm{TDOA}}}} - \Delta {\varphi _{{\rm{TDOA}}}}, {\hat \varphi _{{\rm{TDOA}}}}{\rm{ \!+\! }}$ $\Delta {\varphi _{{\rm{TDOA}}}}],\theta \! \in \! {B_\theta },r \! \in \! [\!{\hat r_{{\rm{TDOA}}}} \!-\! \Delta {r_{{\rm{TDOA}}}},{\hat r_{{\rm{TDOA}}}} \!+\! \Delta {r_{{\rm{TDOA}}}}\! ]\} $ . 之后进入细定位阶段, 根据B_s对集群质心坐标集合B进行筛选, 仅对收缩空间B_s内的集群质心使用式(13)进行可控响应功率计算, 之后根据式(14)取功率最大的质心位置作为声源.

3 仿真实验与结果分析

本文仿真实验中, SRP-PHAT算法的全网格搜索方法和搜索空间聚类方法的网格尺寸分别取径向距离的精度为20 cm、10 cm和5 cm, 搜索空间收缩聚类算法取5 cm的细粒度网格尺寸, 俯仰角和方向角的精度为1°, 设置∆φ=10、∆θ=10、∆r=1, 若最终估计声源位置的球坐标在该设定范围内, 则标记为正确定位.

表1以信噪比为30 dB, 混响时间为0 s的情况, 对3种算法进行了对比, 根据表1可以看出: 本文提出的SSSC算法搜索的格子数相比SSC算法减少了96%, 从而使得计算时间比SSC算法减少了97%, 具有较好的实时性, 且定位准确率达95%.

图6和图7是对3种算法在球坐标系的3个维度上分别进行准确率判断, 其中, FGS和SSC算法选取与搜索空间收缩聚类算法相同网格粒度的数据作为展示. 可以看出, 由于PHAT加权函数对混响的抑制作用, 3个算法的准确率在仅有混响的情况下几乎不变. SSC算法无法较好的估计出声源的径向距离, 改进后的算法通过引入TDOA算法计算出准确的径向距离, 从而弥补了SSC算法的缺点. 在低信噪比的情况下, FGS和SSC算法的准确率降低, 而改进后的算法在3个维度上的准确率都保持在80%以上, 这是因为其在TDOA定位阶段进行了离群值校正, 从而提高了鲁棒性.

4 结论

本文分析了基于TDOA的声源定位算法和SRP-PHAT声源定位算法, 在二者的基础上提出了基于SRP-PHAT的搜索空间收缩聚类优化算法, 本文选取了不同的网格步长、信噪比和混响时间, 对全空间搜索、搜索空间聚类和搜索空间聚类算法进行了对比分析, 从实验结果可以看出搜索空间聚类算法减少了需要进行功率计算的网格数, 从而减少了运行时间, 具有较好的实时性, 并且在球坐标的3个维度都有着高准确率, 在实际应用中有一定的价值.