CN106443587A - 一种高分辨率的快速反卷积声源成像算法 - Google Patents

Abstract

本发明公开了一种高分辨率的快速反卷积声源成像算法，其特征是首先在反卷积声源成像算法的点扩展函数矩阵构造过程中，利用点扩展函数的近似空间平移不变性，计算声源计算平面中心位置处声源的点扩展函数，然后通过对中心位置处的点扩展函数向上和向下循环移位的方法构造出点扩展函数矩阵，避免了对全部点扩展函数的计算，降低计算量，提高了计算速度和效率；其次在声源源强能量分布的反卷积重构过程中，利用声源的空间稀疏先验，结合压缩感知理论，通过正交匹配追踪算法来实现声源源强能量分布的快速稀疏反卷积重构，降低了迭代次数，提高计算效率和分辨率。本发明具有高计算效率和空间分辨率，能更快更好地识别与定位声源在空间中的位置。

Description

一种高分辨率的快速反卷积声源成像算法

技术领域

本发明主要运用于噪声源的定位与分析领域中，是一种高分辨率的快速反卷积声源成像算法。

背景技术

传统波束形成技术基于延时求和、互谱成像函数等方法，对各传感器接收信号进行后处理，使声源计算平面上对应真实声源的聚焦点的输出量被加强，而其他聚焦点的输出量被衰减，从而识别声源。但其不仅会在真实声源位置输出具有一定宽度的“主瓣”，还会在非声源位置输出“旁瓣”，为了有效消除主瓣宽度和旁瓣干扰的影响，发展出了反卷积声源成像方法，即DAMAS，其空间分辨率明显优于延时求和、互谱成像等传统波束形成技术。但是其由于需要构建完整的点扩展函数矩阵，而且需要反复进行迭代来解方程组，所以有着计算量大，耗时长的局限，计算效率很低。为了提高其计算效率，近年来，发展出了多种基于快速傅里叶变换的反卷积声源成像方法，即FFT-DAMAS，以下两种最具代表性：第一种是DAMAS的扩展方法DAMAS2，其假设点扩展函数只取决于观测点与声源点间的相对位置，而与具体位置无关，具有近似空间平移不变性，利用快速傅里叶变换将声源源强能量分布与点扩展函数间的卷积转化为波数域的乘积，然后利用循环迭代来进行求解；第二种是基于快速傅里叶变换的非负最小二乘反卷积声源成像方法FFT-NNLS，其基本思想是在互谱成像波束形成输出结果、点扩展函数、声源分布之间建立差函数，然后利用快速傅里叶变换和点扩展函数的空间平移不变性来极小化差函数。这两种方法与最初的反卷积声源成像方法相比，具有收敛快、计算效率高、声源定位准确等优点，所以被广泛使用。

但DAMAS2方法和FFT-NNLS方法在计算效率和分辨率方面还有进一步提升的空间。首先DAMAS2方法和FFT-NNLS方法均为迭代算法，为了得到比较精确的结果，通常需要进行上千次的循环迭代，并且由于傅里叶变换引入了卷绕误差，为减小卷绕误差对计算结果的影响，通常需要将测量数据点数补零到原来的四倍，导致计算规模急剧增大，所以尽管DAMAS2方法和FFT-NNLS方法在计算效率上优于最初的反卷积声源成像方法，但是在计算效率和分辨率上还是有所不足。

发明内容

本发明是为避免上述现有技术所存在的不足之处，提出一种高分辨率的快速反卷积声源成像算法，利用点扩展函数矩阵近似空间平移不变性，提高计算速度和效率，通过正交匹配追踪算法来实现反卷积过程，提高计算结果的分辨率，降低所需的迭代次数，避免DAMAS2方法和FFT-NNLS方法中傅里叶变换带来的卷绕误差，从而获得更高的计算效率、更高的计算精度和更高的分辨率。

本发明解决技术问题采用的技术方案是：

本发明高分辨率的快速反卷积声源成像算法的特点是按如下步骤进行：

步骤a、在K个声源辐射形成的声场中呈阵列布置M个传感器，形成测量面W，采集获得各个传感器声压数据，声源的个数K小于传感器的数量M，所述传感器为传声器或质点振速传感器；

