CN107478325A

CN107478325A - 一种非稳态平面声源的自由场还原方法

Info

Publication number: CN107478325A
Application number: CN201710798515.4A
Authority: CN
Inventors: 耿林
Original assignee: Anhui University
Current assignee: Anhui University
Priority date: 2017-09-07
Filing date: 2017-09-07
Publication date: 2017-12-15
Anticipated expiration: 2037-09-07
Also published as: CN107478325B

Abstract

本发明公开了一种非稳态平面声源的自由场还原方法，将两个测量平面上的时域声压信号进行二维空间傅里叶变换获得相应面上各波数处的声压时域波数谱；利用两个测量平面上的声压时域波数谱和已知的声压与声压之间的时域脉冲响应函数，将非稳态向外声场和向内声场分离出来；以目标声源表面反射系数为边界条件，建立非稳态向内声场和散射声场之间的关系，并计算出非稳态散射声场；从分离的非稳态向外声场中去除非稳态散射声场，还原出目标声源自由场条件下在测量平面上辐射的时域声压信号。本发明在时域上实现了非稳态散射声的消除，还原出目标声源自由场条件下在测量平面上辐射的时域声压，为实际声场环境下精确分析目标声源的时变辐射特性和振动特性提供了前处理手段。

Description

一种非稳态平面声源的自由场还原方法

技术领域

本发明涉及物理专业中噪声类领域，具体涉及一种非稳态平面声源的自由场还原方法。

背景技术

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

到目前为止，国内外学者已提出多种自由场还原方法：Hu等推导了基于二维空间傅里叶变换的平面近场声全息的完备计算公式，它包含了干扰声和散射声的影响，从而实现了目标声场的准确重建；Hald等和Fernandez-Grande等分别提出了基于双面声压测量和单面声压与振速测量的统计最优近场声全息的自由场还原技术，并分别实现了内空间真实辐射声强的重建和在近似表面刚性边界条件下消除散射声影响；Braikia等提出了基于球面波叠加法的自由场还原技术，并结合镜像原理消除了散射声的影响，实现了小内空间的声场重建；Langrenne等提出了基于边界元法的自由场还原技术，在近似表面刚性边界条件下还原出目标声源自由场条件下的辐射声场；Bi等提出了基于等效源法的自由场还原技术，实现了任意阻抗表面的自由场还原和导纳边界下的腔体内主动声源识别。在上述自由场还原方法中，还原的对象都是稳态声源，无法用于非稳态声源。非稳态声场的还原还需要在时域上消除干扰声和散射声的影响，还原出目标声源随时间变化的声学特性。

为此，Zhang等在2012年提出了基于双面声压测量的时域声场分离方法，在时域和空间域上消除了干扰声源的影响。为了减少计算时间，Bi等在2014年提出了基于声压和加速度测量的实时声场分离方法。为了避免有限差分引起的误差，并提高分离精度，Bi等提出了基于单层声压和振速测量的非稳态声场分离方法。

但是这些时域分离方法仅仅能够去除测量平面另一侧干扰声源产生的非稳态声场，无法消除干扰声在目标声源表面产生的非稳态散射声场。

发明内容

为了消除非稳态散射声场，本发明提供了一种采用双面声压测量的非稳态平面声源的自由场还原方法，以期获得还原非稳态目标声源自由场条件下辐射的非稳态声场的能力，进而为实际声场环境下精确分析目标声源的时变辐射特性和振动特性提供前处理手段。

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

一种非稳态平面声源的自由场还原方法，步骤如下：

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

步骤2、对于测量平面H₁和测量平面H₂上的时域声压信号p(x,y,z_H1,t)和p(x,y,z_H2,t)分别按照式(1)和式(2)进行二维空间傅里叶变换获得测量平面H₁和测量平面H₂上的声压时域波数谱P(k_x,k_y,z_H1,t)和P(k_x,k_y,z_H2,t)：

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

步骤3、构建测量平面H₁上声压时域波数谱P(k_x,k_y,z_H1,t)、测量平面H₂上声压时域波数谱P(k_x,k_y,z_H2,t)、已知的时域脉冲响应函数G(k_x,k_y,Δz₂₁,t)和非稳态向内声场在测量平面H₂上的声压时域波数谱P_in(k_x,k_y,z_H2,t)之间的关系如式(3)所示：

