CN102494755A - 一种非稳态声场分离方法 - Google Patents

Abstract

本发明公开了一种非稳态声场分离方法，其特征是对于同时存在目标声源M_o和干扰声源M_d，并且均辐射非稳态的声压信号的被测声场，在被测声场中布置一个与测量面S₁平行、且相隔距离为Δz的辅助测量面S₂，同步测量辅助测量面S₂与测量面S₁上的声压时域信号；并将声压时域信号分别进行二维空间傅里叶变换获得其声压波数谱；利用声压波数谱和时域脉冲响应函数，并通过解卷运算分离出只有目标声源M_o在测量面S₁上所辐射的声压波数谱；对分离出的声压波数谱进行二维空间反傅里叶变换，最终获得测量面S₁上只有目标声源M_o所辐射的声压时域信号。本发明方法既可以实现声场在空间域的分离，也可以实现声场在时域的分离，特别适用于非自由声场环境下非稳态声场的分离。

Description

一种非稳态声场分离方法

技术领域

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

背景技术

在实际工程中，目标声源位于测量面的一侧，而在测量面的另一侧往往存在干扰声源，这些干扰声源所产生的声场影响了对目标声源所辐射声场的准确测量，因此需要采用一定的声场分离方法将干扰声源的影响从测量结果中分离出来。到目前为止，国内外学者已提出多种声场分离方法，这些方法可大致分为五类：一是基于空间傅里叶变换法(SFT)的声场分离技术：G.V.Frisk等在1980年首次提出采用SFT法间接测量海洋底面的反射系数，并建立了基于SFT法的双面声场分离理论；M.Tamura于1990年将G.V.Frisk等提出的方法用于测量斜入射时材料的反射系数；M.T.Cheng等对M.Tamura提出的方法进行了推广，建立了柱面坐标下的双测量面声场分离公式，并用于实现散射声场的分离。二是基于统计最优近场声全息(SONAH)的声场分离技术：J.Hald在对SONAH研究的基础上，提出了基于双全息面声压测量的统计最优声场分离技术；F.Jacobsen等在J.Hald提出的方法的基础上，提出了基于声压和速度测量的统计最优声场分离技术。三是基于球面波叠加的声场分离方法：1956年，J.Pachner采用球面波叠加法实现了任意波场中行波和驻波声场的分离；G.Weinreich等在1980年对J.Pachner提出的方法作了进一步改进，建立了基于双球面测量的声场分离理论。四是基于边界元法(BEM)的声场分离技术：C.Langrenne等在2007年提出一种基于边界元法的双面声场分离方法；随后，E.G.Williams等在2008年提出一种基于逆边界元法和声压速度测量(Cauchy数据)的声场分离方法。五是基于等效源法的声场分离技术：C.X.Bi提出的基于等效源法的声场分离技术，适用于任意形状测量面，且计算稳定性好、计算精度高。但是，上述声场分离方法都仅仅考虑了稳态声场，仅实现了声场在空间域的分离。当目标声源和干扰声源辐射的声场是非稳态时，既需要实现声场在空间域的分离，也需要实现声场在时域的分离，因此上述声场分离方法将不再适用。

发明内容

本发明是为避免上述现有技术所存在的不足之处，提供一种非稳态声场分离方法，既实现声场在空间域的分离，又实现声场在时域的分离。

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

本发明非稳态声场分离方法的特点是对于同时存在目标声源M_o和干扰声源M_d，目标声源M_o和干扰声源M_d均辐射非稳态的声压信号的被测声场，在所述被测声场中布置一个与测量面S₁平行、且相隔距离为Δz的辅助测量面S₂，同步测量辅助测量面S₂与测量面S₁上的声压时域信号；将测量面S₁和辅助测量面S₂上的声压时域信号分别进行二维空间傅里叶变换获得其声压波数谱；利用所述声压波数谱和时域脉冲响应函数，并通过解卷运算分离出只有目标声源M_o在测量面S₁上所辐射的声压波数谱；对分离出的声压波数谱进行二维空间反傅里叶变换，从而获得测量面S₁上只有目标声源M_o所辐射的时域信号。

本发明方法的特点按如下步骤进行：

