CN111965599A - 二维动态网格压缩波束形成的声源识别方法 - Google Patents

Abstract

本发明涉及声波识别领域，适用于任意的平面传声器阵列，能够克服基不匹配的问题，获取超高精度的声源成像，具体涉及一种二维动态网格压缩波束形成的声源识别方法，包括以下步骤：1)建立压缩感知问题模型；2)构造替代函数；3)求解替代函数；4)通过步骤2)与步骤3)反复循环，构造和迭代缩减替代函数，使网格坐标θ和

逐渐收敛至真实声源所在的位置，源强分布矢量q也在这个过程中得到愈加准确的估计。本发明的优点在于，能够克服传统压缩波束形成存在的基不匹配问题，对不在网格点上的声源也能够准确定位，其性能显著优于传统压缩波束形成，满足了实际应用中获取超高精度声源识别成像的需求。

Description

二维动态网格压缩波束形成的声源识别方法

技术领域

本发明涉及声波识别领域，具体涉及一种二维动态网格压缩波束形成的声源识别方法。

背景技术

波束形成作为一种信号处理技术，通过对阵列传声器测量声信号进行空间滤波从而实现声源波达方向(Direction-of-arrival,DOA)的估计和源强量化，其中基于平面传声器阵列的二维波束形成得到了最为广泛的应用。压缩波束形成是近年来提出并得到迅速发展的一种新型波束形成方法。其基于压缩感知理论和主要声源通常为稀疏分布的事实，建立起阵列传声器测量声压信号和预设声源网格点处未知源强之间的欠定线性方程组，并利用l₁范数最小化、正交匹配追踪(Orthogonal matching pursuit,OMP)等稀疏促进算法求解上述方程组从而获得声源DOA和源强的准确估计，即便是在传声器数量很少、测量干扰较大的情况下也能获得准确的声源识别结果。然而在实际应用中，DOA几乎不会与预设的网格点重合，即出现基不匹配问题。基不匹配会劣化DOA估计性能，同时使源强量化出现误差。细化网格点虽能一定程度上缓解基不匹配缺陷，但以加重计算量为代价，且过密网格点会增大感知矩阵的列相干性反而可能劣化声源识别性能。因此，实际应用中传统压缩波束形成不能满足超高精度的声源识别要求。

发明内容

本发明的目的是针对现有技术对应的不足，提供一种二维动态网格压缩波束形成的声源识别方法，适用于任意的平面传声器阵列，能够克服基不匹配的问题，获取超高精度的声源成像。

本发明的目的是采用下述方案实现的：一种二维动态网格压缩波束形成的声源识别方法，包括以下步骤：

1)建立压缩感知问题模型：

1-1)以由若干传声器组成的平面阵列的中心为原点建立三维坐标系，用于声源识别试验；

1-2)在三维坐标系中，设置任意声源i的入射方向与z轴的夹角为θ_i，声源i的入射方向在x-y平面上的投影与x轴的夹角为

且0°≤θ_i≤90°，

得到声源i的声源坐标

用于刻画声源i的波达方向；

1-3)在三维坐标系中，所述平面阵列中的所有传声器均在x-y平面上，设置任一传声器的坐标为(x,y)，则声源i到坐标为(x,y)的任一传声器之间的传递函数G(x,y)用下列公式(1)表示：

式中，

为复数单位，

λ为声源i的声波波长，e为自然常数，π为圆周率，

为复数集；

1-4)设置声源i在三维坐标系的原点产生的声压复幅值为q_i，用于表示声源i的源强，所述平面阵列中的任一传声器测得声源i的声压信号p(x,y)用公式(2)表示：

式中，q_i为声源i在坐标原点产生的声压复幅值，L为声源的总个数，G(x,y)为声源i到任一坐标为(x,y)的传声器之间的传递函数；

1-5)设置在平面阵列中的所有传声器的x轴坐标构成的M维矢量为

设置在平面阵列中的所有传声器的y轴坐标构成的M维矢量为

将上述在平面阵列中的所有传声器的坐标构成的M维矢量代入公式(1)和公式(2)中，得到由声源i到所有传声器之间的传递函数构成的列向量为

构造出所述平面阵列中的所有传声器测得声源i的声压信号p(x,y)的M维矢量来表示真实声源产生的声压信号

1-6)由步骤1-5)可得出所述平面阵列中的所有传声器测得声源i的声压信号p(x,y)的M维矢量来表示真实声源产生的声压信号p^★，用下列公式(3)表示：

