CN109932682B - 二维多快拍无网格压缩波束形成声源识别方法 - Google Patents

Abstract

本发明公开了一种二维多快拍无网格压缩波束形成声源识别方法，包括以下步骤：步骤1、获取测量声压矩阵P^★；步骤2、迭代求解下式：

步骤3、估计声源DOA；步骤4、估计声源强度。本发明的技术效果是：由于采用IPM计算ANM的计算复杂度为O((N_u+AB×L+L²)²(AB+L)^2.5),本发明的计算复杂度为O((AB+L)³)，二者相比，本发明明显降低。本发明能准确估计声源DOA和量化声源强度，且现有技术相比，计算效率显著提高。

Description

二维多快拍无网格压缩波束形成声源识别方法

技术领域

本发明属于声场识别技术领域。

背景技术

基于平面传声器阵列测量的压缩波束形成是实现声源二维波达(Direction-of-arrival,DOA)方向估计和强度量化的有效途径。无网格压缩波束形成基于的原子范数最小化(Atomic Norm Minimization，ANM)为凸优化问题，可转化为半正定规划后利用现成的CVX工具箱中的SDPT3求解器求解，该求解器采用内点方法(Interior Point Method,IPM)，该方法适用于小维度矩阵问题，对大维度矩阵问题则比较耗时。二维多快拍无网格压缩波束形成中，平面阵列传声器数目多及快拍数目多的因素使其形成大维度矩阵，IPM耗时问题较为明显，因此，有必要在此基础上探索一种提高计算效率的新方法。

发明内容

针对现有技术存在的问题，本发明所要解决的技术问题就是提供一种二维多快拍无网格压缩波束形成声源识别方法，它既能准确估计声源DOA和量化声源强度，还能减少运算时间，提高计算效率。

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

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

测量声压矩阵

为：

P^★＝P+N

测量声压矩阵P^★通过传声器阵列测量得到，

为噪声干扰，

为复数集，A为矩形传声器阵列的行数，B为矩形传声器阵列的列数，L为快拍次数；

步骤2、重构声源在阵列传声器处产生的声压P；

步骤1)、建立重构P的数学模型

式中，

是该式计算结果，||·||_F表示Frobenius范数，P为待求声压矩阵，Z为辅助矩阵，τ为规则化参数：

T_b(·)为二重Toeplitz算子，T_b(u)为AB×AB维的Hermitian二重Toeplitz矩阵，将任意给定向量u映射为A×A维的块Toeplitz型Hermitian矩阵：

T_b(u)中，每个块T_a(0≤a≤A-1)都是一个B×B维的Toeplitz矩阵：

u、E为辅助量，N_u＝(A-1)(2B-1)+B，E为Hermitian矩阵，“tr(·)”表示求矩阵的迹，上标“^H”表示共轭转置，“≥0”表示半正定；

步骤2)、求解P

引入Hermitian拉格朗日乘子矩阵

和惩罚参数ρ，q为迭代次数索引，I₁和I₂分别为L×L和AB×AB维的单位矩阵，初始化Z⁰＝Λ⁰＝0，引入如下矩阵划分：

进行q+1次迭代时具有如下变量更新：

上述迭代式中：

“diag(·)”表示形成以括号内向量为对角元素的对角矩阵，

表示Kronecker乘积，

为实数集，M^-1表示矩阵M的逆矩阵，

为T_b(·)的伴随算子；对任意给定矩阵

Θ_k为第k个对角线的元素均为1其它元素均为0的基本Toeplitz矩阵，

是(A-1,B-1)的半空间，

意味着k₁＝0时，k₂∈{0,1,…,B-1}，k₁∈{1,2,…,A-1}时，k₂∈{-(B-1),-(B-2),…,B-1}；

当连续两次迭代的u间的相对变化量和P间的相对变化量均小于10^-4，或者最大迭代次数被完成时，迭代终止；

步骤3、估计声源DOA；

步骤4、估计声源强度。

本发明的技术效果是：

根据文献Z.Yang,J.Li,P.Stoica,L.Xie.Sparse methods for direction-of-arrival estimation[J].arXiv Preprint arXiv:1609.09596[cs.IT],2017:1-65.(杨在，李剑，P.Stoica，谢利华，“用于波达方向估计的稀疏方法[J]”.arXiv预印arXiv:1609.09596[cs.IT],2017:1-65.)中记载：基于IPM的SDPT3求解器的计算复杂度可写为

n₁为求解变量的数目，n₂×n₂为正半定矩阵的维度；由于本发明仅需要计算特征值分解，其计算复杂度可写为

