利用三分法进行多圆联动搜索的声源定位算法

覃泓铭 王 玫，2 罗丽燕* 韦金泉

1(桂林电子科技大学认知无线电与信息处理教育部重点实验室 广西 桂林 541004) 2(桂林理工大学信息科学与工程学院 广西 桂林 541007)

0 引 言

声源定位技术由于广泛应用于智能交通[1]、安防监控[2]、电话会议[3]、人机交互[4]、导航[5]等而备受关注。传统的声源定位方法可分为三类：基于可控波束形成的方法、基于高分辨率谱估计的方法和基于到达时间差的方法。其中,基于可控波束形成的定位方法[6-7]需要做全局搜索,计算量大,并且依赖于声源信号和噪声信号的先验知识。基于高分辨率谱估计的定位方法[8-9]需要在空间中搜索目标声源位置,这种方法的定位精度高但计算复杂度非常高。相比前两种定位技术,基于到达时间差(TDOA)的定位方法[10-11]计算复杂度较低、定位速度较快性好且硬件成本低,因此成为声源定位领域中的研究热点。

基于TDOA的定位方法的中心思想是求解由几何关系构成的方程组,求解方程组的复杂度会随着接收节点/基站数量的增加而加大,并且设备误差与时延估计误差可能会导致方程组无解,因此基于该方法的声源定位问题往往被表述为一个非线性优化问题。已经存在很多方法用于解决优化问题,例如最大似然法、非线性最小二乘法和加权最小二乘法,但这些方法具有计算复杂度高、需要估计初始位置的缺点。Chan算法、Fang算法、Taylor级数展开法与最小二乘法的结合[12-14]解决了确定初始位置与声源位置最优估计的问题,但仍存在计算复杂度高、对初始位置敏感的问题。为了解决这一问题,相关研究者从几何分析的角度出发,提出了利用圆的几何关系搜索声源位置的定位算法[15-17]。文献[15]提出了基于收缩圆的定位方法，其思想是以恒定的步长收缩所有圆的半径,直到相交区域的面积达到较小的阈值。文献[16]中提出一种基于GCC-PHAT与收缩圆的优化方法,相比于收缩圆方法,定位精度与速度得到了改善。文献[17]中提出了多圆联动搜索目标的定位算法,其思想是：首先计算初始搜索半径,根据TDOA的测量值计算每个接收节点相对于搜索半径的联动半径,并根据搜索步长改变搜索半径进行搜索,找到具有最大步长的最佳半径,然后减小步长进行搜索以获得更准确的估计位置。该方法将传统的多维搜索转换为一维搜索,并且初始值的计算简单,降低了计算复杂度。相比于传统的搜索算法,该方法有效提高了搜索速度和定位精度,但该方法通过逐步提高步长的精度来提升定位精度,在搜索过程中存在重复搜索的问题。本文提出一种利用三分法进行多圆联动搜索的声源定位算法，优化了文献[17]中减小搜索步长进行重复搜索的问题。此外,本文对于初始值的设定进行了改进,减少了搜索迭代次数。较于文献[17],本文算法不仅减少了迭代次数,还降低了每次迭代搜索时的计算复杂度。实验结果表明,本文方法在提高定位速度的同时,有效提升了定位精度。

1 基于TDOA的声源定位方法

基于TDOA的声源定位基本思想[18]是：通过测量声源到各个接收节点的到达时间差，得到至少两个距离差，通过距离差可以获得双曲线方程，求解双曲线的交点即可得到声源的位置。考虑如图1所示的二维声源定位系统，假设声源位置为S(x,y)，接收节点位置为Ni(xi,yi)，声源到各个接收节点的距离为ri，到各接收节点之间的距离差为dij，以节点0为参考节点，声源到各传感节点与参考节点的时间差为τi0，c为声速，则有：

(1)

di0=r0-ri=cτi0

(2)

(i=0,1,2,…,n-1)

(3)

图1 声源定位系统示意图

声源位置就是由式(4)构成的方程组的解,该方程组是基于双曲线的非线性方程组,在理想条件下,通过联立方程组即可求出解,即双曲线的交点。但由于设备误差、时延估计误差的存在,这个交集往往是空集,因此求解该方程组的重点是如何采用搜索算法估计声源位置。声源定位中,两种广泛使用的声源定位算法是Chan算法和Taylor算法。Chan算法将如上所述的非线性方程组转化为线性方程组,然后采用最小二乘法(LS)获得最佳解。Taylor算法将非线性方程组通过Taylor级数展开的方法整理成线性方程组,选择一个初始值后进行迭代搜索,直到坐标精度满足搜索的终止条件，即可获得最佳解。

1.1 Chan算法

