CN112147571B - 基于正则正交匹配追踪和蝙蝠算法的声源方位角估计方法 - Google Patents

Abstract

本发明属于声源定位领域，具体涉及一种基于正则正交匹配追踪和蝙蝠算法的声源方位角估计方法。包括如下步骤：步骤(1)：用软件模拟阵列采集的声音信号；步骤(2)：声音信号预处理；步骤(3)：建立基于压缩感知理论的声源方位角估计稀疏表达模型；步骤(4)：应用正则正交匹配追踪算法与蝙蝠算求解稀疏表达模型，获得具有最大声源强度的网格节点的位置，该位置即代表目标声源相对于传声器阵列的方位角。本发明结合了正则正交匹配追踪和蝙蝠算法，提高了正则正交匹配追踪算法的原子寻优过程的求解效率，与传统方位角估计算法相比可以用更少的声音信号数据得到更优的精度，与现有的基于压缩感知原理的方位角估计算法相比提升了求解效率。

Description

基于正则正交匹配追踪和蝙蝠算法的声源方位角估计方法

技术领域

本发明属于声源定位领域，具体涉及一种基于正则正交匹配追踪和蝙蝠算法的声源方位角估计方法。

背景技术

通常声音信号属于宽带信号，目前有四种方法解决宽带信号的声源方位角估计问题：基于子空间、基于最大似然、基于时延估计这三种传统理论的方法和基于信号稀疏分解理论的算法。

基于子空间的算法运算量较小，但是当信噪比较低时宽带信号的时频聚集性不好，特别是对于多个声源，传统时频分析方法的交叉项对算法结果有较大的影响。

基于最大似然的算法有较高的分辨性能，但是算法复杂度高、运算量大，而且不易于实现。

基于时延估计的算法原理简单，易于实现，但是这种算法的精度依赖于时延差的测量精度，而且这种算法只适用于近场的声源，声源的距离越远，算法的误差越大。

信号的稀疏分解重构理论，即压缩感知理论，在诸多研究领域都有应用，基于该理论的声源方向角估计技术不受奈奎斯特采样定理的限制，它与传统技术相比所需的声音信号数据比较少，因此数据处理效率有了很大的提高，而且它的估计结果准确度也高于传统方法。

目前的基于压缩感知理论的方位角估计算法求解方式中，有一种正交匹配追踪算法，该算法的方位角分辨性能还有待提升，而且该算法的计算量较大，不能实时给出；还有一种正则正交匹配追踪算法，与正交匹配追踪算法相比，它的估计结果更接近于全局最优解，所以准确度更高，但是它要花费更长的迭代计算时间。

发明内容

本发明的目的在于提供一种基于正则正交匹配追踪与蝙蝠算法的声源方位角估计方法，解决了在声源方位角估计过程中声音信号可用数据过少、求解精度与求解效率不能兼得的问题。

实现本发明目的的技术解决方案为：一种基于正则正交匹配追踪和蝙蝠算法的声源方位角估计方法，包括如下步骤：

步骤(1)：用软件模拟阵列采集的声音信号；

步骤(2)：声音信号预处理；

步骤(3)：建立基于压缩感知理论的声源方位角估计稀疏表达模型；

步骤(4)：应用正则正交匹配追踪算法与蝙蝠算求解稀疏表达模型，获得具有最大声源强度的网格节点的位置，该位置即代表目标声源相对于传声器阵列的方位角。

进一步的，所述步骤(1)“用软件模拟阵列采集的声音信号”具体为：

在MATLAB模型中添加圆形传声器阵列的阵列流型矩阵，用函数生成多声源声音信号，该信号与阵列流型矩阵结合可模拟阵列采集到的声音信号，获得阵列采集到的声音信号在时域的时域测量数据。

进一步的，所述步骤(2)“声音信号预处理”具体为：

对步骤(1)采集的时域测量数据进行离散傅里叶变换处理，获得声音信号在频域的测量数据；选择所需的频段和频率点个数，获取所述频段内各个频率点的测量数据。

进一步的，所述步骤(3)“建立基于压缩感知理论的声源方位角估计稀疏表达模型”包括对二维声源平面的网格划分和测量矩阵的建立，具体步骤如下：

步骤(3-1)：假设声源点与传声器阵列处在同一平面，对圆形传声器阵列所对应的圆进行等角度的网格划分，以每个网格节点为可能的声源方位角，共划分了3600×1＝3600个网格，每两个相邻的网格节点所代表的方位角相差0.1°；

