CN110488223A - 一种声源定位方法 - Google Patents

Abstract

一种声源定位方法属于声源定位技术领域。包括以下步骤：建立均匀圆阵麦克风阵列，并以此构建空间坐标系确定各麦克风的坐标，获取麦克风采集的声源信号；根据各麦克风提供的声源信号，利用广义互相关估计出声源信号到达各麦克风之间的时间差；由参考阵元与所得的时间差共同构建各麦克风接收信号所对应的矢量形式，最后通过计算谱函数确定DOA估计值。本发明的方法及系统，克服了声源信号间的相互干扰，提高了定位精度，而且本发明有较高的抗混响性，能在强混响情况下准确估计出声源的位置。

Description

一种声源定位方法

技术领域

本发明属于声源定位技术领域，特别涉及到阵列信号处理，具体为一种基于麦克风阵列进行的声源定位方法。

背景技术

近年来，基于麦克风阵列的声源定位技术在工业检测、智能会议室、军事勘察等领域都有广泛的应用。声源定位的基本原理是将一组麦克风按照一定的空间几何结构排列，接收到声源信号后，用特定的算法确定声源在麦克风坐标系中的方位角和俯仰角，即声源信号的DOA。基于麦克风阵列的声源定位方法主要有：基于最大输出功率的可控波束形成的定位方法、基于时延(TDOA，Time Difference ofArrival)估计的声源定位方法和基于高分辨率的空间谱估计方法。基于最大输出功率的可控波束形成的定位方法对麦克风阵列接收到的语音信号进行滤波、加权求和，然后直接控制麦克风指向使波束有最大输出功率的方向，因其计算复杂度高不可能被用于实时处理系统；基于TDOA的声源定位方法是通过估计声源到麦克风之间的时延差并根据麦克风的位置来估计声源位置，只适用于单声源，对于多声源其定位效果不理想；且房间的混响会影响定位精度。基于高分辨率的空间谱估计方法的精度和分辨率得到了很大的提高，但在实际应用中定位效果不佳。目前用于实际的声源定位系统存在定位精度低等问题，难以快速的确定声源位置，因此，在声源定位背景下的两个关键问题为：麦克风阵列模型的构建和声源定位算法的改进，这一直是研究的重点。

总体来说，基于麦克风阵列的声源定位系统设计已经取得了很大的成就，但各种理论算法尚不够完善,在应用配置上仍有难点需要去进一步深入研究。阵型的选择和声源的定位算法的研究是麦克风阵列进行声源定位的关键技术。本专利进一步研究声源定位算法，设计硬件平台使其能够满足声源定位的要求，完善系统使其适合更多的应用场景。

发明内容

本发明的目的是提供一种科学合理，适用性强，定位可视化的声源定位系统，使得能够实时、准确地计算出被测物产生的声场，确定出声源的DOA值。

本发明的目的是由以下技术方案来实现的：麦克风阵列模型的建立与声源定位算法的改进，它包括的内容有：

1)麦克风阵列模型的建立

本发明采用均匀圆阵，将M个完全相同的麦克风均匀分布在平面X-Y上一个半径为r的圆周上，为了计算方便在圆心处加上一个阵元，并以此阵元作为参考阵元。以整个阵列模型的中心位置即参考阵元为空间坐标系的圆心，坐标系的X轴为参考阵元与第一个阵元之间的连线，建立如图1所示的声源定位系统的空间坐标系。均匀圆阵模型如图1所示。

第m个阵元与X轴之间的夹角为：

在空间坐标系中第m个阵元的坐标为：

远场中入射声源的单位向量可表示为：

因此第m个阵元与参考阵元之间的延时可表示为：

(1)-(4)式中M和r的意义如前文表述，式(4)中c为声速。阵列的响应矩阵为：

A(ω,τ)＝[e^jωτ1,e^jωτ2,...,e^jωτM]^T (5)

式(5)中w＝2πf为载频，λ为信号的波长。将式(4)代入式(5)，阵列的响应矩阵可化简为：