Chan算法[19]是非递归双曲线方程组解法,具有解析表达式解。其主要的特点为在测量误差服从理想高斯分布时定位精度高,并且可以通过增加基站数量来提高算法精度。该算法的核心思想如下：

将式(1)展开并将其线性化得：

(x0-xi)x+(y0-yi)y=ki+r0di0

(4)

(5)

由此可得到矩阵表达式：

AX=F

(6)

(7)

由此可求出声源的估计位置：

X=(ATA)-1ATF

(8)

1.2 Taylor算法

Taylor级数展开法[20]是一种需要初始估计位置的递归算法,在每一次递归中通过求解TDOA测量误差的局部线性最小二乘(LS)解来改进估计位置。该算法的核心思想如下：

将式(3)根据Taylor级数展开得：

ht=Gtδ+ε

(9)

(10)

式中：Rm，1表示由声源到传感节点m与传感节点1的时延差计算出来的距离差；Rm-R1为声源到传感节点m与到传感节点1的真实距离差。

(11)

其加权最小二乘解为：

(12)

式中：Q为TDOA测量值的协方差矩阵。初始迭代时,令x=x0，y=y0,在下一次迭代中，则有：

x(1)=x0+Δxy(1)=y0+Δy

(13)

重复以上过程,直到Δx、Δy足够小,并满足一定的阈值ε,使得：

|Δx|+|Δy|<ε

(14)

经过k次迭代,得到的解(xk,yk)即为声源的估计位置。

2 算法简介

2.1 算法的基本原理

算法基于TDOA测量值,通过构建几何关系将非线性问题转化为线性问题,将多维搜索转化为一维搜索。算法的基本原理可由下描述：

假设声源位置为S(x,y),参考节点位置为N0(x0,y0),其他参考节点位置为Ni(xi,yi),声源到各个接收节点的距离为ri,声源到各接收节点与到参考节点的时间差为τi0,则有：

(15)

(16)

易知,式(3)可由式(15)、式(16)联立得到。对于式(15)、式(16)的求解,考虑如图2所示的定位系统,假设估计的声源到参考节点的距离为Rk,到其他接收节点的距离为Ri,则有：

(17)

(18)

联立式(17)、式(18)可得：

(19)

令：

(20)

则式(19)的解为：

(21)

(22)

图2 接收节点部署示意图

由方程可知,若改变参考圆半径Rk,其他圆的半径Ri也会随着改变,即不同的Rk对应不同的解(估计声源位置)。如图3所示(★：声源位置；●：接收节点；-‥-：Rk为初始值时的几何关系；---：Rk改变一次时的几何关系；—：Rk改变N次时的几何关系)，不断地搜索Rk,直到Rk=r0,此时的估计声源位置即为真实声源位置。

图3 多圆联动搜索的基本思想

在搜索Rk的过程中,初始值的设定以及搜索算法直接影响到算法收敛速度,即声源搜索速度。文献[17]根据参考圆B0与各个参考圆Bi相交时,方程具有解的条件设定了Rk的初始搜索半径,并以初始搜索半径为最小值,通过改变搜索步长来寻找最优半径,最后根据最优半径估计出声源位置。由于只设定了初始搜索半径的最小值,该搜索过程为由小往大的单向搜索。为了进一步加快搜索速度,本文根据接收节点位置的几何关系设定了Rk的初始搜索半径的最小值与最大值,即搜索半径区间,并同时改变最小值与最大值进行双向搜索,不断地缩小Rk的搜索区间。此外,文献[17]首先采用较大的搜索步长来估计声源的粗略位置,然后通过减小步长进行精确搜索。当搜索步长较大时,算法收敛速度较快,但定位精度相应较低,当搜索步长较小时,定位精度得到有效的提高,但重复搜索加长了搜索时间,整个搜索过程速度较慢。本文将三分法运用到双向搜索过程中,改善了文献[17]中重复搜索的问题,进一步加快了搜索速度。同时,利用三分法进行搜索,最终搜索区间的下限值无限逼近于上限值,有效提高了定位精度。

2.2 搜索半径的初始值

设Lij为接收节点i到参考节点j之间的距离,则：

(23)

由声源、节点Ni、节点Nj构成的三角形满足三角形三边关系：

Rk+Rk-cτij≥Lij

(24)

即：

(25)

声源应该位于由4个节点形成的矩形区域中，则Rk应满足：

Rk≤max{Lij}

(26)

所以Rk应满足：

min{(L1j-d1j)/2}≤Rk≤max{Lij}

(27)

以式(27)确定的搜索区间进行搜索,在搜索过程中,每改变一次搜索区间,误差函数的值都会随着改变。当误差函数为最小值时即停止搜索,此时的Rk为全局最优解,即最佳半径。

2.3 搜索的停止条件

根据式(7)，可定义误差函数：

(28)

