CN113312841B - 具有声源稀疏度自适应性的变范数等效源近场声全息算法 - Google Patents

Abstract

本发明公开了一种具有声源稀疏度自适应性的变范数等效源近场声全息算法，是在声源近场辐射区域内布置全息面并测量获得全息面上声压P_H；在目标重建面远离全息面一侧布置等效源面；利用等效源与全息面之间的声压传递矩阵建立声压与等效源间关系；采用带有自适应范数约束惩罚项的迭代正则化算法求解等效源源强，再用求得的源强以及等效源与目标重建面之间的传递矩阵计算获得目标重建面上的声场数据。本发明自适应地改变正则化求解过程中l_p范数惩罚项中的p值，使求解所得等效源源强的稀疏度与声源的实际稀疏度更加一致，对空间稀疏型声源、空间连续型声源或源强稀疏度介于充分稀疏与连续之间的声源皆能够实现精确求解。

Description

具有声源稀疏度自适应性的变范数等效源近场声全息算法

技术领域

本发明涉及噪声源的识别领域，具体涉及一种具有声源稀疏度自适应性的变范数等效源近场声全息算法。

背景技术

近场声全息技术是一种强大的声源识别和声场可视化工具，全息变换算法是近场声全息的核心部分，经过多年研究，目前已经发展出适用于平面、柱面和球面等可分离坐标的空间傅立叶变换法和统计最优算法，以及适用于任意形状声源的边界元法和等效源法等多种近场声全息算法。其中等效源法由于其可以用于任意形状物体的辐射声场重建和预测，并避免了BEM中特征频率处解的非唯一性问题，克服了计算边界处声学量的奇异积分问题且具有较高的计算精度等，在工程中具有广泛的应用前景。

而基于等效源法的近场声全息在计算精度上仍有一定的进步空间。欲获得精确的重建结果，等效源源强的求解过程尤为关键，该过程在数学上属于具有不适定性的逆问题，故通常采用正则化方法稳定求解过程。在常规的等效源法中，通常采用最小l₂范数约束惩罚项的Tikhonov正则化方法，由于在l₂范数约束准则下往往会导致一个光滑解，从而获得相对连续的等效源强解，故在处理空间连续型声源时常规等效源法有着较好的重建效果，而在处理空间稀疏型声源时的效果不尽如人意。针对基于常规等效源法的近场声全息的缺点，后续有学者将压缩感知理论引入到等效源法中，实现了基于压缩等效源法的近场声全息。该方法在对源强求解时采用l₁范数约束惩罚项，使所得结果总会导致一个稀疏解，故该方法在处理空间稀疏型声源时有较高的重建精度，但在处理空间连续型声源时却无法获得更高重建精度的近场声全息计算结果。然而在实际应用中，不仅有空间连续型声源和空间稀疏型声源，还存在大量稀疏度介于两者之间的声源，由于上述两种方法均采用固定范数，没有根据声源实际稀疏度的不同调整正则化方法中采用的惩罚项范数，故无法很好地匹配该类声源的实际稀疏度以获得精确的重建结果，由此可见上述两种方法均存在求解过程灵活性差且普适性低的局限性。

发明内容

本发明是为避免上述现有技术所存在的不足，提供一种对具有不同稀疏度的各种声源皆能灵活地实现精确求解的具有声源稀疏度自适应性的变范数等效源近场声全息算法。

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

本发明具有声源稀疏度自适应性的变范数等效源近场声全息算法的特点是按如下步骤进行：

步骤a、于声源的近场区域内布置全息面并通过测量获得所述全息面上的测量声压：

在由目标声源产生的声场中，于近场辐射区域内布置全息面H，所述全息面H为任意形状曲面，在所述全息面H上划分测量网格点，通过传感器测量获得各测量网格点处的测量声压P_H；

步骤b、在声源与全息面H之间布置目标重建面T，在目标重建面T远离全息面H的一侧布置等效源面S_e，所述等效源面S_e与目标重建面T之间的垂直距离为d_h，布置在所述等效源面S_e上的等效源的点数不大于全息面H上测量网格点的点数；

步骤c、建立各等效源与全息面H上测量声压P_H之间的关系如式(1)：

P_H＝G_HpQ (1)

G_Hp为等效源面S_e上各等效源与全息面H上各测量网格点之间的声压传递矩阵；

Q为等效源源强向量；

步骤d、以h为步长，将l_p线性划分为m个数值，则有l_p＝{ξ_i，i＝1，2，...，m}，并将m个l_p的数值分别代入

形成求解等效源源强向量Q的l_p范数正则化过程；l_p是数值范围为0到2的实数；

步骤e、按如下过程采用迭代重加权算法将最小l_p问题转化为最小l₂问题进行求解，是指将l_p范数正则化过程转化为使用最小l₂范数作惩罚项的标准Tikhonov正则化过程，求解等效源源强Q：

第一步：将采用常规等效源法求解所得的等效源源强向量作为等效源源强向量初始解Q⁰，令初始迭代次数为1，设定最大迭代次数为t_max，用t表征迭代次数；