步骤b、将声源计算平面离散成一网格面，所述网格面为聚焦面T，在所述聚焦面T中包含有N个网格点，每个网格点为聚焦点，由式(1)计算出第n个聚焦点的互谱成像波束形成输出量b(r_n)：

其中，为导向矢量，i为虚数单位，k为声波波数，k＝2πf/c，π为圆周率，f为声源频率，c为声速，为聚焦面T上第n个的聚焦点到测量面W上所有传感器的距离组成的列向量，H为共轭转置；

U为测量声压数据的互谱矩阵， 表示测量面W上所有传感器坐标组成的列向量，表示所有传感器所接受到的声压组成的列向量，为抑制传感器阵列自噪声的干扰，将互谱矩阵U的主对角线元素全部置零，即对互谱矩阵U去自谱；

利用式(1)对每个聚焦点进行计算，获得所有聚焦点的互谱成像波束形成输出量组成的列向量 表示聚焦面T上所有聚焦点坐标组成的列向量；

步骤c、由式(2)计算出聚焦面T的中心点上的声源源强能量与聚焦面T上第n个聚焦点的互谱成像波束形成输出量b(r_n)之间的点扩展函数psf(r_n/r_z)，其中r_z表示聚焦面T的中心点的坐标，并且当聚焦面T上聚焦点个数N为奇数时，z＝(N+1)/2，为偶数时，z＝N/2：

式(2)中，为传递函数组成的列向量， 为聚焦面T中心点到所有传感器的距离组成的列向量；

利用式(2)对每个聚焦点进行计算，得到聚焦面T中心点上的声源源强能量与聚焦面T上所有聚焦点的互谱成像波束形成输出量之间的点扩展函数组成的列向量

步骤d、通过如下循环移位的方式构建所有聚焦点上的声源源强能量与所有聚焦点的互谱成像波束形成输出量之间的点扩展函数组成的矩阵 表示聚焦面T上所有聚焦点坐标组成的列向量：

第n₁个聚焦点上的声源源强能量与聚焦面T上所有聚焦点的互谱成像波束形成输出量之间的点扩展函数组成的列向量为点扩展函数矩阵中的第n₁列；

当n₁＜z时，是通过将中的每个元素向上循环平移(z-n₁)位得到；

当n₁＞z时，是通过将中的每个元素向下循环平移(n₁-z)位得到；

通过上述方式构建出点扩展函数矩阵中除第z列以外的其他N-1列，并且已知点扩展函数矩阵中的第z列从而得出完整的点扩展函数矩阵

步骤e、利用点扩展函数矩阵和互谱成像波束形成输出量建立式(3)所示的矩阵方程：

式(3)中，为所有聚焦点上的声源源强能量组成的列向量；

步骤f、运用正交匹配追踪算法对式(3)进行求解，从而获得声源源强能量分布实现声源定位。

本发明高分辨率的快速反卷积声源成像算法的特点也在于：所述步骤f中运用正交匹配追踪算法对式(3)进行求解是按如下步骤进行：

第一步：初始化残差r₀＝b，残差下降率r₀′＝0，索引集Λ₀＝φ，选取的点扩展函数矩阵列向量组成的矩阵Φ₀＝φ，令初始迭代次数为j＝1；

第二步：通过式(4)找到中的列向量与残差向量内积最大值索引位置λ_j：

第三步：按如下方式更新索引集Λ_j和Φ_j：

将第二步得到λ_j放入第j-1次迭代得到的索引集Λ_j-1中，再将中的第λ_j列放入第j-1次迭代得到的矩阵Φ_j-1中，并且将的第λ_j列全部赋值为零；

第四步：根据最小二乘法求解出声源源强能量x_j为：x_j＝arg min||Φ_jx-b||²；

第五步：根据求解的x_j值计算新的残差值r_j为：r_j＝b-Φ_jx_j；则残差下降率r_j′为：r_j′＝||r_j-1||²-||r_j||²；其中，r_j-1为第j-1次迭代得到的残差；

第六步：若r′_j-1≤r_j′，则将j增加1作为新的迭代次数，并且返回到第二步继续迭代；否则停止迭代，进入第七步，并且取出迭代过程中所得到的索引集Λ_j和x_j；

第七步：根据第六步得到的索引集Λ_j和x_j获得声源源强能量分布列向量是将声源源强能量分布列向量在索引位置Λ_j处的项全部赋值为x_j中所对应的值，剩余项全部为零。