步骤a、同步测量测量面S₁和辅助测量面S₂上的声压时域信号

在由目标声源M_o和干扰声源M_d构成的非稳态被测声场中，位于目标声源M_o和干扰声源M_d之间有测量面S₁，在测量面S₁与干扰声源M_d之间设置一个与测量面S₁平行、且相隔距离为Δz的辅助测量面S₂；在测量面S₁与辅助测量面S₂上分别均匀分布M个测量网格点，同步测量辅助测量面S₂与测量面S₁上各个网格点处的声压时域信号；所述Δz的取值大于零，且不大于测量网格点的间隔。

步骤b、对测量面S₁与辅助测量面S₂上的声压时域信号分别进行二维空间傅里叶变换获得声压波数谱，变换过程为：

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

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

在式(11)和式(12)中：

j表示虚数单位；k_x、k_y分别为x、y方向的波数分量；

z₁、z₂分别为测量面S₁、辅助测量面S₂的z方向的空间坐标；

p(x，y，z₁，t)为测量面S₁上的声压时域信号；

p(x，y，z₂，t)为辅助测量面S₂上的声压时域信号；

P(k_x，k_y，z₁，t)为测量面S₁上的声压波数谱；

P(k_x，k_y，z₂，t)为辅助测量面S₂上的声压波数谱；

步骤c、建立测量面S₁与辅助测量面S₂上的声压波数谱、时域脉冲响应函数、目标声源M_o在测量面S₁上所辐射的声压波数谱之间的传递关系

P(k_x，k_y，z₁，t)-P(k_x，k_y，z₂，t)*h(k_x，k_y，Δz，t)

＝P_o(k_x，k_y，z₁，t)*[δ(t)-h(k_x，k_y，Δz，t)*h(k_x，k_y，Δz，t)]    (13)

式(13)中

“*”表示卷积运算；

δ(t)为已知的Dirac delta函数；

h(k_x，k_y，Δz，t)为已知的时域脉冲响应函数；

P_o(k_x，k_y，z₁，t)为目标声源M_o在测量面S₁上所辐射的声压波数谱，为待求量；

步骤d、通过解卷运算分离出只有目标声源M_o在测量面S₁上所辐射的声压波数谱P_o(k_x，k_y，z₁，t)，其解卷过程如下：

P_o＝B^-1A                                                           (14)

式(14)中

A＝[A(k_x，k_y，t₁) A(k_x，k_y，t₂) Λ A(k_x，k_y，t_N)]^T                 (15)

P_o＝[P_o(k_x，k_y，z₁，t₁) P_o(k_x，k_y，z₁，t₂) Λ P_o(k_x，k_y，z₁，t_N)]^T (16)

$B=\left[\begin{array}{cccc}B(k_x,k_y,t_1)& 0& \mathrm{\Λ}& 0\\ B(k_x,k_y,t_2)& B(k_x,k_y,t_1)& O& M\\ M& O& O& 0\\ B(k_x,k_y,t_N)& \mathrm{\Λ}& B(k_x,k_y,t_2)& B(k_x,k_y,t_1)\end{array}\right]---\left(17\right)$

在式(15)至式(17)中，t_n(n＝1，2，Λ N)表示各离散时间点，N表示时间离散点数；

在式(15)中，A(k_x，k_y，t_n)表示当t＝t_n时A(k_x，k_y，t)的取值，其中

A(k_x，k_y，t)＝P(k_x，k_y，z₁，t)-P(k_x，k_y，z₂，t)*h(k_x，k_y，Δz，t)  (18)

在式(16)中，P_o(k_x，k_y，z₁，t_n)表示当t＝t_n时P_o(k_x，k_y，z₁，t)的取值；

在式17中，B(k_x，k_y，t_n)表示当t＝t_n时B(k_x，k_y，t)的取值，其中

B(k_x，k_y，t)＝δ(t)-h(k_x，k_y，Δz，t)*h(k_x，k_y，Δz，t)    (19)

步骤e、对由步骤d所分离出的波数谱P_o(k_x，k_y，z₁，t)进行二维空间反傅里叶变换，获得测量面S₁上只有目标声源M_o所辐射的时域信号，其变换过程如下：

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