第二步：利用向量Q^t-1作为后验加权系数构建加权矩阵W如式(2)：

所述加权矩阵W是利用向量Q^t-1作为主对角元素构建的对角阵，所述向量Q^t-1为第t-1次迭代计算获得的等效源源强向量；对于首次迭代计算过程，向量Q^t-1即为等效源源强向量初始解Q⁰；

第三步：利用所述加权矩阵W的逆矩阵W^-1对向量Q^t加权，并构造加权后向量Q^t的惩罚项为||W^-1Q^t||，形成如式(3)的极小化过程以获得等效源源强向量Q^t的正则解：

min{||P_H-G_HpQ^t||²+λ||W^-1Q^t||} (3)

λ为正则化参数，正则化参数λ采用L曲线法选择获得；

第四步：通过预条件处理令W^-1Q^t≡X^t，G_HpW≡A^t，将式(3)所示的对等效源源强向量Q^t求解的正则化过程转化为如式(4)对中间变量X^t求解的标准Tikhonov正则化过程，并采用奇异值分解法对式(4)求解以获得中间变量X^t的正则解：

min{||P_H-A^tX^t||²+λ||X^t||} (4)

经过第t次迭代获得等效源源强向量Q^t如式(5)：

Q^t≡WX^t (5)

第五步：若t≤t_max则将迭代次数赋值为t+1，并转入第二步继续进行迭代；若t>t_max则终止迭代，并令Q^t为等效源源强向量Q；

步骤f、通过K折交叉验证法选出最匹配当前等效源源强稀疏度的l_p-optimal：

所述K折交叉验证法是指：将测量声压P_H均分为K份轮流作验证集，剩余的K-1份作训练集代入步骤e中用于计算；通过m次K折交叉验证法获得由m个l_p值分别作范数约束惩罚项进行正则化过程求解时所得的平均误差，取其中最小平均误差对应的l_p值为l_p-optimal，形成如式(6)所示的以等效源源强向量Q为求解对象的l_p-optimal范数正则化过程：

步骤g、输入完整的测量声压P_H，重复步骤e，对所述l_p-optimal范数正则化过程求解，令第五步中求得的等效源源强解Q为输出的最终源强解Q_f；

步骤h、利用P_T＝G_TpQ_f计算获得目标重建面T上的声压P_T；利用V_T＝G_TvQ_f计算获得目标重建面T上的法向振速V_T，其中，G_Tp为各等效源与目标重建面之间的声压传递矩阵，G_Tv为各等效源与目标重建面T之间的振速传递矩阵。

本发明具有声源稀疏度自适应性的变范数等效源近场声全息算法的特点也在于：所述d_h的取值为0.03m。

本发明具有声源稀疏度自适应性的变范数等效源近场声全息算法的特点也在于：所述步长h＝0.1，m＝21，l_p＝{0，0.1，0.2，...，2}。

本发明具有声源稀疏度自适应性的变范数等效源近场声全息算法的特点也在于：所述K折交叉验证法中K＝8。

本发明具有声源稀疏度自适应性的变范数等效源近场声全息算法的特点也在于：所述最大迭代次数t_max不超过10次。

本发明具有声源稀疏度自适应性的变范数等效源近场声全息算法的特点也在于：设置所述最大迭代次数t_max为5次。

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

1、本发明方法通过采用带有自适应范数约束惩罚项的迭代正则化算法，使求解过程能够根据声源稀疏度的不同进行自适应调节，故无需以满足一定的源强稀疏度条件为前提即可对任意稀疏度的等效源源强实现精确求解，求解过程灵活并具有普适的应用范围；

2、本发明方法在整个适宜于近场声全息运用的频段内，对于空间连续型声源能够达到与基于常规等效源的近场声全息相同的精度水平；对于空间非连续型声源(即空间稀疏型声源与源强稀疏度介于充分稀疏与连续之间的声源)，相比基于常规等效源法与基于压缩等效源法的近场声全息，本方法皆具有精度优势；

3、本发明方法性能稳定，具备优良的鲁棒性和抗干扰能力，在不同的信噪比条件下的计算精度均优于基于常规等效源法与压缩等效源法的近场声全息。

附图说明

图1为本发明采用的等效源计算模型；

图2为l_p变范数中不同的p值会导致不同稀疏度解的几何说明；

图3a为单极子声源频率为2000Hz时，目标重建面上的理论声压分布；

图3b为单极子声源频率为2000Hz时，本发明方法计算所得的声压分布；

图3c为单极子声源频率为2000Hz时，压缩等效源法计算所得的声压分布；

图3d为单极子声源频率为2000Hz时，常规等效源计算所得的声压分布；

图4为以单极子为目标声源时，不同频率下本发明方法与常规等效源法、压缩等效源法重建误差的比较；

图5a为0.5m×0.5m简支铝板声源频率为2000Hz时，目标重建面上的理论声压分布；

图5b为0.5m×0.5m简支铝板声源频率为2000Hz时，本发明方法计算所得的声压分布；