式中，q_i为声源i在坐标原点产生的声压复幅值，L为声源的总个数，

为声源i到所有传声器之间的传递函数；

1-7)考虑测量噪声干扰，通过添加一定SNR的高斯白噪声来模拟真实情况：

所述平面阵列中的所有传声器实际接收到的声压信号p为真实声源产生的声压信号p^★与噪声信号n之和，用公式(4)如下所示：

p＝p^★+n (4)

式中，

为噪声信号，且满足信噪比SNR＝20lg(||p^★||₂/||n||₂)；

1-8)将声源i的坐标中θ_i的变化方向定义为θ轴，

的变化方向定义为

轴，建立以θ轴和

轴为坐标轴的二维坐标系，将声源i坐标的可行域沿θ轴和

轴两个方向离散为N个网格点，假设实际声源的波达方向均与网格点重合，建立压缩感知问题模型，用于求解源强分布矢量q，从而将声源识别问题转变为求解如下的欠定线性方程组，得到声源i的声源坐标

及其强度，用公式(5)如下所示：

p＝Aq (5)

式中，

为所有固定网格点到所有传声器之间的传递函数构成的感知矩阵，

为源强分布矢量，

为复数集；

1-9)所述源强分布矢量q的估计值是通过求解如下的l₁范数最小化问题获得，采用公式(6)表示：

式中，ξ为对噪声信号n的容差，通常取ξ＝||n||₂，p为在平面阵列中的所有传声器实际接收到的声源信号，

为所有固定网格点到所有传声器之间的传递函数构成的感知矩阵，

为源强分布矢量，

为复数集；

2)构造替代函数：

2-1)所述二维坐标系中的坐标

为动态网格坐标，建立关于所述动态网格坐标的函数，用于表示采用动态网格坐标构造的感知矩阵

由于采用的是动态网格坐标，感知矩阵不再是常矩阵，所述感知矩阵

用下列公式(7)所示：

式中，

为实数集；

2-2)将公式(7)代入步骤1-8)中的公式(6)中，得到公式(8)如下所示：

式中，ξ为对噪声信号n的容差，通常取ξ＝||n||₂，p为在平面阵列中的所有传声器实际接收到的声源信号，

为所有动态变化的网格点到所有传声器之间的传递函数构成的感知矩阵，

为源强分布矢量，

为复数集，

为实数集；

2-3)在公式(8)的基础上引入对数求和罚函数来代替l₁范数最小化，用于促进解的稀疏性，对公式(8)进行优化后，将公式(8)改写为公式(9)所示：

式中，ξ为对噪声信号n的容差，通常取ξ＝||n||₂，p为在平面阵列中的所有传声器实际接收到的声源信号，

为所有动态变化的网格点到所有传声器之间的传递函数构成的感知矩阵，

为源强分布矢量，

为复数集，

为实数集，∈为正参数，作为迭代进程的控制参数，且用于保证对数函数被正确定义；

2-4)将公式(9)中的∈由初值随着迭代过程逐渐缩减至0，使公式(9)的全局最优解收敛至真实解附近，增加正则化惩罚项

得到下列公式：

式中，p为在平面阵列中的所有传声器实际接收到的声源信号，

为采用动态的网格坐标和非常矩阵的感知矩阵构造的关于声源i的声源坐标

的函数，

为源强分布矢量，

为复数集，

为实数集，λ为正则化参数，其大小会影响重建精度和解的稀疏度；

2-5)用

和λ^(t)分别表示第t次迭代后获得的源强分布矢量和λ，在第t+1次，构造一个关于

的替代函数

用公式(11)如下所示：

式中，

表示定义，且满足

(当且仅当

时等号成立)；

2-6)将公式(11)代入公式(10)，改写为替代函数形式，用公式(12)表示:

式中，p为在平面阵列中的所有传声器实际接收到的声源信号，

为所有动态变化的网格点到所有传声器之间的传递函数构成的感知矩阵，

为源强分布矢量，

为复数集，

为实数集，λ为正则化参数，其大小会影响重建精度和解的稀疏度；

2-8)优化最小化框架通过公式(12)所示的替代函数来间接实现对公式(9)所示的目标函数的优化；

3)求解替代函数：

3-1)为求解公式(12)，首先通过移除常数项将公式(12)简化，如公式(13)所示：

式中，p为在平面阵列中的所有传声器实际接收到的声源信号，

为所有动态变化的网格点到所有传声器之间的传递函数构成的感知矩阵，