2)声源定位算法的改进

经过多通道数据采集系统获得的声源信号是几组数字信号，直接对所得的声源信号求协方差矩阵，再计算空间谱函数并进行峰值搜索时会出现对称双峰现象，定位系统分辨率较低。针对以上缺陷，本发明提出利用所得的信号求阵列信号的方向矩阵，并构建各阵元接收信号所对应的矢量形式，通过谱函数计算结果寻求峰值，再确定DOA估计值的方法。

基于麦克风阵列的改进MUSIC算法分为两个步骤：第一步为时延估计，用于计算各个阵元与参考阵元之间的时间差，根据时间差构建方向矩阵。第二步为声源信号的DOA估计，首先利用阵列的方向矩阵与参考阵元输出信号求得各阵元对应的矢量信号，并以此构建阵列的协方差矩阵，再划分信号子空间与噪声子空间得出谱函数，最后求得DOA估计值。改进MUSIC算法的系统框图如图2所示。

改进的MUSIC算法主要采用的广义互相关算法求延时，主要是通过求两信号之间的互功率谱，并在频域内给予一定的加权来抑制噪声和反射的影响，再反变换到时域，从而得到两个信号之间的互相关函数，该互相关函数的峰值位置为两个信号之间的相对时延。

假设麦克风接收的信号的表达式如式(7)和式(8)：

假设麦克风接收的信号的表达式如式(7)和式(8)：

x₁(n)＝a₁S(n-τ₁)+w₁(n) (7)

x₂(n)＝a₂S(n-τ₂)+w₂(n) (8)

式(1)中：a₁、a₂为声波从声源到麦克风之间的传播衰减系数，它们的取值范围为0～1；S(n)表示声源信号；S(...)表示在理想状态下，麦克风接收的声源信号；τ₁、τ₂分别代表声源与麦克风1(参考阵元)和麦克风2(圆周上任意一个麦克风)之间的传播时间；w₁(n)、w₂(n)分别代表麦克风1和2接收信号中的高斯白噪声，且声源信号S(n)与噪声信号w₁(n)、w₂(n)两者之间是互不相关的。

利用阵列麦克风1(参考阵元)与麦克风2的接收信号之间的互相关函数来求广义互相关函数，则两路信号之间的延时就是互相关函数的峰值处。首先使用FIR滤波器对接收的声源信号进行滤波处理，则求得的互相关函数如下：

R₁₂(τ)＝E{x₁(n)x₂(n-τ)} (9)

将式(7)和(8)代入(9)求得：

R₁₂(τ)＝a₁a₂E{S(n-τ₁)S(n-τ₂-τ)}+a₁E{S(n-τ₁)w₂(n-τ)}+a₂E{S(n-τ₂-τ)w₁(n)}+E{ω₁(n)w₂(n-τ)}(10)

w₁(n)、w₂(n)是互不相关的噪声信号，且S(n)与w₁(n)、w₂(n)也互不相关，因此可将式(10)化简如下：

R₁₂(τ)＝a₁a₂E{S(n-τ₁)S(n-τ₂-τ)}＝a₁a₂R_s(τ-(τ₁-τ₂)) (11)

式(9)-(11)中：τ表示麦克风1与麦克风2之间的时间差，即时迟；x₂(n-τ)表示麦克风2相对于麦克风1接收的信号；E{...}表示期望值；R_s表示理想状态下，麦克风接收信号的的互相关函数。

利用麦克风1和麦克风2的互相关函数可以很容易的确定当τ＝τ₁-τ₂时，R₁₂(τ)达到最大值，因此R₁₂(τ)的峰值处对应的τ值即是两个麦克风之间的延时。在短时间内的信号处理中，可能会导致R₁₂(τ)的峰值不够突出，从而导致所求的麦克风之间的时迟τ值的精度下降。因此可以用对滤波之后的信号进行加权处理来制止噪声信号以及混响的不利因素的影响。对x₁(t)和x₂(t)进行FFT变换可得X₁(w)和X₂(w)，则两路输入信号的广义互相关函数为R₁₂(τ)可以表示为：