W1(k_x,k_y,t)＝P_in(k_x,k_y,z_H2,t)*W2(k_x,k_y,t) (3)

式(3)中，

W1(k_x,k_y,t)＝P(k_x,k_y,z_H1,t)*G(k_x,k_y,Δz₂₁,t)-P(k_x,k_y,z_H2,t) (4)

W2(k_x,k_y,t)＝G(k_x,k_y,Δz₂₁,t)*G(k_x,k_y,Δz₂₁,t)-δ(t) (5)

“*”表示两个时间函数的卷积运算，Δz₂₁＝z_H2-z_H1，δ(t)为Dirac函数；

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

对于每个波数对(k_x,k_y)，式(6)可以写成式(7)的矩阵形式：

W1＝W2P_in (7)

式(7)中，

W1＝[W1(k_x,k_y,t₁) W1(k_x,k_y,t₂) … W1(k_x,k_y,t_N)]^T (8)

P_in＝[P_in(k_x,k_y,z_H2,t₁) P_in(k_x,k_y,z_H2,t₂) … P_in(k_x,k_y,z_H2,t_N)]^T (9)

“T”表示矩阵的转置；

式(7)中，W1和W2是已知的，可以通过式(11)奇异值分解求逆解得P_in

P_in＝[W2]⁺W1 (11)

式(11)中，“+”表示矩阵的伪逆；

最终获得非稳态向内声场在测量平面H₂上的声压时域波数谱P_in(k_x,k_y,z_H2,t)；非稳态向外声场在测量平面H₁上的声压时域波数谱P_out(k_x,k_y,z_H1,t)按照式(12)计算：

P_out(k_x,k_y,z_H1,t)＝P(k_x,k_y,z_H1,t)-P_in(k_x,k_y,z_H2,t)*G(k_x,k_y,Δz₂₁,t) (12)

步骤4、以目标声源S_o表面反射系数为边界条件，构建非稳态散射声场在测量平面H₁上的散射声压时域波数谱P_s(k_x,k_y,z_H1,t)与非稳态向内声场在测量平面H₂上的声压时域波数谱P_in(k_x,k_y,z_H2,t)之间的关系如式(13)所示：

P_s(k_x,k_y,z_H1,t)＝P_in(k_x,k_y,z_H2,t)*G(k_x,k_y,Δz₂₀,t)*R(k_x,k_y,t)*G(k_x,k_y,Δz₁₀,t) (13)

式(13)中，R(k_x,k_y,t)为目标声源表面材料的复平面波反射系数时域波数谱，Δz₂₀＝z_H2-z₀，Δz₁₀＝z_H1-z₀，z₀表示目标声源S_o表面所在位置；

将式(12)与式(13)相减，还原出目标声源S_o自由场条件下在测量平面H₁上的声压时域波数谱P_f(k_x,k_y,z_H1,t)如式(14)所示：

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

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

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

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

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

(1)本发明方法直接在时域上进行自由场还原计算。

(2)本发明方法可以进一步消除干扰声在目标声源表面产生的非稳态散射声场，彻底还原出目标声源自由场条件下辐射的非稳态声场，因而更适用于精确分析实际声场环境下目标声源的时变辐射特性和振动特性。

附图说明

图1为本发明非稳态平面声源的自由场还原方法示意图；

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

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

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

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

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

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

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

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

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

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

图4(e)为t＝4.69ms时刻利用时域声场分离技术分离的声压场p_e；

图4(f)为t＝6.99ms时刻利用时域声场分离技术分离的声压场p_e；

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

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

图5为测量平面H₁上本发明方法和时域声场分离技术的相位评价因子T₁值和幅值评价因子T₂值。

具体实施方式

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

本发明一种非稳态平面声源的自由场还原方法，其特征是在目标声源和干扰声源之间布置两个测量平面，同步采集两测量平面上各网格点处的时域声压信号；将两个测量平面上的时域声压信号进行二维空间傅里叶变换获得相应面上各波数处的声压时域波数谱；然后利用两个测量平面上的声压时域波数谱和已知的声压与声压之间的时域脉冲响应函数，将非稳态向外声场(包含目标声源自由场条件下辐射的非稳态声场和干扰声在目标声源表面产生的非稳态散射声场)和向内声场(测量面另一侧干扰声源辐射的非稳态声场)分离出来；再以目标声源表面反射系数为边界条件，建立非稳态向内声场和散射声场之间的关系，从而计算出非稳态散射声场；最后从分离的向外声场中去除非稳态散射声场，从而彻底还原出目标声源自由场条件下在测量平面上辐射的时域声压信号。本发明方法直接在时域上实现了非稳态散射声的消除，进一步还原出目标声源自由场条件下在测量平面上辐射的时域声压，可为实际声场环境下精确分析目标声源的时变辐射特性和振动特性提供前处理手段。

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

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