n₃×n₃为进行特征值分解的矩阵的维度；在二维多快拍无网格压缩波束形成中，n₁＝N_u+AB×L+L²，n₂＝AB+L，n₃＝AB+L。

所以采用IPM计算ANM的计算复杂度为O((N_u+AB×L+L²)²(AB+L)^2.5)；本发明给出的计算复杂度为O((AB+L)³)，二者相比，本发明明显降低，相应地，计算效率显著提高。

本发明能准确估计声源DOA和量化声源强度，且现有技术相比，计算效率显著提高。

附图说明

本发明的附图说明如下：

图1为传声器阵列测量布局；

图2为声源频率为2000Hz、3000Hz和4000Hz时的仿真重构结果对比图；

图2(a)、(c)、(e)对应IPM；图2(b)、(d)、(f)对应本发明；

图3为试验布局图；

图4为声源频率为2000Hz、3000Hz和4000Hz时的实验声场重构结果对比图；

图4(a)、(c)、(e)对应IPM；图4(b)、(d)、(f)对应本发明。

具体实施方式

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

本发明包括以下步骤：

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

二维无网格压缩波束形成声源识别是利用矩形传声器阵列测量声信号。如图1所示的传声器阵列测量布局，符号“●”表示传声器，a＝0,1,…,A-1，b＝0,1,…,B-1分别为x、y维的传声器索引，Δx、Δy分别为x、y维的传声器间隔，θ_i、φ_i分别为i号声源DOA的仰角和方位角(0°≤θ≤90°、0°≤φ≤360°)。记

为各快拍下i号声源的强度(声源在(0,0)号传声器处产生的声压)组成的行向量，l＝1,2,…,L为快拍索引，s_i,l为第l快拍下i号声源的强度，

为复数集。

假设声源辐射平面声波，波长为λ，t_1i≡sinθ_i cosφ_iΔx/λ，t_2i≡sinθ_i sinφ_iΔy/λ，则各快拍下声源在(a,b)号传声器处产生声压组成的行向量

可表示为：

式(1)中，k为声源总数，

为虚数单位。

构建矩阵：

列向量

标量

行向量

其中，上标“^T”表示转置，符号

表示Kronecker乘积，“||·||₂”表示

范数，

为正实数集，||ψ_i||₂＝1，与式(1)对应，有：

存在噪声干扰

时，测量声压矩阵

可表示为：

P^★＝P+N (3)

在仿真试验中，添加噪声为独立同分布高斯白噪声，信噪比(Signal-to-noiseRatio,SNR)定义为SNR＝20log₁₀(||P||_F/||N||_F)，由此可确定||N||_F＝||P||_F10^-SNR20，其中，“||·||_F”表示Frobenius范数。

试验测试时，

由传声器阵列测量得到。

步骤2、重构声源在阵列传声器处产生的声压P

步骤1)、建立重构P的数学模型二维多快拍无网格压缩波束形成后处理的第一步是滤除测量声压P^★中的噪声、重构声源在阵列传声器处产生的声压P，该步骤通过对声源施加稀疏约束来实现。无网格设置时，所有D(θ_i,φ_i)＝d(t₁,t₂)ψ为i号声源到阵列传声器声压的传递矩阵，

称为原子集合。用P的原子范数来量度声源稀疏性，定义式为：

式(4)中，“inf”表示下确界。上述P的重构问题可写为：

式(5)中，ε为噪声控制参数，通常取为||N||_F。式(5)可转化为如下半正定规划：

式(6)中，

表示计算结果，P为待求声压矩阵；u、E为辅助量，N_u＝(A-1)(2B-1)+B，E为Hermitian矩阵，“tr(·)”表示求矩阵的迹，上标“^H”表示共轭转置，“≥0”表示半正定，T_b(·)为二重Toeplitz算子，T_b(u)为AB×AB维的Hermitian二重Toeplitz矩阵，将任意给定向量

映射为A×A维的块Toeplitz型Hermitian矩阵：

式(7)中，每个块T_a(0≤a≤A-1)都是一个B×B维的Toeplitz矩阵：

式(6)表示的原子范数最小化的等价半正定规划为标准凸优化问题，可利用CVX工具箱中的SDPT3求解器求解式(6)，该求解器采用内点方法(IPM，S.Boyd,L.Vandenberghe.Convex Optimization[M].Cambridge:Cambridge University Press,2004:1-684.(S.Boyd,L.Vandenberghe.凸优化[M].剑桥:剑桥大学出版社,2004:1-684.))，适用于小维度矩阵问题，对大维度矩阵问题则比较耗时。

步骤2)、使用本发明求解式(6)，重构声源声压P、获得

