使用麦克风阵列对声源定位的方法

摘要

提供一种使用麦克风阵列对声源定位的方法，三个麦克风构成正三角形，该方法包括：建立坐标系，其原点与正三角形的重心重合，第一麦克风位于坐标系的纵轴上；将正三角形的重心与正三角形的三个顶点连接并延长，将全圆周分为6个相等的区间；计算声源分别到达第一至第三麦克风的时间差；将时间差乘以音速，得到声源到第二麦克风的距离与声源到第一麦克风的距离之间的第一距离差以及声源到第三麦克风的距离与声源到第一麦克风的距离之间的第二距离差；根据第一和第二距离差的符号以及第一、第二距离差的绝对值大小比较关系，初步确定声源相对于正三角形重心的方位角Φ所落入的区间；采用0.618法对声源与正三角形重心的距离R以及方位角Φ逼近优化点。

说明书

技术领域

本发明涉及声源的定位，更具体地讲，涉及一种使用麦克风(MIC)阵列来对声源定位的方法。

背景技术

阵列信号处理已经广泛应用于通信、雷达、声纳、医学和航天航空等诸多领域。近年来，随着反恐战的开展，各国已经研制了各种探测爆炸点或枪炮射击点的位置的设备。在实际应用中，由于视觉定位往往不易实现，从而发展了一些声学定位的系统，例如，美国的回力棒系统、以色列拉斐尔公司的反狙击手声探侧系统，加拿大麦克唐纳·迪特维利公司和加拿大防务研究和发展委员会也共同开发了一种名为“雪貂”(Ferret)的小型武器探测和定位系统。Microsoft公司开发的新的Vista操作系统集成了对MIC阵列的支持，Intel公司的HDAudoi规范可以实现16个MIC、32KHz采样，这些为声波信号阵列处理提供了良好的支持。

目前，使用麦克风阵列来对声源进行定位的方法总体来说可以分为两类：定向波速形成和时间延迟估计。定向波速形成是基于最大输出功率可控波束方法的一种形式，该方法对麦克风阵列接收到的语音信号进行滤波延迟并加权求和形成特定指向的波束，然后计算各个指向上的波束能量，其中搜索空间中功率最大的方向就是对声源位置的一个最大近似。在基于时间延迟估计的方法中，由于声源发出的信号到达各个麦克风存在时间差，该方法根据各个通道信号间的相关关系对这个时间延迟进行估计，从而来对声源位置进行估计。

定向波速形成的方法为了搜索能够使波束功率最大的位置，需要对整个空间进行扫描，需要进行大量计算，因此很难满足实时的要求。

基于时间延迟估计的方法计算速度快，能达到实时的要求。这类方法可以适用于如同语音这样的非平稳宽带信源，但是在只有单个信源时能较好地工作，定位的精度与到达时间差(TDOA)的估计精度有关。

发明内容

本发明的目的在于提供一种使用麦克风阵列对声源定位的方法，该方法能够快速对声源定位。

根据本发明的一方面，提供一种使用麦克风阵列对声源定位的方法，所述麦克风阵列为构成正三角形的三个麦克风，所述方法包括：建立坐标系，所述坐标系的原点与所述正三角形的重心重合，三个麦克风中的第一麦克风位于坐标系的纵轴上；将正三角形的重心与正三角形的三个顶点连接并延长，从而将以所述正三角形的重心为圆心的全圆周分为6个相等的区间；计算声源分别到达三个麦克风中的第一麦克风、第二麦克风、第三麦克风的时间差；将计算的时间差乘以音速，得到声源到第二麦克风的距离与声源到第一麦克风的距离之间的第一距离差以及声源到第三麦克风的距离与声源到第一麦克风的距离之间的第二距离差；根据第一距离差和第二距离差的符号以及第一距离差、第二距离差的绝对值大小比较关系，初步确定声源相对于正三角形重心的方位角Φ所落入的区间；采用0.618法对声源与正三角形重心的距离R以及方位角Φ逼近优化点，获得最终估计的距离R和方位角Φ。