与已有技术相比，本发明有益效果体现在：

1、本发明方法由于使用了正交匹配追踪算法来解卷积，而且避免了傅里叶变换的卷绕误差，所以在整个适宜于反卷积算法的频段内，其分辨率与计算精度均明显优于DAMAS2方法和FFT-NNLS方法。

2、本发明方法中利用点扩展函数的近似空间平移不变性，仅计算声源计算平面中心位置处声源的点扩展函数，通过循环移位方法构造点扩展函数矩阵，而且通过正交匹配追踪算法来解卷积其迭代次数少，因此本发明与DAMAS2方法和FFT-NNLS方法相比，具有有着更高的计算效率。

3、本发明方法具备优良的鲁棒性和抗干扰能力，在低信噪比条件下声源的定位精度明显优于DAMAS2方法和FFT-NNLS方法，性能稳定。

附图说明

图1为本发明声源识别布局示意图；

图2a、图2b和图2c分别是在声源频率为1600Hz时，DAMAS2方法、FFT-NNLS方法以及本发明方法的声源定位效果；

图3a、图3b和图3c分别是在声源频率为2000Hz时，DAMAS2方法、FFT-NNLS方法以及本发明方法的声源定位效果；

图4a、图4b和图4c分别是在声源频率为2400Hz时，DAMAS2方法、FFT-NNLS方法以及本发明方法的声源定位效果；

图5为声源频率为2400Hz时，本发明方法、DAMAS2方法和FFT-NNLS方法的计算时间随聚焦点数目的变化曲线；

图6a、图6b和图6c分别是在信噪比0dB的低信噪比条件下，声源频率为1600Hz时，DAMAS2方法、FFT-NNLS方法以及本发明方法的声源定位效果；

图7a、图7b和图7c分别是在信噪比0dB的低信噪比条件下，声源频率为2000Hz时，DAMAS2方法、FFT-NNLS方法以及本发明方法的声源定位效果；

图8a、图8b和图8c分别是在信噪比0dB的低信噪比条件下，声源频率为2400Hz时，DAMAS2方法、FFT-NNLS方法以及本发明方法的声源定位效果；

具体实施方式

本发明高分辨率的快速反卷积声源成像算法是按如下步骤进行：

步骤a、如图1所示，在K个声源辐射形成的声场中呈阵列布置M个传感器，形成测量面W，采集获得各个传感器声压数据，声源的个数K小于传感器的数量M，传感器为传声器或质点振速传感器。

步骤b、将声源计算平面离散成一网格面，网格面为聚焦面T，在聚焦面T中包含有N个网格点，每个网格点为聚焦点，由式(1)计算出第n个聚焦点的互谱成像波束形成输出量b(r_n)：

式(1)中，为导向矢量，i为虚数单位，k为声波波数，k＝2πf/c，π为圆周率，f为声源频率，c为声速，为聚焦面T上第n个的聚焦点到测量面W上所有传感器的距离组成的列向量，H为共轭转置。

U为测量声压数据的互谱矩阵， 表示测量面W上所有传感器坐标组成的列向量，表示所有传感器所接受到的声压组成的列向量，为抑制传感器阵列自噪声的干扰，将互谱矩阵U的主对角线元素全部置零，即对互谱矩阵U去自谱。

利用式(1)对每个聚焦点进行计算，获得所有聚焦点的互谱成像波束形成输出量组成的列向量 表示聚焦面T上所有聚焦点坐标组成的列向量。

步骤c、由式(2)计算出聚焦面T的中心点上的声源源强能量与聚焦面T上第n个聚焦点的互谱成像波束形成输出量b(r_n)之间的点扩展函数psf(r_n/r_z)，其中r_z表示聚焦面T的中心点的坐标，并且当聚焦面T上聚焦点个数N为奇数时，z＝(N+1)/2，为偶数时，z＝N/2：

式(2)中，为传递函数组成的列向量， 为聚焦面T中心点到所有传感器的距离组成的列向量。

利用式(2)对每个聚焦点进行计算，得到聚焦面T中心点上的声源源强能量与聚焦面T上所有聚焦点的互谱成像波束形成输出量之间的点扩展函数组成的列向量

