CN103728013B - 噪声源识别方法 - Google Patents

Abstract

本发明旨在提供一种噪声源识别方法，包括以下步骤：A、在测量声场中设置测量面，测量面上均匀设置质点振速传感器，测量该测量面上法向质点振速；B、在声源内部设定虚源面，并在虚源面上连续分布设置等效源，所述等效源为单极子声源；C、建立等效源与测量面上法向质点振速之间的传递关系；D、建立等效源与声源面上声压与法向振速之间的传递关系；E、根据上述传递关系求解重建面声压值。本发明噪声源识别方法克服现有技术成本高、计算效率与精度低的缺陷，采用基于Fourier波叠加法结合连续分布等效源与快速FFT算法计算积分，具有较高的计算精度和稳定性。

Description

噪声源识别方法

技术领域

本发明涉及噪声源识别方法领域，具体涉及一种采用振速测量和Fourier波叠加法的噪声源识别方法。

背景技术

近场声全息是一种有效的声源识别技术，现有的基于近场声全息技术的噪声源识别方法有（1）基于空间Fourier变换的声源识别技术。Williams等建立了基于基于平面和柱面空间Fourier变换的声源识别方法。Lee等又将该方法推广到球面坐标系下，建立了基于球面空间声场变换的声识别方法，与平面和柱面声场变换不同的时，球面NAH不能借助快速Fourier变换（FFT）进行计算。这类噪声源识别方法的特点是实现容易，计算速度快，但测量面形状受限，只能是平面、柱面和球面等，而且由于采用Fourier变换，孔径问题所引起窗效应和卷绕误差会大大降低声源识别精度。（2）基于边界元的声源识别方法。对于形状复杂的声源，边界元法较之于上一类方法更加有效。该方法通过离散声源表面，建立声源表面上的声压、振速与测量面上声压之间的关系，从而获得源表面上的声压和法向振速，进而获得空间任意点上的声学量。这一方法对测量面和声源面没有具体要求，且二者之间也没有任何关系，但在计算过程需要处理奇异积分，随着频率的增大，所需测量点增大多,计算精度和效率会因此大大降低。（3）基于等效源法的声源识别方法。Helmholtz方程最小均方误差法（HELS）和波叠加法是两个重要的基于等效源的声源识别技术。这两种方法均具有测量面和声源面适应性好，无需处理奇异积分，测量点较少，计算效率较高等。然而HELS对长宽高比例较大的声源重建效果较差，尽管波叠加法能够较好地识别长宽高比例较大的声源，然而目前已有的基于波叠加法的声源识别方法主要通过增加等效源面上配置的离散等效源数目来改进重建精度，当在声源频率较高的情况下，则需要更多的等效源数目才能保证重建精度，这就需要增加的测量点，从而导致测量成本的增加。

发明内容

本发明旨在提供一种噪声源识别方法，该噪声源识别方法克服现有技术成本高、计算效率与精度低的缺陷，采用基于Fourier波叠加法结合连续分布等效源与快速FFT算法计算积分，具有较高的计算精度和稳定性。

本发明的技术方案如下：一种噪声源识别方法，包括以下步骤：

A、在测量声场中设置测量面，测量面上呈网格式设置质点振速传感器，所述的质点振速传感器测量位于测量面法向上的质点振速；

B、在声源内部设定虚源面，并在虚源面上连续分布设置等效源，将等效源设置为单极子声源；

C、采用波叠加法，建立等效源与测量面法向上质点振速之间的传递关系，计算等效源强函数，获得源强展开系数列向量；

D、根据等效源源强展开系数列向量，建立测量面法向振速与声源面上声压与法向振速之间的传递关系；

E、根据步骤C、D中的传递关系确定重建面声压和法向振速幅值。

所述步骤C的具体步骤如下：

根据波叠加法，声源外声场的辐射声压是由分布在声源内部的等效源加权叠加而成，测量面上场点P处的声压和法向质点振速为：