为源强分布矢量，

为复数集，

为实数集，λ为正则化参数，其大小会影响重建精度和解的稀疏度；

3-2)D^(t)的表达式如公式(14)所示：

式中，

表示第t次迭代后获得的源强分布矢量，∈为正参数。

3-3)固定θ和

后，通过拉格朗日乘子法可以得到公式(13)关于源强分布矢量q的优化解

如公式(15)所示：

式中，p为在平面阵列中的所有传声器实际接收到的声源信号，

为所有动态变化的网格点到所有传声器之间的传递函数构成的感知矩阵，λ^(t)表示第t次迭代后的正则化参数。

将公式(15)代入公式(13)，消除一个未知量，进一步将公式(13)简化为公式(16)：

式中，p为在平面阵列中的所有传声器实际接收到的声源信号，

为所有动态变化的网格点到所有传声器之间的传递函数构成的感知矩阵，λ^(t)表示第t次迭代后的正则化参数。

3-4)采用梯度下降法求解公式(16)，从而获得新一轮的网格坐标

和

每个网格点的更新如公式(17)所示：

式中，α为初始步长，所述初始步长α的范围为10²～10⁴；

3-5)步骤3-4)的公式(17)中的

的计算公式如公式(18)、公式(19)所示：

3-6)采用梯度下降法，对公式(17)中的初始步长α进行修正，以确保更新后的

和

满足公式(20)：一般按照0.1倍的变化规律对公式(17)中的初始步长α进行修正。

4)通过步骤2)与步骤3)反复循环，构造和迭代缩减替代函数，使网格坐标θ和

逐渐收敛至真实声源所在的位置，源强分布矢量q也在这个过程中得到愈加准确的估计：

4-1)将

和

代入公式(15)即可得到对源强分布矢量q的新一轮估计：

式中，

和

为第t+1次迭代后的网格坐标；

4-2)设置λ的初始值λ⁰，根据公式(20)对下一轮迭代中λ的取值进行更新，用于在迭代过程中保证重构精度和稀疏度之间的平衡：

式中，δ为缩放常数，所述缩放常数δ取值范围为0.05～1，M为传声器的数量，p为在平面阵列中的所有传声器实际接收到的声源信号，

为第t+1次迭代后的网格点到所有传声器之间的传递函数构成的感知矩阵，

为第t+1次迭代后的源强分布矢量；

4-3)设σ为收敛误差，定义

为第t+1次迭代的收敛误差，设置收敛阈值和允许最大迭代次数；

4-4)当收敛误差σ＜收敛阈值或迭代次数≥允许最大迭代次数时，停止构造和迭代缩减替代函数，收敛误差σ收敛到0附近，网格坐标θ和

收敛至真实声源所在的位置，源强分布矢量q得到准确的估计。

还可以在步骤4)中，设置预设阈值κ，删除源强分布矢量q中小于预设阈值κ和其对应的网格点坐标，用于实时动态减少感知矩阵

的维度。

本发明的优点在于，能够克服传统压缩波束形成存在的基不匹配问题，对不在网格点上的声源也能够准确定位，其性能显著优于传统压缩波束形成，满足了实际应用中获取超高精度声源识别成像的需求。

附图说明

图1为本发明的流程图；

图2为本发明的平面阵列波束形成的坐标系示意图；

图3为声源识别仿真流程示意图；

图4为3000Hz时不同声源布置下的仿真识别成像图；

图5为3000Hz时扬声器声源的试验识别成像图。

具体实施方式

如图1至图5所示，一种二维动态网格压缩波束形成的声源识别方法，包括以下步骤：

1)建立压缩感知问题模型：

1-1)以由若干传声器组成的平面阵列的中心为原点建立三维坐标系，用于声源识别试验；

1-2)在三维坐标系中，设置任意声源i的入射方向与z轴的夹角为θ_i，声源i的入射方向在x-y平面上的投影与x轴的夹角为

且0°≤θ_i≤90°，

得到声源i的声源坐标

用于刻画声源i的波达方向；

1-3)在三维坐标系中，所述平面阵列中的所有传声器均在x-y平面上，设置任一传声器的坐标为(x,y)，则声源i到坐标为(x,y)的任一传声器之间的传递函数G(x,y)用下列公式(1)表示：

式中，

为复数单位，

λ为声源i的声波波长，e为自然常数，π为圆周率，

为复数集；