本发明非稳态声场分离方法的特点也在于：

所述各网格点上的声压时域信号是采用双传声器阵列在两测量面上一次快照获得。

所述的测量面S₁和辅助测量面S₂均为平面。

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

理论模型：

在测量面S₁的两侧分别布置着目标声源M_o和干扰声源M_d，在测量面S₁与干扰声源M_d之间设置一个与测量面S₁平行、且相隔距离为Δz的辅助测量面S₂。

根据声波的叠加原理，在测量面S₁上所测得的声压等于目标声源M_o在测量面S₁上所辐射的声压与干扰声源M_d在测量面S₁所辐射的声压之和，即

p(x，y，z₁，t)＝p_o(x，y，z₁，t)+p_d(x，y，z₁，t)            (21)

式(21)中，下标“o”表示目标声源M_o，下标“d”表示干扰声源M_d。

同理，在辅助测量面S₂上所测得的声压等于目标声源M_o在辅助测量面S₂上所辐射的声压与干扰声源M_d在辅助测量面S₂所辐射的声压之和，即

p(x，y，z₂，t)＝p_o(x，y，z₂，t)+p_d(x，y，z₂，t)            (21)

定义声压p(x，y，z，t)关于x、y的二维空间傅里叶变换为

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

P(k_x，k_y，z，t)表示声压波数谱。

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

P(k_x，k_y，z₁，t)＝P_o(k_x，k_y，z₁，t)+P_d(k_x，k_y，z₁，t)     (24)

P(k_x，k_y，z₂，t)＝P_o(k_x，k_y，z₂，t)+P_d(k_x，k_y，z₂，t)     (25)

由声压在时域波数域的传播方程可知

P_o(k_x，k_y，z₂，t)＝P_o(k_x，k_y，z₁，t)*h(k_x，k_y，Δz，t)    (26)

P_d(k_x，k_y，z₁，t)＝P_d(k_x，k_y，z₂，t)*h(k_x，k_y，Δz，t)    (27)

式(26)和(27)中，h(k_x，k_y，Δz，t)为时域脉冲响应函数，其表达式为

$h(k_x,k_y,\mathrm{\Δz},t)=\mathrm{\δ}(t-\mathrm{\Δz}/c)-$

$\mathrm{\Δz}\sqrt{{k_x}^2+{k_y}^2}\frac{J_1\left(c\sqrt{{k_x}^2+{k_y}^2}\sqrt{t^2-{\mathrm{\Δz}}^2/c^2}\right)}{\sqrt{t^2-{\mathrm{\Δz}}^2/c^2}}H(t-\mathrm{\Δz}/c)---\left(28\right)$

式(28)中，c为声速，J₁表示第一类贝塞尔函数，阶数为1，H(t)表示Heaviside函数。

将式(26)代入式(25)中得

P(k_x，k_y，z₂，t)＝P_o(k_x，k_y，z₁，t)*h(k_x，k_y，Δz，t)+P_d(k_x，k_y，z₂，t)    (29)

将式(27)代入式(24)中得

P(k_x，k_y，z₁，t)＝P_o(k_x，k_y，z₁，t)+P_d(k_x，k_y，z₂，t)*h(k_x，k_y，Δz，t)    (30)

将式(29)两边分别与时域脉冲响应函数h(k_x，k_y，Δz，t)进行卷积，然后与式(30)相减可进一步得到

P(k_x，k_y，z₁，t)-P(k_x，k_y，z₂，t)*h(k_x，k_y，Δz，t)

＝P_o(k_x，k_y，z₁，t)-P_o(k_x，k_y，z₁，t)*h(k_x，k_y，Δz，t)*h(k_x，k_y，Δz，t)

＝P_o(k_x，k_y，z₁，t)*[δ(t)-h(k_x，k_y，Δz，t)*h(k_x，k_y，Δz，t)]            (31)

式(31)中，P(k_x，k_y，z₁，t)和P(k_x，k_y，z₂，t)可分别通过对测得的声压p(x，y，z₁，t)和p(x，y，z₂，t)进行二维空间傅里叶变换得到，δ(t)和h(k_x，k_y，Δz，t)均为已知函数，通过对式(31)进行解卷运算，求得目标声源M_o在测量面S₁上所辐射的声压波数谱P_o(k_x，k_y，z₁，t)。

