CN109490833B - 一种改进型传播矩阵的gcc逆模型快速声源识别方法 - Google Patents

Abstract

本发明公开了一种改进型传播矩阵的GCC逆模型快速声源识别方法，它包括：骤1、建立GCC逆模型；步骤2、改进型传播矩阵求解：

M为麦克风总数，i代表第i号麦克风，k代表第k号麦克风，Δt_ki,n为声音从网格点n传播到麦克风对(i，k)的时间差，Δt_ki,l为声音从网格点l传播到麦克风对(i，k)的时间差，w是声源的圆频率，w_max为上限频率；步骤3、GCC逆模型求解，重构声源强度分布并成像。本发明解决了现有传播矩阵对时间阈值的依赖性和对实验配置的单一适应性，能自动随实验配置相应变化，表现出有很好的自适应性。

Description

一种改进型传播矩阵的GCC逆模型快速声源识别方法

技术领域

本发明属于声场识别技术领域，具体涉及一种GCC逆模型快速声源识别方法。

背景技术

基于麦克风阵列信号处理的时域波束形成在宽带声源识别领域占据不可缺少的地位，常用的方法之一是GCC逆模型算法。该算法在GCC成像结果和理论噪声图像之间建立线性逆模型，通过最小化两者之差获取高空间分辨率的声源识别结果，GCC逆模型算法的关键是求取用于建立逆模型的传播矩阵。现有方法中，该矩阵的求解过程中需要人为设定时间阈值。由于该时间阈值对声源识别结果相当敏感，且该阈值随不同的实验配置而变化，进而使得求解的传播矩阵不具有自适应性，因此实际应用中对最佳时间阈值的选择存在困难。

发明内容

针对现有技术存在的问题，本发明所要解决的技术问题就是提供一种改进型传播矩阵的GCC逆模型快速声源识别方法，它能避免矩阵求解过程中设置时间阈值，自适应不同实验配置，通过该矩阵建立GCC逆模型，在时域下快速获得高分辨率的宽带噪声图像。

本发明所要解决的技术问题是通过这样的技术方案实现的，它包括以下步骤：

步骤1、建立GCC逆模型

采用Combo阵列波束形成的坐标系，x中任意位置声源的相对能量辐射到b中的每一个位置，存在一个传播矩阵A∈R_N×N，则：b＝Ax

列向量b∈R_N×1为N个网格点的计算得到的相对能量输出，x∈R_N×1为点声源在N个网格点上的真实相对能量分布；A的矩阵形式为：

矩阵元素a_n,l中，n、l为网格上任意两个点，N为网格点总数；

在真实声源数量少于网格点数和存在噪声干扰的条件下，则该式的稀疏解为：

x是稀疏解，σ为干扰的二范数，||·||₁、||·||₂分别是矩阵或向量的一范数和二范数，s.t.代表约束条件；

步骤2、改进型传播矩阵求解

式中，M为麦克风总数，i代表第i号麦克风，k代表第k号麦克风，Δt_ki,n为声音从网格点n传播到麦克风对(i，k)的时间差，Δt_ki,l为声音从网格点l传播到麦克风对(i，k)的时间差，w是声源的圆频率，w_max为上限频率；

步骤3、GCC逆模型求解，重构声源强度分布并成像。

本发明的技术效果是：

本方法发明在等强度源、不等强度源、和相干源等条件下均可清晰化声源图像，准确定位声源位置，有效抑旁瓣鬼影的产生，改进型传播矩阵在GCC逆模型声源识别中性能优越；本发明的改进型传播矩阵解决了现有传播矩阵对时间阈值的依赖性和对实验配置的单一适应性，能自动随实验配置相应变化，表现出有很好的自适应性。

附图说明

本发明的附图说明如下：

图1为Combo阵列波束构成的坐标系；

图2为现有传播矩阵的时间阈值影响声源识别的测试图；

图3为现有传播矩阵时间阈值对不同实验配置适应性的测试图；

图4为仿真改进型传播矩阵对不同声源模式识别的效果图；

图5为验证试验的布局图；

图6为不同位置处扬声器声源的试验识别成像图。

具体实施方式

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

本发明包括以下步骤：

步骤1、建立GCC逆模型

本实施例采用图1所示Combo阵列波束形成的坐标系，图1中，“*”所在的平面为网格平面，“+”是声源位置，声源处于网格点上，原点位于阵列中心，用实心·表示，空心o代表传声器位置，麦克风总数M为36，假设声源位于网格点上，N为网格点总数，n＝1,2,,N为网格点索引，网格点l处有一单极子源，任意网格点n处得到GCC波束形成输出：

式(1)中，i代表第i号麦克风，k代表第k号麦克风，Δt_ki,n为声音从网格点n传播到麦克风对(i，k)的时间差，w是声源的圆频率，C_ik(w)为麦克风对接收声压的互功率谱，j为虚数单位，W(w)是频率加权函数。本发明基于PHAT加权推导，故选用PHAT加权：

