CN103592022A - 采用声压和质点加速度测量的实时声场分离方法 - Google Patents

Abstract

本发明公开了一种采用声压和质点加速度测量的实时声场分离方法，其特征是在目标声源和干扰声源之间布置两个相互平行的测量平面H和测量平面H₁，同步采集两个测量平面上的声压时域信号；利用两个测量平面上的声压时域信号，通过有限差分法获得测量平面H上的质点加速度时域信号；再利用测量平面H上的声压时域波数谱和质点加速度时域波数谱、已知的时域脉冲响应函数，实时分离出目标声源单独在测量平面H上所辐射的声压时域波数谱，进而获得目标声源单独在测量平面H上所辐射的声压时域信号。本发明方法无需进行任何求逆和正则化运算处理，计算速度快、稳定性高；本发明方法具备实时分离声场的能力，可用于噪声干扰环境下现场分析目标声源的时变辐射特性。

Description

采用声压和质点加速度测量的实时声场分离方法

技术领域

本发明涉及物理专业中噪声类领域声场分离方法。

背景技术

在实际工程中，目标声源位于测量平面的一侧，而在测量平面的另一侧往往存在干扰声源，这些干扰声源所产生的声场影响了对目标声源所辐射声场的准确测量，因此需要采用一定的声场分离方法将干扰声源的影响从测量结果中分离出来。到目前为止，国内外学者已提出多种声场分离方法，这些方法可大致分为五类：一是基于空间傅里叶变换法（SFT）的声场分离技术：G.V.Frisk等在1980年首次提出采用SFT法间接测量海洋底面的反射系数，并建立了基于SFT法的双面声场分离理论；M.Tamura于1990年将G.V.Frisk等提出的方法用于测量斜入射时材料的反射系数；M.T.Cheng等对M.Tamura提出的方法进行了推广，建立了柱面坐标下的双测量声场分离式，并用于实现散射声场的分离。二是基于球面波叠加的声场分离方法：1956年，J.Pachner采用球面波叠加法实现了任意波场中行波和驻波声场的分离；G.Weinreich等在1980年对J.Pachner提出的方法作了进一步改进，建立了基于双球面测量的声场分离理论；Y.Braikia等在2013年分别采用双球面声压测量的分离方法和单球面声压振速测量的分离方法实现了小空间范围内的声源识别。三是基于统计最优近场声全息（SONAH）的声场分离技术：J.Hald在对SONAH研究的基础上，提出了基于双全息面声压测量的统计最优声场分离技术；F.Jacobsen等在J.Hald提出的方法的基础上，提出了基于声压振速测量的统计最优声场分离技术；E.Fernandez-Grande等在2011年提出了基于双面振速测量的统计最优声场分离技术。四是基于边界元法（BEM）的声场分离技术：C.Langrenne等在2007年提出一种基于BEM的双面声场分离方法；随后，E.G.Williams等在2008年提出一种基于BEM和声压振速测量(Cauchy数据)的声场分离方法。五是基于等效源法（ESM）的声场分离技术：C.X.Bi提出的基于ESM的声场分离技术，适用于任意形状的测量面，且计算稳定性好、计算精度高；E.Fernandez-Grande等在2012年提出了基于ESM的双面振速测量的分离技术和单面声压振速测量的分离技术。但是，上述声场分离方法都仅仅实现了单一频率或频带下的声场分离。若要实现任意时刻下的声场分离，上述声场分离方法将不再适用。为实现任意时刻下的声场分离，X.Z.Zhang等在2012年提出一种基于双平面声压测量的时域声场分离方法。在该分离方法中，每个时刻目标声源所辐射的声压信号需要通过解卷运算得到，由于该解卷运算过程具有不适定性，需要采用奇异值分解和正则化来稳定求解过程。而正则化的采用使得该分离方法耗时长，且分离精度很大程度上依赖于正则化参数的合理选取。

发明内容

为实现任意时刻下的声场分离，本发明提供一种计算速度快、稳定性高的采用声压和质点加速度测量的实时声场分离方法，以期获得实时分离声场的能力，进而用于噪声干扰环境下现场分析目标声源的时变辐射特性。

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