1-4)设置声源i在三维坐标系的原点产生的声压复幅值为q_i，用于表示声源i的源强，所述平面阵列中的任一传声器测得声源i的声压信号p(x,y)用公式(2)表示：

式中，q_i为声源i在坐标原点产生的声压复幅值，L为声源的总个数，G(x,y)为声源i到任一坐标为(x,y)的传声器之间的传递函数；

1-5)设置在平面阵列中的所有传声器的x轴坐标构成的M维矢量为

设置在平面阵列中的所有传声器的y轴坐标构成的M维矢量为

将上述在平面阵列中的所有传声器的坐标构成的M维矢量代入公式(1)和公式(2)中，得到由声源i到所有传声器之间的传递函数构成的列向量为

构造出所述平面阵列中的所有传声器测得声源i的声压信号p(x,y)的M维矢量来表示真实声源产生的声压信号

1-6)由步骤1-5)可得出所述平面阵列中的所有传声器测得声源i的声压信号p(x,y)的M维矢量来表示真实声源产生的声压信号p^★，用下列公式(3)表示：

式中，q_i为声源i在坐标原点产生的声压复幅值，L为声源的总个数，

为声源i到所有传声器之间的传递函数；

1-7)考虑测量噪声干扰，通过添加一定SNR的高斯白噪声来模拟真实情况：

所述平面阵列中的所有传声器实际接收到的声压信号p为真实声源产生的声压信号p^★与噪声信号n之和，用公式(4)如下所示：

p＝p^★+n (4)

式中，

为噪声信号，且满足信噪比SNR＝20lg(||p^★||₂/||n||₂)；

1-8)将声源i的坐标中θ_i的变化方向定义为θ轴，

的变化方向定义为

轴，建立以θ轴和

轴为坐标轴的二维坐标系，将声源i坐标的可行域沿θ轴和

轴两个方向离散为N个网格点，假设实际声源的波达方向均与网格点重合，建立压缩感知问题模型，用于求解源强分布矢量q，从而将声源识别问题转变为求解如下的欠定线性方程组，得到声源i的声源坐标

及其强度，用公式(5)如下所示：

p＝Aq (5)

式中，

为所有固定网格点到所有传声器之间的传递函数构成的感知矩阵，

为源强分布矢量，

为复数集；

1-9)所述源强分布矢量q的估计值是通过求解如下的l₁范数最小化问题获得，采用公式(6)表示：

式中，ξ为对噪声信号n的容差，通常取ξ＝||n||₂，p为在平面阵列中的所有传声器实际接收到的声源信号，

为所有固定网格点到所有传声器之间的传递函数构成的感知矩阵，

为源强分布矢量，

为复数集；

2)构造替代函数：

2-1)所述二维坐标系中的坐标

为动态网格坐标，建立关于所述动态网格坐标的函数，用于表示采用动态网格坐标构造的感知矩阵

由于采用的是动态网格坐标，感知矩阵不再是常矩阵，所述感知矩阵

用下列公式(7)所示：

式中，

为实数集；

2-2)将公式(7)代入步骤1-8)中的公式(6)中，得到公式(8)如下所示：

式中，ξ为对噪声信号n的容差，通常取ξ＝||n||₂，p为在平面阵列中的所有传声器实际接收到的声源信号，

为所有动态变化的网格点到所有传声器之间的传递函数构成的感知矩阵，

为源强分布矢量，

为复数集，

为实数集；

2-3)在公式(8)的基础上引入对数求和罚函数来代替l₁范数最小化，用于促进解的稀疏性，对公式(8)进行优化后，将公式(8)改写为公式(9)所示：

式中，ξ为对噪声信号n的容差，通常取ξ＝||n||₂，p为在平面阵列中的所有传声器实际接收到的声源信号，

为所有动态变化的网格点到所有传声器之间的传递函数构成的感知矩阵，

为源强分布矢量，

为复数集，

为实数集，∈为正参数，作为迭代进程的控制参数，且用于保证对数函数被正确定义；

2-4)将公式(9)中的∈由初值随着迭代过程逐渐缩减至0，使公式(9)的全局最优解收敛至真实解附近，增加正则化惩罚项

得到下列公式：

式中，p为在平面阵列中的所有传声器实际接收到的声源信号，

为采用动态的网格坐标和非常矩阵的感知矩阵构造的关于声源i的声源坐标

的函数，

为源强分布矢量，

为复数集，

为实数集，λ为正则化参数，其大小会影响重建精度和解的稀疏度；

2-5)用