根据本发明的另一方面，提供一种使用麦克风阵列对声源定位的方法，所述麦克风阵列为构成正三角形的三个麦克风，所述方法包括：建立坐标系，所述坐标系的原点与所述正三角形的重心重合，三个麦克风中的第一麦克风位于坐标系的纵轴上；计算声源分别到达三个麦克风中的第一麦克风、第二麦克风、第三麦克风的时间差；将计算的时间差乘以音速，得到声源分别到达第一麦克风、第二麦克风、第三麦克风的距离差；根据三角形边角关系建立关于声源到各个麦克风的距离差的联立方程组；对联立方程组中的任意两个方程进行联立求解，得到声源与正三角形重心的距离R以及声源相对于正三角形重心的方位角Φ。

附图说明

通过结合附图，从下面的实施例的描述中，本发明这些和/或其它方面及优点将会变得清楚，并且更易于理解，其中：

图1是根据本发明的采用正三角形麦克风阵列对声源定位的示意图；

图2是根据本发明第一实施例的采用正三角形麦克风阵列对声源定位的方法的流程图；

图3是根据本发明的采用时间延迟特性进行方位分割的示意图；

图4是根据本发明第二实施例的采用正三角形麦克风阵列对声源定位的方法的流程图；

具体实施方式

本发明针对基于时间延迟定位声源的方法进行了改进和创新。在本发明中，采用了准L1相关技术提高时间延迟估计的速度和抗野值干扰能力，进一步采用时间延迟特性进行方位分割以提高计算速度和去除局部极值点，同时可采用0.618法加快逼近速度，从而可以用普通微型计算机在数十毫秒级解决声源定位问题，满足实用性的速度和精度要求。

以下，参照附图来详细说明本发明的原理以及实施例。在本发明中，采用由三个麦克风构成的麦克风阵列来对声源定位，三个麦克风被布置为正三角形，即，三个麦克风分别位于正三角形的三个顶点上。

图1是根据本发明的采用正三角形麦克风阵列对声源定位的示意图，图2是根据本发明第一实施例的采用正三角形麦克风阵列对声源定位的方法的流程图。

参照图1和图2，在步骤201，建立坐标系，坐标系的原点O(0，0)与第一麦克风a、第二麦克风b、第三麦克风c构成的正三角形的重心重合，三个麦克风中的一个麦克风(例如第一麦克风a)位于坐标系的纵轴上。

在步骤202，计算声源p分别到达第一麦克风a、第二麦克风b、第三麦克风c的时间差。

在步骤203，将计算的时间差乘以声音在空气中传播的速度(音速)，得到声源p分别到达第一麦克风a、第二麦克风b、第三麦克风c的距离差。

在计算所述时间差的步骤中，可采用准L1算法估计所述时间差。普通相关函数可以认为是2范数下的测度，L1相关是1范数下的测度，L1相关与L2相关有几乎相同的各种性质。但是，L1相关比L2相关容易计算，且抗干扰能力更强。

两个麦克风处的两个声音信号x(n)和y(n)的L1互相关函数定义如下：

$\mathrm{Rxy}\left(m\right)=\underset{n=1}{\overset{N}{\Sigma }}((\mathrm{abs}\left(x\left(n\right)\right)+\mathrm{abs}\left(y\left(n\right)\right))/2)-\underset{n=1}{\overset{N}{\Sigma }}\left(\mathrm{abs}(x\left(n\right)-y(n+m))\right)$

其中，n表示声音信号的采样时间点，N表示声音信号的总采样点数，m表示时间偏移量。m是采样间隔的整数倍。

由于上式的第一项对平稳随机过程是常数，因此选择最大的互相关值Rxy(m)，等效于选择最小的第二项，这就是准L1互相关函数。

对于本发明，对声音信号进行采样，已知两个麦克风处的两个声音信号x(n)和y(n)的波形基本相同但时间延迟不同，针对上述互相关函数改变m值，使得两个声音信号x(n)和y(n)的互相关值Rxy(m)最大，即，调整x(n)与y(n)错位对齐，此时的m值即是两个声音信号x(n)和y(n)的时间差。另外，本领域技术人员应该理解，在上述互相关函数中，abs()表示绝对值函数，x(n)和y(n)表示声音信号的幅度。