本发明采用声压和质点加速度测量的实时声场分离方法的特点是按如下步骤进行：

步骤a、在目标声源M_o和干扰声源M_d之间布置相互平行的测量平面H和测量平面H₁；所述目标声源M_o和干扰声源M_d均可辐射任意线性声场信号；在所述测量平面H和测量平面H₁上分别均匀分布有M个测量网格点，测量平面H和测量平面H₁上的测量网格大小和测量网格点位置相同；以目标声源M_o的中心为坐标原点建立笛卡尔坐标系，所述笛卡尔坐标系中的xoy坐标平面平行于测量平面H，z轴方向垂直于测量平面H；测量平面H上的测量网格点坐标为(x,y,z_H)，测量平面H₁上的测量网格点坐标为(x,y,z_H1)；同步采集测量平面H上各测量网格点(x,y,z_H)处的声压时域信号p(x,y,z_H,t)和测量平面H₁上各测量网格点(x,y,z_H1)处的声压时域信号p(x,y,z_H1,t)；

步骤b、按照式(1)进行有限差分获得测量平面H上各测量网格点(x,y,z_H)处的质点加速度时域信号a(x,y,z_H,t)：

$a(x,y,z_H,t)=-\frac{p(x,y,z_{H1},t)-p(x,y,z_H,t)}{{\mathrm{\ρ}}_0(z_{H1}-z_H)}---\left(1\right)$

式(1)中，t为时间；ρ₀为介质密度；

步骤c、对于测量平面H上的声压时域信号p(x,y,z_H,t)按照式(2)进行二维空间傅里叶变换获得声压时域波数谱P(k_x,k_y,z_H,t)，对测量平面H上的质点加速度时域信号a(x,y,z_H,t)按照式(3)进行二维空间傅里叶变换获得质点加速度时域波数谱A(k_x,k_y,z_H,t)，

$P(k_x,k_y,z_H,t)={\∫}_{-\∞}^{\∞}{\∫}_{-\∞}^{\∞}p(x,y,z_H,t)e^{j(k_{x}x+k_{y}y)}\mathrm{dxdy}---\left(2\right)$

$A(k_x,k_y,z_H,t)={\∫}_{-\∞}^{\∞}{\∫}_{-\∞}^{\∞}a(x,y,z_H,t)e^{j(k_{x}x+k_{y}y)}\mathrm{dxdy}---\left(3\right)$

在式(2)和式(3)中，j表示虚数单位；k_x、k_y分别为x、y方向的波数；

步骤d、构建声压时域波数谱P(k_x,k_y,z_H,t)、质点加速度时域波数谱A(k_x,k_y,z_H,t)、已知的时域脉冲响应函数h(k_x,k_y,0,t)和目标声源M_o单独在测量平面H上所辐射的声压时域波数谱P_o(k_x,k_y,z_H,t)之间的关系如式(4)所示：

P_o(k_x,k_y,z_H,t)=0.5[P(k_x,k_y,z_H,t)+A(k_x,k_y,z_H,t)*h(k_x,k_y,0,t)]  (4)

在式(4)中，“*”表示两个时间函数的卷积运算；

将式(4)中的时间t离散为t_n=(n-1)Δt，其中Δt为采样时间间隔，n=1,2,...,N，N为采样点总数，式(4)的离散形式如式(5)所示：

$P_o(k_x,k_y,z_H,t_n)=0.5[P(k_x,k_y,z_H,t_n)+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}A(k_x,k_y,z_H,t_i)h(k_x,k_y,0,t_{n-i+1})]---\left(5\right)$

步骤e、对于所述声压时域波数谱P_o(k_x,k_y,z_H,t_n)按照式(6)进行二维空间傅里叶反变换，分离出每个时刻下目标声源M_o单独在测量平面H上所辐射的声压时域信号p_o(x,y,z_H,t_n)：

