CN109798975A - 一种采用声压和质点加速度测量的非稳态平面声源的自由场实时还原方法 - Google Patents

Abstract

本发明公开了一种采用声压和质点加速度测量的非稳态平面声源的自由场实时还原方法，将测量平面H和辅助测量平面H₁上的时域声压进行有限差分获得平面H上的时域质点加速度，并进行二维空间傅里叶变换获得平面H上的声压和质点加速度时域波数谱；再利用平面H上的声压时域波数谱和质点加速度时域波数谱、已知的声压与质点加速度和声压与声压之间的时域脉冲响应函数以及目标声源表面反射系数，推导自由场还原公式，实时还原出目标声源自由场条件下在平面H上所辐射的时变声压信号。本发明方法在时域上消除了非稳态散射声场的影响，无需进行任何求逆和正则化运算处理，具备实时还原声场的能力，可用于实际声场环境下现场分析目标声源的时变辐射特性和振动特性。

Description

一种采用声压和质点加速度测量的非稳态平面声源的自由场 实时还原方法

技术领域

本发明涉及物理专业中噪声类领域，具体涉及一种采用声压和质点加速度测量的非稳态平面声源的自由场实时还原方法。

背景技术

实际声场环境是非自由场，它不仅包含目标声源的辐射声，而且还会包含来自测量面另一侧的干扰声以及干扰声在目标声源表面产生的散射声。除了干扰声所产生的声场能够影响目标声源辐射的声场外，散射声也能够影响对目标声源辐射声场的准确测量，因此需要一种自由场还原方法将干扰声和散射声的影响从测量结果中消除。

到目前为止，国内外学者已提出多种自由场还原方法：Hald等提出了基于统计最优近场声全息(SONAH)的自由场还原方法，首先将向外声场和向内声场进行分离，然后通过目标声源表面边界条件计算得到散射声场，最后从向外声场中去除散射声场，彻底还原出平面目标声源在自由场条件下辐射的声场。Hu等利用二维空间傅里叶变换和平面声源的表面反射系数，推导了非自由场中基于双面声压测量的平面近场声全息的完备计算公式，从而实现了噪声环境中平面目标声源的准确声场重建。为了适应于球形声源，Braikia等提出了基于球面波叠加法(SWSM)的双面声压测量的自由场还原技术，通过镜像原理进一步去除了半球测量面外部干扰声在目标声源表面上散射声的影响，从而实现了小内空间的声场重建。毕传兴等将目标声源表面导纳引入Braikia的方法，避免了刚性表面假设，并去除背侧干扰声在目标声源表面产生的散射声。Langrenne等提出了基于边界元法(BEM)的双面声压测量的自由场还原方法，将任意外形目标声源在自由场条件下辐射的声场从混合声场中彻底还原出来，但会带入奇异值积分、计算效率低和解非唯一性等问题。为了高效地实现任意形状声源的自由场还原，Bi等提出了基于等效源法的自由场还原方法。然而，上述这些自由场还原方法都是基于双面声压测量，两测量面之间的距离会影响还原精度，甚至不合理的距离会导致错误的还原结果。为了避免双测量面之间距离的影响，Fernandez-Grande等将单面声压-振速测量和SONAH结合，提出了基于单面测量的自由场还原方法，实现了目标声源的声场还原。Braikia等也将单面声压-振速测量引入基于SWSM的自由场还原方法，实现了非消声室环境下目标声源的识别。Langrenne等提出了基于单面声压-振速测量的BEM的自由场还原方法，在封闭空间里还原了任意外形目标声源在自由场条件下的声压场。