步骤d、通过如下循环移位的方式构建所有聚焦点上的声源源强能量与所有聚焦点的互谱成像波束形成输出量之间的点扩展函数组成的矩阵 表示聚焦面T上所有聚焦点坐标组成的列向量：

第n₁个聚焦点上的声源源强能量与聚焦面T上所有聚焦点的互谱成像波束形成输出量之间的点扩展函数组成的列向量为点扩展函数矩阵中的第n₁列；

当n₁＜z时，是通过将中的每个元素向上循环平移(z-n₁)位得到；

当n₁＞z时，是通过将中的每个元素向下循环平移(n₁-z)位得到；

通过上述方式构建出点扩展函数矩阵中除第z列以外的其他N-1列，并且已知点扩展函数矩阵中的第z列从而得出完整的点扩展函数矩阵

点扩展函数矩阵中，n和n₁都表示任意一个聚焦点，它们取值没有关联，它们的范围都是1到N，而且n对应点扩展函数矩阵的行，n₁对应点扩展函数矩阵的列，例如第n₁个聚焦点上的声源源强能量与第n个聚焦点的互谱成像波束形成输出量b(r_n)之间的点扩展函数为点扩展函数矩阵第n行第n₁列的数。

本发明方法运用循环移位的方式来构建点扩展函数矩阵相比于通过公式直接算出点扩展函数矩阵其大大缩短了计算所用的时间，提高了本方法的计算效率。

步骤e、利用点扩展函数矩阵和互谱成像波束形成输出量建立式(3)所示的矩阵方程：

式(3)中，为所有聚焦点上的声源源强能量组成的列向量。

步骤f、运用正交匹配追踪算法对式(3)进行求解，从而获得声源源强能量分布实现声源定位。

本发明方法运用正交匹配追踪算法来进行声源源强能量的重构，避免了由于傅里叶变换引入的卷绕误差，使获得的声源源强能量更加准确，从而提高计算结果的分辨率和精度。

具体实施中，运用正交匹配追踪算法对式(3)进行求解是按如下步骤进行：

第一步：初始化残差r₀＝b，残差下降率r₀′＝0，索引集Λ₀＝φ，选取的点扩展函数矩阵列向量组成的矩阵Φ₀＝φ，令初始迭代次数为j＝1。

第二步：通过式(4)找到中的列向量与残差向量内积最大值索引位置λ_j：

第三步：按如下方式更新索引集Λ_j和Φ_j：

将第二步得到λ_j放入第j-1次迭代得到的索引集Λ_j-1中，再将中的第λ_j列放入第j-1次迭代得到的矩阵Φ_j-1中，并且将的第λ列全部赋值为零。

第四步：根据最小二乘法求解出声源源强能量x_j为：x_j＝arg min||Φ_jx-b||²。

第五步：根据求解的x_j值计算新的残差值r_j为：r_j＝b-Φ_jx_j；则残差下降率r_j′为：r_j′＝||r_j-1||²-||r_j||²；其中，r_j-1为第j-1次迭代得到的残差。

第六步：若r′_j-1≤r_j′，则将j增加1作为新的迭代次数，并且返回到第二步继续迭代；否则停止迭代，进入第七步，并且取出迭代过程中所得到的索引集Λ_j和x_j。

第七步：根据第六步得到的索引集Λ_j和x_j获得声源源强能量分布列向量是将声源源强能量分布列向量在索引位置Λ_j处的项全部赋值为x_j中所对应的值，剩余项全部为零。

本发明方法与DAMAS2方法和FFT-NNLS方法相比，具有高分辨率、高的计算效率和稳定性强等特点，能更快、更好地识别与定位声源在声源计算平面上的位置。

理论模型

本发明方法在声源定位的过程中仍采用反卷积声源成像的原理，首先找到互谱成像波束形成输出量与声源源强能量之间的传递关系，然后通过解方程求出声源源强能量，来实现声源的定位，其过程如下：

图1中测量面W位于xoy平面，黑色圆点表示阵列传感器，设总共有M个传感器，而且传感器个数要大于实际声源个数，传感器为传声器或质点振速传感器；将声源计算平面离散成一网格面，网格面为聚焦面T，在聚焦面T中包含有N个网格点，每个网格点为聚焦点。