和λ^(t)分别表示第t次迭代后获得的源强分布矢量和λ，在第t+1次，构造一个关于

的替代函数

用公式(11)如下所示：

式中，

表示定义，且满足

(当且仅当

时等号成立)；

2-6)将公式(11)代入公式(10)，改写为替代函数形式，用公式(12)表示:

式中，p为在平面阵列中的所有传声器实际接收到的声源信号，

为所有动态变化的网格点到所有传声器之间的传递函数构成的感知矩阵，

为源强分布矢量，

为复数集，

为实数集，λ为正则化参数，其大小会影响重建精度和解的稀疏度；

2-8)优化最小化框架通过公式(12)所示的替代函数来间接实现对公式(9)所示的目标函数的优化；

3)求解替代函数：

3-1)为求解公式(12)，首先通过移除常数项将公式(12)简化，如公式(13)所示：

式中，p为在平面阵列中的所有传声器实际接收到的声源信号，

为所有动态变化的网格点到所有传声器之间的传递函数构成的感知矩阵，

为源强分布矢量，

为复数集，

为实数集，λ为正则化参数，其大小会影响重建精度和解的稀疏度；

3-2)D^(t)的表达式如公式(14)所示：

式中，

表示第t次迭代后获得的源强分布矢量，∈为正参数。

3-3)固定θ和

后，通过拉格朗日乘子法可以得到公式(13)关于源强分布矢量q的优化解

如公式(15)所示：

式中，p为在平面阵列中的所有传声器实际接收到的声源信号，

为所有动态变化的网格点到所有传声器之间的传递函数构成的感知矩阵，λ^(t)表示第t次迭代后的正则化参数。

将公式(15)代入公式(13)，消除一个未知量，进一步将公式(13)简化为公式(16)：

式中，p为在平面阵列中的所有传声器实际接收到的声源信号，

为所有动态变化的网格点到所有传声器之间的传递函数构成的感知矩阵，λ^(t)表示第t次迭代后的正则化参数。

3-4)采用梯度下降法求解公式(16)，从而获得新一轮的网格坐标

和

每个网格点的更新如公式(17)所示：

式中，α为初始步长，所述初始步长α的范围为10²～10⁴，本实施例中，所述初始步长α＝10³。

3-5)步骤3-4)的公式(17)中的

的计算公式如公式(18)、公式(19)所示：

3-6)采用梯度下降法，本实施例中，按照0.1倍的变化规律对公式(17)中的初始步长α进行修正，以确保更新后的

和

满足公式(20)：

4)通过步骤2)与步骤3)反复循环，构造和迭代缩减替代函数，使网格坐标θ和

逐渐收敛至真实声源所在的位置，源强分布矢量q也在这个过程中得到愈加准确的估计：

4-1)将

和

代入公式(15)即可得到对源强分布矢量q的新一轮估计：

式中，

和

为第t+1次迭代后的网格坐标；

4-2)设置λ的初始值λ⁰，本实施例中，λ⁰＝10³，根据公式(20)对下一轮迭代中λ的取值进行更新，用于在迭代过程中保证重构精度和稀疏度之间的平衡：

式中，δ为缩放常数，所述缩放常数δ取值范围为0.05～1，M为传声器的数量，p为在平面阵列中的所有传声器实际接收到的声源信号，

为第t+1次迭代后的网格点到所有传声器之间的传递函数构成的感知矩阵，

为第t+1次迭代后的源强分布矢量；

4-3)设σ为收敛误差，定义

为第t+1次迭代的收敛误差，设置收敛阈值和允许最大迭代次数，本实施例中，设置收敛阈值为10^-5，允许最大迭代次数为200；

4-4)当收敛误差σ＜收敛阈值或迭代次数≥允许最大迭代次数时，本实施例中，即当收敛误差σ＜10^-5或迭代次数≥200时，停止构造和迭代缩减替代函数，收敛误差σ收敛到0附近，网格坐标θ和

收敛至真实声源所在的位置，源强分布矢量q得到准确的估计。

还可以在步骤4)中，设置预设阈值κ，删除源强分布矢量q中小于预设阈值κ和其对应的网格点坐标，用于实时动态减少感知矩阵

的维度。

为验证建立本方法的准确性、对比探究本方法在性能上的提升，本实施例中，基于MATLAB软件编程如图3所示的流程进行由正向声场模拟和反向声源识别两个部分组成的声源识别仿真模拟实验，其中，正向声场模拟是用程序模拟了包括了器具搭建、假设理论声源以及传声器测量的步骤，从而得到了理论里的传声器信号；而反向声源识别是通过建立压缩感知问题模型，构造替代函数并不断迭代求解替代函数，得到最终收敛的解来对区域进行成像。