p(x,y,z_H1,t)＝p_out(x,y,z_H1,t)+p_in(x,y,z_H1,t) (16)

p(x,y,z_H2,t)＝p_out(x,y,z_H2,t)+p_in(x,y,z_H2,t) (17)

式(16)和(17)中,

p_out(x,y,z_H1,t)＝p_f(x,y,z_H1,t)+p_s(x,y,z_H1,t) (18)

p_out(x,y,z_H2,t)＝p_f(x,y,z_H2,t)+p_s(x,y,z_H2,t) (19)

p_out(x,y,z_H1,t)和p_out(x,y,z_H2,t)表示非稳态向外声场分别在测量平面H₁和H₂上测量网格点(x,y,z_H1)和(x,y,z_H2)处的时域声压，p_in(x,y,z_H1,t)和p_in(x,y,z_H2,t)分别表示干扰声源S_d在测量平面H₁和H₂上测量网格点(x,y,z_H1)和(x,y,z_H2)处辐射的时域声压，p_f(x,y,z_H1,t)和p_f(x,y,z_H2,t)表示目标声源S_o自由场条件下分别在测量平面H₁和H₂上测量网格点(x,y,z_H1)和(x,y,z_H2)处辐射的时域声压，下表“s”表示干扰声在目标声源S_o表面产生的非稳态散射声场,p_s(x,y,z_H1,t)和p_s(x,y,z_H2,t)分别表示在测量平面H₁和H₂上测量网格点(x,y,z_H1)和(x,y,z_H2)处的时域散射声压；

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

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

对式(16)和(17)进行关于x、y的二维空间傅里叶变换，则得到测量平面H₁上各测量网格点(x,y,z_H1)处的声压时域波数谱P(k_x,k_y,z_H1,t)和测量平面H₂上各测量网格点(x,y,z_H2)处的声压时域波数谱P(k_x,k_y,z_H2,t)

P(k_x,k_y,z_H1,t)＝P_f(k_x,k_y,z_H1,t)+P_s(k_x,k_y,z_H1,t)+P_in(k_x,k_y,z_H1,t) (21)

P(k_x,k_y,z_H2,t)＝P_f(k_x,k_y,z_H2,t)+P_s(k_x,k_y,z_H2,t)+P_in(k_x,k_y,z_H2,t) (22)

式(21)和(22)中，P_f(k_x,k_y,z_H1,t)和P_f(k_x,k_y,z_H2,t)表示目标声源S_o自由场条件下分别在测量平面H₁和H₂上的声压时域波数谱，P_in(k_x,k_y,z_H1,t)和P_in(k_x,k_y,z_H2,t)分别表示非稳态向内声场在测量平面H₁和H₂上的声压时域波数谱，P_s(k_x,k_y,z_H1,t)和P_s(k_x,k_y,z_H2,t)分别表示在测量平面H₁和H₂上的散射声压时域波数谱；

步骤2、由声压在时域波数域的关系可知：

P_f(k_x,k_y,z_H2,t)＝P_f(k_x,k_y,z_H1,t)*G(k_x,k_y,Δz₂₁,t) (23)

P_s(k_x,k_y,z_H2,t)＝P_s(k_x,k_y,z_H1,t)*G(k_x,k_y,Δz₂₁,t) (24)

P_in(k_x,k_y,z_H1,t)＝P_in(k_x,k_y,z_H2,t)*G(k_x,k_y,Δz₂₁,t) (25)

其中，“*”表示两个时间函数的卷积运算，Δz₂₁＝z_H2-z_H1，G(k_x,k_y,Δz₂₁,t)为声压与声压之间的时域脉冲响应函数，其表达式为：

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

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

P(k_x,k_y,z_H1,t)＝P_f(k_x,k_y,z_H1,t)+P_s(k_x,k_y,z_H1,t)+P_in(k_x,k_y,z_H2,t)*G(k_x,k_y,Δz₂₁,t) (27)