解卷过程：

定义

A(k_x，k_y，t)＝P(k_x，k_y，z₁，t)-P(k_x，k_y，z₂，t)*h(k_x，k_y，Δz，t)    (32)

B(k_x，k_y，t)＝δ(t)-h(k_x，k_y，Δz，t)*h(k_x，k_y，Δz，t)              (33)

则式(31)重写成为：

A(k_x，k_y，t)＝P_o(k_x，k_y，z₁，t)*B(k_x，k_y，t)                         (34)

在实际计算过程中，需要将时间t离散为t₁，t₂，Λ t_N。根据离散卷积方程，式(34)的离散形式可表示为：

$A(k_x,k_y,t_n)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{P}_o(k_x,k_y,z_1,t_i)B(k_x,k_y,t_{n+1-i}),n=1,\mathrm{\Λ},N.---\left(35\right)$

对于每一个波数对(k_x，k_y)，可将式(35)写成矩阵形式

A＝BP_o                                                               (36)

式(36)中，

A＝[A(k_x，k_y，t₁) A(k_x，k_y，t₂) Λ A(k_x，k_y，t_N)]^T                   (37)

P_o＝[P_o(k_x，k_y，z₁，t₁) P_o(k_x，k_y，z₁，t₂) Λ P_o(k_x，k_y，z₁，t_N)]^T   (38)

$B=\left[\begin{array}{cccc}B(k_x,k_y,t_1)& 0& \mathrm{\Λ}& 0\\ B(k_x,k_y,t_2)& B(k_x,k_y,t_1)& O& M\\ M& O& O& 0\\ B(k_x,k_y,t_N)& \mathrm{\Λ}& B(k_x,k_y,t_2)& B(k_x,k_y,t_1)\end{array}\right]---\left(39\right)$

对式(36)直接求逆可求得

P_o＝B^-1A                                                             (40)

一旦求得目标声源M_o在测量面S₁上所辐射的声压波数谱P_o(k_x，k_y，z₁，t)，再对其进行二维空间反傅里叶变换，即可获得测量面S₁上只有目标声源M_o所辐射的时域信号，其变换过程为：

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

在上述求逆过程中，为保证求解过程的稳定性，采用Tikhonov正则化方法，正则化参数用L曲线法来进行选取。

通过上述方法，既实现了信号在空间域的分离，也实现了信号在时域的分离。

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

1、本发明考虑了声场在时域内的分离，因而与传统的方法相比，本发明方法特别适用于对非稳态声场进行分离。

2、本发明方法也可以用于对稳态声场进行分析。

附图说明

图1为本发明双测量面非稳态声场分离方法示意图；

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

图2(b)为本发明方法测量面网格划分示意图：图中测量面S₁上的测量点A和测量点C点分别正对着目标声源M_o和干扰声源M_d，测量点B位于测量面S₁的中心，测量点D靠近测量面S₁的边缘；

图3(a)为本发明方法测量点A处声压时域信号图：图中虚线表示测量面S₁上的直接测量值(包含目标声源M_o和干扰声源M_d所辐射的声压)，实线表示目标声源M_o在测量面S₁上所辐射的声压理论值，点线表示采用本发明方法分离出的声压值；

图3(b)为本发明方法测量点B处声压时域信号图：图中虚线表示测量面S₁上的直接测量值(包含目标声源M_o和干扰声源M_d所辐射的声压)，实线表示目标声源M_o在测量面S₁上所辐射的声压理论值，点线表示采用本发明方法分离出的声压值；

图3(c)为本发明方法测量点C处声压时域信号图：图中虚线表示测量面S₁上的直接测量值(包含目标声源M_o和干扰声源M_d所辐射的声压)，实线表示目标声源M_o在测量面S₁上所辐射的声压理论值，点线表示采用本发明方法分离出的声压值；

图3(d)为本发明方法测量点D处声压时域信号图：图中虚线表示测量面S₁上的直接测量值(包含目标声源M_o和干扰声源M_d所辐射的声压)，实线表示目标声源M_o在测量面S₁上所辐射的声压理论值，点线表示采用本发明方法分离出的声压值；