$p_o(x,y,z_H,t_n)=\frac{1}{{\left(2\mathrm{\π}\right)}^2\mathrm{}}{\∫}_{-\∞}^{\∞}{\∫}_{-\∞}^{\∞}{P}_o(k_x,k_y,z_H,t_n)e^{-j(k_{x}x+k_{y}y)}{\mathrm{dk}}_x{\mathrm{dk}}_y---\left(6\right).$

本发明采用声压和质点加速度测量的实时声场分离方法的特点也在于：

所述测量平面H和测量平面H₁上各测量网格点的声压时域信号p(x,y,z_H,t)和p(x,y,z_H1,t)是采用声压传感器阵列一次快照测量获得。

所述干扰声源M_d为噪声源、反射源或散射源。

理论模型：

在辐射任意线性声场信号的目标声源M_o和辐射任意线性声场信号的干扰声源M_d之间分别布置两个相互平行的测量平面H和测量平面H₁。根据Euler公式，测量平面H上的质点加速度时域信号a(x,y,z_H,t)与测量平面H上所测得的声压时域信号p(x,y,z_H,t)的偏导数存在成正比的关系，如式(7)所示：

$a(x,y,z_H,t)=-\frac{1}{{\mathrm{\ρ}}_0}\frac{\∂p(x,y,z_H,t)}{\∂z}---\left(7\right)$

在式(7)中，ρ₀为介质密度。

利用测量平面H上的声压时域信号p(x,y,z_H,t)和测量平面H₁上的声压时域信号p(x,y,z_H1,t)，按照式(8)进行有限差分获得测量平面H上的声压时域信号p(x,y,z_H,t)的偏导数：

$\frac{\∂p(x,y,z_H,t)}{\∂z}=\frac{p(x,y,z_H,t)-p(x,y,z_{H1},t)}{z_H-z_{H1}}---\left(8\right)$

联合式(7)和(8)获得测量平面H上的质点加速度时域信号a(x,y,z_H,t)，即：

$a(x,y,z_H,t)=-\frac{p(x,y,z_H,t)-p(x,y,z_{H1},t)}{{\mathrm{\ρ}}_0(z_H-z_{H1})}---\left(9\right)$

根据声波的叠加原理，在测量平面H上所测得的声压时域信号p(x,y,z_H,t)等于目标声源M_o在测量平面H上所辐射的声压时域信号p_o(x,y,z_H,t)与干扰声源M_d在测量平面H所辐射的声压时域信号p_d(x,y,z_H,t)之和，如式(10)所示：

p(x,y,z_H,t)=p_o(x,y,z_H,t)+p_d(x,y,z_H,t)  (10)

由于质点加速度是矢量，在测量平面H上所获得的质点加速度时域信号a(x,y,z_H,t)应等于目标声源M_o在测量平面H上所辐射的质点加速度时域信号a_o(x,y,z_H,t)与干扰声源M_d在测量平面H所辐射的质点加速度时域信号a_d(x,y,z_H,t)之差，如式(11)所示：

a(x,y,z_H,t)=a_o(x,y,z_H,t)-a_d(x,y,z_H,t)  (11)

定义函数f(x,y,z,t)关于x、y的二维空间傅里叶变换为：

$F(k_x,k_y,z,t)={\∫}_{-\∞}^{\∞}{\∫}_{-\∞}^{\∞}f(x,y,z,t)e^{j(k_{x}x+k_{y}y)}\mathrm{dxdy}---\left(12\right)$

式(12)中，F(k_x,k_y,z,t)表示f(x,y,z,t)的时域波数谱。

分别对式(10)和式(11)进行关于x、y的二维空间傅里叶变换可得：

P(k_x,k_y,z_H,t)=P_o(k_x,k_y,z_H,t)+P_d(k_x,k_y,z_H,t)  (13)

A(k_x,k_y,z_H,t)=A_o(k_x,k_y,z_H,t)-A_d(k_x,k_y,z_H,t)  (14)

由声压和质点加速度在时域波数域的关系可知：

P_o(k_x,k_y,z_H,t)=A_o(k_x,k_y,z_H,t)*h(k_x,k_y,0,t)(15)

P_d(k_x,k_y,z_H,t)=A_d(k_x,k_y,z_H,t)*h(k_x,k_y,0,t)  (16)