步骤(3-2)：假设未知的声源信号z，z为一3600×1的列向量，z的元素的值代表了声源在对应的网格节点处的强度；

步骤(3-3)：基于自由场格林函数的Helmoltz方程，建立对应频率点f_j的测量矩阵C(f_j)；

步骤(3-4)：建立频段内各频率点下的声源方位角稀疏表达模型：

y(f_j)＝C(f_j)z(f_j)+e(f_j)，j＝1,2,…,W

其中y(f_j)为传声器阵列的测量值，z(f_j)为频率点f_j下的声源原始信号，e(f_j)为噪声信号，W为频率点个数；

步骤(3-5)：对声源方位角稀疏表达模型中的y(f_j)进行SVD分解：

y(f_j)_SV＝y(f_j)VD_K，z(f_j)_SV＝z(f_j)VD_K，e(f_j)_SV＝e(f_j)VD_K

其中D_K＝[I_K 0]^T，V为对测量值y(f_j)奇异值分解后的酉矩阵，I_K为K×K阶单位矩阵，K为声源信号的稀疏度，假设K＝3；

步骤(3-6)：将稀疏表达模型再次变形为其中/>和/>由y_SV,z_SV,e_SV的每一列首尾相连得到；由N个矩阵C组成对角线得到新矩阵/>N与稀疏度K相关，当K为3时N为9。

进一步的，所述步骤(4)“应用正则正交匹配追踪算法与蝙蝠算求解稀疏表达模型，获得具有最大声源强度的网格节点的位置，该位置即代表目标声源相对于传声器阵列的方位角”具体包括如下步骤：

步骤(4-1)：给迭代索引集合Λ_t、索引对应列向量组成的矩阵A_t以及重构声音信号的残差r_t赋初始值；

步骤(4-2)：用蝙蝠算法求解传感矩阵A的列向量与残差初始值的内积u的H个较大值或所有非零值，H为K的正整数倍，倍数自拟，这些值对应的列序号j组成集合J；K为信号的稀疏度，在数值上等于声源的个数；

步骤(4-3)：从集合J中得到子集合J₀，集合J₀均满足：|u(i)|≤2|u(j)|,i,j∈J₀选择所有子集J₀中具有最大能量的K个J₀；

步骤(4-4)：求Λ_t和A_t，Λ_t＝Λ_t-1∪J₀,A_t＝A_t-1∪a_j，a_j为J₀内的所有列向量；

步骤(4-5)：求y＝A_tθ_t的最小二乘解：y为观测向量；/>为信号稀疏表示系数的估计值；

步骤(4-6)：求残差r_t的更新值，

步骤(4-7)：如果迭代次数t小于信号稀疏度K，则返回步骤(4-2)；如果t大于K或||Λ_t||₀≥2K或残差r_t＝0则停止迭代进入下一步；

步骤(4-8)：重构所得的值为最后一次迭代所得/>用稀疏矩阵Ψ可得重构信号 中K个较大的值对应的方位角即为假设的K个声源的方位角估计值。

进一步的，所述步骤(4-2)“用蝙蝠算法求解传感矩阵A的列向量与残差初始值的内积u的K个较大值或所有非零值”具体包括如下步骤：

步骤(4-2-1)：设定初始参数，设定脉冲频率范围[f_min,f_max]，设定最大迭代次数M，设定初始响度C_k ^t，设定终止阙值ε_t，在360°的空域内随机产生N个蝙蝠的初始位置x_k(k＝1，2，…N)，N的值根据传感矩阵A的列向量个数与最大迭代次数M自拟，x_k与传感矩阵A中的列向量所对应的角度θ_k一一对应；

步骤(4-2-2)：根据蝙蝠个体的脉冲频率f和响度C_k ^t找出当前最佳的H个蝙蝠的位置，记为x_k ^*，(k＝1,2,...K),H为K的正整数倍，倍数自拟，K为信号的稀疏度；

步骤(4-2-3)：种群迭代，蝙蝠个体的位置、频率和速度按照下式进行迭代更新，

式中rand为一随机数，x^*为当前最佳的H个位置中的最优蝙蝠位置；

步骤(4-2-4)：求局部新解，如果rand<r_k，则从最佳解的集合中选一个，其中r_k为发射脉冲的频度，为预设值，在该最佳解附近生成一个局部解如下式

x_new＝x_old+ε×C^t

其中的ε是[-1，1]中的任意值，C^t是所有蝙蝠在此次迭代中的平均响度；