测量面W上所有传感器接收的声压为传递矩阵和所有聚焦点的声源强度的乘积：

式(2-1)中，表示测量面W上所有传感器坐标组成的列向量，表示聚焦面T上所有聚焦点坐标组成的列向量，表示所有传感器所接受到的声压组成的列向量，为传递函数组成的矩阵，表示所有聚焦点上的声源强度组成的列向量。

式(2-1)中：

式(2-2)中，|r_m,n|表示第m个传感器和第n个聚焦点之间的距离，i表示虚数单位，k表示声波波数，k＝2πf/c，π为圆周率，f为声源频率，c为声速。

第n个聚焦点阵列响应输出为：

式(2-3)中，为导向矢量， 表示第n个聚焦点到所有传感器的距离组成的列向量，H表示共轭转置。

第n个聚焦点的互谱成像波束形成输出量为：

b(r_n)＝y(r_n)y(r_n)^H (2-4)

将式(2-3)代入式(2-4)得：

记式(2-5)中间两项为U：

式(2-6)中，U表示测量声压数据的互谱矩阵，为抑制传感器阵列自噪声的干扰，将互谱矩阵U的主对角线元素全部置零，即对互谱矩阵U去自谱，可提高声源识别精度，将去自谱后的式(2-6)代入式(2-5)，相应的变为得：

将式(2-1)代入式(2-6)得：

式(2-8)中：

假设N个聚焦点上的声源为不相干声源，则式(2-9)中非主对角线元素相对于主对角元素可以忽略不计，所以式(2-8)可以简化为：

式(2-10)中，为传递函数矩阵的第n列，q(r_n)表示第n个聚焦点上的声源强度，声源强度平方|q(r_n)|²为声源源强能量。

将式(2-10)代入式(2-7)，由于式(2-10)和式(2-7)中n的取值没有关联，为了避免混淆，把式(2-10)中的n变为n₁，得：

定义下式：

式(2-12)中，表示第n₁个聚焦点上的声源源强能量与聚焦面T上第n个聚焦点的互谱成像波束形成输出量b(r_n)之间的点扩展函数。

将式(2-12)代入式(2-11)得：

为简化表示，将式(2-13)写成矩阵形式：

式(2-14)中，为所有聚焦点上的声源源强能量组成的列向量，表示所有聚焦点的互谱成像波束形成输出量组成的列向量，表示聚焦面T上所有聚焦点坐标组成的列向量，表示所有聚焦点上的声源源强能量与聚焦面T上所有聚焦点的互谱成像波束形成输出量之间的点扩展函数组成的矩阵。

运用式(2-12)得到聚焦面T中心点上的声源源强能量与聚焦面T上所有聚焦点的互谱成像波束形成输出量之间的点扩展函数的列向量其中r_z表示聚焦面T的中心点的坐标，并且当聚焦面T上聚焦点个数N为奇数时，z＝(N+1)/2，为偶数时，z＝N/2；然后通过如下循环移位来构建点扩展函数矩阵

第n₁个聚焦点上的声源源强能量与聚焦面T上所有聚焦点的互谱成像波束形成输出量之间的点扩展函数组成的列向量为点扩展函数矩阵中的第n₁列，当n₁＜z时，是通过将中的每个元素向上循环平移(z-n₁)位得到，当n₁＞z时，是通过将中的每个元素向下循环平移(n₁-z)位得到；通过上述方式构建出点扩展函数矩阵中除第z列以外的其他N-1列，并且已知点扩展函数矩阵中的第z列从而得出完整的点扩展函数矩阵

运用式(2-7)得到所有聚焦点的互谱成像波束形成输出量组成的列向量

运用正交匹配追踪算法对式(2-14)进行求解是按如下步骤进行：

第一步：初始化残差r₀＝b，残差下降率r₀′＝0，索引集Λ₀＝φ，选取的点扩展函数矩阵列向量组成的矩阵Φ₀＝φ，令初始迭代次数为j＝1；

第二步：通过式(2-15)找到中的列向量与残差向量内积最大值索引位置λ_j：

第三步：按如下方式更新索引集Λ_j和Φ_j：

将第二步得到λ_j放入第j-1次迭代得到的索引集Λ_j-1中，再将中的第λ_j列放入第j-1次迭代得到的矩阵Φ_j-1中，并且将的第λ_j列全部赋值为零；