式(15)和式(16)中，h(k_x,k_y,0,t)为声压和质点加速度之间的时域脉冲响应函数，其表达式为：

$h(k_x,k_y,0,t)={\mathrm{\ρ}}_0{\mathrm{cJ}}_0\left(\sqrt{c^2(k_x^2+k_y^2)t}\right)H\left(t\right)---\left(17\right)$

在式(17)中，c为声速，J₀为阶数为0的第一类贝塞尔函数，H(t)为Heaviside函数。

将式(16)代入式(13)中得到：

P(k_x,k_y,z_H,t)=P_o(k_x,k_y,z_H,t)+A_d(k_x,k_y,z_H,t)*h(k_x,k_y,0,t)  (18)

将式(14)两边分别与h(k_x,k_y,0,t)做卷积得到：

A(k_x,k_y,z_H,t)*h(k_x,k_y,0,t)=A_o(k_x,k_y,z_H,t)*h(k_x,k_y,0,t)-A_d(k_x,k_y,z_H,t)*h(k_x,k_y,0,t)  (19)

将式(15)代入式(19)中得到：

A(k_x,k_y,z_H,t)*h(k_x,k_y,0,t)=P_o(k_x,k_y,z_H,t)-A_d(k_x,k_y,z_H,t)*h(k_x,k_y,0,t)  (20)

将式(18)与式(20)相加可进一步得到：

P_o(k_x,k_y,z_H,t)=0.5[P(k_x,k_y,z_H,t)+A(k_x,k_y,z_H,t)*h(k_x,k_y,0,t)]  (21)

式(21)中，P(k_x,k_y,z_H,t)和A(k_x,k_y,z_H,t)可分别通过对测得的声压时域信号p(x,y,z_H,t)和通过有限差分获得的质点加速度时域信号a(x,y,z_H,t)进行二维空间傅里叶变换得到，h(k_x,k_y,0,t)为已知函数。利用式(21)按如下过程可实时分离出每个时刻下目标声源M_o单独在测量平面H上所辐射的声压时域波数谱P_o(k_x,k_y,z_H,t)：

将式(21)中的时间t离散为t_n=(n-1)Δt，其中Δt为采样时间间隔，n=1,2,...,N，N为采样点总数。

取n=1时，分离得到：

$P_o(k_x,k_y,z_H,t_1)=0.5[P(k_x,k_y,z_H,t_1)+\underset{i=1}{\overset{1}{\mathrm{\Σ}}}A(k_x,k_y,z_H,t_1)h(k_x,k_y,0,t_{1-i+1})]---\left(22\right)$

取n=2时，分离得到：

$P_o(k_x,k_y,z_H,t_2)=0.5[P(k_x,k_y,z_H,t_2)+\underset{i=1}{\overset{2}{\mathrm{\Σ}}}A(k_x,k_y,z_H,t_i)h(k_x,k_y,0,t_{2-i+1})]---\left(23\right)$

……

取n=N时，分离得到：

$P_o(k_x,k_y,z_H,t_N)=0.5[P(k_x,k_y,z_H,t_N)+\underset{i=1}{\overset{N}{\mathrm{\Σ}}}A(k_x,k_y,z_H,t_i)h(k_x,k_y,0,t_{N-i+1})]---\left(24\right)$

由式(22)到式(24)的分离过程可知，要分离出t_n时刻的P_o(k_x,k_y,z_H,t_n)，只需要利用t_n时刻的P(k_x,k_y,z_H,t_n)和t_i(i=1,2,...,n)时刻的A(k_x,k_y,z_H,t_i)，而P(k_x,k_y,z_H,t_n)和A(k_x,k_y,z_H,t_i)可通过对测得的声压时域信号p(x,y,z_H,t)和通过有限差分获得的质点加速度时域信号a(x,y,z_H,t)进行二维空间傅里叶变换得到。因此，一旦测得某个时刻下测量平面H和测量平面H₁上的声压时域信号，即可分离出该时刻的声压时域波数谱，实现声场的实时分离。