图5c为0.5m×0.5m简支铝板声源频率为2000Hz时，压缩等效源法计算所得的声压分布；

图5d为0.5m×0.5m简支铝板声源频率为2000Hz时，常规等效源计算所得的声压分布；

图6为以0.5m×0.5m简支铝板为目标声源时，不同频率下本发明方法与常规等效源法、压缩等效源法重建误差的比较；

图7a为0.2m×0.2m简支铝板声源频率为2000Hz时，目标重建面上的理论声压分布；

图7b为0.2m×0.2m简支铝板声源频率为2000Hz时，本发明方法计算所得的声压分布；

图7c为0.2m×0.2m简支铝板声源频率为2000Hz时，压缩等效源法计算所得的声压分布；

图7d为0.2m×0.2m简支铝板声源频率为2000Hz时，常规等效源计算所得的声压分布；

图8为以0.2m×0.2m简支铝板为目标声源时，不同频率下本发明方法与常规等效源法、压缩等效源法重建误差的比较；

图9为以0.2m×0.2m简支铝板为目标声源时，不同信噪比下本发明方法与常规等效源法、压缩等效源法重建误差的比较。

具体实施方式

本实施例中具有声源稀疏度自适应性的变范数等效源近场声全息算法按如下步骤进行：

步骤a、于声源的近场区域内布置全息面并通过测量获得全息面上的测量声压：

在由目标声源产生的声场中，于近场辐射区域内布置全息面H，全息面H为任意形状曲面，在全息面H上划分测量网格点，通过传感器测量获得各测量网格点处的测量声压P_H。

步骤b、在声源与全息面H之间布置目标重建面T，在目标重建面T远离全息面H的一侧布置等效源面S_e，等效源面S_e与目标重建面T之间的垂直距离为d_h，d_h的取值不宜过大或过小，太大将增强声压传递矩阵列向量间的相关性，使矩阵条件数增加从而影响求解精度；太小会产生奇异性，本实施例取值为d_h＝0.03m；布置在等效源面S_e上的等效源的点数不大于全息面H上测量网格点的点数，以保证建立的不是欠定方程组以增加求解稳定性；等效源为单极子、偶极子或四极子，等效源为点源、面源或体源。

步骤c、建立各等效源与全息面H上测量声压P_H之间的关系如式(1)：

P_H＝G_HpQ (1)

G_Hp为等效源面S_e上各等效源与全息面H上各测量网格点之间的声压传递矩阵；

Q为等效源源强向量。

步骤d、以h为步长，将l_p线性划分为m个数值，则有l_p＝{ξ_i,i＝1,2,...,m}，并将m个l_p的数值分别代入

形成求解等效源源强向量Q的l_p范数正则化过程；l_p是数值范围为0到2的实数；步长h的取值不宜过大或过小，过大会导致错失迭代计算中将解收敛至更高精度的l_p值；过小会导致增加无效的计算量，降低求解效率。本实施例中取为h＝0.1，故有m＝21，l_p＝{0,0.1,0.2,...,2}。

步骤e、按如下过程采用迭代重加权算法将最小l_p问题转化为最小l₂问题进行求解，是指将l_p范数正则化过程转化为使用最小l₂范数作惩罚项的标准Tikhonov正则化过程，求解等效源源强Q。

第一步：将采用常规等效源法求解所得的等效源源强向量作为等效源源强向量初始解Q⁰，令初始迭代次数为1，设定最大迭代次数为t_max，用t表征迭代次数。

本算法收敛速度较快，通常只需3-5次迭代计算即可收敛到高精度的等效源源强解，3-5次迭代后若继续增加迭代次数求解精度将基本保持不变，因此本算法在实施过程中通常t_max取3-5次，优选为5次，最大不超过10次；常规等效源法是指采用以等效源源强向量Q的l₂范数||Q||²作惩罚项的标准Tikhonov正则化方法求解等效源源强向量Q，目前已广泛应用于近场声全息计算中的等效源法。

第二步：利用向量Q^t-1作为后验加权系数构建加权矩阵W如式(2)：

加权矩阵W是利用向量Q^t-1作为主对角元素构建的对角阵，向量Q^t-1为第t-1次迭代计算获得的等效源源强向量；对于首次迭代计算过程，向量Q^t-1即为等效源源强向量初始解Q⁰。

第三步：利用加权矩阵W的逆矩阵W^-1对向量Q^t加权，并构造加权后向量Q^t的惩罚项为||W^-1Q^t||，形成如式(3)的极小化过程以获得等效源源强向量Q^t的正则解：

min{||P_H-G_HpQ^t||²+λ||W^-1Q^t||} (3)

λ为正则化参数，正则化参数λ采用L曲线法选择获得。

第四步：通过预条件处理令W^-1Q^t≡X^t，G_HpW≡A^t，将式(3)所示的对等效源源强向量Q^t求解的正则化过程转化为如式(4)对中间变量X^t求解的标准Tikhonov正则化过程，并采用奇异值分解法对式(4)求解以获得中间变量X^t的正则解：