式(6)中ψ₁₂(ω)是广义互相关加权函数；为麦克风2接收的声源信号的FFT变换的共轭。ω为声源信号的载频；e的上标j表示虚数单位。选择不同的加权函数就得到不同的延时算法，常用的加权函数如表1所示。

表1常用的广义加权函数

根据MATLAB仿真数据可以得出PHAT在噪声环境下具有良好的估计结果，故本系统选择PHAT方法。

麦克风接收的模拟信号经过高速信号采集板转换成数字信号，参考阵元接收信号模型为：x₀(n)，其余各阵元接收信号模型为：x₁(n)、x₂(n)...x_M(n)。利用步骤1-5述方法求出各阵元与参考阵元的时迟值分别为：τ₁、τ₂...τ_M，则麦克则麦克风阵列在空间中一个声源的方向矩阵可表示为：

A(ω,τ)＝[e^jωτ1,e^jωτ2,...,e^jωτM]^T (13)

如果空间信号源频率为f，则(13)可表示为：

A(f,τ)＝[e^j2πfτ1,e^j2πfτ2,...,e^j2πfτM]^T (14)

则麦克风阵列接收信号的矢量形式为：

X(t)＝A(f,τ)x₀(n) (15)

其协方差矩阵估计值为：

R＝X(n)X^H(n)/N (16)

式(9)中：N表示总的采集数据个数，X^H(n)表示X(n)的共轭转置。

对协方差矩阵进行特征值分解，且协方差矩阵有M个不同的特征值，并将特征值进行升序排列：

λ1≥λ2≥…≥λM (17)

与声源信号个数K相等的最大特征值对应的特征矢量所构成的空间看作为信号子空间，把剩下的M-K个特征值对应的特征矢量所构成的空间看作噪声子空间，其中，在实际应用场合，采用信息论方法、平滑秩法、矩阵分解法等估计出声源个数K。

式(12)中：U_S表示信号子空间，U_N表示噪声子空间，表示U_S、U_N的共轭转置，Σ_S为K个较大特征值构成的对角矩阵，而Σ_N为M-K个较小特征值构成的对角矩阵。

扫描θ(0°≤θ≤90°)和求得谱函数值：

根据式(14)求得谱峰值：

Pmax＝max(P_MUSIC) (20)

求得峰值P_max所对应的θ和值，θ和值即为声源的在空间坐标系中的俯仰角θ和方位角

假设均匀圆阵的半径为r＝0.15m，阵元个数为M＝9，快拍数N＝1024，信号源数为P＝1，在信噪比为20dB，构建声源信号为S(t)＝sin(2*pi*f*t)，信号频率为f＝2000Hz，模拟构建各阵元接收输出的信号，单频信号不同方位的DOA估计结果如表2所示。

表2单频信号不同方位的DOA仿真结果(单位:度)

宽带信号模拟实验条件与上述相同，声源信号是用频率分别为f₁＝1000Hz、f₂＝1200Hz、f₃＝1500Hz、f₄＝1800Hz、f₅＝2000Hz构成的混频信号，模拟构建各阵元接收输出的信号，宽带信号不同方位的DOA估计结果如表3所示。

表3宽带信号不同方位的DOA仿真结果(单位:度)

附图说明

下面结合附图和实施例对本发明进一步说明。

图1为本发明一种声源定位方法的均匀圆阵模型图。

图2为本发明一种声源定位方法的改进MUSIC算法的系统框图。

图3为本发明一种声源定位方法的硬件部分组成图。

图4为本发明一种声源定位方法的麦克风阵列结构图。

图5为本发明一种声源定位方法的目标声源方位角为仰角为θ＝12°的实验测量结果。

图6为本发明一种声源定位方法的目标声源方位角为仰角为θ＝48°的实验测量结果。

图7为本发明一种声源定位方法的目标声源定位可视化图。

具体实施方式