第四步：根据最小二乘法求解出声源源强能量x_j为：x_j＝arg min||Φ_jx-b||²；

第五步：根据求解的x_j值计算新的残差值r_j为：r_j＝b-Φ_jx_j；则残差下降率r_j′为：r_j′＝||r_j-1||²-||r_j||²；其中，r_j-1为第j-1次迭代得到的残差；

第六步：若r′_j-1≤r_j′，则将j增加1作为新的迭代次数，并且返回到第二步继续迭代；否则停止迭代，进入第七步，并且取出迭代过程中所得到的索引集Λ_j和x_j；

第七步：根据第六步得到的索引集Λ_j和x_j获得声源源强能量分布列向量是将声源源强能量分布列向量在索引位置Λ_j处的项全部赋值为x_j中所对应的值，剩余项全部为零。

仿真1：

验证本发明方法相对于DAMAS2方法和FFT-NNLS方法具有更高的分辨率与定位精度。

仿真中声源为两个位于Z＝1m的平面T上的具有相等强度的两个不相干的点声源，声源坐标分别为：(-0.15m,0,1m)、(0.1m,0,1m)。

测量面W位于Z＝0m平面内，W的尺寸为1mx1m，W上均匀分布11x11个测量网格点，因此测量面W上的测量点之间的间距为0.1m。

聚焦面T位于Z＝1m平面内，T的尺寸为1mx1m，T上均匀分布41x41个聚焦点，因此聚焦面T上的聚焦点之间的间距为0.025m。

为使仿真与实际实施中存在测量噪声的情况更加一致，测量面W声压添加了高斯白噪声，信噪比为20dB。

分别采用本发明方法、DAMAS2方法和FFT-NNLS方法对聚焦面T上的声源进行定位，并对比本发明方法与DAMAS2方法和FFT-NNLS方法定位的效果，以验证本发明方法的优势。采用DAMAS2方法和FFT-NNLS方法计算时，迭代次数为1000次；本发明算法计算时，根据判断条件自动停止迭代。

图2a、图2b、图2c、图3a、图3b、图3c以及图4a、图4b和图4c分别给出了声源频率为1600Hz、2000Hz、2400Hz时DAMAS2方法、FFT-NNLS方法和本发明方法的声源定位效果，其中，图2a、图3a、图4a为DAMAS2方法的声源定位效果，图2b、图3b、图4b为FFT-NNLS方法的声源定位效果，图2c、图3c、图4c为本发明方法的声源定位效果，图中的两个“+”表示声源实际的位置。

从图2a、图3a和图4a可以看出，DAMAS2方法在频率为1600Hz时定位偏差很大而且旁瓣也很大，频率为2000Hz时，定位仍有偏差而且旁瓣比较大，频率为2400Hz时，虽然定位比较准确但是还是有旁瓣；从图2b、图3b和图4b可以看出，FFT-NNLS方法在三个不同频率下的定位比较准确但旁瓣都很大；从图2c、图3c和图4c可以看出，本发明方法在三个不同频率下的定位非常准确而且旁瓣很小。由图2a、图2b和图2c，图3a、图3b和图3c，图4a、图4b和图4c对比可以看出，在不同频率下，本发明方法相比于DAMAS2方法和FFT-NNLS方法有着更准确的定位和更小的旁瓣。

上述分析表明本发明方法相对于DAMAS2方法和FFT-NNLS方法能够更准确的定位出声源的位置，也就是说本发明方法具有更高的分辨率与定位精度。

仿真2

验证本发明方法相对于DAMAS2方法和FFT-NNLS方法具有更高计算效率。

在频率为2400Hz时，把聚焦面T分别划分为289、729、1369、2209、3249个聚焦点，保持聚焦点间距为0.025m，其余仿真参数与仿真1保持不变，分析本发明方法、DAMAS2方法和FFT-NNLS方法计算所需的时间，从而验证本发明方法在计算效率方面的优越性。