在上述基于双面测量和单面测量的自由场还原方法中，还原的对象都是稳态声源，无法用于非稳态声源。非稳态声场的还原还需要在时域上消除干扰声和散射声的影响，还原出目标声源随时间变化的声学特性。为此，Zhang等在2012年提出了基于双面声压测量的时域声场分离方法，在时域和空间域上消除了干扰声源的影响。为了减少计算时间，Bi等在2014年提出了基于声压和加速度测量的实时声场分离方法。为了避免有限差分引起的误差，并提高分离精度，Bi等提出了基于单层声压和振速测量的非稳态声场分离方法。但是这些时域分离方法仅仅能够去除测量平面另一侧干扰声源产生的非稳态声场，无法消除干扰声在目标声源表面产生的非稳态散射声场。为实现非稳态散射声场的消除，Geng等提出了一种基于双面声压测量的非稳态声场还原方法。在该非稳态还原方法中，目标声源在自由场中辐射的时域声压与双测量面上的混合时域声压不能够直接建立关系，需要先进行非稳态声场分离，再去除非稳态散射声场，无法实现实时还原；该还原方法仍然需要复杂的迭代解卷计算，且其过程是病态的逆向问题，需要Tikhonov正则化和奇异值分解来获得合适的解，导致还原方法计算量大，且还原精度很大程度上依赖于正则化参数的合理选取；双测量面之间的距离对还原精度也有影响。

发明内容

为实现任意时刻下非稳态散射声场的消除，本发明提供一种采用声压和质点加速度测量的非稳态平面声源的自由场实时还原方法。该方法无需任何求逆和正则化运算处理，且不存在解卷过程，因而实现简单，稳定性高；同时该方法具备实时还原非稳态声场的能力，可用于实际声场环境下现场分析目标声源的时变辐射特性和振动特性。

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

一种采用声压和质点加速度测量的非稳态平面声源的自由场实时还原方法，步骤如下：

步骤1、在目标声源S_o和干扰声源S_d之间布置相互平行的测量平面H和辅助测量平面H₁；所述的目标声源S_o和干扰声源S_d均可辐射任意线性非稳态声场信号；在所述测量平面H和辅助测量平面H₁上分别均匀布置M个测量网格点，测量平面H和辅助测量平面H₁的测量网格大小和网格点位置相同；以目标声源S_o为研究对象，将声场分为非稳态向内声场和向外声场，其中非稳态向内声场来源于测量平面另一侧的干扰声源S_d，非稳态向外声场包含目标声源S_o自由场条件下辐射的非稳态声场和干扰声在目标声源表面产生的非稳态散射声场；建立笛卡尔坐标系o(x,y,z)，所述笛卡尔坐标系中的xoy坐标平面平行于测量平面H，z轴方向垂直于测量平面H；测量平面H上的测量网格点坐标为(x,y,z_H)，辅助测量平面H₁的测量网格坐标为(x,y,z_H1)；同步采集t时刻测量平面H上各测量网格点(x,y,z_H)处的时域声压信号p(x,y,z_H,t)和辅助测量平面H₁上各测量网格点(x,y,z_H1)处的时域声压信号p(x,y,z_H1,t)；

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

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

步骤3、对于测量平面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)，

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

步骤4、构建测量平面H上声压时域波数谱P(k_x,k_y,z_H,t)与质点加速度时域波数谱A(k_x,k_y,z_H,t)、已知的时域响应函数T₁(k_x,k_y,t)与T₂(k_x,k_y,t)和目标声源S_o自由场条件下在测量平面H上的声压时域波数谱P_f(k_x,k_y,z_H,t)之间的关系如式(4)所示：

式(4)中，

T₁(k_x,k_y,t)＝δ(t)-G_pp(k_x,k_y,Δz_H0,t)*R(k_x,k_y,t)*G_pp(k_x,k_y,Δz_H0,t) (5)

T₂(k_x,k_y,t)＝G_pa(k_x,k_y,0,t)*[δ(t)+G_pp(k_x,k_y,Δz_H0,t)*R(k_x,k_y,t)*G_pp(k_x,k_y,Δz_H0,t)] (6)