式(2)中，|·|表示取模。

将N个网格点的相对能量输出排列为一个列向量b∈R_N×1，而点声源在N个网格点上的真实相对能量分布为x∈R_N×1(这里是假设每个网格点是一个点声源，如果该点不存在声源，x中对应元素值为0，也就是稀疏约束)，x中任一位置的相对能量可辐射到b中的每一个位置，可将该系统考虑成一个线性系统，则存在一个传播模型矩阵(简称传播矩阵)A∈R_N×N，使下式成立：

b＝Ax (3)

考虑到真实声源数量少于网格点数和存在噪声干扰，则(3)式的解x是稀疏的，即可变换到下式：

式(4)中，x是稀疏解，σ为干扰的二范数，||·||₁、||·||₂分别是矩阵或向量的一范数和二范数，s.t.代表约束条件。

式(3)和式(4)为GCC逆模型。

式(4)是一个L1范数最小化问题，可采用贪婪算法和迭代收缩算法求解，但其关键是传播矩阵A的求取。现有方法中，A只与声速、扫描网格点的位置和麦克风的位置有关，A的矩阵形式为：

式(5)中，a_n,l由下式求得：

式(6)中，ε为时间阈值，但现有方法未给出ε的选取方法。

通过大量仿真结果，总结出本实施例实验配置对应的最优ε_op，选取不同ε＝βε_op，β＝0,0.5,1,1.5,2获取声源识别结果，如图2所示，采用不同阈值的识别结果表现出很大差异性，说明该值对声源识别结果相当敏感，严重时会造成识别混乱。

本实施例仿真了ε_op对不同实验配置的适用性，如图3所示：图3a显示了ε_op(识别距离1m)对于1m识别距离的识别结果，图3b展示了ε_op(识别距离1m)对于2m识别距离的识别结果，图3c展示了ε_op(识别距离2m)对于2m识别距离的识别结果，图3d展示了ε_op(识别距离2m)对于1m识别距离的识别结果。从图3看出：单个ε_op仅适应对应的实验配置，而对其他实验配置达不到最优效果，进而导致所求传播矩阵不具备自适应性。

为解决上述问题，本发明采用以下步骤。

步骤2、改进型传播矩阵求解

式(3)是一个线性系统，能量映射关系可以采用下式表示：

为了与现有的传播矩阵A区别，改进型传播矩阵用

表示。

本发明采用PHAT加权，消除互谱的幅值信息，仅保留相位信息。因此这里得到的波束形成输出b只能代表点声源的相对能量，由此反算出的真实点源分布x也只能代表源的相对能量。由于系统是线性的，b可以看成是每个声源通过

的映射关系到各个网格点的能量叠加。网格索引l处的声源相对能量x_l通过

的第l列映射关系

辐射到网格点的相对能量为b的第l列[b_1,l b_2,l … b_n,l…b_N,l]^T。因此，将

的每一列视为相应声源相对能量辐射到各个网格点的辐射系数。假设系统只有单独一个声源可算出

相应一列的辐射系数a_n,l，可用公式a_n,l＝b_n,l/x_l计算。

考虑空间l处有一单极子声源，则第i、k麦克风对的互功率谱可表示为复指数形式：

C_ik(w)＝|C_ik(w)|exp(-jwΔt_ki,l) (8)

Δt_ki,l为声音从网格点n传播到麦克风对(i，k)的时间差。

将式(8)带入式(1)得l处声源辐射到网格点n处的相对能量为:

在实际数据处理中，由于采样频率限制，存在一个上限频率w_max。通过欧拉变换，式(9)中被积函数的实部是一个关于w的偶函数，虚部是一个关于w的奇函数。根据奇函数、偶函数在对称区间的积分特性，式(9)可以写成：

当声源与网格点处于同一点位置时，辐射能量不会衰减，那么网格点b_l,l上的相对能量可代表相同位置处声源相对能量x_l，则可求得：

式(11)中，M为麦克风总数，i代表第i号麦克风，k代表第k号麦克风，Δt_ki,n为声音从网格点n传播到麦克风对(i，k)的时间差，Δt_ki,l为声音从网格点l传播到麦克风对(i，k)的时间差，w是声源的圆频率，w_max为上限频率。

步骤3、GCC逆模型求解

为了验证式(11)所表达的改进型传播矩阵的性能，通常采用贪婪算法和迭代收缩算法求解，贪婪算法代表OMP，迭代收缩算法代表LS1。本实施例采用贪婪算法的代表OMP对GCC逆模型进行求解，重构声源强度分布并成像。

该算法求解过程见“正交匹配追踪应用于反褶积方法的声源映射反问题”，Thomas，美国声学学会，第138卷第6期，第3678～3685页，2015年记录了OMP的求解步骤：