步骤(4-2-5)：再次生成rand，评价得到的局部新解，如果rand<C_k，则这个解符合要求，增大r_k，减小C_k，开始下一次迭代，响度C_k和脉冲数率r_k的更新如下式

其中，α，γ是预设常数，0<α<1，γ＞0；

步骤(4-2-6)：排列并找到蝙蝠的H个最佳位置，对应可能的方位角估计值；

步骤(4-2-7)：若达到最大迭代次数M，保存此时的x_k ^*，(k＝1,2,...H)；

步骤(4-2-8)：输出此时的x_k ^*。

本发明与现有技术相比，其显著优点在于：

本发明结合了正则正交匹配追踪算法和蝙蝠算法，提高了正则正交匹配追踪算法的原子寻优过程的求解效率，该算法与传统算法相比可以用更少的声音信号数据得到更高的求解精度，与正交匹配追踪算法还有正则正交匹配追踪算法相比提升了求解效率。

附图说明

图1为本发明的基于蝙蝠算法的正则正交匹配追踪算法流程图。

图2为本发明的蝙蝠算法流程图。

图3为本发明的圆形传声器阵列示意图。

具体实施方式

下面结合附图对本发明作进一步详细描述。

下面详细描述本发明的实施例，所述实施例是示例性的，旨在用于解释本发明，而不能理解为对本发明的限制。声源定位过程具体如下：

步骤1：在MATLAB模型中添加圆形传声器阵列的阵列流型矩阵，用函数生成多声源声音信号，该信号与阵列流型矩阵结合可模拟阵列采集到的声音信号，采样率设置为44.1kHz/16bt，采样时长设置为5s，获得阵列采集到的声音信号在时域的测量数据。

步骤2：对步骤1采集的时域测量数据进行离散傅里叶变换处理，获得声音信号在频域的测量数据。

步骤3：选择所需的频段和频率点个数，本示例中设置频段为300-1000Hz，获取这段频段内各个频率点的测量数据。

步骤4：建立基于压缩感知理论的声源方位角估计稀疏表达模型，包括对二维声源平面的网格划分和测量矩阵的建立，具体包括以下步骤：

步骤4.1；假设声源点与传声器阵列处在同一平面，对圆形传声器阵列所对应的圆进行等角度的网格划分，以每个网格节点为可能的声源方位角，共划分了3600×1＝3600个网格，每两个相邻的网格节点所代表的方位角相差0.1°。

步骤4.2：假设未知的声源信号z，z为一3600×1的列向量，z的元素的值代表了声源在对应的网格节点处的强度。

步骤4.3：基于自由场格林函数的Helmoltz方程，建立对应频率点f_j的测量矩阵C(f_j)。

步骤4.4：建立频段内各频率点下的声源方位角稀疏表达模型：

y(f_j)＝C(f_j)z(f_j)+e(f_j)，j＝1,2,…,W

其中y(f_j)为传声器阵列的测量值，z(f_j)为频率点f_j下的声源原始信号，e(f_j)为噪声信号，W为频率点个数。

步骤4.5：对声源方位角稀疏表达模型中的y(f_j)进行SVD分解：

y(f_j)_SV＝y(f_j)VD_K，z(f_j)_SV＝z(f_j)VD_K，e(f_j)_SV＝e(f_j)VD_K

其中D_K＝[I_K0]^T，V为对测量值y(f_j)奇异值分解后的酉矩阵，I_K为K×K阶单位矩阵，K为声源信号的稀疏度，本实例中设K＝3。

步骤4.6：将稀疏表达模型再次变形为其中/>和/>由y_SV,z_SV,e_SV的每一列首尾相连得到；由N个矩阵C组成对角线得到新矩阵/>N与稀疏度K相关，本实例中值为9.

步骤5：应用基于蝙蝠算法改进的正则正交匹配追踪算法对变形后的稀疏表达模型进行求解，获得不同频率点下每个网格节点处的声源强度，具体过程为：

(1)对重构声音信号的残差r_t、迭代索引集合Λ_t、索引对应列向量组成的矩阵A_t分别进行初始化；

(2)用蝙蝠算法求解传感矩阵A的列向量与残差初始值的所有内积u中的H个较大值或所有非零值，H为K的正整数倍，倍数自拟，K为信号的稀疏度，其步骤如下：