G_pp(k_x,k_y,Δz_H0,t)和G_pa(k_x,k_y,0,t)表示已知的时域脉冲响应函数，R(k_x,k_y,t)表示已知的目标声源表面材料的复平面波反射系数时域波数谱，z₀表示目标声源S_o表面所在位置，“*”表示两个时间函数的卷积运算，Δz_H0＝z_H-z₀，δ(t)为Dirac函数；

将式(4)中的时间t离散为t_i＝(i-1)Δt，其中Δt为采样时间间隔，i＝1,2,...,I，I为采样点总数，则式(4)的离散形式如式(7)所示：

步骤5、对于所述声压时域波数谱P_f(k_x,k_y,z_H1,t)按照式(8)进行二维空间傅里叶反变换，还原出每个时刻目标声源S_o自由场条件下在测量平面H上辐射的时域声压信号p_f(x,y,z_H,t)；

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

所述的目标声源的尺寸相对于声波波长要大，非稳态散射声场影响显著。

所述目标声源和干扰声源辐射的声场为非稳态，在时域内进行实时还原。

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

(1)本发明方法直接在时域上消除干扰声在目标声源表面产生的非稳态散射声场，还原过程是一个简单的正向卷积叠加过程，无需任何求逆和正则化运算处理，因而计算时间少，具有较高的稳定性。

(2)本发明方法具有实时还原非稳态声场的能力，因而更适用于实际声场环境下现场分析目标声源的时变辐射特性和振动特性。

附图说明

图1为本发明采用声压和质点加速度测量的非稳态平面声源的自由场实时还原方法示意图；

图2为本发明方法声源、测量平面位置分布图和测量平面网格划分示意图；

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

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

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

图3(d)为本发明方法测量平面H上点D处的时域声压信号图；

图4(a)为本发明方法t＝3.40ms时刻目标声源S_o自由场条件下在测量平面H上辐射的理论声压场p_m；

图4(b)为本发明方法t＝3.40ms时刻的混合声压场p_c；

图4(c)为t＝3.40ms时刻利用实时声场分离方法分离的声压场p_e；

图4(d)为t＝3.40ms时刻采用本发明方法还原的声压场p_r；

图4(e)为本发明方法t＝6.88ms时刻目标声源S_o自由场条件下在测量平面H上辐射的理论声压场p_m；

图4(f)为本发明方法t＝6.88ms时刻的混合声压场p_c；

图4(g)为t＝6.88ms时刻利用实时声场分离方法分离的声压场p_e；

图4(h)为t＝6.88ms时刻采用本发明方法还原的声压场p_r；

图5为测量平面H上本发明方法和实时声场分离方法的相位评价因子E_p值和幅值评价因子E_a值。

具体实施方式

下面结合附图以及具体实施方式进一步说明本发明。

本发明一种采用声压和质点加速度测量的非稳态平面声源的自由场实时还原方法，其特征是，在目标声源和干扰声源之间布置两个相互平行的测量平面H和辅助测量平面H₁，同步采集两测量平面上各网格点的时域声压信号，并通过有限差分技术获得测量平面H上的质点加速度时域信号；将测量平面H上的时域声压和质点加速度信号进行二维空间傅里叶变换获得平面H上各波数处的声压时域波数谱和质点加速度时域波数谱；再利用测量平面H上的声压时域波数谱和质点加速度时域波数谱、已知的声压与质点加速度之间的时域脉冲响应函数、已知的声压与声压之间的时域脉冲响应函数和目标声源表面反射系数，同时去除测量平面另一侧干扰声源辐射的非稳态声场和干扰声源在目标声源表面产生的非稳态散射声场，实时还原出目标声源自由场条件下在测量平面H上所辐射的时变声压信号。本发明方法在时域上消除了非稳态散射声场，无需进行任何求逆和正则化运算处理，并具备实时还原非稳态声场的能力，可用于实际声场环境下现场分析目标声源的时变辐射特性和振动特性。

参见图1，按如下步骤还原出目标声源S_o自由场条件下在测量平面H₁上辐射的时域声压信号：