在求解步骤中采用l₁范数最小化和OMP两种传统压缩波束形成方法与本发明进行比较：三种算法对不同位置声源进行识别的成像结果如图4所示，l₁范数最小化、OMP、动态网格压缩波束形成分别对应图4中的(a)和(d)成像图、(b)和(e)成像图、(c)和(f)成像图。

仿真采用阵列直径为0.65m的36通道非规则扇形轮阵列识别阵列前方半空间的观测区域并建立如图2所示的坐标系，模拟声场中存在4个不等强度同频率声源的情况对其进行声源识别，源强从左到右依次被设置为2Pa、

1Pa、2Pa，声源频率均为3000Hz，在测量得到的声压信号中添加SNR为30dB的高斯白噪声。

传统压缩波束形成θ和

两个方向以5°的间隔进行网格划分，动态网格压缩波束形成θ和

两个方向以18°的间隔进行初始网格划分。

设置图4中的(a)、(b)、(c)成像图中的声源坐标为(65°,110°)、(45°,160°)、(30°,210°)、(45°,245°)对应其中传统压缩波束形成基匹配的情况，设置图4中的(d)、(e)、(f)成像图中的声源坐标为(63°,107.5°)、(44°,162°)、(26.5°,213°)、(47°,243°)对应其中传统压缩波束形成基不匹配的情况。

计算结果参照最大输出值进行dB缩放，显示动态范围设置为20dB，同时每幅图上方标出以基准声压2×10^-5Pa为参考的最大输出值。

在动态网格压缩波束形成的仿真中预设κ＝0.5Pa。由图4中的(a)、(d)成像图可见，当声源未落在网格点上，即发生了基不匹配时，l₁范数最小化的识别结果中声源能量被分散到了若干相邻的网格点上，不仅无法准确定位真实声源位置，并且对源强的量化出现了较大的偏差；同样对比图4中的(b)、(e)成像图可见，OMP对声源的定位也出现了偏差，只能识别到真实声源附近的网格点，不能准确识别声源位置。而由图4中的(c)、(f)成像图可见，动态网格压缩波束形成在基不匹配情况下同样能够准确识别真实声源所在位置，综上所述，本发明的性能相较于传统压缩波束形成具有显著优势。

为检验仿真结论的正确性，在半消声室内进行试验验证，将稳态信号激励的扬声器作为声源，采用

公司、半径为0.65m、集成4958型传声器的36通道非规则扇形轮阵列采样声压信号。

传统压缩波束形成θ和

两个方向以5°的间隔进行网格划分，动态网格压缩波束形成θ和

两个方向以18°的间隔进行初始网格划分。

2个扬声器和2个由于地面反射形成的镜像声源之间的位置关系分别用坐标(1.25,0,3m),(1.25,-2.4,3m),(-1.25,0,3m),(-1.25,-2.4,3m)表示，对应的DOA分别为(22.62°,0°),(42.05°,297.51°),(22.62°,180°),(42.05°,242.49°)。

各传声器接收的声压信号经PULSE 3560D型数据采集系统同时采集并传输到PULSE LABSHOP中进行频谱分析，采样频率为16384Hz，频率分辨率为4Hz。进一步，采用MATLAB编制的动态网格压缩波束形成程序和传统压缩波束形成程序分别对频谱成分进行识别成像。

图5为3000Hz时扬声器声源的识别成像图，l₁范数最小化、OMP、动态网格压缩波束形成分别对应图5中的(a)、(b)、(c)成像图，与图4呈现出一致的规律，从而证明了本发明方法切实可行。

以上所述仅为本发明的优选实施例，并不用于限制本发明，本领域的技术人员在不脱离本发明的精神的前提提下，对本发明进行的改动均落入本发明的保护范围。

Claims (5)

1.一种二维动态网格压缩波束形成的声源识别方法，其特征在于，包括以下步骤：

1)建立压缩感知问题模型：

1-1)以由若干传声器组成的平面阵列的中心为原点建立三维坐标系，用于声源识别试验；

1-2)在三维坐标系中，设置任意声源i的入射方向与z轴的夹角为θ_i，声源i的入射方向在x-y平面上的投影与x轴的夹角为

且0°≤θ_i≤90°，