图3所示表明本发明一种声源定位方法由麦克风阵列、数据采集卡系统设计和上位机系统三部分构成。首先，由麦克风阵列将声源信号转化为电信号并对目标声源信号进行放大处理，采用数据采集卡进行多通道数据同步采集，最后由上位机系统中基于改进MUSIC算法获取声源的位置信息，并通过上位机系统的外围设备的显示器输出。麦克风阵列与声数据采集卡系统的硬件部分均用导线相互连接。

1)阵列结构设计

为了尽量减小对声场的影响，本设计选取9个MAX9814型的驻极体咪头麦克风和半径为15cm的均匀圆阵搭建声源定位系统，能够极大的减小声源定位所收录的噪声，如图4所示。

2)数据采集

针对同步数据采集时存在读取数据误差大等一些问题，本发明设计基于LabVIEW的多个通道数据采集系统，最终实现麦克风阵列数据同步采集。根据市场调查与麦克风阵列的需求分析，本系统采用PCI-6221多功能数据采集卡，采样频率为250KS/s，采样精度可达到16位，满足系统多通道数据同步采集的需求。

3)发明方法的测试结果

为了测试本发明的改进MUSIC算法在声源定位应用中的有效性，利用该系统对2KHz的单频信号进行一系列的声源定位测试实验，扬声器放置在声源定位系统前方3.5m处，为了确保实验的正确性，每个位置采集30组数据进行分析。声源的方位角和仰角分别为θ＝12°的实验测量结果如图5所示，声源的方位角和仰角分别为θ＝48°实验测量的结果如图6所示，单频信号的不同方位的DOA估计记录在表2中。基于麦克风阵列的声源定位系统可以将视频图像与空间声场生成的图像有效的相融合，使定位可视化，声源定位图像如图7所示。

表4单频信号的不同方位的DOA估计(单位:度)

4)测试结论

根据实际的实验的测量数据分析，信号的波达方向在两个极端的位置(0°≤θ≤10°和80°≤θ≤90°)时所测量结果的误差比较大，主要误差在于声源的仰角，误差在5°左右，仰角在32°≤θ≤72°范围内时，仰角测量结果比较稳定，误差在2°左右。声源的波方位角的误差较为稳定，最大误差为4°，在方位角为180°左右时，误差较小。根据定位系统的测量结果可以得出，改进的MUSIC算法在基于麦克风阵列的声源定位系统应用中，可以有效的解决双峰问题，同时提高了系统的稳定性，使测量结果更加准确。

本发明的具体实施方式仅为一个实例，并非穷举，本领域技术人员根据本发明的启示所做的任何显而易见的改动，都属于本发明权利保护的范围。

Claims (4)

1.一种声源定位方法，其特征在于，包括以下步骤：

(1)麦克风阵列设计：在声源定位空间坐标系中，将M个完全相同的麦克风均匀分布在平面X-Y上一个半径为r的圆周上，在圆心处加上一个阵元，并以此阵元作为参考阵元，建立均匀圆形麦克风阵列模型；

(2)改进声源定位算法，实现步骤如下：

步骤1：假设麦克风接收的信号的表达式如式(1)和式(2)：

x₁(n)＝a₁S(n-τ₁)+w₁(n) (1)

x₂(n)＝a₂S(n-τ₂)+w₂(n) (2)

式(1)中：a₁、a₂为声波从声源到麦克风之间的传播衰减系数，它们的取值范围为0～1；S(n)表示声源信号；S(...)表示在理想状态下，麦克风接收的声源信号；τ₁、τ₂分别代表声源与麦克风1即参考阵元和麦克风2即圆周上任意一个麦克风之间的传播时间；w₁(n)、w₂(n)分别代表麦克风1和2接收信号中的高斯白噪声，且声源信号S(n)与噪声信号w₁(n)、w₂(n)两者之间是互不相关的；

步骤2：使用FIR滤波器对接收的声源信号进行滤波处理，则求得的互相关函数如下：

R₁₂(τ)＝E{x₁(n)x₂(n-τ)} (3)