对分离出的时域波数谱P_o(k_x,k_y,z_H,t_n)按照式(25)进行二维空间傅里叶反变换，最终获得每个时刻下目标声源M_o单独在测量平面H上所辐射的声压时域信号p_o(x,y,z_H,t_n)

$P_o(x,y,z_H,t_n)=\frac{1}{{\left(2\mathrm{\π}\right)}^2}{\∫}_{-\∞}^{\∞}{\∫}_{-\∞}^{\∞}{P}_o(k_x,k_y,z_H,t_n)e^{-j(k_{x}x+k_{y}y)}d{k}_{x}d{k}_y---\left(25\right)$

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

1、本发明方法的分离过程是一个简单的正向求解过程，不需要进行任何求逆和正则化运算处理，因而计算速度快、稳定性高。

2、本发明方法具备实时分离声场的能力，因而更适用于噪声干扰环境下现场分析目标声源的时变辐射特性。

附图说明

图1为本发明采用声压和质点加速度测量的实时声场分离方法示意图；

图2(a)为本发明方法的实验声源、测量平面位置分布图；

图2(b)为本发明方法的实验测量平面网格划分示意图；

图3(a)为本发明方法的实验测量点A处声压时域信号图；

图3(b)为本发明方法的实验测量点B处声压时域信号图；

图3(c)为本发明方法的实验测量点C处声压时域信号图；

图4(a)为本发明方法实验所获得的相位评价因子Ep的数值分布图；

图4(b)为本发明方法实验所获得的幅值评价因子Ea的数值分布图；

图5(a)为本发明方法t=9.38ms时刻目标声源M_o单独在测量平面H上所辐射的测量声压场p_m；

图5(b)为本发明方法t=10.78ms时刻目标声源M_o单独在测量平面H上所辐射的测量声压场p_m；

图5(c)为本发明方法t=9.38ms时刻的混合声压场p_c；

图5(d)为本发明方法t=10.78ms时刻的混合声压场p_c；

图5(e)为t=9.38ms时刻采用本发明方法分离出的声压场p_e；

图5(f)为t=10.78ms时刻采用本发明方法分离出的声压场p_e。

具体实施方式

参见图1，本实施例采用声压和质点加速度测量的实时声场分离方法按如下步骤分离出目标声源M_o单独在测量平面H上所辐射的声压时域信号：

步骤a、在目标声源M_o和干扰声源M_d之间布置相互平行的测量平面H和测量平面H₁；所述目标声源M_o和干扰声源M_d均可辐射任意线性声场信号；在所述测量平面H和测量平面H₁上分别均匀分布有M个测量网格点，测量平面H和测量平面H₁上的测量网格大小和测量网格点位置相同；以目标声源M_o的中心为坐标原点建立笛卡尔坐标系，所述笛卡尔坐标系中的xoy坐标平面平行于测量平面H，z轴方向垂直于测量平面H；测量平面H上的测量网格点坐标为(x,y,z_H)，测量平面H₁上的测量网格点坐标为(x,y,z_H1)；同步采集测量平面H上各测量网格点(x,y,z_H)处的声压时域信号p(x,y,z_H,t)和测量平面H₁上各测量网格点(x,y,z_H1)处的声压时域信号p(x,y,z_H1,t)。

步骤b、按照式(26)进行有限差分获得测量平面H上各测量网格点(x,y,z_H)处的质点加速度时域信号a(x,y,z_H,t)：

$a(x,y,z_H,t)=-\frac{p(x,y,z_{H1},t)-p(x,y,z_H,t)}{{\mathrm{\ρ}}_0(z_{H1}-z_H)}---\left(26\right)$

式(26)中，t为时间；ρ₀为介质密度。

步骤c、对于测量平面H上的声压时域信号p(x,y,z_H,t)按照式(27)进行二维空间傅里叶变换获得声压时域波数谱P(k_x,k_y,z_H,t)，对测量平面H上的质点加速度时域信号a(x,y,z_H,t)按照式(28)进行二维空间傅里叶变换获得质点加速度时域波数谱A(k_x,k_y,z_H,t)，