图5表示在同一运行环境下，三种算法所需要的时间随聚焦点个数增加的变化情况。其中a、b和c曲线分别表示FFT-NNLS方法、DAMAS2方法和本发明方法的计算时间随聚焦点数目的变化曲线。从图中可以看出随着聚焦点数目的增加，三种算法计算时间都有所上升，而且可以发现DAMAS2方法和FFT-NNLS方法所需时间的曲线始终处于本发明方法所需时间的曲线的上面，DAMAS2方法和FFT-NNLS方法计算所需的时间始终多于本发明方法计算所需要的时间，并且可以看出，当聚焦面T上的聚焦点数量较多时DAMAS2方法和FFT-NNLS方法计算所需的时间远远多于本发明方法计算所需要的时间，例如：当聚焦点数目为1369时，DAMAS2方法所需时间是本发明方法所需时间的9倍，FFT-NNLS方法所需时间是本发明方法所需时间的26倍之多。

上述分析表明本发明方法相对于DAMAS2方法和FFT-NNLS方法具有更高计算效率。

仿真3

验证本发明方法相比于DAMAS2方法和FFT-NNLS方法在低信噪比条件下的优越性。

仿真条件基本与仿真1中相同，仅将测量面W声压的信噪比降低为0dB。也是得到1600Hz、2000Hz、2400Hz三个不同频率下，本发明方法、DAMAS2方法和FFT-NNLS方法声源定位效果。

图6a、图6b、图6c、图7a、图7b、图7c以及图8a、图8b和图8c分别给出了在信噪比为0dB的情况下声源频率为1600Hz、2000Hz、2400Hz时DAMAS2方法、FFT-NNLS方法和本发明方法的声源定位效果，其中，图6a、图7a、图8a为DAMAS2方法在信噪比为0dB的情况下的声源定位效果，图6b、图7b、图8b为FFT-NNLS方法在信噪比为0dB的情况下的声源定位效果，图6c、图7c、图8c为本发明方法在信噪比为0dB的情况下的声源定位效果，图中的两个“+”表示声源实际的位置。

由图6a、图7a、图8a、图6b、图7b和图8b可以看出，DAMAS2方法和FFT-NNLS方法在信噪比为0dB的情况下定位会有偏差而且旁瓣比较大；由图6c、图7c和图8c可以看出，本发明方法在信噪比为0dB的情况下的定位还是非常准确而且旁瓣也很小。而且本仿真的结果与仿真1的结果相比，可以看出DAMAS2方法和FFT-NNLS方法当信噪比降低时，其定位的效果会变差，但信噪比降低对本发明方法的定位效果基本没有影响。

上述分析表明本发明方法比DAMAS2方法和FFT-NNLS方法在低信噪比的情况下定位更精确，而且本发明方法比DAMAS2方法和FFT-NNLS方法的性能更加稳定。因此本发明方法相比于DAMAS2方法和FFT-NNLS方法在低信噪比条件下具有优越性。

Claims (2)

1.一种高分辨率的快速反卷积声源成像算法，其特征是按如下步骤进行：

步骤a、在K个声源辐射形成的声场中呈阵列布置M个传感器，形成测量面W，采集获得各个传感器声压数据，声源的个数K小于传感器的数量M，所述传感器为传声器或质点振速传感器；

步骤b、将声源计算平面离散成一网格面，所述网格面为聚焦面T，在所述聚焦面T中包含有N个网格点，每个网格点为聚焦点，由式(1)计算出第n个聚焦点的互谱成像波束形成输出量b(r_n)：

$b\left(r_n\right)=\frac{1}{M^2-M}v{\left(r_{\stackrel{\→}{m},n}\right)}^{H}Uv\left(r_{\stackrel{\→}{m},n}\right)---\left(1\right)$

其中，为导向矢量，i为虚数单位，k为声波波数，k＝2πf/c，π为圆周率，f为声源频率，c为声速，为聚焦面T上第n个的聚焦点到测量面W上所有传感器的距离组成的列向量，H为共轭转置；

U为测量声压数据的互谱矩阵， 表示测量面W上所有传感器坐标组成的列向量，表示所有传感器所接受到的声压组成的列向量，为抑制传感器阵列自噪声的干扰，将互谱矩阵U的主对角线元素全部置零，即对互谱矩阵U去自谱；

利用式(1)对每个聚焦点进行计算，获得所有聚焦点的互谱成像波束形成输出量组成的列向量 表示聚焦面T上所有聚焦点坐标组成的列向量；