步骤1、按如图1所示，在目标声源S_o和干扰声源S_d之间布置相互平行的测量平面H和辅助测量平面H₁；所述的目标声源S_o和干扰声源S_d均可辐射任意线性非稳态声场信号；在所述测量平面H和辅助测量平面H₁上分别均匀布置M个测量网格点，测量平面H和辅助测量平面H₁的测量网格大小和网格点位置相同；以目标声源S_o为研究对象，将声场分为非稳态向内声场和向外声场，其中非稳态向内声场来源于测量平面另一侧的干扰声源S_d，非稳态向外声场包含目标声源S_o自由场条件下辐射的非稳态声场和干扰声在目标声源表面产生的非稳态散射声场；建立笛卡尔坐标系o(x,y,z)，所述笛卡尔坐标系中的xoy坐标平面平行于测量平面H，z轴方向垂直于测量平面H；测量平面H上的测量网格点坐标为(x,y,z_H)，辅助测量平面H₁的测量网格坐标为(x,y,z_H1)；同步采集t时刻测量平面H上各测量网格点(x,y,z_H)处的时域声压信号p(x,y,z_H,t)和辅助测量平面H₁上各测量网格点(x,y,z_H1)处的时域声压信号p(x,y,z_H1,t)；

步骤2、根据Euler公式，测量平面H上的质点加速度时域信号a(x,y,z_H,t)与测量平面H上所测得的声压时域信号p(x,y,z_H,t)的偏导数存在成正比的关系，如式(9)所示：

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

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

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

步骤3、根据声波叠加原理，t时刻测量平面H上各测量网格点(x,y,z_H)处的时域声压信号p(x,y,z_H,t)和时域质点加速度信号a(x,y,z_H,t)可以表示为：

p(x,y,z_H,t)＝p_f(x,y,z_H,t)+p_s(x,y,z_H,t)+p_in(x,y,z_H,t) (12)

a(x,y,z_H,t)＝a_f(x,y,z_H,t)+a_s(x,y,z_H,t)-a_in(x,y,z_H,t) (13)

式(12)和(13)中,p_in(x,y,z_H,t)和a_in(x,y,z_H,t)分别表示干扰声源S_d在测量平面H上测量网格点(x,y,z_H)处辐射的时域声压和质点加速度，p_f(x,y,z_H,t)和a_f(x,y,z_H,t)分别表示目标声源S_o自由场条件下在测量平面H上测量网格点(x,y,z_H)处辐射的时域声压和质点加速度，下表“s”表示干扰声在目标声源S_o表面产生的非稳态散射声场,p_s(x,y,z_H,t)和a_s(x,y,z_H,t)分别表示在测量平面H上测量网格点(x,y,z_H)处的时域散射声压和质点加速度。

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

式(14)中，F(k_x,k_y,z,t)表示f(x,y,z,t)的时域波数谱，k_x、k_y分别为x、y方向的波数；

对式(12)和(13)进行关于x、y的二维空间傅里叶变换，则得到测量平面H上测量网格点(x,y,z_H)处辐射的声压时域波数谱P(k_x,k_y,z_H,t)和质点加速度时域波数谱A(k_x,k_y,z_H,t)；

P(k_x,k_y,z_H,t)＝P_f(k_x,k_y,z_H,t)+P_s(k_x,k_y,z_H,t)+P_in(k_x,k_y,z_H,t) (15)

A(k_x,k_y,z_H,t)＝A_f(k_x,k_y,z_H,t)+A_s(k_x,k_y,z_H,t)-A_in(k_x,k_y,z_H,t) (16)

式(15)和(16)中，P_f(k_x,k_y,z_H,t)和A_f(k_x,k_y,z_H,t)表示目标声源S_o自由场条件下分别在测量平面H上的声压时域波数谱和质点加速度时域波数谱，P_in(k_x,k_y,z_H,t)和A_in(k_x,k_y,z_H,t)分别表示非稳态向内声场在测量平面H上的声压时域波数谱和质点加速度时域波数谱，P_s(k_x,k_y,z_H,t)和A_s(k_x,k_y,z_H,t)分别表示在测量平面H上的散射声压时域波数谱和质点加速度时域波数谱；

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