互相关比对求时间差的方法(即，准L1算法)虽然可以应用所有采样点的延迟信息，但因为数字化，所以比对只是在整采样点处进行，而最佳延迟量一般不在整采样点处。对此问题可以用提高采样率的方法提高估计精度。采样率提高，导致采样间隔较密，可以更准确地估计延迟时间差。但采样率受硬件限制，存储设备也受硬件限制。因此，在没有提高采样率的情况下，在本发明中还可以采用三点二次曲线插值的方法估计两个采样点之间的值，再进行互相关法估计，从而提高了估计精度。

如图1所示，正三角形的重心为坐标系的原点O(0，0)。正三角形的边长为L，声源p到各个麦克风a、b、c的距离差为：

D1＝R2-R1(1)

D2＝R3-R1(2)

D3＝R3-R2(3)

声源p的极坐标为(R，Φ)，R是声源p与正三角形的重心的距离，Φ是声源p相对于正三角形的重心的方位角，R和Φ是待求未知数。

三个麦克风a、b、c的顶点极坐标分别为麦克风a：(r，π/2)，麦克风b：(r，7π/6)，麦克风c：(r，-π/6)。

另外，参照图2，各个麦克风a、b、c与正三角形的重心的距离r有下述关系：

在步骤204，根据三角形边角关系建立关于声源p到各个麦克风a、b、c的距离差的联立方程组。

参照图2，根据三角形边角关系可以得出：

$R1=\sqrt{R^2+r^2-2\times R\times r\times \cos \left(\Phi \right)}---\left(5\right)$

$R2=\sqrt{R^2+r^2-2\times R\times r\times \cos (\frac{2\pi }{3}-\Phi )}---\left(6\right)$

$R3=\sqrt{R^2+r^2-2\times R\times r\times \cos (\Phi +\frac{2\pi }{3})}---\left(7\right)$

将方程(4)-(7)代入方程(1)-(3)，可以得到声源p到各个麦克风a、b、c的距离差的联立方程组如下：

$D1=\sqrt{R^2+r^2-2\times R\times r\times \cos (\frac{2\pi }{3}-\Phi )}-\sqrt{R^2+r^2-2\times R\times r\times \cos \left(\Phi \right)}---\left(8\right)$

$D2=\sqrt{R^2+r^2-2\times R\times r\times \cos (\Phi +\frac{2\pi }{3})}-\sqrt{R^2+r^2-2\times R\times r\times \cos \left(\Phi \right)}---\left(9\right)$

$D3=\sqrt{R^2+r^2-2\times R\times r\times \cos (\Phi +\frac{2\pi }{3})}-\sqrt{R^2+r^2-2\times R\times r\times \cos (\frac{2\pi }{3}-\Phi )}---\left(10\right)$

在步骤205，对联立方程组中的任意两个方程进行联立求解，得到声源p的极坐标，从而完成对声源S的定位。

例如，对三个方程(8)-(10)中的任意两个方程进行联立求解，均可以解出方位角Φ和距离R的值，得到声源p的极坐标(R，Φ)。

由于实际测量误差，还可以进一步采用方程(8)(9)、(9)(10)、(8)(10)分别联立求解，分别求出方位角Φ和距离R的三组解，再取三组解的平均值作为方位角Φ和距离R的估计值。

图3是根据本发明的采用时间延迟特性进行方位分割的示意图，图4是根据本发明第二实施例的采用正三角形麦克风阵列对声源定位的方法的流程图。

在步骤401，建立坐标系，坐标系的原点O(0，0)与第一麦克风a、第二麦克风b、第三麦克风c构成的正三角形的重心重合，三个麦克风中的一个麦克风(例如第一麦克风a)位于坐标系的纵轴上。

为了防止逼近算法进入局部最小值，考虑时间差的符号和大小，在步骤402，将正三角形的重心与正三角形的三个顶点连接并延长，从而将以正三角形的重心为圆心的全圆周分为6个相等的区间1-6，如图3所示。