步骤c、由式(2)计算出聚焦面T的中心点上的声源源强能量与聚焦面T上第n个聚焦点的互谱成像波束形成输出量b(r_n)之间的点扩展函数psf(r_n/r_z)，其中r_z表示聚焦面T的中心点的坐标，并且当聚焦面T上聚焦点个数N为奇数时，z＝(N+1)/2，为偶数时，z＝N/2：

$psf(r_n/r_z)=\frac{1}{M^2-M}v{\left(r_{\stackrel{\→}{m},n}\right)}^{H}g\left(r_{\stackrel{\→}{m},z}\right)g{\left(r_{\stackrel{\→}{m},z}\right)}^{H}v\left(r_{\stackrel{\→}{m},n}\right)---\left(2\right)$

式(2)中，为传递函数组成的列向量， 为聚焦面T中心点到所有传感器的距离组成的列向量；

利用式(2)对每个聚焦点进行计算，得到聚焦面T中心点上的声源源强能量与聚焦面T上所有聚焦点的互谱成像波束形成输出量之间的点扩展函数组成的列向量

步骤d、通过如下循环移位的方式构建所有聚焦点上的声源源强能量与所有聚焦点的互谱成像波束形成输出量之间的点扩展函数组成的矩阵 表示聚焦面T上所有聚焦点坐标组成的列向量：

第n₁个聚焦点上的声源源强能量与聚焦面T上所有聚焦点的互谱成像波束形成输出量之间的点扩展函数组成的列向量为点扩展函数矩阵中的第n₁列；

当n₁＜z时，是通过将中的每个元素向上循环平移(z-n₁)位得到；

当n₁＞z时，是通过将中的每个元素向下循环平移(n₁-z)位得到；

通过上述方式构建出点扩展函数矩阵中除第z列以外的其他N-1列，并且已知点扩展函数矩阵中的第z列从而得出完整的点扩展函数矩阵

步骤e、利用点扩展函数矩阵和互谱成像波束形成输出量建立式(3)所示的矩阵方程：

$b\left(r_{\stackrel{\→}{n}}\right)=psf(r_{\stackrel{\→}{n}}/r_{{\stackrel{\→}{n}}_1})x\left(r_{{\stackrel{\→}{n}}_1}\right)---\left(3\right)$

式(3)中，为所有聚焦点上的声源源强能量组成的列向量；

步骤f、运用正交匹配追踪算法对式(3)进行求解，从而获得声源源强能量分布实现声源定位。

2.根据权利要求1所述的高分辨率的快速反卷积声源成像算法，其特征是：所述步骤f中运用正交匹配追踪算法对式(3)进行求解是按如下步骤进行：

第一步：初始化残差r₀＝b，残差下降率r₀′＝0，索引集Λ₀＝φ，选取的点扩展函数矩阵列向量组成的矩阵Φ₀＝φ，令初始迭代次数为j＝1；

第二步：通过式(4)找到中的列向量与残差向量内积最大值索引位置λ_j：

${\mathrm{\&lambda;}}_j=\arg {\max }_{n_1=1,\mathrm{...},N}|<r_{j-1},psf(r_{\stackrel{\&RightArrow;}{n}}/r_{n_1})>|---\left(4\right)$

第三步：按如下方式更新索引集Λ_j和Φ_j：

将第二步得到λ_j放入第j-1次迭代得到的索引集Λ_j-1中，再将中的第λ_j列放入第j-1次迭代得到的矩阵Φ_j-1中，并且将的第λ_j列全部赋值为零；

第四步：根据最小二乘法求解出声源源强能量x_j为：x_j＝argmin||Φ_jx-b||²；

第五步：根据求解的x_j值计算新的残差值r_j为：r_j＝b-Φ_jx_j；则残差下降率r_j′为：r_j′＝||r_j-1||²-||r_j||²；其中，r_j-1为第j-1次迭代得到的残差；

第六步：若r′_j-1≤r′_j，则将j增加1作为新的迭代次数，并且返回到第二步继续迭代；否则停止迭代，进入第七步，并且取出迭代过程中所得到的索引集Λ_j和x_j；

第七步：根据第六步得到的索引集Λ_j和x_j获得声源源强能量分布列向量是将声源源强能量分布列向量在索引位置Λ_j处的项全部赋值为x_j中所对应的值，剩余项全部为零。