图4(a)为本发明方法所获得的时域评价因子T₁的数值分布图，图中等高线值为0.99；

图4(b)为本发明方法所获得的时域评价因子T₂的数值分布图，图中等高线值为0.05；

图5(a)为本发明方法t₁＝4.2ms时刻目标声源M_o在测量面S₁上所辐射的理论声压场p_t；

图5(b)为本发明方法t₁＝4.2ms时刻直接测量的声压场p_m(包含目标声源M_o和干扰声源M_d所辐射的声压场)；

图5(c)为t₁＝4.2ms时刻采用本发明方法分离出的声压场p_c；

图5(d)为本发明方法t₂＝5.0ms时刻目标声源M_o在测量面S₁上所辐射的理论声压场p_t；

图5(e)为本发明方法t₂＝5.0ms时刻直接测量的声压场p_m(包含目标声源M_o和干扰声源M_d所辐射的声压场)；

图5(f)为t₂＝5.0ms时刻采用本发明方法分离出的声压场p_c；

图6为本发明方法所获得的空间误差衡量标准

的时间变化图。

具体实施方式

参见图1，本实施例中，在测量面S₁的两侧分别布置着目标声源M_o和干扰声源M_d，干扰声源M_d可以为噪声源、反射或散射源。在测量面S₁与干扰声源M_d之间设置一与测量面S₁平行、且相隔距离为Δz的辅助测量面S₂。在测量面S₁与辅助测量面S₂上分别均匀分布有M个测量网格点。Δz的取值大于零，且不大于测量网格点的间隔。

具体实施步骤为：

a、采用双传声器阵列一次快照获得测量面S₁与辅助测量面S₂上的声压时域信号；

b、对测量面S₁和辅助测量面S₂上的声压时域信号分别进行二维空间傅里叶变换获得其声压波数谱，其变换过程如下：

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

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

在式(42)和式(43)中：

j表示虚数单位；k_x、k_y分别为x、y方向的波数分量；

z₁、z₂分别为测量面S₁、辅助测量面S₂的z方向的空间坐标；

p(x，y，z₁，t)为测量面S₁上的声压时域信号；

p(x，y，z₂，t)为辅助测量面S₂上的声压时域信号；

P(k_x，k_y，z₁，t)为测量面S₁上的声压波数谱；

P(k_x，k_y，z₂，t)为辅助测量面S₂上的声压波数谱；

c、建立测量面S₁与辅助测量面S₂上的声压波数谱、时域脉冲响应函数、目标声源M_o在测量面S₁上所辐射的声压波数谱之间的传递关系

P(k_x，k_y，z₁，t)-P(k_x，k_y，z₂，t)*h(k_x，k_y，Δz，t)

＝P_o(k_x，k_y，z₁，t)*[δ(t)-h(k_x，k_y，Δz，t)*h(k_x，k_y，Δz，t)]    (44)

在式(44)中：

“*”表示卷积运算；

δ(t)为已知的Dirac delta函数；

h(k_x，k_y，Δz，t)为已知的时域脉冲响应函数；

P_o(k_x，k_y，z₁，t)为目标声源M_o在测量面S₁上所辐射的声压波数谱，为待求量；

d、通过解卷运算分离出只有目标声源M_o在测量面S₁上所辐射的声压波数谱P_o(k_x，k_y，z₁，t)，其解卷过程如下：

P_o＝B^-1A                                                           (45)

式(45)中

A＝[A(k_x，k_y，t₁) A(k_x，k_y，t₂) Λ A(k_x，k_y，t_N)]^T                 (46)

P_o＝[P_o(k_x，k_y，z₁，t₁) P_o(k_x，k_y，z₁，t₂) Λ P_o(k_x，k_y，z₁，t_N)]^T (47)

$B=\left[\begin{array}{cccc}B(k_x,k_y,t_1)& 0& \mathrm{\Λ}& 0\\ B(k_x,k_y,t_2)& B(k_x,k_y,t_1)& O& M\\ M& O& O& 0\\ B(k_x,k_y,t_N)& \mathrm{\Λ}& B(k_x,k_y,t_2)& B(k_x,k_y,t_1)\end{array}\right]---\left(48\right)$