P_f(k_x,k_y,z_H,t)＝A_f(k_x,k_y,z_H,t)*G_pa(k_x,k_y,0,t) (17)

P_s(k_x,k_y,z_H,t)＝A_s(k_x,k_y,z_H,t)*G_pa(k_x,k_y,0,t) (18)

P_in(k_x,k_y,z_H,t)＝A_in(k_x,k_y,z_H,t)*G_pa(k_x,k_y,0,t) (19)

其中，“*”表示两个时间函数的卷积运算，G_pa(k_x,k_y,0,t)为声压与质点加速度之间的时域脉冲响应函数，其表达式为：

式(20)中，c为声音的传播速度，J₀为阶数为0的第一类贝塞尔函数，H(t)为Heaviside函数；

联立式(15)-(19)，可以解得：

由声压在时域波数域的关系可知：

P_in(k_x,k_y,z₀,t)＝P_in(k_x,k_y,z_H,t)*G_pp(k_x,k_y,Δz_H0,t) (23)

P_s(k_x,k_y,z_H,t)＝P_s(k_x,k_y,z₀,t)*G_pp(k_x,k_y,Δz_H0,t) (24)

式(23)和(24)中，z₀表示目标声源S_o所在平面在z轴上的位置，Δz_H0＝z_H-z₀；G_pp(k_x,k_y,Δz_H0,t)为声压与声压之间的时域脉冲响应函数，其表达式为：

式(25)中，c为声音的传播速度，J₁为阶数为1的第一类贝塞尔函数，δ(t)为Dirac函数；

由于散射声是干扰声在目标声源表面作用产生的，将目标声源S_o表面反射系数作为边界条件，构建目标声源表面上的散射声压时域波数谱P_s(k_x,k_y,z₀,t)与非稳态向内声场在目标声源S_o表面的声压时域波数谱P_in(k_x,k_y,z₀,t)之间的关系如式(26)所示：

P_s(k_x,k_y,z₀,t)＝P_in(k_x,k_y,z₀,t)*R(k_x,k_y,t) (26)

式(26)，R(k_x,k_y,t)为目标声源表面材料的复平面波的反射系数时域波数谱。

将式(26)代入式(24)中得到：

P_s(k_x,k_y,z_H,t)＝P_in(k_x,k_y,z₀,t)*R(k_x,k_y,t)*G_pp(k_x,k_y,Δz_H0,t) (27)

将式(23)代入式(27)中得到：

P_s(k_x,k_y,z_H,t)＝P_in(k_x,k_y,z_H,t)*G_pp(k_x,k_y,Δz_H0,t)*R(k_x,k_y,t)*G_pp(k_x,k_y,Δz_H0,t) (28)

将式(28)代入式(21)中得到：

将式(22)代入式(29),则可以获得目标声源S_o自由场条件下在测量平面H上的声压时域波数谱P_f(k_x,k_y,z_H,t)的还原公式：

式(30)中，

T₁(k_x,k_y,t)＝δ(t)-G_pp(k_x,k_y,Δz_H0,t)*R(k_x,k_y,t)*G_pp(k_x,k_y,Δz_H0,t) (31)

T₂(k_x,k_y,t)＝G_pa(k_x,k_y,0,t)*[δ(t)+G_pp(k_x,k_y,Δz_H0,t)*R(k_x,k_y,t)*G_pp(k_x,k_y,Δz_H0,t)] (32)

利用式(30)按如下过程可实时还原出每个时刻下目标声源S_o自由场条件下在测量平面H上的声压时域波数谱P_f(k_x,k_y,z_H,t)：

将式(30)中的时间t离散为t_i＝(i-1)Δt，其中Δt为采样时间间隔，i＝1,2,...,I，I为采样点总数。

取i＝1时，还原得到：

取i＝2时，分离得到：

……

