CN103941229A - 一种局部近场声全息法的移动噪声源识别方法 - Google Patents

Abstract

本发明属于噪声类声场重建和可视化领域，具体涉及一种局部近场声全息法的移动噪声源识别方法。本发明包括：测量测量面H上声压；对测量得到的声压进行多普勒效应去除，得到测量面H上不含多普勒效应的声压；测量面H上的声压补零扩展得到测量面H+上声压；计算测量面H+上的声压与重构面S上声压之间的传递矩阵；求解得到重构面S上的声压。本发明采用快速局部近场声全息法为基础，与传统法相比，本发明具有计算简单、计算时间短及计算效率高等优点。

Description

一种局部近场声全息法的移动噪声源识别方法

技术领域

本发明属于噪声类声场重建和可视化领域，具体涉及一种局部近场声全息法的移动噪声源识别方法。

背景技术

局部近场声全息技术是一种重要的声源定位和声场可视化技术，实际测量时该技术只需采用较小测量孔径就可以对大尺寸声源进行声场重建，能够有效减小窗效应和卷绕误差，但在对移动大尺寸声源进行声场重构时该方法就无法处理，为此人们将该技术与移动框架声全息技术结合来解决大尺寸移动声源的声场重构问题。现有技术包括：（1）基于统计最优近场声全息技术的移动声源识别方法（《哈尔滨工程大学学报》（2010年31卷7期））。该方法结合统计最优法和移动框架技术各自特点，可以在局部测量的条件下对移动声源进行声场重建，该方法的缺点是受到矩阵计算的影响，计算时间较长，计算过程复杂计算量大；（2）基于改进统计最优近场声全息技术的移动声源识别方法（《Journal of Sound and Vibration》（2012年331期））。该方法是在统计最优近场声全息技术的基础上，对该方法进行了改进，通过与已有方法的对比发现，该方法能够获得更好的重构精度和更小的重构误差，但是该方法仍然无法避免大规模矩阵计算的问题，计算时间仍然较长。

发明内容

本发明的目的在于提供一种适用于大型移动声源、计算过程简单、计算时间短的局部近场声全息法的移动噪声源识别方法。

本发明的目的是这样实现的：

（1）测量测量面H上声压：

位于声源近场有测量面H，在测量面上呈网格式分布有测量点，相邻网格点之间的距离小于半个波长；移动声源运动方向所在平面为重构面S，移动声源运动距离为L，在重构面S上分布有网格点，网格点的分布与测量面H上的测量点相同；测量所述测量面H上各网格点处的声压幅值和相位信息获得测量面H上含有多普勒效应的声压；

（2）对测量得到的声压进行多普勒效应去除，得到测量面H上不含多普勒效应的声压：

$p_H\left(H\right)=p_H^{\′}\left(H\right)R{(1-\mathrm{Ma}\cos \mathrm{\θ})}^2{e}^{\mathrm{jkR}}\frac{e^{-\mathrm{jkr}}}{r},$ 其中

$R=\frac{\mathrm{Ma}\·[x_j-x_i]}{1-{\mathrm{Ma}}^2}+\frac{\sqrt{{(x_j-x_i)}^2+(1-{\mathrm{Ma}}^2)[{(y_j-y_i)}^2+d^2]}}{1-{\mathrm{Ma}}^2},$ cosθ＝(x_j-x_i)/R、

$r=\sqrt{{(x_j-x_i)}^2+{(y_j-y_i)}^2+d^2},$ Ma=v/c；

pH(H)为测量面H上不含多普勒效应的声压、p'H(H)为测量面H上测量得到含有多普勒效应的声压、r为测量面H上网格点与重构面S上网格点之间的距离、cosθ为测量面H上网格点和重构面S上网格点之间的距离与声源运动方向夹角的余弦、R为声波从移动声源位置到测量面H上网格点的实际传播距离、c为介质中声速、v为移动声源的运动速度、Ma为移动声源的运动马赫数、k为波数、d为测量面H与重构面S之间的距离、(x_i,y_i)为测量面H上网格点的坐标、(x_j,y_j)为重构面S上网格点的坐标，测量面H与重构面S平行，测量面H和重构面S的法向均为z方向；