(a)设定初始参数，设定脉冲频率范围[f_min,f_max]，设定最大迭代次数M，设定初始响度C_k ^t，设定终止阙值ε_t，在360°的空域内随机产生N个蝙蝠的初始位置x_k(k＝1，2，…N)，N的值根据传感矩阵A的列向量个数与最大迭代次数M自拟，x_k与传感矩阵A中的列向量所对应的角度θ_k一一对应；

(b)根据蝙蝠个体的脉冲频度和响度找出当前最佳的H个蝙蝠的位置，记为x_k ^*，(k＝1,2...H)；

(c)种群迭代，蝙蝠个体的位置、频率和速度按照下式

进行迭代更新，rand为一随机数，x^*为当前最佳的H个位置中的最优蝙蝠位置；

(d)求局部新解，如果rand<r_k(r_k为发射脉冲的频度，为预设值)，则从最佳解的集合中选一个，在最佳解附近生成一个局部解如下式

x_new＝x_old+ε×C^t

其中的ε是[-1，1]中的任意值，C^t是所有蝙蝠在此次迭代中的平均响度；

(e)再次生成rand,评价得到的解，如果rand<C_k，则这个解符合要求，增大r_k，减小C_k，开始下一次迭代，响度C_k和脉冲数率r_k的更新如下式

其中，α，γ是预设常数，0<α<1，γ＞0；

(f)排列并找到蝙蝠的H个最佳位置，对应可能的方位角估计值；

(g)若达到最大迭代次数M，保存此时的x_k ^*，(k＝1,2...H)；

(h)输出此时的x_k ^*，对应方位角估计的所有可能解。将这些值对应的列序号j构成集合J；

(3)从集合J中得到集合J₀，满足：|u(i)|≤2|u(j)|,for all i,j∈J₀，选择所有满足要求的子集中具有最大能量的K个J₀；

(4)令Λ_t＝Λ_t-1∪J₀,A_t＝A_t-1∪a_j，a_j为J₀内的所有列向量；

(5)求y＝A_tθ_t的最小二乘解：y为观测向量，为信号稀疏表示系数估计；

(6)更新残差

(7)如果迭代次数t小于信号稀疏度K，则返回步骤二，如果t大于K或||Λ_t||₀≥2K，或残差r_t＝0则停止迭代进入步骤八；

(8)重构所得的值为最后一次迭代所得/>利用稀疏矩阵Ψ可得重构信号本实例中设置的声源个数为3，则重构信号/>的元素中最大的3个数值对应该系统输出的声源方位角的估计值。

Claims (5)

1.一种基于正则正交匹配追踪和蝙蝠算法的声源方位角估计方法，其特征在于，包括如下步骤：

步骤(1)：用软件模拟阵列采集的声音信号；

步骤(2)：声音信号预处理；

步骤(3)：建立基于压缩感知理论的声源方位角估计稀疏表达模型；

步骤(4)：应用正则正交匹配追踪算法与蝙蝠算求解稀疏表达模型，获得具有最大声源强度的网格节点的位置，该位置即代表目标声源相对于传声器阵列的方位角：

步骤(4-1)：给迭代索引集合Λ_t、索引对应列向量组成的矩阵A_t以及重构声音信号的残差r_t赋初始值；

步骤(4-2)：用蝙蝠算法求解传感矩阵A的列向量与残差初始值的内积u的H个较大值或所有非零值，H为K的正整数倍，倍数自拟，这些值对应的列序号j组成集合J；K为信号的稀疏度，在数值上等于声源的个数；

步骤(4-3)：从集合J中得到子集合J₀，集合J₀均满足：|u(i)|≤2|u(j)|,i,j∈J₀选择所有子集J₀中具有最大能量的K个J₀；

步骤(4-4)：求Λ_t和A_t，Λ_t＝Λ_t-1∪J₀,A_t＝A_t-1∪a_j，a_j为J₀内的所有列向量；

步骤(4-5)：求y＝A_tθ_t的最小二乘解：y为观测向量；/>为信号稀疏表示系数的估计值；

步骤(4-6)：求残差r_t的更新值，

步骤(4-7)：如果迭代次数t小于信号稀疏度K，则返回步骤(4-2)；如果t大于K或||Λ_t||₀≥2K或残差r_t＝0则停止迭代进入下一步；

步骤(4-8)：重构所得的值为最后一次迭代所得/>用稀疏矩阵Ψ可得重构信号/> 中K个较大的值对应的方位角即为假设的K个声源的方位角估计值。

2.根据权利要求1所述的方法，其特征在于，所述步骤(1)“用软件模拟阵列采集的声音信号”具体为：