取i＝I时，分离得到：

由式(33)到式(35)的还原过程可知，要分离出t_i时刻的P_f(k_x,k_y,z_H,t_i)，只需要利用t_q(q＝1,2,...,i)时刻的P(k_x,k_y,z_H,t_q)和A(k_x,k_y,z_H,t_q)，而P(k_x,k_y,z_H,t_q)和A(k_x,k_y,z_H,t_q)可通过对测得的声压时域信号p(x,y,z_H,t)和通过有限差分获得的质点加速度时域信号a(x,y,z_H,t)进行二维空间傅里叶变换得到。因此，一旦测得某个时刻下测量平面H和辅助测量平面H₁上的声压时域信号，即可还原出该时刻目标声源S_o在自由场条件下辐射的声压时域波数谱，实现非稳态声场的实时还原。

步骤5、对还原的声压时域波数谱P_f(k_x,k_y,z_H,t)按照式(36)进行二维空间傅里叶反变换，最终获得每个时刻下目标声源S_o自由场条件下在测量平面H上辐射的时域声压信号p_f(x,y,z_H,t)

所述测量平面H和辅助测量平面H₁上各测量网格点的时域声压信号p(x,y,z_H,t)和p(x,y,z_H1,t)是采用声压传感器阵列一次快照测量获得。目标声源的尺寸相对于声波波长要大，非稳态散射声场影响显著。所述目标声源和干扰声源辐射的声场为非稳态，在时域内进行实时还原。

方法的检验：

在测量平面H的一侧布置大尺寸平面声源作为目标声源S_o，在另一侧布置点源作为干扰声源S_d。采用本发明方法将目标声源S_o自由场条件下在测量平面H上辐射的时域声压信号还原出来，并与其理论的时域声压信号进行比较。

本实施例中，选取半径0.2m的固定在无限大障板上的圆形活塞为目标声源S_o，活塞产生一个高斯调制正弦脉冲信号，其表达式为：

式(37)中，频率f₁＝1000Hz。

选取两个点源作为干扰声源S_d，其辐射信号的表达式都为：

s(t)＝0.3sin(2πf₂t)e^-50t (38)

式(38)中，频率f₂＝800Hz。

目标声源S_o、干扰声源S_d、测量平面H和辅助测量平面H₁在笛卡尔坐标系o(x,y,z)中的位置关系如图2所示，目标声源S_o的圆形活塞中心位于(0.25m,0.25m,0m)处，干扰声源S_d的两个点源分别位于(0.25m,0.25m,0.1m)和(0.35m,0.15m,0.1m)处，目标声源S_o所在平面为z₀＝0m，测量平面H和H₁分别位于z_H＝0.04m和z_H1＝0.06m。测量平面H和H₁上均匀分布11×11个测量点，相邻点在x和y方向上的间隔都是0.05m，如图2所示。时域信号采样频率为25.6kHz，采样时间为10ms，且仿真信号中加入信噪比30dB的高斯白噪声。仿真中，固定在无限大障板上的圆形活塞所产生的理论时域声压通过脉冲响应法计算获得，同时假设目标声源S_o表面为刚性，因此干扰声源S_d在目标声源表面产生的时域散射声压可以利用镜像原理计算获得，刚性特性使得表面反射系数时域波数谱R(k_x,k_y,t)取为δ(t)。