$p\left(P\right)={\∫}_{S^{\′}}\mathrm{\σ}\left(Q\right)G(P,Q)d{S}^{\′},P\∈H,Q\∈S^{\′}---\left(1\right);$

$v\left(P\right)={\∫}_{S^{\′}}\mathrm{\σ}\left(Q\right)\frac{\∂G(P,Q)}{{\∂n}_P}d{S}^{\′},P\∈H,Q\∈S^{\′}---\left(2\right);$

其中为虚数单位，d为场点P和源点Q之间的距离，k₀=ω/c为波数，ω为激励的圆频率，c为声速，n_P为P点的法向单位向量；

球坐标系下场点P和源点Q的坐标分别为和

虚源面上连续分布的等效源的源强密度函数为对其进行双向傅里叶展开，得

将上式代入式（1）和（2）可得

对式（4）和（5）分别进行无穷级数截断至有限项，得到

$p\left(P\right)=\underset{m=-m_c}{\overset{+m_c}{\mathrm{\Σ}}}\underset{n=-n_c}{\overset{+n_c}{\mathrm{\Σ}}}{b}_{\mathrm{mn}}{K}_{\mathrm{mn}}---\left(6\right);$

$v\left(P\right)=\underset{m=-m_c}{\overset{+m_c}{\mathrm{\Σ}}}\underset{n=-n_c}{\overset{+n_c}{\mathrm{\Σ}}}{b}_{\mathrm{mn}}{K}_{\mathrm{mn}}^{\′}---\left(7\right);$

其中

将式(8)和(9)离散后,采用快速傅里叶进行计算；

由式（7）可得测量面上各测量点的质点振速可表示为

$-\mathrm{i\ρck}{v}_H\left(P_i\right)=\underset{m=-m_c}{\overset{+m_c}{\mathrm{\Σ}}}\underset{n=-n_c}{\overset{+n_c}{\mathrm{\Σ}}}{b}_{\mathrm{mn}}{K}_{\mathrm{mn}}^{\′i},i=\mathrm{1,2},...,M.---\left(10\right)$

写成矩阵形式为：

$V_H=\frac{1}{-\mathrm{i\ρck}}{K}_H^{\′}b---\left(11\right)$

式中V_H为测量面H上M个测量点处的质点振速幅值构成的列向量；ρ为空气密度；b为等效源源强展开系数组成的列向量，K′_H为等效源源强展开系数与测量面H法向振速之间的传递矩阵；

将K′_H写成分块矩阵的形式，即

$K_H^{\′}={\left[\begin{array}{cccc}{K}_H^{\′1}& K_H^{\′2}& ...& K_H^{\′M}\end{array}\right]}^T,$ 各个分块元素为行向量，即

$K_H^{\′i}=[K_{-m_c,-n_c}^{\′i},K_{-m_c,-n_c+1}^{\′i},...,K_{-m_c,n_c}^{\′i},K_{-m_c+1,-n_c}^{\′i}...,K_{m_c,n_c}^{\′i}],$

$b=[b_{-m_c,-n_c},b_{-m_c,-n_c+1},...,b_{-m_c,n_c},b_{-m_c+1,-n_c}...,b_{m_c,n_c}{]}^T;$

将式（6）写成矩阵形式，即

p(P)=Kb(12)

其中 $K=[K_{-m_c,-n_c},K_{-m_c,-n_c+1},...,K_{-m_c,n_c},K_{-m_c+1,-n_c}...,K_{m_c,n_c}]$

由（11）式可知，当矩阵K′_H的阶数满足M≥N,其中N=(2m_c+1)×(2n_c+1)

即测量点的个数不小于等效源源强展开系数的数目时，则可以通过奇异值分解唯一确定源强展开系数列向量，即

b=(K′_H)⁺V_H(13)

式中“+”表示广义逆。

所述步骤D中等效源源强展开系数与声源面上声压与法向振速之间的传递关系如下：

p_S=K_Sb（14）；

$V_S=\frac{1}{-\mathrm{i\ρck}}{K}_S^{\′}b---\left(15\right);$