将式(1)和(2)代入(3)求得：

R₁₂(τ)＝a₁a₂E{S(n-τ₁)S(n-τ₂-τ)}+a₁E{S(n-τ₁)w₂(n-τ)}+a₂E{S(n-τ₂-τ)w₁(n)}+E{ω₁(n)w₂(n-τ)}(4)

w₁(n)、w₂(n)是互不相关的噪声信号，且S(n)与w₁(n)、w₂(n)也互不相关，因此将式(4)化简如下：

R₁₂(τ)＝a₁a₂E{S(n-τ₁)S(n-τ₂-τ)}＝a₁a₂R_s(τ-(τ₁-τ₂)) (5)

式(3)-(5)中：τ表示麦克风1与麦克风2之间的时间差，即时迟；x₂(n-τ)表示麦克风2相对于麦克风1接收的信号；E{...}表示期望值；R_s表示理想状态下，麦克风接收信号的的互相关函数；

对滤波之后的信号进行加权处理来制止噪声信号以及混响的不利因素的影响；对x₁(n)和x₂(n)进行FFT变换得X₁(ω)和X₂(ω)，则两路输入信号的广义互相关函数为R₁₂(τ)表示为：

式(6)中ψ₁₂(ω)是广义互相关加权函数；为麦克风2接收的声源信号的FFT变换的共轭；ω为声源信号的载频；e的上标j表示虚数单位；

步骤4：选择加权函数确定所采用的延时算法，本方法选择PHAT方法；

步骤5：麦克风接收的模拟信号经过高速信号采集板转换成数字信号，参考阵元接收信号模型为：x₀(n)，其余各阵元接收信号模型为：x₁(n)、x₂(n)...x_M(n)；

求出各阵元与参考阵元的时迟值分别为：τ₁、τ₂...τ_M，则麦克风阵列在空间中一个声源的方向矩阵为：

式(7)中：T为矩阵转置操作；

步骤6：求出接收信号的矢量形式，利用空间信号源频率为f，则式(7)化简为：

则麦克风阵列接收信号的矢量形式为：

X(t)＝A(f,τ)x₀(n) (9)

步骤7：其协方差矩阵估计值为：

R＝X(n)X^H(n)/N (10)

式(9)中：N表示总的采集数据个数，X^H(n)表示X(n)的共轭转置；

步骤8：对协方差矩阵进行特征值分解，且协方差矩阵有M个不同的特征值，并将特征值进行升序排列：

λ1≥λ2≥…≥λM (11)

步骤9：与声源信号个数K相等的最大特征值对应的特征矢量所构成的空间看作为信号子空间，把剩下的M-K个特征值对应的特征矢量所构成的空间看作噪声子空间；

式(12)中：U_S表示信号子空间，U_N表示噪声子空间，表示U_S、U_N的共轭转置，Σ_S为K个较大特征值构成的对角矩阵，而Σ_N为M-K个较小特征值构成的对角矩阵；

步骤10：扫描θ(0°≤θ≤90°)和求得谱函数值：

式(13)中是声源在空间坐标系中位置为的阵列响应：

式(14)中m＝0,1,…M-1，r为半径，λ为声源信号的波长；

步骤11：根据式(14)求得谱峰值：

Pmax＝max(P_MUSIC) (15)

步骤12：求得峰值P_max所对应的θ和值，θ和值即为声源的在空间坐标系中的俯仰角θ和方位角

2.按照权利要求1所述的一种声源定位方法，其特征在于：所述麦克风阵列架构是用硬质材料做成的正八边形构架，根据内接圆构建均匀圆阵麦克风阵列。

3.按照权利要求1所述的一种声源定位方法，其特征在于：麦克风阵列中所使用的麦克风是全指向型驻极体麦克风。

4.按照权利要求1所述的一种声源定位方法，其特征在于：选择多功能数据采集卡，并设计基于LabVIEW的多个通道数据采集系统，最终实现麦克风阵列数据同步采集。