得到声源i的声源坐标

用于刻画声源i的波达方向；

1-3)在三维坐标系中，所述平面阵列中的所有传声器均在x-y平面上，设置任一传声器的坐标为(x,y)，则声源i到坐标为(x,y)的任一传声器之间的传递函数G(x,y)用下列公式(1)表示：

式中，

为复数单位，

λ为声源i的声波波长，e为自然常数，π为圆周率，

为复数集；

1-4)设置声源i在三维坐标系的原点产生的声压复幅值为q_i，用于表示声源i的源强，所述平面阵列中的任一传声器测得声源i的声压信号p(x,y)用公式(2)表示：

式中，q_i为声源i在坐标原点产生的声压复幅值，L为声源的总个数，G(x,y)为声源i到任一坐标为(x,y)的传声器之间的传递函数；

1-5)设置在平面阵列中的所有传声器的x轴坐标构成的M维矢量为

设置在平面阵列中的所有传声器的y轴坐标构成的M维矢量为

将上述在平面阵列中的所有传声器的坐标构成的M维矢量代入公式(1)和公式(2)中，得到由声源i到所有传声器之间的传递函数构成的列向量为

构造出所述平面阵列中的所有传声器测得声源i的声压信号p(x,y)的M维矢量来表示真实声源产生的声压信号

1-6)由步骤1-5)可得出所述平面阵列中的所有传声器测得声源i的声压信号p(x,y)的M维矢量来表示真实声源产生的声压信号p^★，用下列公式(3)表示：

式中，q_i为声源i在坐标原点产生的声压复幅值，L为声源的总个数，

为声源i到所有传声器之间的传递函数；

1-7)所述平面阵列中的所有传声器实际接收到的声压信号p为真实声源产生的声压信号p^★与噪声信号n之和，用公式(4)如下所示：

p＝p^★+n (4)

式中，

为噪声信号，且满足信噪比SNR＝20lg(||p^★||₂/||n||₂)；

1-8)将声源i的坐标中θ_i的变化方向定义为θ轴，

的变化方向定义为

轴，建立以θ轴和

轴为坐标轴的二维坐标系，将声源i坐标的可行域沿θ轴和

轴两个方向离散为N个网格点，假设实际声源的波达方向均与网格点重合，建立压缩感知问题模型，用于求解源强分布矢量q，得到声源i的声源坐标

及其强度，用公式(5)如下所示：

p＝Aq (5)

式中，

为所有固定网格点到所有传声器之间的传递函数构成的感知矩阵，

为源强分布矢量，

为复数集；

1-9)所述源强分布矢量q的估计值是通过求解如下的l₁范数最小化问题获得，采用公式(6)表示：

式中，ξ为对噪声信号n的容差，通常取ξ＝||n||₂，p为在平面阵列中的所有传声器实际接收到的声源信号，

为所有固定网格点到所有传声器之间的传递函数构成的感知矩阵，

为源强分布矢量，

为复数集；

2)构造替代函数：

2-1)所述二维坐标系中的坐标

为动态网格坐标，建立关于所述动态网格坐标的函数，用于表示采用动态网格坐标构造的感知矩阵

所述感知矩阵

用下列公式(7)所示：

式中，

为实数集；

2-2)将公式(7)代入步骤1-8)中的公式(6)中，得到公式(8)如下所示：

式中，ξ为对噪声信号n的容差，通常取ξ＝||n||₂，p为在平面阵列中的所有传声器实际接收到的声源信号，

为所有动态变化的网格点到所有传声器之间的传递函数构成的感知矩阵，

为源强分布矢量，

为复数集，

为实数集；

2-3)在公式(8)的基础上引入对数求和罚函数来代替l₁范数最小化，用于促进解的稀疏性，对公式(8)进行优化后，将公式(8)改写为公式(9)所示：

式中，ξ为对噪声信号n的容差，通常取ξ＝||n||₂，p为在平面阵列中的所有传声器实际接收到的声源信号，

为所有动态变化的网格点到所有传声器之间的传递函数构成的感知矩阵，

为源强分布矢量，

为复数集，

为实数集，∈为正参数，作为迭代进程的控制参数，且用于保证对数函数被正确定义；

2-4)将公式(9)中的∈由初值随着迭代过程逐渐缩减至0，使公式(9)的全局最优解收敛至真实解附近，增加正则化惩罚项

得到下列公式：

式中，p为在平面阵列中的所有传声器实际接收到的声源信号，