（3）测量面H上的声压补零扩展得到测量面H+上声压：

$p_E(H+)=\left\{\begin{array}{cc}{p}_H\left(H\right)& (x_i,y_i)\∈H\\ 0& (x_i,y_i)\∉H\end{array}\right.,$ p_E(H+)＝D·p_H(H+)，

其中D＝diag[D₁₁,…,D_MM]， $D_{\mathrm{ii}}=\left\{\begin{array}{cc}1& (x_i,y_i)\∈H\\ 0& (x_i,y_i)\∉H\end{array}\right.,$ p_E(H+)=[p_E1，…，p_EM]^T，

p_H(H+)=[p_H1，…，p_HM]^T，p_E(H+)为测量面H+上由p_H(H)补零后得到的声压、p_H(H+)为测量面H+上的声压、D为采样算子、D_ii为矩阵D中对角线上的值、M为测量面H+上测量网格的点数、T为矩阵转置；

（4）计算测量面H+上的声压与重构面S上声压之间的传递矩阵：

$W_D={\mathrm{DB}}_{k_c}{F}^{-1}{G}_{D}F$

其中， $G_D=g^{j{k}_{z}d},B_{k_c}=F^{-1}{K}_{k_c}F,k_z=\left\{\begin{array}{cc}\sqrt{k^2-k_x^2-k_y^2}& k_x^2+k_y^2\&le;k^2\\ -j\sqrt{k_x^2+k_y^2-k^2}& k_x^2+k_y^2>k^2\end{array}\right.,$ $L_{k_c}=\mathrm{diag}[L_{11},\&CenterDot;\&CenterDot;\&CenterDot;,L_{\mathrm{MM}}],L_{\mathrm{ii}}=\left\{\begin{array}{cc}1& \sqrt{k_x^2+k_y^2}<k_c\\ 0& \sqrt{k_x^2+k_y^2}>k_c\end{array}\right.;$

W_D为测量面H+上由p_H(H)补零后得到的声压p_E(H+)与重构面S上声压的传递矩阵、G_D为测量面H+上声压p_H(H+)与重构面S上声压的格林函数、F为二维空间傅里叶变换因子、F^-1为二维空间傅里叶变换因子的逆变换、为波数域内带限算子、k_c为截止波数、为低通滤波器、L_ii为矩阵中对角线上的值；

（5）求解得到重构面S上的声压：

$p_S={(\mathrm{\&alpha;I}+W_D^H{W}_D)}^{-1}{W}_D^H{p}_E(H+)$

其中p_S为重构面S上的声压、α为正则化参数、I为单位对角矩阵、为W_D的共轭转置矩阵、为矩阵的逆矩阵。

本发明的有益效果在于：

本发明采用快速局部近场声全息法为基础，与传统法相比，本发明具有计算简单、计算时间短及计算效率高等优点。

附图说明

图1为采用快速局部近场声全息法的移动噪声源识别方法示意图；

图2为频率为500Hz，声源运动速度为3m/s时，全息面H上测量的声压分布；

图3为频率为500kHz，声源运动速度为3m/s时，全息面H上测量的声压分布与声源静置时相同频率下的声压分布比较；

图4为频率为500kHz，声源运动速度为3m/s时，重构面S上真实声压分布；

图5a为频率为500kHz，声源运动速度为3m/s时，采用常规局部近场声全息法重建的重构面S上声压分布；

图5b为频率为500kHz，声源运动速度为3m/s时，采用本发明方法重建的重构面S上声压分布；

图6a为频率为500kHz，声源运动速度为10m/s时，采用常规局部近场声全息法重建的重构面S上声压分布；

图6b为频率为500kHz，声源运动速度为10m/s时，采用本发明方法重建的重构面S上声压分布；

图7为不同频率时，采用本发明方法得到的声压的均方根误差率。

具体实施方式

本发明方法是测量位于移动声源近场的有限孔径测量面H上的声压，对测得的声压进行多普勒去除，然后对其进行补零扩展，计算测量面与重构面之间的传递矩阵，最终实现声源表面声场的重建。

a、测量测量面H上声压

在移动声源的辐射声场中，位于声源近场有测量面H，在测量面上呈网格式分布有测量点，相邻网格点之间的距离小于半个波长；移动声源运动方向所在平面为重构面S，移动声源运动距离为L，在重构面S上分布有网格点，网格点的分布与测量面H上的测量点相同；测量所述测量面H上各网格点处的声压幅值和相位信息获得测量面H上含有多普勒效应的声压。

b、对测量得到的声压进行多普勒效应去除，得到测量面H上不含多普勒效应的声压

$p_H\left(H\right)=p_H^{\&prime;}\left(H\right)R{(1-\mathrm{Ma}\cos \mathrm{\&theta;})}^2{e}^{\mathrm{jkR}}\frac{e^{-\mathrm{jkr}}}{r},$ 其中

$R=\frac{\mathrm{Ma}\&CenterDot;[x_j-x_i]}{1-{\mathrm{Ma}}^2}+\frac{\sqrt{{(x_j-x_i)}^2+(1-{\mathrm{Ma}}^2)[{(y_j-y_i)}^2+d^2]}}{1-{\mathrm{Ma}}^2},$

cosθ＝(xj-xi)/R、

$r=\sqrt{{(x_j-x_i)}^2+{(y_j-y_i)}^2+d^2},$

Ma=v/c；

p_H(H)为测量面H上不含多普勒效应的声压、

p'_H(H)为测量面H上测量得到含有多普勒效应的声压、

r为测量面H上网格点与重构面S上网格点之间的距离、

cosθ为测量面H上网格点和重构面S上网格点之间的距离与声源运动方向夹角的余弦、

R为声波从移动声源位置到测量面H上网格点的实际传播距离、

c为介质中声速、  v为移动声源的运动速度、

Ma为移动声源的运动马赫数、k为波数、

d为测量面H与重构面S之间的距离、

(x_i,y_i)为测量面H上网格点的坐标、(x_j,y_j)为重构面S上网格点的坐标，测量面H与重构面S平行，测量面H和重构面S的法向均为z方向；

c、测量面H上的声压补零扩展得到测量面H+上声压

$p_E(H+)=\left\{\begin{array}{cc}{p}_H\left(H\right)& (x_i,y_i)\&Element;H\\ 0& (x_i,y_i)\&NotElement;H\end{array}\right.$

p_E(H+)＝D·p_H(H+)，其中

D＝diag[D₁₁,…,D_MM]

$D_{\mathrm{ii}}=\left\{\begin{array}{cc}1& (x_i,y_i)\&Element;H\\ 0& (x_i,y_i)\&NotElement;H\end{array}\right.$

p_E(H+)=[p_E1，…，p_EM]^T

p_H(H+)=[p_H1，…，p_HM]^T

p_E(H+)为测量面H+上由p_H(H)补零后得到的声压、

p_H(H+)为测量面H+上的声压、

D为采样算子、D_ii为矩阵D中对角线上的值、

M为测量面H+上测量网格的点数、T为矩阵转置；

d、计算测量面H+上的声压与重构面S上声压之间的传递矩阵

$W_D={\mathrm{DB}}_{k_c}{F}^{-1}{G}_{D}F$ 其中

$G_D=e^{j{k}_{z}d},$

$B_{k_c}=F^{-1}{L}_{k_c}F,$

$k_z=\left\{\begin{array}{cc}\sqrt{k^2-k_x^2-k_y^2}& k_x^2+k_y^2\&le;k^2\\ -j\sqrt{k_x^2+k_y^2-k^2}& k_x^2+k_y^2>k^2\end{array}\right.,$

$L_{k_c}=\mathrm{diag}[L_{11},\&CenterDot;\&CenterDot;\&CenterDot;,L_{\mathrm{MM}}],L_{\mathrm{ii}}=\left\{\begin{array}{cc}1& \sqrt{k_x^2+k_y^2}<k_c\\ 0& \sqrt{k_x^2+k_y^2}>k_c\end{array}\right.;$

W_D为测量面H+上由p_H(H)补零后得到的声压p_E(H+)与重构面S上声压的传递矩阵、

G_D为测量面H+上声压p_H(H+)与重构面S上声压的格林函数、

F为二维空间傅里叶变换因子、F^-1为二维空间傅里叶变换因子的逆变换、

为波数域内带限算子、k_c为截止波数、

为低通滤波器、L_ii为矩阵中对角线上的值；

e、求解得到重构面S上的声压

$p_S={(\mathrm{\&alpha;I}+W_D^H{W}_D)}^{-1}{W}_D^H{p}_E(H+),$ 其中

p_S为重构面S上的声压、

α为正则化参数、I为单位对角矩阵、

为W_D的共轭转置矩阵、为矩阵的逆矩阵。

参见图1，本实施例中，移动声源以匀速v由位置（I）运动到位置（III），运动距离为L，移动声源运动方向所在平面为重构面S，在重构面上分布有网格点，重构面S上网格点的坐标为(x_j,y_j)。在移动声源的辐射声场中，位于声源近场有测量面H，在测量面上呈网格式分布有测量点，相邻网格点之间的距离小于半个波长，(x_i,y_i)为测量面H上网格点的坐标，测量面H与重构面S平行，测量面H与重构面S距离为d。当移动声源运动到（II）位置即x_j＝L/2位置时，在测量面H上采用传声器阵列快照法同时采样获得测量面上的声压，即

$p_H^{\&prime;}\left(H\right)=\frac{1}{4\mathrm{\&pi;}}\frac{p_H[t-R/c}{R{(1-\mathrm{Ma}\cos \mathrm{\&theta;})}^2}---\left(1\right)$

式中p_H为测量面H上不含多普勒效应的声压、p'_H(H)为测量面H上含多普勒效应的声压、c为介质中声速、Ma为移动声源的运动马赫数Ma=v/c、R为声波从移动声源位置到测量面H上网格点的实际传播距离、cosθ为测量面H上网格点和重构面S上网格点之间的距离与声源运动方向夹角的余弦，R和cosθ分别定义为

$R=\frac{\mathrm{Ma}\&CenterDot;[x_j-x_i]}{1-{\mathrm{Ma}}^2}+\frac{\sqrt{{(x_j-x_i)}^2+(1-{\mathrm{Ma}}^2)[{(y_j-y_i)}^2+d^2]}}{1-{\mathrm{Ma}}^2}---\left(2\right)$

cosθ＝(x_j-x_i)/R   （3）

由运动声源理论可知，在已知声源与观测点的位置以及相对运动速度的条件下，可对测量信号进行时域上压缩或展宽的波形进行还原，消除测量信号多普勒效应。对测量得到的声压进行多普勒效应去除，得到测量面H上不含多普勒效应的声压p_H(H)为

$p_H\left(H\right)=p_H^{\&prime;}\left(H\right)R{(1-\mathrm{Ma}\cos \mathrm{\&theta;})}^2{e}^{\mathrm{jkR}}\frac{e^{-\mathrm{jkr}}}{r}---\left(4\right)$

式中r为测量面H上网格点与重构面S上网格点之间的距离，r定义为

$r=\sqrt{{(x_j-x_i)}^2+{(y_j-y_i)}^2+d^2}---\left(5\right)$

将测量面H上去除多普勒效应的声压p_H(H)经补零扩展得到测量面H+上声压p_E(H+)。可以建立测量面H+上补零扩展得到的声压p_E(H+)与测量面H上的声压p_H(H)之间的关系为

$p_E(H+)=\left\{\begin{array}{cc}{p}_H\left(H\right)& (x_i,y_i)\&Element;H\\ 0& (x_i,y_i)\&NotElement;H\end{array}\right.---\left(6\right)$

式中p_E(H+)=[p_E1，…，p_EM]^T，p_H(H+)=[p_H1，…，p_HM]^T，M为测量面H+上测量网格的点数。式（6）还可以表达为

p_E(H+)＝D·p_H(H+)   （7）

式中p_H(H+)为测量面H+上的声压，D＝diag[D₁₁,…,D_MM]，D定义为采样因子

$D_{\mathrm{ii}}=\left\{\begin{array}{cc}1& (x_i,y_i)\&Element;H\\ 0& (x_i,y_i)\&NotElement;H\end{array}\right.---\left(8\right)$

由理想流体媒质中小振幅声波的波动方程，可以得到不依赖于时间变量的稳态声场Helmholtz方程：

▽²p(x,y,z)+k²p(x,y,z)＝0  （9）

式中：p(x,y,z)为空间点(x,y,z)的复声压，k＝ω/c＝2π/λ为流体介质中自由场波数，ω为声波的角频率，λ为特征波长。由格林函数公式可以得到式(9)的解，即测量面H+上的声压p_H(H+)与重构面S上声压p_S的关系为

p_H(H+)＝F^-1G_DFp_S   （10）

式中p_S=[p_S1，…，p_SM]^T，F为二维空间傅里叶变换因子，F^-1为二维空间傅里叶变换因子的逆变换，G_D为测量面H+上声压p_H(H+)与重构面上声压的格林函数，式中k_z为

$k_z=\left\{\begin{array}{cc}\sqrt{k^2-k_x^2-k_y^2}& k_x^2+k_y^2\&le;k^2\\ -j\sqrt{k_x^2+k_y^2-k^2}& k_x^2+k_y^2>k^2\end{array}\right.---\left(11\right)$

由于声压数据具有波数域带限性质，可以由波数域内带限算子对全息面上声压进行带限处理，以消除截止波数k_c以外的高波数成分

$p_H(H+)=B_{k_c}\&CenterDot;p_H(H+)$

式中带限算子 $B_{k_c}=F^{-1}{L}_{k_c}F,$ 为低通滤波器， $L_{k_c}=\mathrm{diag}[L_{11},\&CenterDot;\&CenterDot;\&CenterDot;,L_{\mathrm{MM}}].$ 根据辐射源理论，可以表示为

$L_{\mathrm{ii}}=\left\{\begin{array}{cc}1& \sqrt{k_x^2+k_y^2}<k_c\\ 0& \sqrt{k_x^2+k_y^2}>k_c\end{array}\right.---\left(13\right)$

截止波数k_c的选择可以定义为

$k_c=\left\{\begin{array}{cc}{k}_0& k_0\&le;\min (k_{x\max },k_{y\max })\\ \min (k_{x\max },k_{y\max })& k_0>\min (k_{x\max },k_{y\max })\end{array}\right.---\left(14\right)$

式中k_xmax＝π/Δx,k_ymax＝π/Δy,Δx和Δy分别为测量面H上x方向和y方向的网格点间隔。

$k_0=\sqrt{k^2+{\left(\frac{\ln \left(10\mathrm{SNR}\right)}{20(z_H-z_S)}\right)}^2}---\left(15\right)$

式中SNR为信噪比。把式（12）带入到式（7）中，可以得到测量面H+上由p_H(H)补零后得到的声压p_E(H+)与测量面H+上的声压p_H(H+)之间的传递关系

$p_E(H+)=D\&CenterDot;B_{k_c}\&CenterDot;p_H(H+)---\left(16\right)$

再把式（10）带入式（16）可以得到测量面H+上由p_H(H)补零得到的声压p_E(H+)与重构面S上的声压p_S的传递关系

$\begin{array}{c}{p}_E(H+)=D\&CenterDot;B_{k_c}\&CenterDot;p_H(H+)\\ =D{B}_{k_c}{F}^{-1}{G}_{D}F{p}_s=W_D{p}_S\end{array}---\left(17\right)$

式中W_D为测量面H+上由p_H(H)补零后得到的声压p_E(H+)与重构面S上声压之间的传递矩阵。

根据上式建立的传递关系，对式（17）进行求逆，得到重构面S上的声压

$p_S={(\mathrm{\&alpha;I}+W_D^H{W}_D)}^{-1}{W}_D^H{p}_E(H+)---\left(18\right)$

其中I为单位对角矩阵，为矩阵的逆矩阵，为W_D的共轭转置矩阵，α为正则化参数。

通过对式（18）的求解，实现了对测量面H上不含多普勒效应的声压的补零扩展，可以获得重构面S上的声压。

与已有技术相比，本发明的有益效果：

1、本发明采用快速局部近场声全息法为基础，与传统法相比，本发明具有计算简单、计算时间短及计算效率高等优点。

参见图1，本实施例中，在移动声源的辐射声场中，位于声源近场有测量面H，在测量面上呈网格式分布有测量点，相邻网格点之间的距离小于半个波长；移动声源运动方向所在平面为重构面S，移动声源运动距离为L,在重构面上分布有网格点，网格点的分布与测量面H上的测量点相同。

具体实施步骤为：

a、当移动声源运动到位置IIx_j＝L/2时，在测量面H上采用传声器阵列快照法同时采样获得测量面H上的声压信息；

b、对测量得到的声压进行多普勒效应去除，得到测量面H上不含多普勒效应的声压

$p_H\left(H\right)=p_H^{\&prime;}\left(H\right)R{(1-\mathrm{Ma}\cos \mathrm{\&theta;})}^2{e}^{\mathrm{jkR}}\frac{e^{-\mathrm{jkr}}}{r},$ 其中

$R=\frac{\mathrm{Ma}\&CenterDot;[x_j-x_i]}{1-{\mathrm{Ma}}^2}+\frac{\sqrt{{(x_j-x_i)}^2+(1-{\mathrm{Ma}}^2)[{(y_j-y_i)}^2+d^2]}}{1-{\mathrm{Ma}}^2},$

cosθ＝(x_j-x_i)/R、

$r=\sqrt{{(x_j-x_i)}^2+{(y_j-y_i)}^2+d^2},$

Ma=v/c；

p_H(H)为测量面H上不含多普勒效应的声压、

p'_H(H)为测量面H上测量得到含有多普勒效应的声压、

r为测量面H上网格点与重构面S上网格点之间的距离、

cosθ为测量面H上网格点和重构面S上网格点之间的距离与声源运动方向夹角的余弦、

R为声波从移动声源位置到测量面H上网格点的实际传播距离、

c为介质中声速、v为移动声源的运动速度、

Ma为移动声源的运动马赫数、k为波数、

d为测量面H与重构面S之间的距离、

(x_i,y_i)为测量面H上网格点的坐标、(x_j,y_j)为重构面S上网格点的坐标，测量面H与重构面S平行，测量面H和重构面S的法向均为z方向；

c、测量面H上的声压补零扩展得到测量面H+上声压

$p_E(H+)=\left\{\begin{array}{cc}{p}_H\left(H\right)& (x_i,y_i)\&Element;H\\ 0& (x_i,y_i)\&NotElement;H\end{array}\right.$

p_E(H+)＝D·p_H(H+)，其中

D＝diag[D₁₁,…,D_MM]

$D_{\mathrm{ii}}=\left\{\begin{array}{cc}1& (x_i,y_i)\&Element;H\\ 0& (x_i,y_i)\&NotElement;H\end{array}\right.$

p_E(H+)=[p_E1，…，p_EM]^T

p_H(H+)=[p_H1，…，p_HM]^T

p_E(H+)为测量面H+上由p_H(H)补零后得到的声压、

p_H(H+)为测量面H+上的声压、

D为采样算子、D_ii为矩阵D中对角线上的值、

M为测量面H+上测量网格的点数、T为矩阵转置；

d、计算测量面H+上的声压与重构面S上声压之间的传递矩阵

$W_D={\mathrm{DB}}_{k_c}{F}^{-1}{G}_{D}F$ 其中

$G_D=g^{j{k}_{z}d},$

$B_{k_c}=F^{-1}{K}_{k_c}F,$

$k_z=\left\{\begin{array}{cc}\sqrt{k^2-k_x^2-k_y^2}& k_x^2+k_y^2\&le;k^2\\ -j\sqrt{k_x^2+k_y^2-k^2}& k_x^2+k_y^2>k^2\end{array}\right.,$

$L_{k_c}=\mathrm{diag}[L_{11},\&CenterDot;\&CenterDot;\&CenterDot;,L_{\mathrm{MM}}],L_{\mathrm{ii}}=\left\{\begin{array}{cc}1& \sqrt{k_x^2+k_y^2}<k_c\\ 0& \sqrt{k_x^2+k_y^2}>k_c\end{array}\right.;$

W_D为测量面H+上由p_H(H)补零后得到的声压p_E(H+)与重构面S上声压的传递矩阵、

G_D为测量面H+上声压p_H(H+)与重构面S上声压的格林函数、

F为二维空间傅里叶变换因子、F^-1为二维空间傅里叶变换因子的逆变换、

为波数域内带限算子、k_c为截止波数、

为低通滤波器、L_ii为矩阵中对角线上的值；

e、求解得到重构面S上的声压

$p_S={(\mathrm{\&alpha;I}+W_D^H{W}_D)}^{-1}{W}_D^H{p}_E(H+),$ 其中

p_S为重构面S上的声压、

α为正则化参数、I为单位对角矩阵、

为W_D的共轭转置矩阵、为矩阵的逆矩阵。

以移动点声源为例进行验证：

声源表面声压和辐射声场由移动点声源产生，采用本发明方法分别实现声源辐射声场的重建，并与其真实值比较。

计算中水中声速为1500m/s，水的密度为1000kg/m³，测量面与重构面的距离为0.1m，测量面H的孔径为8m×4m，测量点数为16×8，对测量面H进行补零扩展，扩展测量面H+的孔径为16m×8m，测量点数为32×16，重构面大小为8m×4m。

图2为频率为500Hz，声源运动速度为3m/s时，全息面H上测量的声压分布；图3为频率为500Hz，声源运动速度为3m/s时，全息面H上测量的声压分布与声源静置时相同频率下的声压分布比较。可以很明显看出由于多普勒效应的影响两者存在不同。

以图2所示的测量声压为输入量，按照常规局部近场声全息法重建的重构面S上声压分布如图5（a）所示，按照本发明方法重建的重构面S上声压分布如图5（b）所示。与图4中重构面S上真实声压分布相比，在声压分布上常规局部近场声全息法重建失效，显示出了现有近场声全息技术的局限性；而采用本发明方法的重建结果与真实值得吻合非常好，成功消除了多普勒效应的影响。

改变移动声源的运动速度，由3m/s增加到10m/s，按照常规局部近场声全息法重建的重构面S上声压分布如图6（a）所示，按照本发明方法重建的重构面S上声压分布如图6（b）所示。与图4中重构面S上真实声压分布相比，采用本发明方法的重建结果与真实值得吻合非常好，成功消除了多普勒效应的影响。

对于图5和图6中的计算结果，常规局部近场声全息法的移动声源识别方法需要8分钟才能得到计算结果，而采用本发明方法进行重构时计算时间为15秒，可以很明显看出本发明方法计算简单，能够以更短时间进行声场重建。

为了定量的分析重建精度，下面对本发明方法的重构误差进行计算，定义重构均方误差率为

$\mathrm{MSE}=\sqrt{\frac{\mathrm{\&Sigma;}{|p_t-p_s|}^2}{\mathrm{\&Sigma;}{\left|p_t\right|}^2}}\&times;100\%---\left(19\right)$

式中，p_t为重构面S上的声压理论值，p_s为本发明方法得到的重构面S上声压的重构值。由式（19）计算可得，采用本发明方法得到的重构误差为6.2%，可见本发明方法能取得很好的重建效果。

图7所示为声源频率在0至2500Hz的频带内，采用常规局部近场声全息法和本发明方法得到的重构结果误差。可以看到随着频率的增大，常规局部近场声全息法的重建误差也随之增大，重建结果更加杂乱无章，完全失效，而本发明方法能够获得精度更高的重构结果，与真实声压值吻合的很好。

Claims (1)

1.一种局部近场声全息法的移动噪声源识别方法，其特征在于：

（1）测量测量面H上声压：

位于声源近场有测量面H，在测量面上呈网格式分布有测量点，相邻网格点之间的距离小于半个波长；移动声源运动方向所在平面为重构面S，移动声源运动距离为L，在重构面S上分布有网格点，网格点的分布与测量面H上的测量点相同；测量所述测量面H上各网格点处的声压幅值和相位信息获得测量面H上含有多普勒效应的声压；

（2）对测量得到的声压进行多普勒效应去除，得到测量面H上不含多普勒效应的声压：

$p_H\left(H\right)=p_H^{\&prime;}\left(H\right)R{(1-\mathrm{Ma}\cos \mathrm{\&theta;})}^2{e}^{\mathrm{jkR}}\frac{e^{-\mathrm{jkr}}}{r},$ 其中

$R=\frac{\mathrm{Ma}\&CenterDot;[x_j-x_i]}{1-{\mathrm{Ma}}^2}+\frac{\sqrt{{(x_j-x_i)}^2+(1-{\mathrm{Ma}}^2)[{(y_j-y_i)}^2+d^2]}}{1-{\mathrm{Ma}}^2},$ cosθ＝(x_j-x_i)/R、

$r=\sqrt{{(x_j-x_i)}^2+{(y_j-y_i)}^2+d^2},$ Ma=v/c；

p_H(H)为测量面H上不含多普勒效应的声压、p'_H(H)为测量面H上测量得到含有多普勒效应的声压、r为测量面H上网格点与重构面S上网格点之间的距离、cosθ为测量面H上网格点和重构面S上网格点之间的距离与声源运动方向夹角的余弦、R为声波从移动声源位置到测量面H上网格点的实际传播距离、c为介质中声速、v为移动声源的运动速度、Ma为移动声源的运动马赫数、k为波数、d为测量面H与重构面S之间的距离、(x_i,y_i)为测量面H上网格点的坐标、(x_j,y_j)为重构面S上网格点的坐标，测量面H与重构面S平行，测量面H和重构面S的法向均为z方向；

（3）测量面H上的声压补零扩展得到测量面H+上声压：

$p_E(H+)=\left\{\begin{array}{cc}{p}_H\left(H\right)& (x_i,y_i)\&Element;H\\ 0& (x_i,y_i)\&NotElement;H\end{array}\right.,$ p_E(H+)＝D·p_H(H+)，

其中D＝diag[D₁₁,…,D_MM]， $D_{\mathrm{ii}}=\left\{\begin{array}{cc}1& (x_i,y_i)\&Element;H\\ 0& (x_i,y_i)\&NotElement;H\end{array}\right.,$ p_E(H+)=[p_E1，…，p_EM]^T，

p_H(H+)=[p_H1，…，p_HM]^T，p_E(H+)为测量面H+上由p_H(H)补零后得到的声压、p_H(H+)为测量面H+上的声压、D为采样算子、D_ii为矩阵D中对角线上的值、M为测量面H+上测量网格的点数、T为矩阵转置；

（4）计算测量面H+上的声压与重构面S上声压之间的传递矩阵：

$W_D={\mathrm{DB}}_{k_c}{F}^{-1}{G}_{D}F$

其中， $G_D=g^{j{k}_{z}d},B_{k_c}=F^{-1}{K}_{k_c}F,k_z=\left\{\begin{array}{cc}\sqrt{k^2-k_x^2-k_y^2}& k_x^2+k_y^2\&le;k^2\\ -j\sqrt{k_x^2+k_y^2-k^2}& k_x^2+k_y^2>k^2\end{array}\right.,$ $L_{k_c}=\mathrm{diag}[L_{11},\&CenterDot;\&CenterDot;\&CenterDot;,L_{\mathrm{MM}}],L_{\mathrm{ii}}=\left\{\begin{array}{cc}1& \sqrt{k_x^2+k_y^2}<k_c\\ 0& \sqrt{k_x^2+k_y^2}>k_c\end{array}\right.;$

W_D为测量面H+上由p_H(H)补零后得到的声压pE(H+)与重构面S上声压的传递矩阵、G_D为测量面H+上声压p_H(H+)与重构面S上声压的格林函数、F为二维空间傅里叶变换因子、F^-1为二维空间傅里叶变换因子的逆变换、为波数域内带限算子、k_c为截止波数、为低通滤波器、L_ii为矩阵中对角线上的值；

（5）求解得到重构面S上的声压：

$p_S={(\mathrm{\&alpha;I}+W_D^H{W}_D)}^{-1}{W}_D^H{p}_E(H+)$

其中p_S为重构面S上的声压、α为正则化参数、I为单位对角矩阵、为W_D的共轭转置矩阵、为矩阵的逆矩阵。