在步骤403，计算圆周上任一声源p分别到达第一麦克风a、第二麦克风b、第三麦克风c的时间差。

在第一实施例中描述的准L1算法以及三点二次曲线插值的方法同样可应用于第二实施例中，这里不再重复描述。

在步骤404，将计算的时间差乘以声音在空气中传播的速度(音速)，得到声源p到第二麦克风b的距离与声源p到第一麦克风a的距离之间的第一距离差D1以及声源p到第三麦克风c的距离与声源p到第一麦克风a的距离之间的第二距离差D2。

在步骤405，根据第一距离差D1和第二距离差D2的符号以及D1、D2的绝对值大小比较关系，初步确定方位角Φ所落入的区间(即，方位角Φ落入区间1-6中的哪个区间中)。

参照图3和图4，假设D1为圆周上任一声源p到麦克风b的距离pb减去声源p到麦克风a的距离pa所得的距离差，即D1＝pb-pa；假设D2为圆周上声源p到麦克风c的距离pc减去声源p到麦克风a的距离pa所得的距离差，即D2＝pc-pa，则有如下表所示的关系：

  区间  Φ  D₁  D₂  性质  1  0°-60°  +  +  |D₁|＜|D₂|  2  60°-120°  -  +  |D₁|＜|D₂|  3  120°-180°  -  -  |D₁|＜|D₂|  4  180°-240°  -  -  |D₁|＞|D₂|  5  240°-300°  +  -  |D₁|＞|D₂|  6  300°-360°  +  +  |D₁|＞|D₂|

因此，根据D1、D2的符号以及D1、D2的绝对值大小比较关系，可以初步确定方位角所落入的区间，从而使每一个取值范围内具有单一的优化点，避免了局部最小值。同时减少估计的值域范围，提高估计速度。

在步骤406，用0.618法对双变量(即，距离R和方位角Φ)逼近优化点，获得最终的估计结果。

0.618法也称为黄金分割法，0.618为黄金分割点。可对二维单调唯一最小点的进行0.618法选优。

在本实施例中，先取距离R的初始值。距离R的初始值在50-400米的范围内选取，优选为200米。如果距离R小于50米，则声源借助视力可被容易确定。如果距离R大于400米，则接收的数据差比已经很小，包含的距离信息较少，不便于识别定位。因为将针对距离R的初始值进行0.618法逼近，因此，距离R的初始值可认为是对声源p定位的测量范围。

然后，固定距离R的初始值，在初步确定的方位角Φ的范围内，采用0.618法缩小优化点所在的区间，并逼近误差最小的方位角Φ，直到两次逼近的方位角Φ的误差变化小于预定阈值(例如相对误差为0.01)为止。接着，固定该逼近的方位角Φ，采用0.618法缩小距离R的范围，直到两次逼近的距离R的误差变化小于另一预定阈值为止，得到逼近的距离R。因此，对方位角Φ及距离R进行逼近，直到二者的误差变化达到指定的阈值以内，得到估计的结果。

具体地讲，对于方位角Φ，如上所述，根据D1、D2的符号以及D1、D2的绝对值大小比较关系，可以初步确定方位角Φ落入区间1-6中的哪个区间中。

接着，可按照下面的方式逼近误差最小的方位角Φ：

(a1)取方位角Φ为确定的区间上限S1乘以0.618；

(a2)根据三角形的边角关系，利用所述方位角Φ以及距离R的初始值求出声源p到各个麦克风的距离R1、R2、R3；

(a3)将计算出的距离R1、R2、R3分别除以音速，可以求得声源p到达各个麦克风的时间；

(a4)计算声源p到达第二麦克风b与声源p到达第一麦克风a的时间差t′_ba1，并计算声源p到达第三麦克风c与声源p到达第一麦克风a的时间差t′_ca1；