$P(k_x,k_y,z_H,t)={\∫}_{-\∞}^{\∞}{\∫}_{-\∞}^{\∞}p(x,y,z_H,t)e^{j(k_{x}x+k_{y}y)}\mathrm{dxdy}---\left(27\right)$

$A(k_x,k_y,z_H,t)={\∫}_{-\∞}^{\∞}{\∫}_{-\∞}^{\∞}a(x,y,z_H,t)e^{j(k_{x}x+k_{y}y)}\mathrm{dxdy}---\left(28\right)$

在式(27)和式(28)中，j表示虚数单位；k_x、k_y分别为x、y方向的波数。

步骤d、构建声压时域波数谱P(k_x,k_y,z_H,t)、质点加速度时域波数谱A(k_x,k_y,z_H,t)、已知的时域脉冲响应函数h(k_x,k_y,0,t)和目标声源M_o单独在测量平面H上所辐射的声压时域波数谱P_o(k_x,k_y,z_H,t)之间的关系如式(29)所示：

P_o(k_x,k_y,z_H,t)=0.5[P(k_x,k_y,z_H,t)+A(k_x,k_y,z_H,t)*h(k_x,k_y,0,t)]  (29)

在式(29)中，“*”表示两个时间函数的卷积运算。

将式(29)中的时间t离散为t_n=(n-1)Δt，其中Δt为采样时间间隔，n=1,2,...,N，N为采样点总数。式(29)的离散形式如式(30)所示

$P_o(k_x,k_y,z_H,t_n)=0.5[P(k_x,k_y,z_H,t_n)+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}A(k_x,k_y,z_H,t_i)h(k_x,k_y,0,t_{n-i+1})]---\left(30\right)$

步骤e、对于所述声压时域波数谱P_o(k_x,k_y,z_H,t_n)按照式(31)进行二维空间傅里叶反变换，分离出每个时刻下目标声源M_o单独在测量平面H上所辐射的声压时域信号p_o(x,y,z_H,t_n)

$p_o(x,y,z_H,t_n)=\frac{1}{{\left(2\mathrm{\π}\right)}^2}{\∫}_{-\∞}^{\∞}{\∫}_{-\∞}^{\∞}{P}_o(k_x,k_y,z_H,t_n)e^{-j(k_{x}x+k_{y}y)}{\mathrm{dk}}_x{\mathrm{dk}}_y---\left(31\right)$

测量平面H和测量平面H₁上各测量网格点的声压时域信号p(x,y,z_H,t)和p(x,y,z_H1,t)是采用声压传感器阵列一次快照测量获得。干扰声源M_d为噪声源、反射源或散射源。

方法的检验：

在测量平面H和测量平面H1的一侧布置音箱S1和音箱S3作为目标声源M_o，在另外一侧布置音箱S2作为干扰声源M_d。采用本发明方法将测量平面H上目标声源M_o单独所辐射的声压时域信号分离出来，并与其测量声压时域信号进行比较。

本实验中，目标声源M_o（包含音箱S1和音箱S3）和干扰声源M_d（音箱S2）都辐射余弦调制信号，其表达式为：

$s\left(t\right)=\cos \left(2\mathrm{\π}{f}_{0}t\right)e^{[-{10}^6{(t-0.005)}^2/2]}---\left(32\right)$

在式(32)中，频率f₀=800Hz。

测量平面H和H1与声源之间的位置关系参见图2(a)，测量平面H位于z=0m的平面上，测量平面H₁位于z=-0.02m的平面上。音箱S1和音箱S3的纸盆中心分别位于（0.1m,0.45m,0.145m）和（0.55m,0.35m,0.14m）处，音箱S2的纸盆中心位于(0.4m,0.45m,-0.14m)处。测量平面H和H1的面积大小都为0.7m×0.7m，其上都均匀地分布15×15个测量点，参见图2(b)。时域信号采样频率为25600Hz，采样点数为512。