为采用动态的网格坐标和非常矩阵的感知矩阵构造的关于声源i的声源坐标

的函数，

为源强分布矢量，

为复数集，

为实数集，λ为正则化参数，其大小会影响重建精度和解的稀疏度；

2-5)用

和λ^(t)分别表示第t次迭代后获得的源强分布矢量和λ，在第t+1次，构造一个关于

的替代函数

用公式(11)如下所示：

式中，

表示定义，且满足

(当且仅当

时等号成立)；

2-6)将公式(11)代入公式(10)，改写为替代函数形式，用公式(12)表示:

式中，p为在平面阵列中的所有传声器实际接收到的声源信号，

为所有动态变化的网格点到所有传声器之间的传递函数构成的感知矩阵，

为源强分布矢量，

为复数集，

为实数集，λ为正则化参数，其大小会影响重建精度和解的稀疏度；

2-8)优化最小化框架通过公式(12)所示的替代函数来间接实现对公式(9)所示的目标函数的优化；

3)求解替代函数：

3-1)为求解公式(12)，首先通过移除常数项将公式(12)简化，如公式(13)所示：

式中，p为在平面阵列中的所有传声器实际接收到的声源信号，

为所有动态变化的网格点到所有传声器之间的传递函数构成的感知矩阵，

为源强分布矢量，

为复数集，

为实数集，λ为正则化参数，其大小会影响重建精度和解的稀疏度；

3-2)D^(t)的表达式如公式(14)所示：

式中，

表示第t次迭代后获得的源强分布矢量，∈为正参数。

3-3)固定θ和

后，通过拉格朗日乘子法可以得到公式(13)关于源强分布矢量q的优化解

如公式(15)所示：

式中，p为在平面阵列中的所有传声器实际接收到的声源信号，

为所有动态变化的网格点到所有传声器之间的传递函数构成的感知矩阵，λ^(t)表示第t次迭代后的正则化参数。

将公式(15)代入公式(13)，消除一个未知量，进一步将公式(13)简化为公式(16)：

式中，p为在平面阵列中的所有传声器实际接收到的声源信号，

为所有动态变化的网格点到所有传声器之间的传递函数构成的感知矩阵，λ^(t)表示第t次迭代后的正则化参数。

3-4)采用梯度下降法求解公式(16)，从而获得新一轮的网格坐标

和

每个网格点的更新如公式(17)所示：

式中，α为初始步长；

3-5)步骤3-4)的公式(17)中的

的计算公式如公式(18)、公式(19)所示：

3-6)采用梯度下降法，对公式(17)中的初始步长α进行修正，以确保更新后的

和

满足公式(20)：

4)通过步骤2)与步骤3)反复循环，构造和迭代缩减替代函数，使网格坐标θ和

逐渐收敛至真实声源所在的位置，源强分布矢量q也在这个过程中得到愈加准确的估计：

4-1)将

和

代入公式(15)即可得到对源强分布矢量q的新一轮估计：

式中，

和

为第t+1次迭代后的网格坐标；

4-2)设置λ的初始值，根据公式(20)对下一轮迭代中λ的取值进行更新，用于在迭代过程中保证重构精度和稀疏度之间的平衡：

式中，δ为缩放常数，M为传声器的数量，p为在平面阵列中的所有传声器实际接收到的声源信号，

为第t+1次迭代后的网格点到所有传声器之间的传递函数构成的感知矩阵，

为第t+1次迭代后的源强分布矢量；

4-3)设σ为收敛误差，定义

为第t+1次迭代的收敛误差，设置收敛阈值和允许最大迭代次数；

4-4)当收敛误差σ＜收敛阈值或迭代次数≥允许最大迭代次数时，停止构造和迭代缩减替代函数，网格坐标θ和

收敛至真实声源所在的位置，源强分布矢量q得到准确的估计。

2.根据权利要求1所述的方法，其特征在于：步骤4)中，设置预设阈值κ，删除源强分布矢量q中小于预设阈值κ和其对应的网格点坐标，用于实时动态减少感知矩阵

的维度。

3.根据权利要求1所述的方法，其特征在于：步骤3-4)中，所述初始步长α的范围为10²～10⁴。

4.根据权利要求1所述的方法，其特征在于：步骤3-6)中，按照0.1倍的变化规律对公式(17)中的初始步长α进行修正。

5.根据权利要求1所述的方法，其特征在于：步骤4-2)中，所述缩放常数δ取值范围为0.05～1。

Images