在式(46)至式(48)中，t_n(n＝1，2，Λ N)表示各离散时间点，N表示时间离散点数；

在式(46)中，A(k_x，k_y，t_n)表示当t＝t_n时A(k_x，k_y，t)的取值，其中

A(k_x，k_y，t)＝P(k_x，k_y，z₁，t)-P(k_x，k_y，z₂，t)*h(k_x，k_y，Δz，t) (49)

在式(47)中，P_o(k_x，k_y，z₁，t_n)表示当t＝t_n时P_o(k_x，k_y，z₁，t)的取值；

在式(48)中，B(k_x，k_y，t_n)表示当t＝t_n时B(k_x，k_y，t)的取值，其中

B(k_x，k_y，t)＝δ(t)-h(k_x，k_y，Δz，t)*h(k_x，k_y，Δz，t)           (50)

e、对由步骤d所分离出的波数谱P_o(k_x，k_y，z₁，t)进行二维空间反傅里叶变换，获得测量面S₁上只有目标声源M_o所辐射的时域信号，其变换过程如下：

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

方法的检验：

在测量面S₁两侧分别布置一个单极子，两个单极子均辐射非稳态信号。采用本发明的声场分离方法将测量面S₁上目标声源M_o所辐射的声压分离出来，并与其理论声压值进行比较。

本实施例中，作为目标声源M_o的单极子辐射一线性频率调制(频率调制范围为[200，1800]Hz)、高斯幅值调制信号；作为干扰声源M_d的单极子辐射一Morlet小波信号，其表达式为

$s\left(t\right)=\cos \left(2\mathrm{\π}{f}_{0}t\right)e^{-t^2/2}---\left(52\right)$

式中，f₀＝800Hz。

测量面S₁与辅助测量面S₂的位置关系参见图2(a)，测量面S₁位于z₁＝0m的平面上，辅助测量面S₂位于z₂＝0.03m的平面上，因此测量面S₁与辅助测量面S₂之间的距离Δz为0.03m。测量面S₁与辅助测量面S₂的大小均为0.648m×0.648m，其上均匀地分布19×19个测量点，参见图2(b)。目标声源M_o位于(0.252，0.324，-0.15)m处，干扰声源M_d位于(0.432，0.324，0.18)m处。时域信号采样频率为34400Hz，采样点数为256。

在本实施例的数据处理过程中，为进一步减少截断误差对分离结果的影响，将测量面矩阵维数由19×19通过补零扩展为37×37。

为检验本发明方法在时域内的声场分离效果，在测量面S₁上选取了四个测量点，即测量点A、测量点B、测量点C和测量点D，其位置分别为A(0.252，0.324，0)m、B(0.324，0.324，0)m、C(0.432，0.324，0)m、D(0.540，0.324，0)m，其中测量点A和测量点C分别正对着目标声源M_o和干扰声源M_d，测量点B位于测量面S₁的中心，测量点D靠近测量面S₁的边缘。参见图3(其中，图3(a)、3(b)、3(c)、3(d)分别对应测量点A、测量点B、测量点C和测量点D，图中虚线表示测量面S₁上的直接测量值(包含目标声源M_o和干扰声源M_d所辐射的声压)，图中实线表示目标声源M_o在测量面S₁上所辐射的声压理论值，图中点线表示采用本发明方法分离出的声压值)，比较图中的虚线和实线可以看出，干扰声源M_d对目标声源M_o在测量面S₁上所辐射的声压造成了较大的干扰；比较图中的实线和点线可以看出，采用本发明方法可以很好地消除干扰声源M_d的影响，从而分离出只有目标声源M_o在测量面S₁上所辐射的声压值。

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

$T_1(x_i,y_j)=\frac{<p_t(x_i,y_j,z_1,t)p_c(x_i,y_j,z_1,t)>}{\sqrt{<p_t^2(x_i,y_j,z_1,t)><p_c^2(x_i,y_j,z_1,t)>}}---\left(53\right)$

$T_2(x_i,y_j)=\frac{|p_t^{\mathrm{rms}}(x_i,y_j,z_1)-p_c^{\mathrm{rms}}(x_i,y_j,z_1)|}{p_t^{\mathrm{rms}}(x_i,y_j,z_1)}---\left(54\right)$