min{||P_H-A^tX^t||²+λ||X^t||} (4)

经过第t次迭代获得等效源源强向量Q^t如式(5)：

Q^t≡WX^t (5)

第五步：若t≤t_max则将迭代次数赋值为t+1，并转入第二步继续进行迭代；若t>t_max则终止迭代，并令Q^t为等效源源强向量Q。

本实施例中每次迭代都通过前一步迭代计算获得的等效源源强向量的估计值作后验加权系数，通过对正则化过程中的惩罚项进行迭代重加权来限制解空间的大小，从而提高计算结果的精度。

步骤f、通过K折交叉验证法选出最匹配当前等效源源强稀疏度的l_p-optimal：

K折交叉验证法是指：将测量声压P_H均分为K份轮流作验证集，剩余的K-1份作训练集代入步骤e中用于计算；通过m次K折交叉验证法获得由m个l_p值分别作范数约束惩罚项进行正则化过程求解时所得的平均误差，取其中最小平均误差对应的l_p值为l_p-optimal，形成如式(6)所示的以等效源源强向量Q为求解对象的l_p-optimal范数正则化过程：

本实施例中，由于输入量测量声压P_H所包含的元素数量并不多，故一般K取为6-8折，K值优选为8，不超过10折。

步骤g、输入完整的测量声压P_H，重复步骤e，对l_p-optimal范数正则化过程求解，令第五步中求得的等效源源强解Q为输出的最终源强解Q_f。

步骤h、利用P_T＝G_TpQ_f计算获得目标重建面T上的声压P_T；利用V_T＝G_TvQ_f计算获得目标重建面T上的法向振速V_T，其中，G_Tp为各等效源与目标重建面之间的声压传递矩阵，G_Tv为各等效源与目标重建面T之间的振速传递矩阵。

理论模型

本发明方法在声场计算过程中仍采用等效源模型，如图1所示，物体辐射的声场可以由置于该辐射体内部的一系列等效源产生的声场线性叠加代替以近似实际声场；但与现有常用的常规等效源法与压缩等效源法不同的是，在等效源源强求解过程中，本发明采用带有自适应范数约束惩罚项的新型迭代正则化模型，并以此来代替常规等效源法中的Tikhonov正则化模型和压缩等效源法中的l₁范数正则化模型，实现等效源源强灵活且精确的求解，具有普适的应用范围。

根据等效源法基本原理，如图1所示模型中声源产生的声场由等效源面S_e上的一系列等效源辐射声场线性叠加来近似。则声场中任意一点r处的声压p(r)由等效源表示如式(11)：

式(11)中，n表示等效源面S_e上第n个等效源，N为等效源面S_e上离散的等效源总数，j为虚数单位，ρ为媒质密度，c为声速，k为声波波数，r_n和q(r_n)分别表示等效源面S_e上第n个等效源的位置坐标和其对应的源强，g(r,r_n)为场点r与r_n之间的格林函数如式(12)：

式(12)中，|r-r_n|为场点r与第n个等效源所在位置r_n之间的距离。

全息面H位于声源的声场中，则全息面H上M个测量网格点处的声压表示为式(13)：

式(13)中，m表示全息面上第m个测量网格点，m＝1，2，…，M，r_m表示全息面H上第m个测量网格点的位置坐标，p_H(r_m)为r_m处的声压，g(r_m,r_n)为全息面H上第m个测量网格点与第n个等效源所在位置r_n之间的格林函数；

将式(13)表达为矩阵的形式如式(1)，式(1)中的P_H表征全息面H上测量声压，此时即为全息面H上的M个测量网格点处的测量声压组成的列向量，式(1)中的Q表征等效源源强向量，此时即为N个等效源源强组成的列向量，式(1)中的G_Hp表征等效源面上各等效源与全息面上各测量网格点之间的声压传递矩阵，此时即为N个等效源到全息面上M个测量网格点处的声压传递矩阵，并有：

Q＝[q(r₁),q(r₂),…,q(r_N)]^T (14)

G_Hp(m,n)＝jρckg(r_m,r_n) (15)

G_Hp(m,n)为声压传递矩阵G_Hp中位于第m行、第n列的矩阵元素。

对于式(1)逆向求解获得等效源源强向量Q如式(16)：

为G_Hp的广义逆矩阵。

欲保证近场声全息计算结果的精确度，等效源源强的求解过程尤为关键。然而，由于等效源源强向量Q的求解过程是逆问题，其具有不适定性，考虑到实际中噪声不可避免，而输入时微小的误差会引起求解结果的巨大变化，故实际中通常采用正则化方法求解源强向量Q以稳定求解过程。

常规等效源法采用Tikhonov正则化方法求解Q，其正则化过程可归结为求解如式(17)所表征的极小化问题：

min{||P_H-G_HpQ||²+λ||Q||²} (17)