在MATLAB模型中添加圆形传声器阵列的阵列流型矩阵，用函数生成多声源声音信号，该信号与阵列流型矩阵结合可模拟阵列采集到的声音信号，获得阵列采集到的声音信号在时域的时域测量数据。

3.根据权利要求2所述的方法，其特征在于，所述步骤(2)“声音信号预处理”具体为：

对步骤(1)采集的时域测量数据进行离散傅里叶变换处理，获得声音信号在频域的测量数据；选择所需的频段和频率点个数，获取所述频段内各个频率点的测量数据。

4.根据权利要求3所述的方法，其特征在于，所述步骤(3)“建立基于压缩感知理论的声源方位角估计稀疏表达模型”包括对二维声源平面的网格划分和测量矩阵的建立，具体步骤如下：

步骤(3-1)：假设声源点与传声器阵列处在同一平面，对圆形传声器阵列所对应的圆进行等角度的网格划分，以每个网格节点为可能的声源方位角，共划分了3600×1＝3600个网格，每两个相邻的网格节点所代表的方位角相差0.1°；

步骤(3-2)：假设未知的声源信号z，z为3600×1的列向量，z的元素的值代表了声源在对应的网格节点处的强度；

步骤(3-3)：基于自由场格林函数的Helmoltz方程，建立对应频率点f_j的测量矩阵C(f_j)；

步骤(3-4)：建立频段内各频率点下的声源方位角稀疏表达模型：

y(f_j)＝C(f_j)z(f_j)+e(f_j)，j＝1,2,…,W

其中y(f_j)为传声器阵列的测量值，z(f_j)为频率点f_j下的声源原始信号，e(f_j)为噪声信号，W为频率点个数；

步骤(3-5)：对声源方位角稀疏表达模型中的y(f_j)进行SVD分解：

y(f_j)_SV＝y(f_j)VD_K，z(f_j)_SV＝z(f_j)VD_K，e(f_j)_SV＝e(f_j)VD_K

其中D_K＝[I_K 0]^T，V为对测量值y(f_j)奇异值分解后的酉矩阵，I_K为K×K阶单位矩阵，K为声源信号的稀疏度，假设K＝3；

步骤(3-6)：将稀疏表达模型再次变形为其中/>和/>由y_SV,z_SV,e_SV的每一列首尾相连得到；由N个矩阵C组成对角线得到新矩阵/>N与稀疏度K相关，当K为3时N为9。

5.根据权利要求4所述的方法，其特征在于，所述步骤(4-2)“用蝙蝠算法求解传感矩阵A的列向量与残差初始值的内积u的H个较大值或所有非零值，H为K的正整数倍”具体包括如下步骤：

步骤(4-2-1)：设定初始参数，设定脉冲频率范围[f_min,f_max]，设定最大迭代次数M，设定初始响度C_k ^t，设定终止阈值ε_t，在360°的空域内随机产生N个蝙蝠的初始位置x_k，k＝1，2，…N，N的值根据传感矩阵A的列向量个数与最大迭代次数M自拟，x_k与传感矩阵A中的列向量所对应的角度θ_k一一对应；

步骤(4-2-2)：根据蝙蝠个体的脉冲频率f和响度C_k ^t找出当前最佳的H个蝙蝠的位置，记为x_k ^*，k＝1,2,…H，H为K的正整数倍，倍数自拟，K为信号的稀疏度；

步骤(4-2-3)：种群迭代，蝙蝠个体的位置、频率和速度按照下式进行迭代更新，

式中rand为一随机数，x^*为当前最佳的H个位置中的最优蝙蝠位置；

步骤(4-2-4)：求局部新解，如果rand<r_k，则从最佳解的集合中选一个，其中r_k为发射脉冲的频度，为预设值，在该最佳解附近生成一个局部解如下式

x_new＝x_old+ε×C^t

其中的ε是[-1，1]中的任意值，C^t是所有蝙蝠在此次迭代中的平均响度；

步骤(4-2-5)：再次生成rand，评价得到的局部新解，如果rand<C_k，则这个解符合要求，增大r_k，减小C_k，开始下一次迭代，响度C_k和脉冲数率r_k的更新如下式

其中，α，γ是预设常数，0<α<1，γ>0；

步骤(4-2-6)：排列并找到蝙蝠的H个最佳位置，对应可能的方位角估计值；

步骤(4-2-7)：若达到最大迭代次数M，保存此时的x_k ^*，k＝1,2,...H；

步骤(4-2-8)：输出此时的x_k ^*。