将式(27)两边分别与时域脉冲响应函数G(k_x,k_y,Δz₂₁,t)进行卷积，则得到：

将式(23)和(24)代入式(28)中得到：

将式(29)与式(22)相减可进一步得到：

定义

W1(k_x,k_y,t)＝P(k_x,k_y,z_H1,t)*G(k_x,k_y,Δz₂₁,t)-P(k_x,k_y,z_H2,t) (31)

W2(k_x,k_y,t)＝G(k_x,k_y,Δz₂₁,t)*G(k_x,k_y,Δz₂₁,t)-δ(t) (32)

则式(30)可以进一步表示为：

W1(k_x,k_y,t)＝P_in(k_x,k_y,z_H2,t)*W2(k_x,k_y,t) (33)

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

对于每个波数对(k_x,k_y)，式(34)可以写成式(35)的矩阵形式：

W1＝W2P_in (35)

式(35)中，

W1＝[W1(k_x,k_y,t₁) W1(k_x,k_y,t₂) … W1(k_x,k_y,t_N)]^T (36)

P_in＝[P_in(k_x,k_y,z_H2,t₁) P_in(k_x,k_y,z_H2,t₂) … P_in(k_x,k_y,z_H2,t_N)]^T (37)

“T”表示矩阵的转置；

式(35)中，W1和W2是已知的，可以通过式(39)奇异值分解求逆解得P_in

P_in＝[W2]⁺W1 (39)

式(39)中，“+”表示矩阵的伪逆；

由于求解声压时域波数谱P_in的过程一般是病态的，因此为了获得合适的解，在每对波数的求解公式(39)中需要采用标准Tikhonov正则化技术，其中正则化参数可以通过广义交叉验证方法来选取。最终获得非稳态向内声场在测量平面H₂上的声压时域波数谱P_in(k_x,k_y,z_H2,t)；

步骤3、由声压在时域波数域的关系可知：

P_s(k_x,k_y,z_H1,t)＝P_s(k_x,k_y,z₀,t)*G(k_x,k_y,Δz₁₀,t) (40)

P_in(k_x,k_y,z₀,t)＝P_in(k_x,k_y,z_H2,t)*G(k_x,k_y,Δz₂₀,t) (41)

式(40)和(41)中，z₀表示目标声源S_o所在平面在z轴上的位置，Δz₁₀＝z_H1-z₀；Δz₂₀＝z_H2-z₀，G(k_x,k_y,Δz₁₀,t)和G(k_x,k_y,Δz₂₀,t)为声压与声压之间的时域脉冲响应函数，其表达式分别为：

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

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

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

将式(41)代入式(44)中得到：

P_s(k_x,k_y,z₀,t)＝P_in(k_x,k_y,z_H2,t)*G(k_x,k_y,Δz₂₀,t)*R(k_x,k_y,t) (45)

将式(45)代入式(40)中得到：

P_s(k_x,k_y,z_H1,t)＝P_in(k_x,k_y,z_H2,t)*G(k_x,k_y,Δz₂₀,t)*R(k_x,k_y,t)*G(k_x,k_y,Δz₁₀,t) (46)

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

将式(47)中P_f(k_x,k_y,z_H1,t)往等式左边移，P(k_x,k_y,z_H1,t)往等式右边移，则可以获得目标声源S_o自由场条件下在测量平面H₁上的声压时域波数谱P_f(k_x,k_y,z_H1,t)的还原公式：

式(48)中，P_in(k_x,k_y,z_H2,t)是通过步骤2求解获得；

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

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

方法的检验：

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

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

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

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

s(t)＝0.2sin(2πf₂t)e^-100t (51)

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

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