式(17)中：||||表示l₂范数，惩罚项||Q||²是一个l₂范数约束，由于l₂范数约束总会导致一个连续解，会获得相对连续的等效源强解，故在处理空间连续型声源时常规等效源法有着较好的重建效果，而在处理空间稀疏型声源时的效果较差。

针对基于常规等效源法的近场声全息的缺点，后续有学者将压缩感知理论引入到等效源法中，实现了一种在稀疏框架下的等效源法，用l₁范数惩罚项代替常规等效源法正则化过程中的l₂范数惩罚项，形成了基于压缩等效源法的近场声全息，其正则化过程可归结为求解如式(18)所表征的极小化问题：

min{||P_H-G_HpQ||²+λ||Q||₁} (18)

针对式(18)所示的极小化问题通常采用CVX工具包进行求解，其中，||||₁表示l₁范数，惩罚项||Q||₁是一个l₁范数约束。由于在正则化求解中使用l₁范数约束时，在满足一定的稀疏条件下总会导致一个稀疏解，故该方法在处理空间稀疏型声源时有较高的重建精度，但在处理空间连续型声源时却无法获得精确的近场声全息重建结果。

上述由式(17)和式(18)所表征的两种极小化问题在其正则化求解过程中均采用固定范数，没有根据声源实际稀疏度的不同实现正则化方法中惩罚项范数的自适应变化，故对于实际应用中大量存在的稀疏度介于充分稀疏与连续之间的声源，无法很好地匹配其实际稀疏度以获得更优的重建效果，由此可见上述两种方法均存在求解过程灵活性差且普适性低的局限性。

考虑到当运用l_p变范数最小化求解时，会随着l_p数值的不同导致其解稀疏程度的不同，因此为解决上述问题以获得更精确的重建结果，本发明方法将l_p变范数引入等效源法并形成了具有稀疏度自适应性的求解过程。

l_p范数可以理解为是某一向量

的p-范数，其定义由式(19)所表征：

其中，||u||_p表示对向量u求p-范数，u_i为向量u中的元素，i表示向量u中第i个元素，p为取值范围由0到2的实数。

压缩等效源法的l₁范数正则化与常规等效源法的l₂范数正则化分别对应着p＝1与p＝2时的情况；不同的p值导致不同稀疏度解的几何说明如图2所示。

以对式(20)所表征的最优化问题求解为例：

min||x||_ps.t.Ax＝b (20)

对式(20)求解的几何意义是在图2中所示直线L(满足Ax＝b)上找一点使得l_p范数单位球的半径最小。显然，由于当0≤p＜1时l_p球是呈凹状的，当球的半径逐渐增加时，球与直线的交点(即式(20)的解)将位于坐标轴上，而这样的一个点是稀疏的。当p＝1时，l_p球呈菱形，在一定条件下同样会导致一个稀疏解。当1＜p≤2，l_p球呈外凸状，当球的半径逐渐增加时与图像的切点必不位于坐标轴上，即此时的解是不稀疏的。

综上，本发明方法将l_p变范数引入等效源法，构建带有稀疏度自适应范数约束惩罚项的新型正则化方法，以此代替Tikhonov正则化来求解等效源源强向量。

针对经本发明方法求解获得的最终源强Q_f，结合等效源与目标重建面之间的声压、振速传递矩阵能够计算获得目标重建面T上的声压、质点振速等声场数据，从而实现声场的重建和预测；其中，目标重建面T上的声压P_T和法向振速V_T按式(21)和式(22)计算获得：

P_T＝G_TpQ_f (21)

V_T＝G_TvQ_f (22)

G_Tp和G_Tv分别为各等效源与目标重建面T之间的声压传递矩阵和振速传递矩阵，并有：

G_Tp(m,n)＝jρckg(r'_m,r_n) (23)

G_Tp(m,n)为联系各等效源与目标重建面T的声压传递矩阵G_Tp中位于第m行、第n列的矩阵元素，G_Tv(m,n)为联系各等效源与目标重建面T的声压传递矩阵G_Tv中位于第m行、第n列的矩阵元素，r'_m为目标重建面T上第m个网格点的位置，g(r'_m,r_n)为目标重建面T上第m个网格点与第n个等效源所在位置r_n之间的格林函数，n为目标重建面T的外法线方向，

表示求法向导数。

仿真检验1：用于验证本方法在处理空间稀疏型声源时相对于常规等效源法和压缩等效源法具有更高的声场重建精度。

声源为两个位于Z＝0m的平面内的单极子声源，声源位置分别为(0m,0m,0m)和(0m,0.3m,0m)；全息面H位于声源上方Z＝0.08m的平面内，H的尺寸为1m×1m，H上均匀分布21×21个测量网格点，因此全息面H上的测量间隔为0.05m；目标重建面T位于Z＝0.03m的平面内，并且目标重建面T的尺寸大小以及其上的网格点的分布与全息面H上测量网格点的分布一致；等效源面S_e布置在目标重建面T远离全息面H的一侧，与目标重建面T之间的垂直距离为d_h为0.03m，因此等效源面S_e与声源面同在Z＝0m的平面内，并且S_e的尺寸大小以及其上的网格点的分布与全息面H上测量网格点的分布一致；为进一步增加计算难度，并检验方法的鲁棒性，对全息面声压添加了高斯白噪声，信噪比为30dB。