(The different steps of OMP can be summarized as follows：sss

(1)Initialize i＝1 and the residual r^(i-1)＝Y；

(2)Set Q 1/4 0 and the initial set of source indices J^(i-1)＝Φ；

(3)Search the scan point index such that:j＝arg max|A^H|r^(i-1)|；

(4)Update the index set:J⁽ⁱ⁾＝J^(i-1)U j；

(5)Compute:P_J(i)＝A_J(i)A

_J(i),where A_J(i)is the propagation matrixrestricted to sources defined by the index set J⁽ⁱ⁾；

(6)Compute the residual:r⁽ⁱ)＝r^(i-1)-P_J(i)r^(i-1)；

(7)i＝i+1.Go back to step(3)until a stop criterion is reached.

The final source vector is given by Q_J＝A

_JY,where A

_J(i)is the Moore-Penrose pseudo-inverse and J is the final index set.)。

在OMP中，迭代搜索波束形成输出b与传播矩阵

每一列元素的内积的最大值，本实施例设置迭代终止条件为残差的L2范数达到0.01。

仿真模拟试验

为验证本发明的改进型传播矩阵在声源识别中的可靠性能，进行仿真模拟，具体流程为：

1、按照图1进行仿真基础设置，在“+”位置假设具有特定强度辐射零均值标准差为1的高斯白噪声点声源，声源平面距离阵列平面1m；

2、根据式(1)和式(2)正向计算GCC波束形成输出；

3、根据式(3)和式(4)建立GCC逆模型；

4、根据式(11)求解传播矩阵；

5、基于OMP分别采用传播矩阵A和改进型传播矩阵

求解该逆模型重构声源强度分布并成像。

这里，采用几种不同声源模式(四个等强度声源(相干或不相干)和四个不等强度声源)，采样频率为4.4kHz，采用16384个时域采样点来计算接收信号的互功率谱。

为便于对比，本专利申请中各图中均参考最大输出值进行dB缩放，显示动态范围均为0～-15dB。该仿真模拟的效果如图4所示：图4a列为传统GCC声源识别结果；图4b列为采用未改进A矩阵(时间阈值ε＝ε_op)通过OMP求解逆模型的声源识别结果；图4c列为采用改进型

矩阵通过OMP求解逆模型的声源识别结果；图4d列为采用未改进A矩阵(时间阈值ε＝0.5ε_op)通过OMP求解逆模型的声源识别结果。

从图4看出：本发明中的改进型传播矩阵适用于等强度声源(相干或不相干)和四个不等强度声源的识别，避免设置时间阈值，可准确定位声源位置，抑制旁瓣鬼影，其声源识别性能与采用最优时间阈值ε_op时相当。于综上所述，本发明的性能优越，避免设置时间阈值，自适应性好。

验证试验

验证试验的布局如图5所示，采用白噪声信号激励的扬声器作为声源，采用Brüel&Kjaer公司36通道combo阵列采样声压信号。各传声器接收的声压信号经PULSE 3560D型数据采集系统同时采集并传输到PULSE LABSHOP中进行频谱分析，得声压信号的互功率谱，采样频率为65536Hz，信号加汉宁窗，采用64段平均、66.7％的重叠率，频率分辨率为8Hz。进一步，采用MATLAB编制的程序计算各聚焦点的输出量并成像。

图6为试验扬声器声源的识别成像图：图6a为传统GCC声源识别试验结果；图6b为采用未改进A矩阵(时间阈值ε＝ε_op)通过OMP求解逆模型的声源识别试验结果；图6c为采用改进型

矩阵通过OMP求解逆模型的声源识别试验结果；图6d为采用未改进A矩阵(时间阈值ε＝0.5ε_op)通过OMP求解逆模型的声源识别试验结果。试验结果与图4保持了完全一致的结果，本发明方法是可靠的。

Claims (1)

1.一种改进型传播矩阵的GCC逆模型快速声源识别方法，其特征是，包括以下步骤：

步骤1、建立GCC逆模型

采用Combo阵列波束形成的坐标系，x中任一点的相对能量辐射到b中的每一个点，存在一个传播矩阵A∈R_N×N，则：b＝Ax

列向量b∈R_N×1为N个网格点计算得到相对能量输出，x∈R_N×1为点声源在N个网格点上的真实相对能量分布；A的矩阵形式为：

矩阵元素a_n,l中，n、l为网格上任意两个点，N为网格点总数；

在真实声源数量少于网格点数和存在噪声干扰的条件下，则该式的稀疏解为：

x是稀疏解，σ为干扰的二范数，||·||₁、||·||₂分别是矩阵或向量的一范数和二范数，s.t.代表约束条件；

步骤2、改进型传播矩阵求解

式中，M为麦克风总数，i代表第i号麦克风，k代表第k号麦克风，Δt_ki,n为声音从网格点n传播到麦克风对(i，k)的时间差，Δt_ki,l为声音从网格点l传播到麦克风对(i，k)的时间差，w是声源的圆频率，w_max为上限频率；

步骤3、GCC逆模型求解，重构声源强度分布并成像。

Images