其中p_S为声源面S上的声压；V_S为声源面S上的法向振速；K_S为等效源源强展开系数与声源面S上声压之间的传递矩阵,其计算方法与K_H的式（8）的计算方法相同；K′_S为等效源源强展开系数与声源面S上法向振速之间的传递矩阵,其计算方法与K′_H的式（9）的计算方法相同；“+”表示广义逆。

所述声源面为平面、柱面或球面；所述测量面为平面或曲面。

本发明噪声源识别方法的优点如下：

1、本发明采用基于Fourier波叠加法，与其他形式的基于等效源的声源识别方法相比，在同等重建精度下，具有测量点少，测量成本低的特点；

2、本发明采用连续分布等效源并采用快速FFT算法计算积分，在频率高时，能够获得更好的计算精度和稳定性。

为体现本发明的有益效果，在实施例中进行了本发明基于Fourier波叠加法与离散波叠加法在声源表面声压和法向振速的重建的对比分析，将附图中的效果进行对比可以看出，本发明方法在重建面上无论是声压重建值还是振速重建值都与理论值非常吻合，而采用离散波叠加法的重建效果与理论值有一定差异，在相同条件下，本发明基于Fourier波叠加法的噪声源识别方法具有成本低、精度高的优点。

附图说明

图1为本发明方法噪声源识别方法的流程图

图2为本实施例方法虚源面与测量面布置示意图

图3为本实施例方法测量面上的振速幅值分布

图4为本实施例方法测量面上的振速相位分布

图5为采用离散波叠加法和本实施例方法重构得到的各重建点声压实部与理论解的比较图

图6为采用离散波叠加法和本实施例方法重构得到的各重建点声压虚部与理论解的比较图

图7为采用离散波叠加法和本实施例方法重构得到的各重建点振速实部与理论解的比较图

图8为采用离散波叠加法和本实施例方法重构得到的各重建点振速虚部与理论值的比较图

图9为不同频率下，采用离散波叠加法和本发明方法的声压重建误差图

图10为不同频率下，采用离散波叠加法和本发明方法的振速重建误差图

图2中各部分名称及标号如下:

1为测量面，2为传声器，3为虚源面，4为声源面。

具体实施方式

下面结合附图与实施例具体说明本发明。

实施例1

如图1所示，本实施例采用振速测量和Fourier波叠加法的噪声源识别方法的步骤如下：

A、在测量声场中设置测量面，测量面上呈网格式设置质点振速传感器，所述的质点振速传感器测量位于测量面法向上的质点振速；

B、在声源内部设定虚源面，并在虚源面上连续分布设置等效源，将等效源设置为单极子声源；

C、采用波叠加法，建立等效源与测量面法向上质点振速之间的传递关系，计算等效源强函数，获得源强展开系数列向量；

D、根据等效源源强展开系数列向量，建立测量面法向振速与声源面上声压与法向振速之间的传递关系；

E、根据步骤C、D中的传递关系确定重建面声压和法向振速幅值。

所述步骤C具体步骤如下：

根据波叠加法，声源外声场的辐射声压是由分布在声源内部的等效源加权叠加而成，测量面上场点P处的声压和法向质点振速为：

$p\left(P\right)={\∫}_{S^{\′}}\mathrm{\σ}\left(Q\right)G(P,Q)d{S}^{\′},P\∈H,Q\∈S^{\′}---\left(1\right);$

$v\left(P\right)={\∫}_{S^{\′}}\mathrm{\σ}\left(Q\right)\frac{\∂G(P,Q)}{{\∂n}_P}d{S}^{\′},P\∈H,Q\∈S^{\′}---\left(2\right);$

其中为虚数单位，d为场点P和源点Q之间的距离，k₀=ω/c为波数，ω为激励的圆频率，c为声速，n_P为P点的法向单位向量；

球坐标系下场点P和源点Q的坐标分别为和

虚源面上连续分布的等效源的源强密度函数为对其进行双向傅里叶展开，得

将上式代入式（1）和（2）可得