为检验本发明方法在时域内的声场分离效果，在测量平面H上选取了三个测量点，即测量点A、测量点B和测量点C，其位置分别为A(0.1m,0.45m,0m),B(0.4m,0.45m,0m)和C(0.55m,0.35m,0m)，其中测量点A和测量点C分别正对着目标声源M_o的音箱S1和音箱S3，测量点B正对着干扰声源M_d的音箱S2。图3(a)、图3(b)和图3(c)分别对应测量点A、测量点B和测量点C，图中曲线a表示目标声源M_o单独所辐射的测量声压时域信号，图中曲线b表示包含目标声源M_o和干扰声源M_d所辐射的混合声压时域信号，图中曲线c表示采用本发明方法分离出的声压时域信号，比较图中的曲线a和曲线b可以看出，干扰声源M_d对目标声源M_o在测量平面H上所辐射的声压时域信号造成了较大的干扰；比较曲线a和曲线c可以看出，采用本发明方法可以很好地消除干扰声源M_d的影响，从而分离出目标声源M_o单独在测量平面H上所辐射的声压时域信号。

为了更加客观地评价本发明方法的分离效果，在此定义了两个评价因子，它们的表达式分别为：

$\mathrm{Ep}(x,y,z)=\frac{<p_m(x,y,z,t)p_e(x,y,z,t)>}{\sqrt{<p_m^2(x,y,z,t)><p_e^2(x,y,z,t)>}}---\left(33\right)$

$\mathrm{Ea}(x,y,z)=\frac{|\sqrt{<p_m^2(x,y,z,t)>}-\sqrt{<p_e^2(x,y,z,t)>}|}{\sqrt{<p_m^2(x,y,z,t)>}}---\left(34\right)$

在式(33)和式(34)中，<>表示求平均值，下标“m”表示目标声源单独辐射的测量声压值，下标“e”表示分离声压值。评价因子Ep是用来衡量目标声源单独辐射的测量声压值和分离声压值之间的相位误差，当Ep的值越靠近1时，相位误差越小。评价因子Ea是用来衡量目标声源单独辐射的测量声压值和分离声压值之间的幅值误差，当Ea的值越靠近0时，幅值误差越小。运用式(33)和(34)分别计算了测量平面H上各个测量点处的Ep和Ea。图4(a)表示所有测量点的Ep值，图中等高线值为0.9。图4(b)表示所有测量点的Ea值，图中等高线值为0.1。参见图4(a)和图4(b)，在大多数测量点处，不论是相位还是幅值，测量声压值和分离声压值都吻合地较好，只是在测量平面边缘处吻合程度略低。

为检验本发明方法在空间域的声场分离效果，选取了两个时刻t=9.38ms和t=10.78ms。图5(a)和图5(b)分别为9.38ms和10.78ms时刻目标声源M_o单独在测量平面H上所辐射的测量声压场p_m，图5(c)和图5(d)分别为9.38ms和10.78ms时刻的混合声压场p_c（包含目标声源M_o和干扰声源M_d所辐射的声压场），图5(e)和图5(f)分别为9.38ms和10.78ms时刻采用本发明方法分离出的声压场p_e。比较图5(a)和图5(c)、图5(b)和图5(d)可以看出，干扰声源M_d对目标声源M_o单独在测量平面H上所辐射的声压场造成了严重的干扰；比较图5(a)和图5(e)、图5(b)和图5(f)可以看出，采用本发明方法可以较好地消除干扰声源M_d的影响，从而获得目标声源M_o单独在测量平面H上所辐射的声压场。

Claims (3)

1.采用声压和质点加速度测量的实时声场分离方法，其特征是按如下步骤进行：

步骤a、在目标声源M_o和干扰声源M_d之间布置相互平行的测量平面H和测量平面H₁；所述目标声源M_o和干扰声源M_d均可辐射任意线性声场信号；在所述测量平面H和测量平面H₁上分别均匀分布有M个测量网格点，测量平面H和测量平面H₁上的测量网格大小和测量网格点位置相同；以目标声源M_o的中心为坐标原点建立笛卡尔坐标系，所述笛卡尔坐标系中的xoy坐标平面平行于测量平面H，z轴方向垂直于测量平面H；测量平面H上的测量网格点坐标为(x,y,z_H)，测量平面H₁上的测量网格点坐标为(x,y,z_H1)；同步采集测量平面H上各测量网格点(x,y,z_H)处的声压时域信号p(x,y,z_H,t)和测量平面H₁上各测量网格点(x,y,z_H1)处的声压时域信号p(x,y,z_H1,t)；

