一种基于声音信号到达时间差的声源定位解析解计算方法及装置
阅读说明：本技术 一种基于声音信号到达时间差的声源定位解析解计算方法及装置 (Sound source positioning analytic solution calculation method and device based on sound signal arrival time difference ) 是由 李帛翰 黄志尧 邱奕臻 冀海峰 王保良 于 2021-08-30 设计创作，主要内容包括：本发明公开了一种基于声音信号到达时间差的声源定位解析解计算方法及装置,装置包括麦克风阵列、前置放大器、数据采集模块以及计算机。本发明采用五个麦克风组成阵列,其构型满足声源位置解算线性方程组的系数矩阵可逆,同时获取五路声音信号,利用PHAT为加权函数的广义互相关算法得到观测信号之间的时间延迟,然后根据时间延迟求解得到声源的三维坐标以及声源发声到麦克风接收信号时间差的解析解。本发明通过特定选取麦克风阵列,实现三维空间中点声源坐标的解析求解,具有准确度高、体积小、功耗低和对环境影响小等特点,适用于卫星定位系统无法覆盖场景,如室内仓库等环境、禁止电磁波产生或周边存在电磁波强烈干扰源场景中的应用。(The invention discloses a sound source positioning analysis solution calculation method and device based on sound signal arrival time difference. The method adopts five microphones to form an array, the configuration of the array meets the condition that a coefficient matrix of a sound source position resolving linear equation set is reversible, five paths of sound signals are obtained simultaneously, time delay among observed signals is obtained by utilizing a generalized cross-correlation algorithm with PHAT as a weighting function, and then an analytic solution of three-dimensional coordinates of a sound source and a time difference from sound production of the sound source to signal receiving of the microphones is obtained according to the time delay solution. The method realizes the analytic solution of the midpoint sound source coordinate in the three-dimensional space by specifically selecting the microphone array, has the characteristics of high accuracy, small volume, low power consumption, small influence on the environment and the like, and is suitable for the application of a satellite positioning system which cannot cover the scene, such as an indoor warehouse and other environments, and the scene in which the generation of electromagnetic waves is forbidden or strong interference sources of the electromagnetic waves exist at the periphery.)
技术领域
本发明涉及声源定位技术，尤其涉及一种基于声音信号到达时间差的声源定位解析解计算方法及装置。
背景技术
随着科学技术的发展，越来越复杂的场景对目标定位技术提出更高的要求。早期的主动式目标定位技术通过主动发射探测信号，如电磁波、超声波等，接收经环境反射回来的探测信号后，分析计算得到目标在空间中的具体位置。与之对应的被动式目标定位技术，由于不需要测量装置主动发出探测信号，相比之下具有体积小、功耗低、隐蔽性强、适用范围广和对环境影响小等优点。因此，被动式定位技术如今被许多学者深入研究，并将成果成功应用在包括国防、安保、仓储、生态环境监测等领域。
声源定位技术是一种典型的被动式定位技术，主要原理是通过声音接收器接收声源发出的声音信号，然后通过一系列的信号处理手段分析信号的相关参数得到声源的位置信息。现有的声源定位技术主要分为基于信号到达时间差的声源定位(TDOA)，基于信号到达时间的声源定位(TOA)，基于信号入射角度的声源定位(DOA)，基于信号到达强度的声源定位(RSSI)等方法。其中，基于TDOA的声源定位方法是被广泛使用的声源定位方法之一，该方法流程中还细分为前后两个步骤：时延估计(时间差估计)和位置解算。其中，时延估计，即通过声音信号计算声源到不同麦克风间的时间延迟，广泛采用的算法为互相关算法；位置解算，即通过时间延迟解算得到声源的位置信息，所采用的主流方法有牛顿迭代法、球面插值法等。然而位置解算主流方法的计算结果如今仍只能得到声源位置的数值解，这一结果会出现迭代后发散以及受迭代初值影响较大的问题。同时，现有的常用方法一般需要数量较多的麦克风构成阵列，在求解数值解时存在计算量大等缺陷。因此，需要寻求一种能够尽可能地减少麦克风数量，花费更少计算代价的，实现求解声源位置解析解的方法。
发明内容
本发明针对基于信号到达时间差的声源定位技术的现状，提供了一种基于声音信号到达时间差的声源定位解析解计算方法及装置，该定位装置中的麦克风阵列构型在三维空间中的定位表现为在任意位置均能够获得解析解。本发明采用五个麦克风组成麦克风阵列采集声音信号，选择广义互相关算法中适合的加权函数进行时延估计，通过时间延迟求解得到声源位置的解析解。该麦克风构型和方法在声源定位方面具有计算量小、准确性高等优势。
本发明的技术方案如下：
本发明一方面提供了一种声源定位麦克风阵列的构型，其特征在于包含五个麦克风，坐标为Mi(xi，yi，zi)，i＝1，2，......，5；令麦克风M1为基准麦克风，其可以为五个麦克风中的任意一个；根据其余四个麦克风声音到达与基准麦克风声音到达的四个时间差，可知三维空间声源坐标有解析解的条件为矩阵A可逆，所述声源定位麦克风阵列的构型中的五个麦克风坐标被配置为满足该条件，其中矩阵A的表达式为：
其中c为声速，所述基准麦克风记为M1(x1，y1，z1)，不失一般性地，其可以为五个麦克风中的任意一个，ti1为声源发出的声音信号被第i个麦克风Mi接收的时间(ti0，即声源发出的声音信号传播到第i个麦克风Mi所需的时间)与基准麦克风M1接收到的声音信号时间(t10，即声源发出的声音信号传播到基准麦克风M1所需的时间)之差，ti1＝ti0-t10。
判断矩阵A是否可逆的一种方法为判断A的行列式是否不等于0：
为保证矩阵A可逆，减少求解矩阵A的行列式以及逆矩阵的过程中的计算量，且获得有效的解析解，麦克风阵列优化构型为：其中的四个麦克风在同一平面内且分别位于一个凹四边形的四个顶点处，第五个麦克风在该平面外任意位置。该构型对应的矩阵A可逆的证明如下：
如图2所示建立坐标系以及设定麦克风坐标。不失一般性地，设麦克风M1、M2、M3、M4在X-Y平面内，即z2＝z3＝z4＝0，点M1坐标为(0，0，0)，麦克风M4落于麦克风M1、M2、M3作为顶点的三角形内部，且不落在三角形的边上，麦克风M1、M2、M3、M4四点构成一凹四边形；麦克风M5在该平面外，即z5≠0；设第i个麦克风的位置向量为：mi＝[xi，yi，zi]T，i＝1，2，......，5；声源位置向量为：p＝[x，y，z]T；声源到麦克风Mi的距离为： i＝1，2，......，5，其中||·||表示向量的2范数(向量的模)。
点M1、M2、M3、M4同处X-Y平面，点M1、M2、M3所构成三角形为典型的凸集，点M4位于该三角形内部，且不落在三角形的边上，根据线性代数和凸集理论，可得：
m4＝αm2+βm3 (1)
x4＝αx2+βx3 (2)
y4＝αy2+βy3 (3)
其中α＞0，β＞0，且α+β＜1。相应地，矩阵A可表示为：
矩阵A的行列式为：
在detA的表达式中，声速c不为0；z5不为0；向量m2和m3不共线，即 因此，在此优化构型下，矩阵A可逆的证明问题就转化为证明t41-at21-βt31≠0。
ct21＝c(t20-t10)＝D2-D1 (6)
ct31＝c(t30-t10)＝D3-D1 (7)
ct41＝c(t30-t10)＝D4-D1 (8)
根据α＞0，β＞0，且α+β＜1，公式(6)(7)(8)带入c(t41-αt21-βt31)：
c(t41-αt21-βt31)＝D4-D1-α(D2-D1)-β(D3-D1)＝D4-αD2-βD3-(1-α-β)D1
＝||m4-p||-α||m2-p||-β||m3-p||-(1-α-β)||m1-p||
＝||m4-p||-||α(m2-p)||-||β(m3-p)||-||-(1-α-β)p|| (9)
根据向量不等式公式可得：
上式中||m4-p||-||α(m2-p)||-||β(m3-p)||-||-(1-α-β)p||＝0的充要条件为向量m2-p、m3-p、m4-p三个向量平行且同向。而在此优化构型下，麦克风M4落于麦克风M1、M2、M3作为顶点的三角形内部，且不落在三角形的边上，因此向量m2-p、m3-p、m4-p必不平行且不同向，取等条件不成立，因此有：
c(t41-αt21-βt31)＝||m4-p||-||α(m2-p)||-||β(m3-p)||-||-(1-α-β)p||＜0 (11)
由于声速c不为0，则t41-αt21-βt31＜0，即t41-αt21-βt31不为0。
综上所述，在所描述的构型下，detA不为0，即矩阵A可逆。
在此基础上进一步地减少求解矩阵A的行列式以及逆矩阵的过程中的计算量，令基准麦克风M1为原点(0，0，0)，麦克风M2坐标为(a，0，0)，麦克风M3坐标为(0，a，0)，麦克风M4坐标为(b，b，0)，麦克风M5坐标为(0，0，a)，为保证点M4位于点M1、M2、M3构成的三角形内，其中a＞2b＞0。该构型下的矩阵A可表示为：
对应的逆矩阵A-1为：
其中detA＝16a2c2(at41-bt31-bt21)，相应的声源坐标P(x，y，z)和t10的解析表达式为
本发明还提供了一种基于声音信号到达时间差的声源定位装置，其包括上述构型的麦克风阵列、前置放大器、数据采集模块以及计算机。每个麦克风接收到的声音信号经前置放大器放大，再由数据采集模块进行声音信号采集以及模数转换，数据采集模块通过数据线与计算机进行通信，由计算机解算麦克风之间的时间延迟后计算得到声源位置的解析解。
本发明还提供的一种基于声音信号到达时间差的声源定位解析解计算方法，其包括以下步骤：
步骤一：采用上述构型的麦克风阵列中的五个麦克风分别接收声源发出的信号，各麦克风的接收信号经前置放大器进行放大，再经数据采集模块进行模数转换和采样，得到声音信号si，i为麦克风编号，i＝1，2，......，5；
步骤二：对声音信号si进行分帧加窗、去均值、归一化以及离散快速傅里叶变换处理，得到分u帧处理后的频域声音信号S′ij，其中j＝1，2，......，u；
步骤三：根据麦克风阵列的构型，选择麦克风M1接收到的声音信号为基准，利用以PHAT为加权函数的广义互相关算法，计算其余四组声音信号相对于麦克风M1接收到的声音信号的时间延迟ti1，i＝2，3，......，5；；
步骤四：使用权利要求2所述的方法，根据时延估计得到的四组声音信号相对于麦克风M1接收到的声音信号的时间延迟ti1，并得到它们与声源位置坐标P(x，y，z)之间的方程；
步骤五：使用权利要求2所述的方法，化简步骤四中的方程，得到关于声源坐标P(x，y，z)和声源发出的声音信号传播到基准麦克风M1所需的时间t10的有解析解的线性方程组，求解得到声源坐标P(x，y，z)和的t10解析表达式。
优选的，步骤二中，分帧操作选取五组声音信号上相同的位置作为起点，以w＝10～30ms为帧长，v为帧移，向后选取u帧，其中v为任意大于0的时间，通常取0＜v＜w。
优选的，步骤二中，在归一化和离散快速傅里叶变换之间，还包括巴特沃斯滤波的步骤。
优选的，所述的步骤三，具体为：
以麦克风M1接收到的声音信号为基准，利用以PHAT为加权函数的广义互相关算法，计算其余四组声音信号与基准信号的相对时间延迟ti1：首先，由S′ij各自与S′1j进行共轭后相乘得到互功率谱函数然后，选取PHAT加权函数对互功率谱加权后进行离散快速傅里叶逆变换得到互相关函数 接着，其余四组声音信号的相对采样点延迟将相对采样点延迟转换为相对时间延迟其中fs为模数转换部分设定的采样频率；最后，舍去ti1j中的粗大误差，取平均值计算得到四组声音信号相对于麦克风M1接收到的声音信号的时间延迟ti1。
优选的，根据五个麦克风坐标Mi(xi，yi，zi)以及时延估计得到的四组声音信号相对于麦克风M1接收到的声音信号的时间延迟ti1，可以得到它们与声源位置坐标P(x，y，z)之间的关系：
其中，c为声速，T为环境的摄氏温度。
所述步骤五为：化简步骤四中的方程，得到声源坐标P(x，y，z)和t10的解析表达式。
首先，将公式(15)～(19)两边取平方并展开，得到：
x2-2xx1+x1 2+y2-2yy1+y1 2+z2-2zz1+z1 2＝c2t10 2 (21)
x2-2xx2+x2 2+y2-2yy2+y2 2+z2-2zz2+z2 2＝c2t21 2+2c2t21t10+c2t10 2 (22)
x2-2xx3+x3 2+y2-2yy3+y3 2+z2-2zz3+z3 2＝c2t31 2+2c2t31t10+c2t10 2 (23)
x2-2xx4+x4 2+y2-2yy4+y4 2+z2-2zz4+z4 2＝c2t41 2+2c2t41t10+c2t10 2 (24)
x2-2xx5+x5 2+y2-2yy5+y5 2+z2-2zz5+z5 2＝c2t51 2+2c2t51t10+c2t10 2 (25)
然后，将公式(22)～(25)两侧分别减去公式(21)并移项化简，得到：
转换为线性方程组Aξ＝σ形式：
其中ξ＝[x，y，z，t10]T，
当满足方程中A可逆，此时可以得到声源坐标P(x，y，z)和t10的解析表达式为：
本发明与现有技术相比具有的有益效果：
1)本发明采用了一种新的麦克风阵列构型。该阵列中的五个麦克风在数量上已经达到计算三维空间定位问题的最小值，构型满足解析解存在的条件，并且在定位表现上在任意位置均能够准确定位。
2)本发明位置解算通过求解线性方程组得到解析解。可以克服传统的数值解需要数量较多的麦克风构成阵列，在求解时存在的计算量大等缺陷。解析解求解方法具有稳定性强、准确性高，计算量小的特点，有利于实时定位。
3)本发明对信号进行分帧加窗操作后，对每一帧的时延估计采用去掉粗大误差并取平均的方式，相较于直接对原信号进行一次时延估计的方式，能够提高算法的稳定性。同时由于声音信号经过时域上固定频率fs采样，单次时延估计的分辨率仅有1/fs，而多帧平均值可以得到分辨率更高的时延估计。
4)本发明在广义互相关算法中，采用PHAT加权函数，在效果上相当于白化滤波，在时域上使得广义互相关函数的峰值更加尖锐。
附图说明
图1是一种基于声音信号到达时间差的声源定位装置的结构示意图；
图2是应用在本发明中的一个麦克风阵列的优化空间构型。
图3是图2所示麦克风阵列的优化空间构型的一个进一步优化构型。
具体实施方式
本发明在基于声音信号到达时间差的声源定位问题的基础上，提出了一种新的麦克风阵列的构型，以及针对该构型的声源定位方法。该方法较传统方法具有稳定性强，准确性高，计算量小的特点。
如图1所示，一种基于声音信号到达时间差的声源定位的装置，包括麦克风阵列、前置放大器、数据采集模块以及计算机。五个麦克风用于接收声音信号，并以能够保证声源位置解析解存在的构型排列，每个麦克风接收到的声音信号经前置放大器放大，通过排线与数据采集模块相连实现麦克风供电以及声音信号读入。数据采集模块与电源模块相连实现声音信号采集以及信号模数转换处理部分的装置供电，并且通过数据线与计算机进行通信。
本发明的麦克风阵列包含五个麦克风，五个麦克风坐标为Mi(xi，yi，zi)，i＝1，2，......，5；令麦克风M1为基准麦克风，其可以为五个麦克风中的任意一个；根据其余四个麦克风声音到达与基准麦克风声音到达的四个时间差，可知三维空间声源坐标有解析解的条件为矩阵A可逆，所述声源定位麦克风阵列的构型中的五个麦克风坐标被配置为满足该条件，其中矩阵A的表达式为：
其中c为声速，所述基准麦克风记为M1(x1，y1，z1)，ti1为声源发出的声音信号被第i个麦克风Mi接收的时间(ti0，即声源发出的声音信号传播到第i个麦克风Mi所需的时间)与基准麦克风M1接收到的声音信号时间(t10，即声源发出的声音信号传播到基准麦克风M1所需的时间)之差，ti1＝ti0-t10。
一种判断矩阵A是否可逆的方法，可以判断A的行列式是否不等于0：
如图2所示麦克风阵列构型的一种优化构型，为减少求解矩阵A的行列式以及逆矩阵的过程中的计算量，且获得有效的解析解，可以令四个麦克风(M1、M2、M3、M4)在同一平面内且分别位于一个凹四边形的四个顶点处，第五个麦克风(M5)在该平面外任意位置。
在此基础上进一步地减少求解矩阵A的行列式以及逆矩阵的过程中的计算量，如图3所示，令麦克风M1为原点(0，0，0)，麦克风M2坐标为(a，0，0)，麦克风M3坐标为(0，a，0)，麦克风M4坐标为(b，b，0)，麦克风M5坐标为(0，0，a)，其中a＞2b＞0。该构型下的矩阵A可逆，且形式为：
对应的逆矩阵为：
其中detA＝16a2c2(at41-bt31-bt21)，此时相应的声源坐标P(x，y，z)和t10的解析表达式为：
本发明方法的实施过程中，首先通过上述麦克风阵列接收声源发出的声音信号通过前置放大器和数据采集模块发送给计算机，得到五组声音信号序列si(i＝1，2，......，5)，然后对序列做分帧加窗、去均值、归一化、巴特沃斯滤波以及离散快速傅里叶变换处理，得到处理后的频域声音信号(i＝1，2，......，5)，接着可以以任一麦克风Mk接收到的声音信号为基准利用PHAT加权的广义互相关算法估计其余四组声音信号的相对时间延迟，在本实施例中，仅以麦克风M1接收到的声音信号为基准进行说明，其余四组声音信号相对于麦克风M1接收到的声音信号的时间延迟ti1(i＝2，3，......，5)，最后根据时延估计结果解算声源坐标的解析解P(x，y，z)。具体操作如下：
1)对计算机接收到的声音信号si(i＝1，2，......，5)进行分帧加窗：
s′ij(n)＝hamming(n)si(N+n+(j-1)×r)(n＝1，2......q)(i＝1，2，......，5)(j＝1，2，......，u)
其中，j代表u帧中的第j帧，N代表分帧起始点，q代表帧长，r代表帧移，hamming(n)代表长度为wlen的汉明窗。分帧操作选取五组声音信号上相同的位置作为起点，时域上以10～30ms为帧长，任意时间长度v为帧移，通常两帧之间存在重叠，向后选取u帧。
2)对分帧加窗后的信号进行去均值、归一化、巴特沃斯滤波以及离散快速傅里叶变换处理：
其中，bw代表巴特沃斯滤波器的时域冲激响应函数。
3)以麦克风M1为基准，计算互功率谱：
4)计算PHAT加权的广义互相关函数：
5)计算时间延迟：
6)计算声速c：
其中，T为已知的环境摄氏温度。
7)根据五个麦克风坐标Mi(xi，yi，zi)以及时延估计得到的四组声音信号相对于麦克风M1接收到的声音信号的时间延迟ti1，可以得到它们与声源位置坐标P(x，y，z)之间的关系：
首先，将公式(15)～(19)两边平方并展开，得到：
x2-2xx1+x1 2+y2-2yy1+y1 2+z2-2zz1+z1 2＝c2t10 2 (21)
x2-2xx2+x2 2+y2-2yy2+y2 2+z2-2zz2+z2 2＝c2t21 2+2c2t21t10+c2t10 2 (22)
x2-2xx3+x3 2+y2-2yy3+y3 2+z2-2zz3+z3 2＝c2t31 2+2c2t31t10+c2t10 2 (23)
x2-2xx4+x4 2+y2-2yy4+y4 2+z2-2zz4+z4 2＝c2t41 2+2c2t41t10+c2t10 2 (24)
x2-2xx5+x5 2+y2-2yy5+y5 2+z2-2zz5+z5 2＝c2t51 2+2c2t51t10+c2t10 2 (25)
然后，将公式(22)～(25)两侧分别减去公式(21)并移项化简，得到：
转换为线性方程组Aξ＝σ形式：
其中ξ＝[x，y，z，t10]T，
当满足方程中A可逆，此时可以得到声源坐标P(x，y，z)和t10的解析表达式为：
本发明的阵列中的五个麦克风在数量上已经达到计算三维空间定位问题的最小值，构型满足解析解存在的条件，并且在定位表现上在任意位置均能够准确定位。通过本发明的麦克风构型和基于声音信号到达时间差的声源位置的解析解求解方法，可以实现三维空间中全方位的声源定位，位置解算得到的解析解又有着准确性高，稳定性强的特点。
- 上一篇：一种医用注射器针头装配设备
- 下一篇：一种基于双麦克风的声音定位方法及其装置