(a5)利用先前通过准L1互相关算法(互相关值)得到的声源p到达第二麦克风b与声源p到达第一麦克风a的时间差t_ba以及声源p到达第三麦克风c与声源p到达第一麦克风a的时间差t_ca，计算误差error1＝|t_ba-t′_ba1|+|t_ca-t′_ca1|；

(a6)取方位角Φ为上限S1乘以(1-0.618)，按照步骤(a2)-(a5)计算误差error2＝|t_ba-t′_ba2|+|t_ca-t′_ca2|；这里，t′_ba2和t′_ca2分别是方位角Φ为上限S1乘以(1-0.618)时计算的声源p到达第二麦克风b与声源p到达第一麦克风a的时间差以及计算的声源p到达第三麦克风c与声源p到达第一麦克风a的时间差；

(a7)如果|error1-error2|大于预定阈值(例如可以是0.01)，并且error1＞error2，则确定方位角Φ的范围在确定的区间的下限S2与上限S1乘以0.618之间；如果|error1-error2|大于所述预定阈值，并且error1＜error2，则确定方位角Φ的范围在上限S1乘以(1-0.618)与上限S1之间；

(a8)根据再次确定的方位角Φ所落入的区间，重复上述步骤(a1)-(a7)，直到|error1-error2|小于所述预定阈值。由于最后所取的两个方位角Φ之间相差很小，因此最后所取的两个方位角Φ之一或二者的平均值为最终估计的方位角Φ。

例如，假设初步确定方位角落入区间1中(0到60°)，则取方位角Φ为60°×0.618，为37.08°。因为在准L1相关算法中已经得到声源p到达各个麦克风的时间差t_ba、t_ca，因此可根据三角形的边角关系，将方位角Φ(37.08°)以及距离R的初始值代入方程(5)-(7)，可求出声源p到各个麦克风的距离R1、R2、R3。将计算出的距离R1、R2、R3分别除以音速，可以求得声源p到达各个麦克风的时间，进而可以求得声源p到达第二麦克风b与声源p到达第一麦克风a的时间差t′_ba1，并求得声源p到达第三麦克风c与声源p到达第一麦克风a的时间差t′_ca1。然后，计算误差error1＝|t_ba-t′_ba1|+|t_ca-t′_ca1|。

接着，取方位角Φ为60°×(1-0.618)，为22.92°。按照以上同样的方式求出误差error2＝|t_ba-t′_ba2|+|t_ca-t′_ca2|。

当|error1-error2|大于预定阈值(例如0.01)时，继续按照0.618法搜索。当error1＞error2时，所求方位角Φ的范围是在0～37.08°之间，取方位角Φ为37.08°×0.618，计算error1，之后取方位角Φ为37.08°×(1-0.618)，计算error2；反之，当error1＜error2时，所求方位角Φ的范围是在22.92°～60°之间，取方位角Φ为22.92°+(60°-22.92°)×0.618，计算error1，之后取方位角Φ为22.92°+(60°-22.92°)×(1-0.618)，计算error2。按照上述相同的方式继续求解。直到|error1-error2|小于所述预定阈值为止，最后所取的两个方位角Φ之一或二者的平均值为最终估计的方位角Φ。

在最终估计了方位角Φ之后，可对距离R进行逼近估计。

具体地讲，距离R的初始值以及零值确定了一区间(0，初始值)，可按照下面的方式逼近误差最小的距离R：

(b1)取距离R为确定的区间的上限S3乘以0.618；

(b2)根据三角形的边角关系，利用最终估计的方位角Φ以及距离R求出声源p到各个麦克风的距离R1、R2、R3；

(b3)将计算出的距离R1、R2、R3分别除以音速，可以求得声源p到达各个麦克风的时间；

(b4)计算声源p到达第二麦克风b与声源p到达第一麦克风a的时间差t′_ba3，并计算声源p到达第三麦克风c与声源p到达第一麦克风a的时间差t′_ca3；

(b5)利用先前通过准L1互相关算法(互相关值)得到的声源p到达第二麦克风b与声源p到达第一麦克风a的时间差t_ba以及声源p到达第三麦克风c与声源p到达第一麦克风a的时间差t_ca，计算误差error3＝|t_ba-t′_ba3|+|t_ca-t′_ca3|；