针对现有的基于IPM的CVX工具箱中的SDPT3求解器，本发明提供一种替代方法，将式(6)重写为：

其中，Z为辅助矩阵，τ为规则化参数；优选

引入Hermitian拉格朗日乘子矩阵

和惩罚参数ρ，本发明迭代求解式(9)。令q为迭代次数索引，I₁和I₂分别为L×L和AB×AB维的单位矩阵，初始化Z⁰＝Λ⁰＝0，引入如下矩阵划分，即把每次计算得到的Λ，Z分成四个小矩阵：

进行q+1次迭代时采用如下变量更新：

上述迭代式中：

“diag(·)”表示形成以括号内向量为对角元素的对角矩阵，

表示Kronecker乘积，

为实数集，M^-1表示矩阵M的逆矩阵，

为T_b(·)的伴随算子，对任意给定矩阵

Θ_k为第k个对角线的元素均为1其它元素均为0的基本Toeplitz矩阵，

是(A-1,B-1)的半空间，

意味着k₁＝0时，k₂∈{0,1,…,B-1}，k₁∈{1,2,…,A-1}时，k₂∈{-(B-1),-(B-2),…,B-1}。

式(15)为将Hermitian矩阵

投影到半正定锥上，可通过特征值分解该矩阵并令所有负特征值为零来实现。

令ρ＝1，当连续两次迭代的u间的相对变化量(||u^q-u^q-1||₂/||u^q-1||₂)和P间的相对变化量(||P^q-P^q-1||_F/||P^q-1||_F)均小于10-4或者最大迭代次数被完成时，迭代终止。选择最大迭代次数为1000。

通过上述迭代得到

并计算出

根据式(7)和式(8)得出

步骤3、估计声源DOA；

该步骤的运算与现有二维多快拍无网格压缩波束形成声源识别方法相同。参考文献Z.Yang,L.Xie,P.Stoica.Vandermonde decomposition of multilevel Toeplitzmatrices with application to multidimensional super-resolution[J].IEEETransactions on Information Theory,2016,62(6):3685-3701.(“多重Toeplitz矩阵的范德蒙德分解及其在多维超分辨率中的应用[J].”，杨在,谢利华,P.Stoica，《IEEE信息理论学报》，2016,62(6):3685-3701)。

获得声源个数

及声源

步骤4、估计声源强度；

由声源

根据步骤1中的公式：t_1i≡sinθ_i cosφ_iΔx/λ，t_2i≡sinθ_i sinφ_iΔy/λ，

得到感知矩阵

记

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

的计算式：

式中，A为感知矩阵，上标“⁺”表示伪逆，

为使用本发明迭代求解所得声压。

仿真模拟试验

为验证建立本发明的准确性和对比其性能的提高，进行声源识别仿真模拟。

仿真设置如下：在特定位置假设具有特定强度辐射特定频率声波的点声源(假设四个声源，仰角依次为45°、20°、70°和45°，方位角依次为90°、200°、160°和270°，各快拍下声源强度的均方根依次为100dB、97dB、94dB和90dB(参考2.0×10^-5Pa))。

仿真具体流程为：

1、根据式(2)及式(3)计算各传声器接收声压信号(用A＝B＝8、Δx＝Δy＝0.035m的传声器阵列测量声信号，信噪比SNR取20dB，快拍总数设为10)；

2、根据步骤2的本方法发明，由原子范数最小化的等价半正定规划，得到

和

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

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

本方法发明与用IPM的结果进行对比，如图2所示，声源辐射声波的频率：图2(a)和图2(b)为2000Hz，图2(c)和图(d)为3000Hz，图2(e)和图2(f)为4000Hz；每幅图中，“○”指示真实声源分布，“*”指示重构声源分布，并为便于对比，各图中均参考最大输出值进行dB缩放，显示动态范围均为0～-20dB，同时，每幅图的上方还列出了以标准声压大小2.0×10^-5Pa为参考的最大输出值。

由图2可以看出：使用IPM的二维多快拍无网格压缩波束形成(图2(a)、图2(c)、图2(e))能准确重构2000Hz、3000Hz和4000Hz时的声源分布。使用本发明的二维多快拍无网格压缩波束形成(图2(b)、图2(d)、图2(f))在三个频率下的重构能达到相同效果，显见，四个声源的DOA均被准确估计、强度均被准确量化。

计算在2.2GHz Intel(R)Core(TM)i5-5200U的CPU上由MATLAB R2014a完成，IPM用时约9s，本发明用时约0.9s，可见本发明的计算耗时仅为IPM的计算耗时的1/10。此外，改变声源位置、声源数目、频率、SNR等参数进行仿真，获得类似的结果。可见采用本发明的二维多快拍无网格压缩波束形成能准确估计声源DOA和量化声源强度；本发明克服了IPM的计算速度慢的缺陷，计算效率显著提高。