式中，<>表示求平均值，下标“t”表示声压理论值，下标“c”表示声压分离值，上标“rms”表示求均方根值。评价因子T₁是用来衡量声压理论值和声压分离值之间的相位误差，当T₁的值越靠近1时，相位误差越小。评价因子T₂是用来衡量声压理论值和声压分离值之间的幅值误差，当T₂的值越靠近0时，幅值误差越小。运用公式(53)、(54)分别计算了测量面S₁上各个测量点处的T₁和T₂。参见图4(a)和图4(b)，在大多数测量点处，不论是相位还是幅值，声压理论值和声压分离值都吻合地较好，只是在测量面边缘处吻合程度略低。

为检验本发明方法在空间域的声场分离效果，选取了两个时刻t₁＝4.2ms和t₂＝5.0ms。参见图5(图5(a)和5(d)分别为t₁和t₂时刻目标声源M_o在测量面S₁上所辐射的理论声压场p_t，图5(b)和5(e)分别为t₁和t₂时刻直接测量的声压场p_m(包含目标声源M_o和干扰声源M_d所辐射的声压场)，图5(c)和5(f)分别为t₁和t₂时刻采用本发明方法分离出的声压场p_c)，比较图5(a)和图5(b)、图5(d)和图5(e)可以看出，干扰声源M_d对目标声源M_o在测量面S₁上所辐射的理论声压场造成了严重的干扰；比较图5(a)和图5(c)、图5(d)和图5(f)可以看出，采样本发明方法可以很好地消除干扰声源M_d的影响，从而还原只有目标声源M_o在测量面S₁上所辐射的声压场。

为了更加客观地评价本发明方法在空间域的分离效果，在此定义一个空间误差衡量标准，其表达式为

$E_{x,y}^r\left(t_n\right)=\frac{\sqrt{<{(p_c(x,y,t_n)-p_t(x,y,t_n))}^2>}}{\sqrt{<p_t^2(x,y,t_n)>}}---\left(55\right)$

该标准反映了目标声源M_o在测量面S₁上所辐射的理论声压场与采用本发明方法分离出的声压场之间的相对误差。运用公式(55)计算了每个时刻的相对误差。参见图6，在大多数时刻，目标声源M_o在测量面S₁上所辐射的理论声压场与采用本发明方法分离出的声压场之间的相对误差都较小；较大的误差出现在信号的两端，这是因为在信号的两端理论声压值较小(即公式(55)的分母较小)，从而造成误差变大。

上述例子表明不论是在时域还是在空间域，采用本发明方法都可以很好地将干扰声源M_d的影响从测量声压中消除，从而分离出目标声源M_o在测量面S₁上所辐射的声压场。

Claims (4)

1.一种非稳态声场分离方法，其特征是对于同时存在目标声源M_o和干扰声源M_d，目标声源M_o和干扰声源M_d均辐射非稳态的声压信号的被测声场，在所述被测声场中布置一个与测量面S₁平行、且相隔距离为Δz的辅助测量面S₂，同步测量辅助测量面S₂与测量面S₁上的声压时域信号；将测量面S₁和辅助测量面S₂上的声压时域信号分别进行二维空间傅里叶变换获得其声压波数谱；利用所述声压波数谱和时域脉冲响应函数，并通过解卷运算分离出只有目标声源M_o在测量面S₁上所辐射的声压波数谱；对分离出的声压波数谱进行二维空间反傅里叶变换，从而获得测量面S₁上只有目标声源M_o所辐射的时域信号。

2.根据权利要求1所述的非稳态声场分离方法，其特征是按如下步骤进行：

步骤a、同步测量测量面S₁和辅助测量面S₂上的声压时域信号

在由目标声源M_o和干扰声源M_d构成的非稳态被测声场中，位于目标声源M_o和干扰声源M_d之间有测量面S₁，在测量面S₁与干扰声源M_d之间设置一个与测量面S₁平行、且相隔距离为Δz的辅助测量面S₂；在测量面S₁与辅助测量面S₂上分别均匀分布M个测量网格点，同步测量辅助测量面S₂与测量面S₁上各个网格点处的声压时域信号；所述Δz的取值大于零，且不大于测量网格点的间隔；

