一种三维多快拍无网格压缩波束形成声源识别方法与流程

本发明属于声场识别技术领域，具体涉及一种三维多快拍无网格压缩波束形成声源识别方法。

背景技术：

压缩波束形成是实现声源波达方向(direction-of-arrival,doa)估计和强度量化的有效途径。传统压缩波束形成假设声源分布在一组预先设置的离散网格点上，每个网格点代表一个观测方向，施加稀疏约束体现为最小化声源分布向量的l1范数或声源分布矩阵的l2,1范数。上述假设不成立时，声源未落在网格点上，结果准确度降低，该问题称作基不匹配。

目前，基于线性传声器阵列测量的一维和基于平面传声器阵列测量的二维无网格压缩波束形成方法能解决该难题，但这些方法仅能实现三维空间的局部区域内声源的doa估计和源强量化，无法实现三维空间内全方位的声源doa估计和源强量化。

基于线性传声器阵列测量的一维无网格压缩波束形成依赖于一重toeplitz矩阵的vandermonde分解，基于平面传声器阵列测量的二维无网格压缩波束形成依赖于二重toeplitz矩阵的vandermonde分解，基于立体传声器阵列测量的三维无网格压缩波束形成则依赖于三重toeplitz矩阵的vandermonde分解。toeplitz矩阵层级越多，其构建越复杂，其vandermonde分解也越难实现，需要多次生成庞大复杂的排序矩阵来调整toeplitz矩阵层级的顺序。因此，针对三维无网格压缩波束形成，不论在构建toeplitz矩阵，还是在寻找原子范数最小化的等价半正定规划形式，以及在求解toeplitz矩阵的vandermonde分解来估计声源doa方面，都是现有的一维和二维无网格压缩波束形成的处理手段无法实现的。

本发明为立体传声器阵列测量探索出一套完整的三维无网格压缩波束形成的技术手段，用于实现三维空间内声源的全方位doa估计和源强量化。

技术实现要素：

针对现有技术存在的问题，本发明所要解决的技术问题就是提供一种三维多快拍无网格压缩波束形成声源识别方法，它能实现三维空间内全方位的doa估计和源强量化，且能克服基不匹配难题。

本发明的构思是：针对长方体传声器阵列测量，以多快拍数据模型为基础，基于原子范数最小化(atomicnormminimization,anm)重构源产生的信号，并引入矩阵束与配对(matrixpencilandpairing,mapp)方法来估计doa，最终量化源强。

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

步骤1、获取测量声压矩阵p^★

测量信号为：

p^★＝p+n

测量声压矩阵p^★通过传声器阵列测量得到，为噪声干扰，为复数集，a为长方体传声器阵列的行数，b为长方体传声器阵列的列数，c为长方体传声器阵列的层数，l为快拍数目；

步骤2、重构声源在阵列传声器处产生的声压p、获取三重toeplitz矩阵建立重构p的数学模型

式中，是该式计算结果，||·||f表示frobenius范数，p为待求声压矩阵，定义三重toeplitz算子tbb(·)将u映射为a×a维的块toeplitz型hermitian矩阵：

每个块又是一个b×b维的块toeplitz矩阵：

每个块又是一个c×c维的toeplitz矩阵：

为辅助量，“tr(·)”表示求矩阵的迹，上标“^h”表示共轭转置，表示半正定；

利用基于内点法的cvx工具箱中的sdpt3求解器求解上式得源声压及

步骤3、估计声源doa；

步骤4、估计声源强度。

本发明的技术效果是：

实现了实现三维空间内声源的全方位doa估计和源强量化。

附图说明

本发明的附图说明如下：

图1本发明传声器阵列测量布局；

图2本发明仿真成像图；

(a)为2000hz；(b)为3000hz；(c)为4000hz。

具体实施方式

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

本发明包括以下步骤：

步骤1、获取测量声压矩阵p^★

如图1所示，测量用传声器阵列的布局中，符号“●”表示传声器，均匀分布在长方体内，a＝0,1,…,a-1、b＝0,1,…,b-1、c＝0,1,…,c-1分别为x、y、z维上的传声器索引，δx、δy、δz分别为x、y、z维上的传声器间隔。i号声源的doa用(θi,φi)表示，其中，θi为仰角，是从z轴正向到原点与i号声源连线的角度，落于[0°,180°]区间，φi为方位角，是从x轴正向到原点与i号声源连线在xy平面投影的角度，落于[0°,360°]区间。用声源在(0,0,0)号传声器处产生的信号作为源强，第l快拍下i号声源的强度记为si,l，各快拍下i号声源的强度组成的行向量记为l为快拍总数，为复数集。各快拍下声源在(a,b,c)号传声器处产生信号组成的行向量可表示为：

式(1)中，i为声源总数，为虚数单位，t1i≡sinθicosφiδx/λ，t2i≡sinθisinφiδy/λ，t3i≡cosθiδz/λ，λ为波长。

记声源在所有传声器处产生信号组成的矩阵为：

构建列向量其中，上标“t”表示转置，符号表示kronecker乘积，与式(1)相对应有：

存在干扰时，测量信号为：

p^★＝p+n(3)

本发明要解决的问题即以p^★和||n||f为输入来准确获得源的doa和强度。