(b6)取距离R为上限S3乘以(1-0.618)，按照步骤(b2)-(b5)计算误差error4＝|t_ba-t′_ba4|+|t_ca-t′_ca4|；这里，t′_ba4和t′_ca4分别是距离R为上限S3乘以(1-0.618)时计算的声源p到达第二麦克风b与声源p到达第一麦克风a的时间差以及计算的声源p到达第三麦克风c与声源p到达第一麦克风a的时间差；

(b7)如果|error3-error4|大于另一预定阈值(例如可以是0.01)，并且error3＞error4，则确定距离R的范围在确定的区间的下限S4与上限S3乘以0.618之间；如果|error3-error4|大于所述另一预定阈值，并且error3＜error4，则确定距离R的范围在上限S3乘以(1-0.618)与上限S3之间；

(b8)根据再次确定的距离R所落入的区间，重复上述步骤(b1)-(b7)，直到|error3-error4|小于所述另一预定阈值。由于最后所取的两个距离R之间相差很小，因此最后所取的两个距离R之一或二者的平均值为最终估计的距离R。

例如，如果距离R的初始值为200米，则取距离R为200×0.618，为123.6米。因为在准L1相关算法中已经得到声源p到达各个麦克风的时间差t_ba、t_ca，因此可根据三角形的边角关系，将距离R(123.6米)以及前面已经求得的最终估计的方位角Φ代入方程(5)-(7)，可求出声源p到各个麦克风的距离R1、R2、R3。将计算出的距离R1、R2、R3分别除以音速，可以求得声源p到达各个麦克风的时间，进而可以求得声源p到达第二麦克风b与声源p到达第一麦克风a的时间差t′_ba3，并求得声源p到达第三麦克风c与声源p到达第一麦克风a的时间差t′_ca3。然后，计算误差error3＝|t_ba-t′_ba3|+|t_ca-t′_ca3|。

接着，取距离R为200×(1-0.618)，为76.4米。按照以上同样的方式求出误差error4＝|t_ba-t′_ba4|+|t_ca-t′_ca4|。

当|error3-error4|大于预定阈值(例如0.01)时，继续按照0.618法搜索。当error3＞error4时，所求距离R的范围是在0～123.6米之间，取距离R为123.6×0.618，计算error3，之后取距离R为123.6×(1-0.618)，计算error4；反之，当error3＜error4时，所求距离R的范围是在76.4～200米之间，取距离R为76.4+(200-76.4)×0.618，计算error3，之后取距离R为76.4+(200-76.4)×(1-0.618)，计算error4。按照上述相同的方式继续求解。直到|error3-error4|小于所述预定阈值为止，最后所取的两个距离R之一或二者的平均值为最终估计的距离R。

此外，由于在实际应用中，距离R的估计不如方位角Φ的估计重要，因此也可以在初步估计方位角Φ之后，采用固定的一些R值(例如，50米、100米、150米、200米、300米)进行估计，找出误差最小的值即可。

0.618法比二分之一法能够更快地得到逼近结果。

根据本发明，分析了延时特性，提出了方位分割的方法，采用0.618法加快逼近速度，从而减少计算的复杂度，实现了快速估计的算法。对延时的估计，提出了准L1相关的算法，减少了乘法运算，提高了估计的精度。另外，采用了插值算法进一步提高延时的估计精度，从而得到了快速准确的算法。

另外，在麦克风布局上，采用正三角形布局对各方向的敏感度一致性好。

虽然针对由3个麦克风构成的正三角形阵列来描述本发明的麦克风定位方法，但是也可采用由4个麦克风构成的正方形阵列来实现本发明的麦克风定位方法，只需要调整相应的几何关系即可。同样，采用正方形的麦克风布局对各方向的敏感度一致性好。

虽然本发明是参照其示例性的实施例被具体描述和显示的，但是本领域的普通技术人员应该理解，在不脱离由权利要求限定的本发明的精神和范围的情况下，可以对其进行形式和细节的各种改变。