对式（4）和（5）分别进行无穷级数截断至有限项，得到

$p\left(P\right)=\underset{m=-m_c}{\overset{+m_c}{\mathrm{\Σ}}}\underset{n=-n_c}{\overset{+n_c}{\mathrm{\Σ}}}{b}_{\mathrm{mn}}{K}_{\mathrm{mn}}---\left(6\right);$

$v\left(P\right)=\underset{m=-m_c}{\overset{+m_c}{\mathrm{\Σ}}}\underset{n=-n_c}{\overset{+n_c}{\mathrm{\Σ}}}{b}_{\mathrm{mn}}{K}_{\mathrm{mn}}^{\′}---\left(7\right);$

其中

将式(8)和(9)离散后,采用快速傅里叶进行计算；

由式（7）可得测量面上各测量点的法向质点振速可表示为

$-\mathrm{i\ρck}{v}_H\left(P_i\right)=\underset{m=-m_c}{\overset{+m_c}{\mathrm{\Σ}}}\underset{n=-n_c}{\overset{+n_c}{\mathrm{\Σ}}}{b}_{\mathrm{mn}}{K}_{\mathrm{mn}}^{\′i},i=\mathrm{1,2},...,M.---\left(10\right);$

写成矩阵形式为：

$V_H=\frac{1}{-\mathrm{i\ρck}}{K}_H^{\′}b---\left(11\right);$

式中V_H为测量面H上M个测量点处的法向质点振速幅值构成的列向量，b为等效源源强展开系数组成的列向量，K′_H为等效源源强展开系数与测量面H法向振速之间的传递矩阵；

将K′_H写成分块矩阵的形式，即

$K_H^{\′}={\left[\begin{array}{cccc}{K}_H^{\′1}& K_H^{\′2}& ...& K_H^{\′M}\end{array}\right]}^T,$ 各个分块元素为行向量，即

$K_H^{\′i}=[K_{-m_c,-n_c}^{\′i},K_{-m_c,-n_c+1}^{\′i},...,K_{-m_c,n_c}^{\′i},K_{-m_c+1,-n_c}^{\′i}...,K_{m_c,n_c}^{\′i}],$

$b=[b_{-m_c,-n_c},b_{-m_c,-n_c+1},...,b_{-m_c,n_c},b_{-m_c+1,-n_c}...,b_{m_c,n_c}{]}^T;$

将式（6）写成矩阵形式，即

p(P)=Kb(12)；

其中 $K=[K_{-m_c,-n_c},K_{-m_c,-n_c+1},...,K_{-m_c,n_c},K_{-m_c+1,-n_c}...,K_{m_c,n_c}]$

由（11）式可知，当矩阵K′_H的阶数满足M≥N,其中N=(2m_c+1)×(2n_c+1)

即测量点的个数不小于等效源源强展开系数的数目时，则可以通过奇异值分解唯一确定源强展开系数列向量，即

b=(K′_H)⁺V_H(13)；

式中“+”表示广义逆。

所述步骤D中等效源源强展开系数与声源面上声压与法向振速之间的传递关系如下：

p_S=K_Sb（14）；

$V_S=\frac{1}{-\mathrm{i\ρck}}{K}_S^{\′}b---\left(15\right);$

其中p_S为声源面S上的声压；V_S为声源面S上的法向振速；ρ为空气密度；K_S为等效源源强展开系数与声源面S上声压之间的传递矩阵,其计算方法与K_H的式（8）的计算方法相同；K′_S为等效源源强展开系数与声源面S上法向振速之间的传递矩阵,其计算方法与K′_H的式（9）的计算方法相同；“+”表示广义逆。

本实施例声源采用脉动球，对于半径为R的脉动球，在场点P处的声压的解析解为

$p\left(P\right)=v\left(\frac{R}{r}\right)\left(\frac{\mathrm{i\ρ}\left(2\mathrm{\πf}\right)R}{1+\mathrm{ikR}}\right)e^{-\mathrm{ik}(r-R)}---\left(16\right);$