仿真中分别采用本发明方法、常规等效源法以及压缩等效源法对100-2000Hz范围内20个频率(频率间隔100Hz)的单极子声源进行声压重建计算，所有计算中最大迭代次数t_max均取5次。

图3a、图3b、图3c和图3d分别给出了单极子声源频率为2000Hz时目标重建面上的理论声压分布图和本发明方法、压缩等效源法以及常规等效源法重建获得的目标重建面上声压分布图。其中图3a为按照单极子声压解析表达式计算得到的目标重建面理论声压；对比图3a、图3b、图3c和图3d可知，采用本发明方法和压缩等效源法获得的目标重建面声压分布与理论值非常一致；本发明方法获得的重建结果中，声压峰值处与声压值较低处皆能够很好地与理论声压值拟合；而常规等效源法获得的重建结果，由于采用l₂范数作惩罚项的Tikhonov正则化，其计算所得的源强解更为连续，等效源源强解有较大误差，故两个声压峰值处的重建声压值与理论值有偏差，并且在声压较低处无法很好地拟合，与理论值有很大差别。

仿真检验1表明，以空间稀疏型声源为目标声源时，本发明方法在分辨率与重建精度方面，相较常规等效源法的重建结果都有很大优势；本方法与压缩等效源法都能够很好地拟合理论声压分布，而本发明方法可以拟合地更精确。为了定量描述重建结果的精确程度，定义声压重建误差E_P为:

E_P＝||P_Theoretical-P_{Reconstructed}||/||P_Theoretical|| (25)

P_Theoretical表示目标重建面上的理论声压，P_{Reconstructed}表示计算获得的重建声压。

按式(25)对本发明方法在100-2000Hz频率范围内的重建误差进行了计算，并与常规等效源法以及压缩等效源法在相同频率范围内的重建误差进行比较，如图4所示，图4中曲线a1为在100-2000Hz频率范围内采用本方法所得的重建误差，曲线b1为在100-2000Hz频率范围内采用常规等效源法所得的重建误差，曲线c1为在100-2000Hz频率范围内采用压缩等效源法所得的重建误差。

图4表明在所有计算频率处，本发明方法的重建误差均小于常规等效源法和压缩等效源法；本发明方法使用由其声源稀疏度自适应过程所选择的l_p-optimal范数约束惩罚项，与压缩等效源法所采用的l₁范数惩罚项相比能更精准地匹配当前声源的实际稀疏度，故能导致一个更为稀疏的源强解。因此相较于压缩等效源法，本发明方法在所有频率处皆仍具有精度优势，能够达到更小的重建误差。

仿真检验1表明，本发明方法在处理空间稀疏型声源时的重建精度方面均优于基于常规等效源法以及压缩等效源法的近场声全息。

仿真检验2：用于验证本方法在处理空间连续型声源时相对于压缩等效源法具有更高的声场重建精度，与常规等效源法的声场重建精度一致。

声源为一块位于Z＝0m平面内的简支铝板声源，板的面积为0.5m×0.5m，厚度为0.003m，在板中心施加大小为5N的激励力；全息面H位于声源上方Z＝0.08m的平面内，H的尺寸为0.5m×0.5m，H上均匀分布21×21个测量网格点，因此全息面H上的测量间隔为0.025m；目标重建面T位于Z＝0.03m的平面内，并且目标重建面T的尺寸大小以及其上的网格点的分布与全息面H上测量网格点的分布一致；等效源面S_e布置在目标重建面T远离全息面H的一侧，与目标重建面T之间的垂直距离d_h为0.03m，因此等效源面S_e与声源面同在Z＝0m的平面内，并且S_e的尺寸大小以及其上的网格点的分布与全息面H上测量网格点的分布一致；为增加计算难度，并检验方法鲁棒性，对全息面声压添加高斯白噪声，信噪比为30dB。

仿真中同样分别采用本发明方法、常规等效源法以及压缩等效源法对100-2000Hz范围内20个频率(频率间隔100Hz)的0.5m×0.5m简支铝板声源进行了声压重建计算，所有计算中最大迭代次数t_max均取5次。

图5a、图5b、图5c和图5d分别给出了0.5m×0.5m简支铝板声源频率为2000Hz时目标重建面上的理论声压分布图和本发明方法、压缩等效源法以及常规等效源法重建获得的目标重建面上声压分布图。其中图5a为按照简支铝板声压解析表达式计算得到的目标重建面理论声压；对比图5a、图5b、图5c和图5d可知，采用本发明方法和常规等效源法获得的目标重建面声压分布能够很好地与理论声压值拟合；而由于压缩等效源法在正则化过程中采用l₁范数惩罚项，其计算所得的源强解较为稀疏，并不匹配当前简支板声源相对连续的稀疏度，故其重建声压分布情况尤其在声压峰值处的数值与理论值有较大偏差。