在仿真试验中，添加噪声为独立同分布高斯白噪声，信噪比(signal-to-noiseratio,snr)定义为snr＝20log10(||p||f/||n||f)，由此可确定||n||f＝||p||f10^-snr/20，其中，“||·||f”表示frobenius范数。

试验测试时，由传声器阵列测量得到。

步骤2、重构声源在阵列传声器处产生的声压p、获取三重toeplitz矩阵令“||·||2”表示l2范数，式(2)可写为：

式(4)所表示信号模型的原子，即构成p的基本单元，为d(t1,t2,t3)ψ，无限势原子集合为：

p的原子范数定义为：

式(6)中，“inf”表示下确界，为正实数集。

用p的原子范数度量源的稀疏性，p的重构问题可写为：

可转化为如下半正定规划：

式(8)中，是该式计算结果，

为辅助量，“tr(·)”表示求矩阵的迹，上标“h”表示共轭转置，表示半正定。定义三重toeplitz算子tbb(·)将u映射为a×a维的块toeplitz型hermitian矩阵：

式(9)中，每个块tα1(0≤α1≤a-1)又是一个b×b维的块toeplitz矩阵：

式(10)中的每个块又是一个c×c维的toeplitz矩阵：

用现成的基于内点法的cvx工具箱中的sdpt3求解器求解式(8)得源声压及

步骤3、估计声源doa

式(8)获得的为半正定三重toeplitz矩阵，其vandermonde分解包含声源doa信息，可基于矩阵束与配对(matrixpencilandpairing,mapp)方法寻找的vandermonde分解，具体步骤如下：

1、特征值分解

式(12)中，为的特征向量构成的酉矩阵，为的特征值构成的对角矩阵。

2、确定声源数目记大于最大特征值的0.01倍(即20db动态范围被考虑)的特征值个数为记为的个较大特征值的平方根构成的对角矩阵，为的个较大特征值对应的特征向量构成的矩阵，令

3、删除y1的后bc行得删除y1的前bc行得计算矩阵束(y1d,y1u)的广义特征值得

4、构建排序矩阵：令为第α个元素为1其它元素为0的列向量，对中的元素重新排序得新的三重toeplitz矩阵，第一重包含b×b个大块，第二重包含a×a个小块，第三重为c×c维矩阵。特征值分解

式(13)中，为的特征向量构成的酉矩阵，为的特征值构成的对角矩阵。记为的个较大特征值的平方根构成的对角矩阵，为的个较大特征值对应的特征向量构成的矩阵，令

5、删除y2的后ac行得删除y2的前ac行得计算矩阵束(y2d,y2u)的广义特征值得

6、令

对中的元素重新排序得新的三重toeplitz矩阵，第一重包含c×c个大块，第二重包含b×b个小块，第三重为a×a维矩阵。特征值分解

式(14)中，为的特征向量构成的酉矩阵，为的特征值构成的对角矩阵。记为的个较大特征值的平方根构成的对角矩阵，为的个较大特征值对应的特征向量构成的矩阵，令

7、删除y3的后ab行得删除y3的前ab行得计算矩阵束(y3d,y3u)的广义特征值得

8.计算下列函数来依次为m＝1,2,…,找到配对的(n,o)：

进而配对和得简记为

9、由得其中，“im”表示取虚部。根据前述式(1)中(t1i,t2i,t3i)与(θi,φi)的关系计算完成源doa估计。

步骤4、估计声源强度

由声源根据步骤1中的公式：t1i≡sinθicosφiδx/λ，t2i≡sinθisinφiδy/λ，t3i≡cosθiδz/λ，

得到感知矩阵

记为量化的各快拍下各声源的强度构成的矩阵，根据式(2)得到声源强度的计算式：

式(16)中，为感知矩阵，上标“+”表示伪逆，为步骤2求得源声压。

仿真模拟试验

为验证建立本发明的准确性、探究本方法发明的可靠性，进行声源识别仿真模拟。

仿真设置如下：在特定位置假设具有特定强度辐射特定频率声波的点声源(具体假设六个声源，doa依次为(45°,90°)、(45°,120°)、(90°,180°)、(120°,180°)、(135°,270°)和(155°,290°)，各快拍下强度的均方根依次为100db、98db、96db、94db、92db和90db(参考2×10^-5pa))。

仿真具体流程为：

1、根据式(2)和式(3)计算各传声器接收声压信号(用a＝b＝c＝7、δx＝δy＝δz＝0.035m的传声器阵列测量声信号，snr取20db，快拍总数设为10)；

2、根据步骤2使用matlab的cvx工具箱中的sdpt3求解器(内点法)求解声压及

3、根据步骤3估计声源doa；

4、根据步骤4估计声源强度。

如图2所示，声源频率为(a)2000hz、(b)3000hz和(c)4000hz时的仿真结果图。每幅图中，“○”指示真实声源分布，“*”指示重构声源分布，并为便于对比，各图中均参考最大输出值进行db缩放，显示动态范围均为0～-20db，同时，每幅图的上方还列出了以标准声压大小2.0×10^-5pa为参考的最大输出值。

由图2看出：本发明的三维无网格压缩波束形成准确估计了每个源的doa并准确量化了每个源的强度。