式中，均匀法向速度v=1m/s，空气密度为ρ=1.2kg/m³，f为声源振动频率，k=2πf/c，计算中声速为c=344m/s；球源半径为0.1m，测量面H为1m×1m的平面，测量面H距声源中心的距离d为0.3m，测量面上均匀地分布11×11个测量点，球形声源表面分布60个重建点（径向和周向间隔都为π/6，不包括在z轴上的那两个点）；

虚源面与测量面布置如图2所示；

本实施例频率为5000Hz；

测量面上的振速幅值分布、振速相位分布分别如图3、4所示；

采用离散波叠加法和本实施例方法重建声压的实部和虚部与理论声压的比较分别如图5、6所示，由图5可以看到，在频率低于3000Hz时，离散波叠加法和本实施例方法都有着良好的重建精度，当频率在大于3000Hz时，离散波叠加法的重建精度开始变差，甚至误差达到300%以上，而离散波叠加法要改善重建精度的话，需要增加离散虚源数目，这就需要增加测量点的数目，引起测量成本的增加；而本实施例方法在测量点数目不变的情况下，且频率在3000-6000Hz之间仍然具有比离散波叠加法更高的重建精度；

采用离散波叠加法和本实施例方法重建振速的实部和虚部与理论振速的比较分别如图7、8所示；

为了定量区分两种方法的重建精度，定义声压重建误差和振速重建误差分别为

${\mathrm{\ϵ}}_p=\sqrt{\underset{i=1}{\overset{L}{\mathrm{\Σ}}}{\left(\right|{\stackrel{\‾}{p}}_i-p_i\left|\right)}^2}/\sqrt{\sqrt{\underset{i=1}{\overset{L}{\mathrm{\Σ}}}{\left|p_i\right|}^2}}\×100(\%)---\left(17\right);$

${\mathrm{\ϵ}}_V=\sqrt{\underset{i=1}{\overset{L}{\mathrm{\Σ}}}{\left(\right|{\stackrel{\‾}{V}}_i-V_i\left|\right)}^2}/\sqrt{\sqrt{\underset{i=1}{\overset{L}{\mathrm{\Σ}}}{\left|V_i\right|}^2}}\×100(\%)---\left(18\right);$

式中，L为声源面的重建点总数，分别为第i个重建点处的重建声压和重建振速，p_i,V_i则为第i个重建点处的理论声压和理论振速。

在不同频率下，采用同样多测量点时，采用本实施例方法的声压重建误差与离散波叠加法声压重建误差的对比如图9所示；

在不同频率下，采用同样多测量点时，采用本实施例方法的振速重建误差与离散波叠加法振速重建误差的对比如图10所示。

Claims (4)

1.一种噪声源识别方法，其特征在于包括以下步骤：

A、在测量声场中设置测量面，测量面上呈网格式设置质点振速传感器，所述的质点振速传感器测量位于测量面法向上的质点振速；

B、在声源内部设定虚源面，并在虚源面上连续分布设置等效源，将等效源设置为单极子声源；

C、采用Fourier波叠加法与快速FFT算法，建立等效源与测量面法向上质点振速之间的传递关系，计算等效源强函数，获得源强展开系数列向量；

D、根据等效源源强展开系数列向量，建立测量面法向振速与声源面上声压与法向振速之间的传递关系；

E、根据步骤C、D中的传递关系确定重建面声压和法向振速幅值。

2.如权利要求1所述的噪声源识别方法，其特征在于：

所述步骤C的具体步骤如下：

根据波叠加法，声源外声场的辐射声压是由分布在声源内部的等效源加权叠加而成，测量面上场点P处的声压和法向质点振速为：

p(P)＝∫_S′σ(Q)G(P,Q)dS′,P∈H,Q∈S′(1)；

$v\left(P\right)={\∫}_{S^{\′}}\mathrm{\σ}\left(Q\right)\frac{\∂G(P,Q)}{\∂n_P}{\mathrm{dS}}^{\′},P\∈H,Q\∈S^{\′}---\left(2\right);$