步骤b、对测量面S₁与辅助测量面S₂上的声压时域信号分别进行二维空间傅里叶变换获得声压波数谱，变换过程为：

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

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

在式(1)和式(2)中：

j表示虚数单位；k_x、k_y分别为x、y方向的波数分量；

z₁、z₂分别为测量面S₁、辅助测量面S₂的z方向的空间坐标；

p(x，y，z₁，t)为测量面S₁上的声压时域信号；

p(x，y，z₂，t)为辅助测量面S₂上的声压时域信号；

P(k_x，k_y，z₁，t)为测量面S₁上的声压波数谱；

P(k_x，k_y，z₂，t)为辅助测量面S₂上的声压波数谱；

步骤c、建立测量面S₁与辅助测量面S₂上的声压波数谱、时域脉冲响应函数、目标声源M_o在测量面S₁上所辐射的声压波数谱之间的传递关系

P(k_x，k_y，z₁，t)-P(k_x，k_y，z₂，t)*h(k_x，k_y，Δz，t)

＝P_o(k_x，k_y，z₁，t)*[δ(t)-h(k_x，k_y，Δz，t)*h(k_x，k_y，Δz，t)]    (3)

式(3)中

“*”表示卷积运算；

δ(t)为已知的Dirac delta函数；

h(k_x，k_y，Δz，t)为已知的时域脉冲响应函数；

P_o(k_x，k_y，z₁，t)为目标声源M_o在测量面S₁上所辐射的声压波数谱，为待求量；

步骤d、通过解卷运算分离出只有目标声源M_o在测量面S₁上所辐射的声压波数谱P_o(k_x，k_y，z₁，t)，其解卷过程如下：

P_o＝B^-1A                                                           (4)

式(4)中

A＝[A(k_x，k_y，t₁) A(k_x，k_y，t₂) Λ A(k_x，k_y，t_N)]^T                 (5)

P_o＝[P_o(k_x，k_y，z₁，t₁) P_o(k_x，k_y，z₁，t₂) Λ P_o(k_x，k_y，z₁，t_N)]^T (6)

$B=\left[\begin{array}{cccc}B(k_x,k_y,t_1)& 0& \mathrm{\&Lambda;}& 0\\ B(k_x,k_y,t_2)& B(k_x,k_y,t_1)& O& M\\ M& O& O& 0\\ B(k_x,k_y,t_N)& \mathrm{\&Lambda;}& B(k_x,k_y,t_2)& B(k_x,k_y,t_1)\end{array}\right]---\left(7\right)$

在式(5)至式(7)中，t_n(n＝1，2，Λ N)表示各离散时间点，N表示时间离散点数；在式(5)中，A(k_x，k_y，t_n)表示当t＝t_n时A(k_x，k_y，t)的取值，其中

A(k_x，k_y，t)＝P(k_x，k_y，z₁，t)-P(k_x，k_y，z₂，t)*h(k_x，k_y，Δz，t)  (8)

在式(6)中，P_o(k_x，k_y，z₁，t_n)表示当t＝t_n时P_o(k_x，k_y，z₁，t)的取值；

在式(7)中，B(k_x，k_y，t_n)表示当t＝t_n时B(k_x，k_y，t)的取值，其中

B(k_x，k_y，t)＝δ(t)-h(k_x，k_y，Δz，t)*h(k_x，k_y，Δz，t)        (9)

步骤e、对由步骤d所分离出的波数谱P_o(k_x，k_y，z₁，t)进行二维空间反傅里叶变换，获得测量面S₁上只有目标声源M_o所辐射的时域信号，其变换过程如下：

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

2.根据权利要求1所述的非稳态声场分离方法，其特征是所述各网格点上的声压时域信号是采用双传声器阵列在两测量面上一次快照获得。

3.根据权利要求1所述的非稳态声场分离方法，其特征是所述测量面S₁和辅助测量面S₂均为平面。

4.根据权利要求1所述的非稳态声场分离方法，其特征是所述干扰声源M_d为噪声源、反射源或散射源。

Images