式中：第一项为真实声源到各接收节点的距离差，第二项为根据Rk估计的声源到各接收节点的距离差。

(29)

在理想情况下，εi(x,y)=0，但由于设备误差、时延估计误差的存在，εi(x,y)≠0，所以当误差满足式(29)时，即停止搜索。

图4 三分法的基本思想

文献[17]中,改变步长进行迭代搜索的算法计算复杂度为O(mn),其中m为改变步长精度的次数。由三分法的基本思想可知,利用三分法进行搜索算法的计算复杂度为O(2log3n),易知O(2log3n)

利用三分法搜索Rk从而搜索误差的最小值,即可获得估计声源的最佳位置,该过程由以下步骤完成：

步骤1计算搜索半径Rk的初始值：

(30)

(31)

步骤2进行第i次迭代搜索(可见算法1)：

(32)

则

(33)

反之

(34)

算法1基于三分法的半径搜索算法

3) fori=1;i<100;i++

9) else

12) end

14) break

15) end

16) end

3 实验及分析

为了验证本文方法的定位性能,采用文献[17]提出的多圆联动搜索目标定位算法与本文算法进行比较。实验选用Behringer ECM 8000全向型麦克风4个、哈曼卡顿aura studio音响1个、M-AUDIO M-TRACK QUAD 4通道外置声卡1个,另外使用1台联想笔记本控制收发信号及数据处理。如图5所示,实验场景较为空旷,其中一面为墙面。受麦克风信号线长限制,麦克风以矩形(6.5 m×6.5 m)阵型进行摆放,坐标分别为N0(0,0)、N1(6.5,0)、N2(6.5,6.5)、N3(0,6.5)。麦克风与音响高度均为1.3 m,声源信号为0至16 k chirp信号,采样频率为32 kHz。选取3个声源位置S1(2,2)、S2(3,4)、S3(5.5,3),在每个声源位置分别采集20组独立数据。每组定位结果如图6-图8所示,平均误差如表1所示。

图5 实验场景

图6 声源位于(2,2)处的误差对比结果

图7 声源位于(3,4)处的误差对比结果

图8 声源位于(5.5,3)处的误差对比结果

表1 两种方法的声源定位平均误差

从表1可以看出,从位置1到位置3的平均误差在增大,原因在于声源越靠近墙面(见图5),麦克风接收到的反射波干扰越大,时延估计的误差则越大,导致定位误差增大,本文暂不讨论时延估计误差对定位的影响。在定位精度方面,本文方法的定位精度与文献[17]方法的定位精度相差不大,但略微优于文献[17]方法。为证明本文方法在计算复杂度方面的优越性,取声源位置为(3,4)时的数据进行计算复杂度的对比,结果如表2、表3所示。

表2 文献[17]方法的Rk搜索过程 单位：m

表3 本文方法的Rk搜索过程 单位：m

从表2可以看出,文献[17]的方法首先以3.08 m为初始半径、1 m为步长进行搜索,找到最佳半径位于[5.80,6.08]区间。然后将步长减小为0.1 m进行搜索,找到最佳半径位于[5.18,5.28]区间。最后将步长减小为0.01 m进行搜索,直到误差满足搜索结束条件。搜索结束时,Rk为5.21 m,对应的估计声源坐标为(3.135,4.170),经过10次迭代搜索找到声源最佳估计位置。而从表3可以看出,本文方法首先确定最佳半径位于[3.080,9.192]区间,然后通过三分法不断改变区间的左右界限,并结合误差判断是否结束搜索。搜索结束时,Rk为5.17 m,对应的估计声源坐标为(3.105,4.146),仅经过5次迭代搜索便可找到声源最佳估计位置。对于文献[17]的方法,通过改变三次步长进行迭代搜索,易知其计算复杂度为O(3n),而对于本文方法,通过将n/3处的元素与2n/3处的元素比较,将搜索范围缩小为原来的n/3,其计算复杂度为O(2log3n)。

综上所述,相比于文献[17]提出的多圆联动搜索声源定位算法,本文算法不仅减小了迭代搜索的次数,并且减小了每次搜索时的计算复杂度。在加快定位速度的同时,有效提升了定位精度。

4 结 语

对于TDOA声源定位系统中的定位估计阶段,定位方程求解速度与定位精度是关键问题。较传统的迭代搜索定位算法,多圆联动搜索声源定位算法有效提高了搜索速度和定位精度。本文提出利用三分法进行多圆联动搜索的声源定位算法,将非线性方程组转化为线性方程组,并在迭代搜索的过程中引入三分法进行搜索,进一步提高了定位求解速度与定位精度。理论分析和实验结果均证明了本文算法的优越性。未来的工作将研究减小时延估计误差的方法并与本文算法相结合,进一步提高该算法的定位精度。