为检验本发明方法在时域内的声场还原效果，在测量平面H上选取了四个测量点，即测量点A、测量点B、测量点C和测量点D，其位置分别为A(0.15m,0.5m,0.04m)、B(0.1m,0.25m,0.04m)、C(0.25m,0.25m,0.04m)和D(0.2m,0m,0.04m)。图3(a)、图3(b)、图3(c)和图3(d)分别对应测量点A、测量点B、测量点C和测量点D，图中曲线a表示目标声源S_o自由场条件下在测量平面H上辐射的理论时域声压信号，图中曲线b表示包含目标声源S_o、干扰声源S_d和干扰声在目标声源表面产生的散射声辐射的混合时域声压信号，图中曲线c表示利用实时声场分离方法分离的时域声压信号，图中曲线d表示采用本发明方法还原的时域声压信号。比较图中的曲线a和曲线b可以看出，干扰声源S_d和散射声对目标声源S_o在测量平面H上辐射的时域声压信号造成了较大的干扰；比较曲线a和曲线c可以看出，利用实时声场分离方法仅能够消除干扰声源S_d的影响，无法消除散射声的影响；比较曲线a和曲线d可以看出，采用本发明方法可以进一步很好地消除散射声的影响，从而彻底地还原出目标声源S_o自由场条件下在测量平面H上辐射的时域声压信号。

为检验本发明方法在空间域上的声场还原效果，选取了两个时刻t＝3.40ms和t＝6.88ms。图4(a)和图4(e)分别为3.40ms和6.88ms时刻目标声源S_o自由场条件下在测量平面H上辐射的理论声压场p_m，图4(b)和图4(f)分别为3.40ms和6.88ms时刻的混合声压场p_c(包含目标声源S_o、干扰声源S_d和干扰声在目标声源表面产生的散射声所辐射的声压场)，图4(c)和图4(g)分别为3.40ms和6.88ms时刻利用实时声场分离方法分离的声压场p_e，图4(d)和图4(h)分别为3.40ms和6.88ms时刻采用本发明方法还原的声压场p_r。比较图4(a)和图4(b)、图4(e)和图4(f)可以看出，干扰声源S_d和散射声对目标声源S_o自由场条件下在测量平面H上的声压空间分布造成了严重的改变；比较图4(a)和图4(c)、图4(e)和图4(g)可以看出，利用实时声场分离方法仅能够去除干扰声源S_d对目标声源声压空间分布的改变，也表明散射声的影响不可以忽略；比较图4(a)和图4(d)、图4(e)和图4(h)可以看出，采用本发明方法可以进一步消除散射声对目标声源声压空间分布的改变，从而彻底还原出每个时刻目标声源S_o自由场条件下在测量平面H上辐射的声压场。

为了量化评价本发明方法的还原效果，相位评价因子E_p和幅值评价因子E_a被引入，其分别定义为：

式(39)和(40)中，p_m(x,y,z,t)表示目标声源S_o自由场条件下在测量平面H上辐射的理论时域声压信号，p_D(x,y,z,t)表示利用实时声场分离方法分离的时域声压信号或采用本发明方法还原的时域声压信号，＜＞表示求平均值。评价因子E_p是用来衡量目标声源自由场条件下辐射的理论声压值和分离声压值或还原声压值之间的相位误差，当E_p的值越靠近1时，相位误差越小。评价因子E_a是用来衡量目标声源自由场条件下辐射的理论声压值和分离声压值或还原声压值之间的幅值误差，当E_a的值越靠近0时，幅值误差越小。运用式(39)和(40)分别计算了实时声场分离方法和本发明方法在测量平面H上各测量点处的E_p和E_a值。在图5中，曲线a表示本发明方法在各测量点处的E_p值，曲线b表示本发明方法在各测量点处的E_a值，曲线c表示实时声场分离方法在各测量点处的E_p值，曲线d表示实时声场分离方法在各测量点处的E_a值。参见图5中曲线c和曲线d，可以看出利用实时声场分离方法分离的声压值与理论声压值存在较大的差别，主要因为实时声场分离方法无法消除非稳态散射声的影响。参见图5中曲线a和曲线b，可以看出在大多数测量点处，E_p值都大于0.95，E_a值都小于0.12，表明本发明方法还原的声压值和理论声压值吻合地很好，本发明方法可以有效地消除非稳态散射声的影响。

上述实施例表明，采用本发明方法不仅可以消除来自测量面另一侧干扰声源的影响，而且可以进一步消除非稳态散射声的影响，从而彻底还原出目标声源S_o自由场条件下在测量平面H上辐射的时域声压信号。