其中i²＝-1为虚数单位，d为场点P和源点Q之间的距离，k₀＝ω/c为波数，ω为激励的圆频率，c为声速，n_P为P点的法向单位向量；

球坐标系下场点P和源点Q的坐标分别为和

虚源面上连续分布的等效源的源强密度函数为对其进行双向傅里叶展开，得

将上式代入式(1)和(2)可得

对式(4)和(5)分别进行无穷级数截断至有限项，得到

$p\left(P\right)=\underset{m=-m_c}{\overset{+m_c}{\Σ}}\underset{n=-n_c}{\overset{+n_c}{\Σ}}{b}_{mn}{K}_{mn}---\left(6\right);$

$v\left(P\right)=\underset{m=-m_c}{\overset{+m_c}{\Σ}}\underset{n=-n_c}{\overset{+n_c}{\Σ}}{b}_{mn}{K}_{mn}^{\′}---\left(7\right);$

其中

将式(8)和(9)离散后,采用快速傅里叶进行计算；

由式(7)可得测量面上各测量点的质点振速可表示为

$-{\mathrm{i\ρckv}}_H\left(P_i\right)=\underset{m=-m_c}{\overset{+m_c}{\Σ}}\underset{n=-n_c}{\overset{+n_c}{\Σ}}{b}_{mn}{K}_{mn}^{\′i},i=1,2,...,M.---\left(10\right)$

写成矩阵形式为：

$V_H=\frac{1}{-i\mathrm{\ρ}ck}{K}_H^{\′}b---\left(11\right)$

式中V_H为测量面H上M个测量点处的质点振速幅值构成的列向量；ρ为空气密度；b为等效源源强展开系数组成的列向量，K′_H为等效源源强展开系数与测量面H法向振速之间的传递矩阵；

将K′_H写成分块矩阵的形式，即

$K_H^{\′}={\left[\begin{array}{cccc}{K}_H^{\′1}& K_H^{\′2}& ...& K_H^{\′M}\end{array}\right]}^T,$ 各个分块元素为行向量，即

$K_H^{\′i}=\[K_{-m_c,-n_c}^{\′i},K_{-m_c,-n_c+1}^{\′i},...,K_{-m_c,n_c}^{\′i},K_{-m_c+1,-n_c}^{\′i}...,K_{m_c,n_c}^{\′i}\],$

$b={\[b_{-m_c,-n_c},b_{-m_c,-n_c+1},...,b_{-m_c,n_c},b_{-m_c+1,-n_c}...,b_{m_c,n_c}\]}^T;$

将式(6)写成矩阵形式，即

p(P)＝Kb(12)

其中 $K=\[K_{-m_c,-n_c},K_{-m_c,-n_c+1},...,K_{-m_c,n_c},K_{-m_c+1,-n_c}...,K_{m_c,n_c}\]$

由(11)式可知，当矩阵K′_H的阶数满足M≥N,其中N＝(2m_c+1)×(2n_c+1)

即测量点的个数不小于等效源源强展开系数的数目时，则可以通过奇异值分解唯一确定源强展开系数列向量，即

b＝(K′_H)⁺V_H(13)

式中“+”表示广义逆。

3.如权利要求2所述的噪声源识别方法，其特征在于：

所述步骤D中等效源源强展开系数与声源面上声压与法向振速之间的传递关系如下：

p_S＝K_Sb(14)；

$V_S=\frac{1}{-i\mathrm{\ρ}ck}{K}_S^{\′}b---\left(15\right);$

其中p_S为声源面S上的声压；V_S为声源面S上的法向振速；K_S为等效源源强展开系数与声源面S上声压之间的传递矩阵,其计算方法与K_H的式(8)的计算方法相同；K′_S为等效源源强展开系数与声源面S上法向振速之间的传递矩阵,其计算方法与K′_H的式(9)的计算方法相同；“+”表示广义逆。

4.如权利要求1～3任意一项所述的噪声源识别方法，其特征在于：所述声源面为平面、柱面或球面；所述测量面为平面或曲面。