为检验本发明方法在时域内的声场还原效果，在测量平面H₁上选取了四个测量点，即测量点A、测量点B、测量点C和测量点D，其位置分别为A(0.2m,0.7m,0.04m)、B(0.1m,0.5m,0.04m)、C(0.65m,0.2m,0.04m)和D(0.35m,0.05m,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＝4.69ms和t＝6.99ms。图4(a)和图4(e)分别为4.69ms和6.99ms时刻目标声源S_o自由场条件下在测量平面H₁上辐射的理论声压场p_m，图4(b)和图4(f)分别为4.69ms和6.99ms时刻的混合声压场p_c(包含目标声源S_o、干扰声源S_d和干扰声在目标声源表面产生的散射声所辐射的声压场)，图4(c)和图4(g)分别为4.69ms和6.99ms时刻利用时域声场分离技术分离的声压场p_e，图4(d)和图4(h)分别为4.69ms和6.99ms时刻采用本发明方法还原的声压场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₁上辐射的声压场。

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

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

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

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

Claims

1.一种非稳态平面声源的自由场还原方法，其特征在于，包括如下步骤；

<mrow> <mi>P</mi> <mrow> <mo>(</mo> <msub> <mi>k</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>k</mi> <mi>y</mi> </msub> <mo>,</mo> <msub> <mi>z</mi> <mrow> <mi>H</mi> <mn>1</mn> </mrow> </msub> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <msubsup> <mo>&Integral;</mo> <mrow> <mo>-</mo> <mi>&infin;</mi> </mrow> <mi>&infin;</mi> </msubsup> <msubsup> <mo>&Integral;</mo> <mrow> <mo>-</mo> <mi>&infin;</mi> </mrow> <mi>&infin;</mi> </msubsup> <mi>p</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <msub> <mi>z</mi> <mrow> <mi>H</mi> <mn>1</mn> </mrow> </msub> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <msup> <mi>e</mi> <mrow> <mi>j</mi> <mrow> <mo>(</mo> <msub> <mi>k</mi> <mi>x</mi> </msub> <mi>x</mi> <mo>+</mo> <msub> <mi>k</mi> <mi>y</mi> </msub> <mi>y</mi> <mo>)</mo> </mrow> </mrow> </msup> <mi>d</mi> <mi>x</mi> <mi>d</mi> <mi>y</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow>

<mrow> <mi>P</mi> <mrow> <mo>(</mo> <msub> <mi>k</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>k</mi> <mi>y</mi> </msub> <mo>,</mo> <msub> <mi>z</mi> <mrow> <mi>H</mi> <mn>2</mn> </mrow> </msub> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <msubsup> <mo>&Integral;</mo> <mrow> <mo>-</mo> <mi>&infin;</mi> </mrow> <mi>&infin;</mi> </msubsup> <msubsup> <mo>&Integral;</mo> <mrow> <mo>-</mo> <mi>&infin;</mi> </mrow> <mi>&infin;</mi> </msubsup> <mi>p</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <msub> <mi>z</mi> <mrow> <mi>H</mi> <mn>2</mn> </mrow> </msub> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <msup> <mi>e</mi> <mrow> <mi>j</mi> <mrow> <mo>(</mo> <msub> <mi>k</mi> <mi>x</mi> </msub> <mi>x</mi> <mo>+</mo> <msub> <mi>k</mi> <mi>y</mi> </msub> <mi>y</mi> <mo>)</mo> </mrow> </mrow> </msup> <mi>d</mi> <mi>x</mi> <mi>d</mi> <mi>y</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </mrow>

W1(k_x,k_y,t)＝P_in(k_x,k_y,z_H2,t)*W2(k_x,k_y,t) (3)

式(3)中，

W1(k_x,k_y,t)＝P(k_x,k_y,z_H1,t)*G(k_x,k_y,Δz₂₁,t)-P(k_x,k_y,z_H2,t) (4)

W2(k_x,k_y,t)＝G(k_x,k_y,Δz₂₁,t)*G(k_x,k_y,Δz₂₁,t)-δ(t) (5)

<mrow> <mi>W</mi> <mn>1</mn> <mrow> <mo>(</mo> <msub> <mi>k</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>k</mi> <mi>y</mi> </msub> <mo>,</mo> <msub> <mi>t</mi> <mi>n</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>n</mi> </munderover> <msub> <mi>P</mi> <mrow> <mi>i</mi> <mi>n</mi> </mrow> </msub> <mrow> <mo>(</mo> <msub> <mi>k</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>k</mi> <mi>y</mi> </msub> <mo>,</mo> <msub> <mi>z</mi> <mrow> <mi>H</mi> <mn>2</mn> </mrow> </msub> <mo>,</mo> <msub> <mi>t</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mi>W</mi> <mn>2</mn> <mrow> <mo>(</mo> <msub> <mi>k</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>k</mi> <mi>y</mi> </msub> <mo>,</mo> <msub> <mi>t</mi> <mrow> <mi>n</mi> <mo>+</mo> <mn>1</mn> <mo>-</mo> <mi>i</mi> </mrow> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </mrow>

对于每个波数对(k_x,k_y)，式(6)可以写成式(7)的矩阵形式：

W1＝W2P_in (7)

式(7)中，

W1＝[W1(k_x,k_y,t₁) W1(k_x,k_y,t₂) … W1(k_x,k_y,t_N)]^T (8)

“T”表示矩阵的转置；

P_in＝[W2]⁺W1 (11)

式(11)中，“+”表示矩阵的伪逆；

<mrow> <mtable> <mtr> <mtd> <mrow> <msub> <mi>P</mi> <mi>f</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>k</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>k</mi> <mi>y</mi> </msub> <mo>,</mo> <msub> <mi>z</mi> <mrow> <mi>H</mi> <mn>1</mn> </mrow> </msub> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <mi>P</mi> <mrow> <mo>(</mo> <msub> <mi>k</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>k</mi> <mi>y</mi> </msub> <mo>,</mo> <msub> <mi>z</mi> <mrow> <mi>H</mi> <mn>1</mn> </mrow> </msub> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>-</mo> <mo>{</mo> <msub> <mi>P</mi> <mrow> <mi>i</mi> <mi>n</mi> </mrow> </msub> <mrow> <mo>(</mo> <msub> <mi>k</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>k</mi> <mi>y</mi> </msub> <mo>,</mo> <msub> <mi>z</mi> <mrow> <mi>H</mi> <mn>2</mn> </mrow> </msub> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>*</mo> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mo>&lsqb;</mo> <mi>G</mi> <mrow> <mo>(</mo> <msub> <mi>k</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>k</mi> <mi>y</mi> </msub> <mo>,</mo> <msub> <mi>&Delta;z</mi> <mn>20</mn> </msub> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>*</mo> <mi>R</mi> <mrow> <mo>(</mo> <msub> <mi>k</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>k</mi> <mi>y</mi> </msub> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>*</mo> <mi>G</mi> <mrow> <mo>(</mo> <msub> <mi>k</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>k</mi> <mi>y</mi> </msub> <mo>,</mo> <msub> <mi>&Delta;z</mi> <mn>10</mn> </msub> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>+</mo> <mi>G</mi> <mrow> <mo>(</mo> <msub> <mi>k</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>k</mi> <mi>y</mi> </msub> <mo>,</mo> <msub> <mi>&Delta;z</mi> <mn>21</mn> </msub> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>&rsqb;</mo> <mo>}</mo> </mrow> </mtd> </mtr> </mtable> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>14</mn> <mo>)</mo> </mrow> </mrow>

<mrow> <msub> <mi>p</mi> <mi>f</mi> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <msub> <mi>z</mi> <mrow> <mi>H</mi> <mn>1</mn> </mrow> </msub> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <msup> <mrow> <mo>(</mo> <mn>2</mn> <mi>&pi;</mi> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mfrac> <msubsup> <mo>&Integral;</mo> <mrow> <mo>-</mo> <mi>&infin;</mi> </mrow> <mi>&infin;</mi> </msubsup> <msubsup> <mo>&Integral;</mo> <mrow> <mo>-</mo> <mi>&infin;</mi> </mrow> <mi>&infin;</mi> </msubsup> <msub> <mi>P</mi> <mi>f</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>k</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>k</mi> <mi>y</mi> </msub> <mo>,</mo> <msub> <mi>z</mi> <mrow> <mi>H</mi> <mn>1</mn> </mrow> </msub> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <msup> <mi>e</mi> <mrow> <mo>-</mo> <mi>j</mi> <mrow> <mo>(</mo> <msub> <mi>k</mi> <mi>x</mi> </msub> <mi>x</mi> <mo>+</mo> <msub> <mi>k</mi> <mi>y</mi> </msub> <mi>y</mi> <mo>)</mo> </mrow> </mrow> </msup> <msub> <mi>dk</mi> <mi>x</mi> </msub> <msub> <mi>dk</mi> <mi>y</mi> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>15</mn> <mo>)</mo> </mrow> <mo>.</mo> </mrow>

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

3.根据权利要求1所述的非稳态平面声源的自由场还原方法，其特征在于：所述的目标声源的尺寸相对于声波波长要大，非稳态散射声场影响显著。

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