试验验证

为检验仿真结论的正确性和使用本发明的二维多快拍无网格压缩波束形成在实际应用中的有效性，进行验证试验。

图3为试验布局图，测量阵列为

公司的64个4958型传声器均匀分布的矩形阵列(A＝B＝8，Δx＝Δy＝0.035m)，声源为稳态白噪声信号激励的2个扬声器，从左至右，笛卡尔坐标依次为(2.24,0,5)m和(-2.24,0,5)m，仰角依次为24.13°和24.13°，方位角依次为0°和180°，由于试验在半消声室内进行，地面存在反射，扬声器关于地面还存在2个镜像声源，笛卡尔坐标依次为(2.24,-2.2,5)m和(-2.24,-2.2,5)m，仰角依次为32.13°和32.13°，方位角依次为315.52°和224.48°。各传声器测得的声压经PULSE 3560D型数据采集***同时采集并传输到PULSE LABSHOP中进行快速傅里叶变换，得声压频谱，采样频率为16384Hz，每个快拍长1s、包含2¹⁴个采样点，10个快拍被采用。进一步，采用IPM和本发明的二维多快拍无网格压缩波束形成后处理声压频谱来重构声源分布，其余设置均与仿真一致。

图4给出了2000Hz、3000Hz和4000Hz频率下采用图1所示64个传声器均匀分布的矩形阵列时的声源分布重构结果，每幅图中，“*”指示重构声源分布，参考最大输出值进行dB缩放，“○”指示真实声源DOA，不包含强度信息。显见，采用IPM(图4(a)、4(c)、4(e))和本发明(图4(b)、4(d)、4(f))的二维多快拍无网格压缩波束形成均准确估计了这些声源的DOA。IPM用时约6s，本发明用时约1s。

试验结论与仿真结论一致，证明仿真结论正确，本发明的计算效率显著更高。

Claims (2)

1.二维多快拍无网格压缩波束形成声源识别方法，其特征是，包括以下步骤：

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

测量声压矩阵P^★∈C^AB×L为：

P^★＝P+N

测量声压矩阵P^★通过传声器阵列测量得到，N∈C^AB×L为噪声干扰，C为复数集，A为矩形传声器阵列的行数，B为矩形传声器阵列的列数，L为快拍次数；

步骤2、重构声源在阵列传声器处产生的声压P；

步骤1)、建立重构P的数学模型

式中，

是该式计算结果，||·||_F表示Frobenius范数，P为待求声压矩阵，Z为辅助矩阵，τ为规则化参数：

T_b(·)为二重Toeplitz算子，T_b(u)为AB×AB维的Hermitian二重Toeplitz矩阵，将任意给定向量u映射为A×A维的块Toeplitz型Hermitian矩阵：

T_b(u)中，每个块T_a，0≤a≤A-1，都是一个B×B维的Toeplitz矩阵：

u、E为辅助量，N_u＝(A-1)(2B-1)+B，E为Hermitian矩阵，“tr(·)”表示求矩阵的迹，上标“H”表示共轭转置，“≥0”表示半正定；

步骤2)、求解P

引入Hermitian拉格朗日乘子矩阵Λ∈C^{(AB+L)×(AB+L)}和惩罚参数ρ，q为迭代次数索引，I₁和I₂分别为L×L和AB×AB维的单位矩阵，初始化Z⁰＝Λ⁰＝0，引入如下矩阵划分：

进行q+1次迭代时具有如下变量更新：

上述迭代式中：

“diag(·)”表示形成以括号内向量为对角元素的对角矩阵，

表示Kronecker乘积，R为实数集，M^-1表示矩阵M的逆矩阵，

为T_b(·)的伴随算子；对任意给定矩阵A∈C^AB×AB，

Θ_k为第k个对角线的元素均为1其它元素均为0的基本Toeplitz矩阵，H是(A-1,B-1)的半空间，(k₁,k₂)∈H意味着k₁＝0时，k₂∈{0,1,…,B-1}，k₁∈{1,2,…,A-1}时，k₂∈{-(B-1),-(B-2),…,B-1}；

当连续两次迭代的u间的相对变化量和P间的相对变化量均小于10^-4，或者最大迭代次数被完成时，迭代终止；

步骤3、估计声源DOA；

步骤4、估计声源强度。

2.根据权利要求1所述的二维多快拍无网格压缩波束形成声源识别方法，其特征是，

所述的最大迭代次数为1000。

Images