仿真检验2表明，以空间连续型声源为目标声源时，本发明方法相较压缩等效源法具有更高的分辨率，重建精度高；并且本发明方法在分辨率与重建精度方面都能够达到与常规等效源法的一样的效果。

为了定量描述重建结果的精确程度，仍按式(25)对本发明方法在100-2000Hz频率范围内的重建误差进行了计算，并与常规等效源法以及压缩等效源法在相同频率范围内的重建误差进行了比较，如图6所示，图6中曲线a2为在100-2000Hz频率范围内采用本方法所得的重建误差，曲线b2为在100-2000Hz频率范围内采用常规等效源法所得的重建误差，曲线c2为在100-2000Hz频率范围内采用压缩等效源法所得的重建误差。

图6表明在所有计算频率处，本发明方法的重建误差均小于压缩等效源法；本发明方法使用由其声源稀疏度自适应过程所挑选的l_p-optimal范数约束惩罚项进行正则化求解，能够求得一个符合当前简支板声源稀疏度相对连续的源强解，故能达到与常规等效源法基本一致的重建效果，并且其重建误差在数值上不超过5％。

仿真检验2表明，本发明方法在处理空间连续型声源时的重建精度方面优于基于压缩等效源法的近场声全息，并能够达到与压缩等效源法的近场声全息一样的精度。

仿真检验3：用于验证本方法在处理稀疏度介于充分稀疏与连续之间的声源时相对于常规等效源法和压缩等效源法皆具有更高的声场重建精度。

仿真条件基本与仿真检验2相同，仅将简支铝板目标声源的尺寸换为0.2m×0.2m。

仿真中同样分别采用本发明方法、常规等效源法以及压缩等效源法对100-2000Hz范围内20个频率(频率间隔100Hz)的0.2m×0.2m简支铝板声源进行了声压重建计算，所有计算中最大迭代次数t_max均取5次。

图7a、图7b、图7c和图7d分别给出了0.2m×0.2m简支铝板声源频率为2000Hz时目标重建面上的理论声压分布图和本发明方法、压缩等效源法以及常规等效源法重建获得的目标重建面上声压分布图。其中图7a为按照简支铝板声压解析表达式计算得到的目标重建面理论声压；对比图7a、图7b、图7c和图7d可知，采用本发明方法获得的目标重建面声压分布与理论值非常一致；而压缩等效源法获得的重建结果，在声压峰值处拟合得较差，常规等效源法获得的重建结果在板边缘处的重建结果不够精确，也与理论值存在偏差。

仿真检验3表明，以源强稀疏度介于充分稀疏与连续之间的声源为目标声源时，本发明方法相对于常规等效源法和压缩等效源法能够更准确地识别声源，更精确地重建出声源的细节信息，具有更高的分辨率和更高的精度。

为了定量描述重建结果的精确程度，仍按式(25)对本发明方法在100-2000Hz频率范围内的重建误差进行了计算，并与常规等效源法以及压缩等效源法在相同频率范围内的重建误差进行了比较，如图8所示，图8中曲线a3为在100-2000Hz频率范围内采用本方法所得的重建误差，曲线b3为在100-2000Hz频率范围内采用常规等效源法所得的重建误差，曲线c3为在100-2000Hz频率范围内采用压缩等效源法所得的重建误差。

图8表明在所有计算频率处，本发明方法的重建误差均小于常规等效源法与压缩等效源法；本发明方法在由其声源稀疏度自适应过程所挑选的l_p-optimal范数约束准则下，能够使所得等效源源强的稀疏度与声源的实际稀疏度更加一致，故能求得一个更为精确的源强解，并且其重建误差在数值上不超过4％。

仿真检验3表明，本发明方法在处理源强稀疏度介于充分稀疏与连续之间的声源时的重建精度方面均优于基于常规等效源法以及压缩等效源法的近场声全息。

仿真检验4：用于验证本方法在不同信噪比条件下的优越性。

仿真对象以及全息面H、目标重建面T与等效源面S_e的布置和网格分布皆与仿真3中相同；仿真中，在1000Hz的频率下，分别采用本发明方法、常规等效源法以及压缩等效源法对10-50dB范围内9个信噪比(信噪比间隔5dB)的简支板声源进行了声压重建计算，所有计算中最大迭代次数t_max均取5次。

为了定量描述重建结果的精确程度，仍按式(25)对本发明方法在10-50dB信噪比范围内的重建误差进行了计算，并与常规等效源法以及压缩等效源法在相同信噪比范围内的重建误差进行了比较，如图9所示，图9中曲线a4为在10-50dB信噪比范围内采用本方法所得的重建误差，曲线b4为在10-50dB信噪比范围内采用常规等效源法所得的重建误差，曲线c4为在10-50dB信噪比范围内采用压缩等效源法所得的重建误差。图9表明在所有信噪比处，本发明方法的重建误差均小于常规等效源法与压缩等效源法。

仿真检验4表明，在不同信噪比条件下，本发明方法的稳定性均优于基于常规等效源法以及压缩等效源法的近场声全息，具备良好的抗干扰能力。