步骤b、按照式(1)进行有限差分获得测量平面H上各测量网格点(x,y,z_H)处的质点加速度时域信号a(x,y,z_H,t)：

$a(x,y,z_H,t)=-\frac{p(x,y,z_{H1},t)-p(x,y,z_H,t)}{{\mathrm{\&rho;}}_0(z_{H1}-z_H)}---\left(1\right)$

式(1)中，t为时间；ρ₀为介质密度；

步骤c、对于测量平面H上的声压时域信号p(x,y,z_H,t)按照式(2)进行二维空间傅里叶变换获得声压时域波数谱P(k_x,k_y,z_H,t)，对测量平面H上的质点加速度时域信号a(x,y,z_H,t)按照式(3)进行二维空间傅里叶变换获得质点加速度时域波数谱A(k_x,k_y,z_H,t)，

$P(k_x,k_y,z_H,t)={\&Integral;}_{-\&infin;}^{\&infin;}{\&Integral;}_{-\&infin;}^{\&infin;}p(x,y,z_H,t)e^{j(k_{x}x+k_{y}y)}\mathrm{dxdy}---\left(2\right)$

$A(k_x,k_y,z_H,t)={\&Integral;}_{-\&infin;}^{\&infin;}{\&Integral;}_{-\&infin;}^{\&infin;}a(x,y,z_H,t)e^{j(k_{x}x+k_{y}y)}\mathrm{dxdy}---\left(3\right)$

在式(2)和式(3)中，j表示虚数单位；k_x、k_y分别为x、y方向的波数；

步骤d、构建声压时域波数谱P(k_x,k_y,z_H,t)、质点加速度时域波数谱A(k_x,k_y,z_H,t)、已知的时域脉冲响应函数h(k_x,k_y,0,t)和目标声源M_o单独在测量平面H上所辐射的声压时域波数谱P_o(k_x,k_y,z_H,t)之间的关系如式(4)所示：

P_o(k_x,k_y,z_H,t)=0.5[P(k_x,k_y,z_H,t)+A(k_x,k_y,z_H,t)*h(k_x,k_y,0,t)]  (4)

在式(4)中，“*”表示两个时间函数的卷积运算；

将式(4)中的时间t离散为t_n=(n-1)Δt，其中Δt为采样时间间隔，n=1,2,...,N，N为采样点总数，式(4)的离散形式如式(5)所示：

$P_o(k_x,k_y,z_H,t_n)=0.5[P(k_x,k_y,z_H,t_n)+\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}A(k_x,k_y,z_H,t_i)h(k_x,k_y,0,t_{n-i+1})]---\left(5\right)$

步骤e、对于所述声压时域波数谱P_o(k_x,k_y,z_H,t_n)按照式(6)进行二维空间傅里叶反变换，分离出每个时刻下目标声源M_o单独在测量平面H上所辐射的声压时域信号p_o(x,y,z_H,t_n)：

$p_o(x,y,z_H,t_n)=\frac{1}{{\left(2\mathrm{\&pi;}\right)}^2}\mathrm{}{\&Integral;}_{-\&infin;}^{\&infin;}{\&Integral;}_{-\&infin;}^{\&infin;}{P}_o(k_x,k_y,z_H,t_n)e^{-j(k_{x}x+k_{y}y)}{\mathrm{dk}}_x{\mathrm{dk}}_y---\left(6\right).$

2.根据权利要求1所述的采用声压和质点加速度测量的实时声场分离方法，其特征是：所述测量平面H和测量平面H₁上各测量网格点的声压时域信号p(x,y,z_H,t)和p(x,y,z_H1,t)是采用声压传感器阵列一次快照测量获得。

3.根据权利要求1所述的采用声压和质点加速度测量的实时声场分离方法，其特征是：所述干扰声源M_d为噪声源、反射源或散射源。

Images