提供以上实施例仅仅是为了描述本发明的目的，而并非要限制本发明的范围。本发明的范围由所附权利要求限定。不脱离本发明的精神和原理而做出的各种等同替换和修改，均应涵盖在本发明的范围之内。

Claims (4)

1.一种采用声压和质点加速度测量的非稳态平面声源的自由场实时还原方法，其特征在于，按如下步骤进行：

步骤1、在目标声源S_o和干扰声源S_d之间布置相互平行的测量平面H和辅助测量平面H₁；所述的目标声源S_o和干扰声源S_d均可辐射任意线性非稳态声场信号；在所述测量平面H和辅助测量平面H₁上分别均匀布置M个测量网格点，测量平面H和辅助测量平面H₁的测量网格大小和网格点位置相同；以目标声源S_o为研究对象，将声场分为非稳态向内声场和向外声场，其中非稳态向内声场来源于测量平面另一侧的干扰声源S_d，非稳态向外声场包含目标声源S_o自由场条件下辐射的非稳态声场和干扰声在目标声源表面产生的非稳态散射声场；建立笛卡尔坐标系o(x,y,z)，所述笛卡尔坐标系中的xoy坐标平面平行于测量平面H，z轴方向垂直于测量平面H；测量平面H上的测量网格点坐标为(x,y,z_H)，辅助测量平面H₁的测量网格坐标为(x,y,z_H1)；同步采集t时刻测量平面H上各测量网格点(x,y,z_H)处的时域声压信号p(x,y,z_H,t)和辅助测量平面H₁上各测量网格点(x,y,z_H1)处的时域声压信号p(x,y,z_H1,t)；

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

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

步骤3、对于测量平面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)，

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

步骤4、构建测量平面H上声压时域波数谱P(k_x,k_y,z_H,t)与质点加速度时域波数谱A(k_x,k_y,z_H,t)、已知的时域响应函数T₁(k_x,k_y,t)与T₂(k_x,k_y,t)和目标声源S_o自由场条件下在测量平面H上的声压时域波数谱P_f(k_x,k_y,z_H,t)之间的关系如式(4)所示：

式(4)中，

T₁(k_x,k_y,t)＝δ(t)-G_pp(k_x,k_y,Δz_H0,t)*R(k_x,k_y,t)*G_pp(k_x,k_y,Δz_H0,t) (5)

T₂(k_x,k_y,t)＝G_pa(k_x,k_y,0,t)*[δ(t)+G_pp(k_x,k_y,Δz_H0,t)*R(k_x,k_y,t)*G_pp(k_x,k_y,Δz_H0,t)] (6)

G_pp(k_x,k_y,Δz_H0,t)和G_pa(k_x,k_y,0,t)表示已知的时域脉冲响应函数，R(k_x,k_y,t)表示已知的目标声源表面材料的复平面波反射系数时域波数谱，z₀表示目标声源S_o表面所在位置，“*”表示两个时间函数的卷积运算，Δz_H0＝z_H-z₀，δ(t)为Dirac函数；

将式(4)中的时间t离散为t_i＝(i-1)Δt，其中Δt为采样时间间隔，i＝1,2,...,I，I为采样点总数，则式(4)的离散形式如式(7)所示：

步骤5、对于所述声压时域波数谱P_f(k_x,k_y,z_H1,t)按照式(8)进行二维空间傅里叶反变换，还原出每个时刻目标声源S_o自由场条件下在测量平面H上辐射的时域声压信号p_f(x,y,z_H,t)；

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

3.根据权利要求1所述的采用声压和质点加速度测量的非稳态平面声源的自由场实时还原方法，其特征在于：所述的目标声源的尺寸相对于声波波长要大，非稳态散射声场影响显著。

4.根据权利要求1所述的采用声压和质点加速度测量的非稳态平面声源的自由场实时还原方法，其特征在于：所述目标声源和干扰声源辐射的声场为非稳态，在时域内进行实时还原。