Claims (6)

1.具有声源稀疏度自适应性的变范数等效源近场声全息方法，其特征是按如下步骤进行：

步骤a、于声源的近场区域内布置全息面并通过测量获得所述全息面上的测量声压：

在由目标声源产生的声场中，于近场辐射区域内布置全息面H，所述全息面H为任意形状曲面，在所述全息面H上划分测量网格点，通过传感器测量获得各测量网格点处的测量声压P_H；

步骤b、在声源与全息面H之间布置目标重建面T，在目标重建面T远离全息面H的一侧布置等效源面S_e，所述等效源面S_e与目标重建面T之间的垂直距离为d_h，布置在所述等效源面S_e上的等效源的点数不大于全息面H上测量网格点的点数；

步骤c、建立各等效源与全息面H上测量声压P_H之间的关系如式(1)：

P_H＝G_HpQ (1)

G_Hp为等效源面S_e上各等效源与全息面H上各测量网格点之间的声压传递矩阵；

Q为等效源源强向量；

步骤d、以h为步长，将

线性划分为m个数值，则有

并将m个

的数值分别代入

形成求解等效源源强向量Q的

范数正则化过程；

是数值范围为0到2的实数；λ为正则化参数，正则化参数λ采用L曲线法选择获得；

步骤e、按如下过程采用迭代重加权算法将最小

问题转化为最小

问题进行求解，是指将

范数正则化过程转化为使用最小

范数作惩罚项的标准Tikhonov正则化过程，求解等效源源强向量Q：

第一步：将采用常规等效源法求解所得的等效源源强向量作为等效源源强向量初始解Q⁰，令初始迭代次数为1，设定最大迭代次数为t_max，用t表征迭代次数；

第二步：利用向量Q^t-1作为后验加权系数构建加权矩阵W如式(2)：

所述加权矩阵W是利用向量Q^t-1作为主对角元素构建的对角阵，所述向量Q^t-1为第t-1次迭代计算获得的等效源源强向量；对于首次迭代计算过程，向量Q^t-1即为等效源源强向量初始解Q⁰；

第三步：利用所述加权矩阵W的逆矩阵W^-1对向量Q^t加权，并构造加权后向量Q^t的惩罚项为||W^-1Q^t||，形成如式(3)的极小化过程以获得等效源源强向量Q^t的正则解：

min{||P_H-G_HpQ^t||²+λ||W^-1Q^t||} (3)

λ为正则化参数，正则化参数λ采用L曲线法选择获得；

第四步：通过预条件处理令W^-1Q^t≡X^t，

将式(3)所示的对等效源源强向量Q^t求解的正则化过程转化为如式(4)对中间变量X^t求解的标准Tikhonov正则化过程，并采用奇异值分解法对式(4)求解以获得中间变量X^t的正则解：

min{||P_H-A^tX^t||²+λ||X^t||} (4)

经过第t次迭代获得等效源源强向量Q^t如式(5)：

Q^t≡WX^t (5)

第五步：若t≤t_max则将迭代次数赋值为t+1，并转入第二步继续进行迭代；若t>t_max则终止迭代，并令Q^t为等效源源强向量Q；

步骤f、通过K折交叉验证法选出最匹配当前等效源源强稀疏度的

所述K折交叉验证法是指：将测量声压P_H均分为K份轮流作验证集，剩余的K-1份作训练集代入步骤e中用于计算；通过m次K折交叉验证法获得由m个

值分别作范数约束惩罚项进行正则化过程求解时所得的平均误差，取其中最小平均误差对应的

值为

形成如式(6)所示的以等效源源强向量Q为求解对象的

范数正则化过程：

步骤g、输入完整的测量声压P_H，重复步骤e，对所述

范数正则化过程求解，令第五步中求得的等效源源强向量Q为输出的最终源强解Q_f；

步骤h、利用P_T＝G_TpQ_f计算获得目标重建面T上的声压P_T；利用V_T＝G_TvQ_f计算获得目标重建面T上的法向振速V_T，其中，G_Tp为各等效源与目标重建面之间的声压传递矩阵，G_Tv为各等效源与目标重建面T之间的振速传递矩阵。

2.根据权利要求1所述的具有声源稀疏度自适应性的变范数等效源近场声全息方法，其特征是：所述d_h的取值为0.03m。

3.根据权利要求1所述的具有声源稀疏度自适应性的变范数等效源近场声全息方法，其特征是：所述步长h＝0.1，m＝21，

4.根据权利要求1所述的具有声源稀疏度自适应性的变范数等效源近场声全息方法，其特征是：所述K折交叉验证法中K＝8。

5.根据权利要求1所述的具有声源稀疏度自适应性的变范数等效源近场声全息方法，其特征是：所述最大迭代次数t_max不超过10次。

6.根据权利要求5所述的具有声源稀疏度自适应性的变范数等效源近场声全息方法，其特征是：设置所述最大迭代次数t_max为5次。

Images