CN112146751A - 一种基于时域等效源法的实时声场分离方法 - Google Patents

Abstract

本发明公开了一种基于时域等效源法的实时声场分离方法，是在目标源和干扰源之间布置一个测量面，同步采集测量面上各测点处的时域声压和时域质点振速；在测量面两侧布置等效源，用于分别模拟目标源声场和干扰源声场；采用线性插值函数对各等效源的时域源强积分进行插值处理，建立任意时刻各测点的声压和质点振速与各等效源源强之间的传递关系；从初始时刻开始依次获得对应时刻测量面两侧各等效源的时域源强积分；根据获得的测量面两侧的各时刻等效源的时域源强积分以及等效源与任意测点之间的传递关系，计算获得目标源声场和干扰源声场，实现实时声场分离。本发明允许测量面为任意形状，能够在混响环境中分离出任意形状目标声源的辐射声场。

Description

一种基于时域等效源法的实时声场分离方法

技术领域

本发明涉及噪声领域声场分离方法，更具体地说是一种基于时域等效源法的实时声场分离方法。

背景技术

在实际工程中，声学测量环境往往不是消声室，而是处于混响环境中，混响环境所产生的干扰声影响了对目标声源所辐射声场的准确测量，因此需要采用声场分离技术将干扰声源的影响从测量结果中分离出去。声场分离技术最初发展于在频域中分离出测量面两侧的简谐声场，其方法主要有：空间傅里叶变换方法、统计最优近场声全息算法、边界元法、等效源法。但是，上述频域声场分离方法都仅仅实现了单一频率下的声场分离。对于随时间变化的瞬态声场，上述声场分离方法将不再适用，而应在每个一瞬时分离出目标声源声场。为实现任意时刻下的声场分离，张小正等于2012年在J.Acoust.Soc.Am.131(3),2180–2189提出了一种基于双平面声压测量的时域声场分离方法，在该方法中，所有时刻的目标声场都通过解卷计算来一次性求解，其中必须使用正则化来改善求解过程中的病态问题。为了避免解卷计算中的病态问题，毕传兴等于2014年在J.Acoust.Soc.Am.135(6),3474–3482提出了一种基于声压和质点加速度测量的时域声场分离方法，该方法在某一时刻的目标声场是由同一时刻的测量声压和之前所有时刻的质点加速度与对应的脉冲响应函数之间卷积的简单叠加得到，因此该方法能实现声场的实时分离。但是，由于质点加速度是通过两个相邻平面上测得声压的有限差分近似得到，而有限差分引起的质点加速度误差会导致错误的分离结果，毕传兴等于2016年在J.Acoust.Soc.Am.139(2),781–789又进一步提出了一种采用单面声压和质点振速测量的时域声场分离方法。在该方法中，声压和质点振速是采用声压和质点振速传感器直接采集得到，因此不存在因有限差分近似而产生的误差。上述的时域声场分离方法均基于平面波展开，这需要测量面为一个平面或两个平行平面。Ma等人于2019年在J.Acoust.Soc.Am.146(1),11–21提出了一种基于球面波展开法的可在混响环境中分离出球形声源辐射的目标声场的时域声场分离方法。该方法需要采用包围目标声源的球形测量面上的声压和法向质点振速来进行声场分离，且法向质点振速是由两个相距很近的球面上声压的有限差分近似。然而，当声源为圆柱形或椭圆柱形等不规则形状时，球形测量面将无法有效分离目标声场，且有限差分近似引起的质点振速的误差也会影响分离精度。

发明内容

本发明是为避免上述现有技术所存在的不足，提供一种基于时域等效源法的实时声场分离方法，允许测量面为任意形状，且能对瞬态声场实现实时分离。

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

本发明基于时域等效源法的实时声场分离方法的特点是：在目标源和干扰源之间布置一个测量面，所述目标源是指目标声源，所述干扰源是指干扰声源，同步采集测量面上各测点处的时域声压和时域质点振速；在测量面两侧布置若干等效源，用于分别模拟目标源声场和干扰源声场；采用线性插值函数对各等效源的时域源强积分进行插值处理，建立任意时刻各测点的声压和质点振速与各等效源源强之间的传递关系；从初始时刻开始依次获得对应时刻测量面两侧各等效源的时域源强积分；根据获得的测量面两侧的各时刻等效源的时域源强积分以及等效源与任意测点之间的传递关系，计算获得目标源声场和干扰源声场，实现实时声场分离；所述测量面为任意形状，能够在混响环境中分离出任意形状目标源的辐射声场。

本发明基于时域等效源法的实时声场分离方法的特点是按如下步骤进行：

步骤a、在目标源和干扰源之间布置一个测量面；在所述测量面上分布M个测点；并有M个传声器一一对应设置在各测点上进行同步测量，获得测量面上各测点处的时域声压和时域质点振速；

步骤b、位于目标源和测量面之间，以及位于干扰源和测量面之间各布置一个虚源面，两个虚源面到测量面的距离相等，在两个虚源面上各分布N个等效源，N不大于M；各等效源的分布与各测点的分布共形；所述等效源为标准点源、面源或体源；测量面上的声场为分布在测量面两侧的N个等效源所产生声场的叠加，测量面上第m个测点处的时域声压p(r_m，t)和时域质点振速v(_rm，t)由式(1)和式(2)所表征：

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

符号“*”表示两个函数之间的卷积；t表示采样时间；

n表示N个等效源中任意一个等效源；m表示M个测点中任意一个测点；

I表示干扰源侧声场，O表示目标源侧声场；

p表示声压，v表示质点振速；r_m表示第m个测点的位置；

r_In表示第n个干扰源侧等效源位置，r_On表示第n个目标源侧等效源位置；

q_In(t)为第n个干扰源侧等效源的时域源强；

q_On(t)为第n个目标源侧等效源的时域源强；

g_p(r_m，r_In，t)是第m个测点与第n个干扰源侧等效源之间的时域声压格林函数；

g_p(r_m，r_On，t)是第m个测点与第n个目标源侧等效源之间的时域声压格林函数；

g_v(r_m，r_In，t)是第m个测点与第n个干扰源侧等效源之间的时域质点振速格林函数；

g_v(r_m，r_On，t)是第m个测点与第n个目标源侧等效源之间的时域质点振速格林函数；

基于时域声压格林函数g_p(r_m，r_n，t)和时域质点振速格林函数g_v(r_m，r_n，t)分别由式(3)和式(4)所示表征，根据按式(3)和式(4)的形式分别获得格林函数g_p(r_m，r_In，t)、g_p(r_m，r_On，t)、g_v(r_m，r_In，t)和g_v(r_m，r_on，t)；

在式(3)和式(4)中：

c是声速；R_mn是第m个测点和第n个等效源之间的距离；x表示测量面的单位法向向量；

ρ是介质密度；δ(t-R_mn/c)表示狄拉克函数；H(t-R_mn/c)表示赫维赛德函数；

通过对时间进行卷积计算，将式(1)变换为式(5)，将式(2)变换为式(6)：

式(5)和式(6)中：

R_Imn为第m个测点和第n个干扰源侧等效源之间的距离；

R_Omn为第m个测点和第n个目标源侧等效源之间的距离；

τ_Imn和τ_Omn均为延迟时刻，并有：τ_Imn＝t-R_Imn/c，τ_Omn＝t-R_Omn/c；

为第n个干扰源侧等效源在τ_Imn时刻的对时间t的时域源强积分；

为第n个干扰源侧等效源在τ_Imn时刻的时域源强；

为第n个目标源侧等效源在τ_Omn时刻的对时间t的时域源强积分；

为第n个目标源侧等效源在τ_Omn时刻的时域源强；

步骤c、采用线性插值函数对各等效源的时域源强积分进行插值处理，第m个测点处在t^j时刻的声压p(r_m，t^j)和质点振速v(r_m，t^j)分别表示为式(7)和式(8)：

式(7)和式(8)中：

τ表示延迟时间；

j表示采样时间离散点数，J为总采样步阶数，j＝1，2，...，J；

k表示延迟时间离散点数，K为总延迟步阶数，k＝1，2，...，K，且：J＝K；

t^j为第j个离散采样步阶，t^j＝t¹+(j-1)Δt，Δt为采样时间间隔，t¹为初始采样步阶；

τ^k为第k个离散延迟步阶，τ^k＝τ¹+(k-1)Δτ，τ¹＝t¹-R_min/c；

Δτ为延迟时间间隔，Δτ＝Δt，R_min为所有等效源和测点之间的最小距离；

为第n个干扰源侧等效源发出声波的时刻，且所述声波在t^j时刻被第m个测点所接收；

为第n个目标源侧等效源发出声波的时刻，且所述声波在t^j时刻被第m个测点所接收；

是第n个干扰源侧等效源在τ^k时刻的对时间t的时域源强积分；

是第n个目标源侧等效源在τ^k时刻的对时间t的时域源强积分；

是以

为参量的线性插值函数；

是以

为参量的线性插值函数；

基于以延迟时间τ为参量的线性插值函数

由式(9)所表征：

根据式(9)所示线性插值函数

的表达式分别获得线性插值函数

和

则所有M个测点在t^j时刻的声压P^j和质点振速V^j分别由式(10)和式(11)所表征：

式(10)和式(11)中：

P^j是由m＝1，2，...，M的p(r_m，t^j)组成的列向量；V^j是由m＝1，2，...，M的v(r_m，t^j)组成的列向量；

是由n＝1，2，...，N的

组成的列向量；

是由n＝1，2，...，N的

组成的列向量；

和

均为M×N的传递矩阵，各传递矩阵的元素分别为：

对于指定的m和n，R_Imn和R_Omn在任意时刻均为常数，则各线性插值函数

和

均具有式(16)所示的线性插值函数

的相同形式：

因此，各传递矩阵

和

均具有由式(17)所表达的传递矩阵G^jk的相同形式：

G^jk＝G^(j-1)(k-1)＝G^(j-2)(k-2)＝…＝G^(j-k+1)(1) (17)

为使G^(j-k+1)(1)中元素不为0，要求满足式(18)：

式(18)中：

τ²为第二个延迟步阶；

为第n个等效源发出声波的时刻，且所述声波在t^j-k+1时刻被第m个测点所接收；

由式(18)推导得式(19)：

式(19)中：R_max是所有等效源与测点之间的最大距离；

令

则将式(10)变换为式(20)、将式(11)变换为式(21)：

式(20)和式(21)中：

步骤d、从初始延迟步阶开始依次获得各延迟步阶测量面两侧各等效源的时域源强积分：

首先，联立式(20)和式(21)，获得式(22)：

式(22)中：

和

分别表示j＝1、k＝1时的传递矩阵

和

和

分别是由n＝1，2，...，N的

和

所组成的列向量；

是第n个干扰源侧等效源在τ^j时刻的对时间t的时域源强积分；

是第n个目标源侧等效源在τ^j时刻的对时间t的时域源强积分；

τ^j为第j个离散延迟步阶；

由式(22)依次获得各延迟步阶的测量面两侧各等效源的时域源强积分，如式(23)所示：

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

由于求等效源时域源强积分的过程属于逆问题，故采用加权最小二乘法稳定求解过程，加权最小二乘解的形式如式(24)：

式(24)中，W为权重对角矩阵；

步骤e、根据获得的测量面两侧的各时刻等效源的时域源强积分以及等效源与任意测点之间的传递关系，计算获得目标源声场，实现实时声场分离如式(25)和式(26)所示：

式(25)和式(26)中：

是由m＝1，2，...，M的p_O(r_m，t^j)组成的列向量；

是由m＝1，2，...，M的v_O(r_m，t^j)组成的列向量；

p_O(r_m，t^j)为目标源侧等效源在t^j时刻，在第m个测点处产生的声压；

v_O(r_m，t^j)为目标源侧等效源在t^j时刻，在第m个测点处产生的质点振速。

本发明基于时域等效源法的实时声场分离方法的特点也在于：所述各测点上的时域声压信号和时域质点振速信号是采用声压和质点振速传感器阵列在测量面上一次测量获得。

本发明基于时域等效源法的实时声场分离方法的特点也在于：所述干扰源为噪声源、反射源或散射源。

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

1、本发明方法中的测量面可以为任意形状面，克服了传统方法只适用于平面和球面等规则形状测量面的缺陷，实现瞬态声场实时分离。

2.本发明方法中测量声压直接用于声场分离，无需进行额外的平面波或球面波展开，因此测量数据处理过程更加简单。

附图说明

图1为本发明中采用平面测量面时目标源、干扰源、测点和等效源的空间分布；

图2为本发明中采用半圆柱测量面时目标源、干扰源、测点和等效源的空间分布；

图3为本发明中采用半球面测量面时目标源、干扰源、测点和等效源的空间分布；

图4为本发明中采用平面测量面时测点A₁处的时域声压信号；

图5为本发明中采用半圆柱面测量面时测点A₁处的时域声压信号；

图6为本发明中采用半球面测量面时测点A₁处的时域声压信号；

图7为本发明中采用平面测量面时测点A₂处的时域声压信号；

图8为本发明中采用半圆柱面测量面时测点A₂处的时域声压信号；

图9为本发明中采用半球面测量面时测点A₂处的时域声压信号；

图10为本发明中在t＝68Δt时刻的声场的空间分布；

图11a、图11b和图11c为本发明中衡量声场分离准确性的量化指标T₁和T₂；

具体实施方式

本实施例中基于时域等效源法的实时声场分离方法是：在目标源和干扰源之间布置一个测量面，目标源是指目标声源，干扰源是指干扰声源，同步采集测量面上各测点处的时域声压和时域质点振速；在测量面两侧布置若干等效源，用于分别模拟目标源声场和干扰源声场；采用线性插值函数对各等效源的时域源强积分进行插值处理，建立任意时刻各测点的声压和质点振速与各等效源源强之间的传递关系；从初始时刻开始依次获得对应时刻测量面两侧各等效源的时域源强积分；根据获得的测量面两侧的各时刻等效源的时域源强积分以及等效源与任意测点之间的传递关系，计算获得目标源声场和干扰源声场，实现实时声场分离；测量面为任意形状，能够在混响环境中分离出任意形状目标源的辐射声场。

具体实施中，基于时域等效源法的实时声场分离方法按如下步骤进行：

步骤a、在目标源和干扰源之间布置一个测量面；在测量面上分布M个测点；并有M个传声器一一对应设置在各测点上进行同步测量，获得测量面上各测点处的时域声压和时域质点振速。

步骤b、位于目标源和测量面之间，以及位于干扰源和测量面之间各布置一个虚源面，两个虚源面到测量面的距离相等，在两个虚源面上各分布N个等效源，N不大于M；各等效源的分布与各测点的分布共形；等效源为标准点源、面源或体源；测量面上的声场为分布在测量面两侧的N个等效源所产生声场的叠加，即：

p(r_m，t)＝p_I(r_m，t)+p_O(r_m，t)

v(r_m，t)＝v_I(r_m，t)+v_O(r_m，t)

其中：

p_I(r_m，t)为干扰源侧等效源在第m个测点处的时域声压；

p_O(r_m，t)为目标源侧等效源在第m个测点处的时域声压；

v_I(r_m，t)为干扰源侧等效源在第m个测点处的时域质点振速；

v_O(r_m，t)为目标源侧等效源在第m个测点处的时域质点振速；

因此有，测量面上第m个测点处的时域声压p(r_m，t)和时域质点振速v(r_m，t)由式(1)和式(2)所表征：

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

符号“*”表示两个函数之间的卷积；t表示采样时间；

n表示N个等效源中任意一个等效源；m表示M个测点中任意一个测点；

I表示干扰源侧声场，O表示目标源侧声场；

p表示声压，v表示质点振速；r_m表示第m个测点的位置；

r_In表示第n个干扰源侧等效源位置，r_On表示第n个目标源侧等效源位置；

q_In(t)为第n个干扰源侧等效源的时域源强；

q_On(t)为第n个目标源侧等效源的时域源强；

g_p(r_m，r_In，t)是第m个测点与第n个干扰源侧等效源之间的时域声压格林函数；

g_p(r_m，r_On，t)是第m个测点与第n个目标源侧等效源之间的时域声压格林函数；

g_v(r_m，r_In，t)是第m个测点与第n个干扰源侧等效源之间的时域质点振速格林函数；

g_v(r_m，r_On，t)是第m个测点与第n个目标源侧等效源之间的时域质点振速格林函数；

基于时域声压格林函数g_p(r_m，r_n，t)和时域质点振速格林函数g_v(r_m，r_n，t)分别由式(3)和式(4)所示表征，根据按式(3)和式(4)的形式分别获得格林函数g_p(r_m，r_In，t)、g_p(r_m，r_On，t)、g_v(r_m，r_In，t)和g_v(r_m，r_On，t)；

在式(3)和式(4)中：

c是声速；R_mn是第m个测点和第n个等效源之间的距离；x表示测量面的单位法向向量；

ρ是介质密度；δ(t-R_mn/c)表示狄拉克函数；H(t-R_mn/c)表示赫维赛德函数；

通过对时间进行卷积计算，将式(1)变换为式(5)，将式(2)变换为式(6)：

式(5)和式(6)中：

R_Imn为第m个测点和第n个干扰源侧等效源之间的距离；

R_Omn为第m个测点和第n个目标源侧等效源之间的距离；

τ_Imn和τ_Omn均为延迟时刻，并有：τ_Imn＝t-R_Imn/c，τ_Omn＝t-R_Omn/c；

为第n个干扰源侧等效源在τ_Imn时刻的对时间t的时域源强积分；

为第n个干扰源侧等效源在τ_Imn时刻的时域源强；

为第n个目标源侧等效源在τ_Omn时刻的对时间t的时域源强积分；

为第n个目标源侧等效源在τ_Omn时刻的时域源强。

步骤c、采用线性插值函数对各等效源的时域源强积分进行插值处理，第m个测点处在t^j时刻的声压p(r_m，t^j)和质点振速v(r_m，t^j)分别表示为式(7)和式(8)：

式(7)和式(8)中：

τ表示延迟时间；

j表示采样时间离散点数，J为总采样步阶数，j＝1，2，...，J；

k表示延迟时间离散点数，K为总延迟步阶数，k＝1，2，...，K，且：J＝K；

t^j为第j个离散采样步阶，t^j＝t¹+(j-1)Δt，Δt为采样时间间隔，t¹为初始采样步阶；

τ^k为第k个离散延迟步阶，τ^k＝τ¹+(k-1)Δτ，τ¹＝t¹-R_min/c；

Δτ为延迟时间间隔，Δτ＝Δt，R_min为所有等效源和测点之间的最小距离；

为第n个干扰源侧等效源发出声波的时刻，且声波在t^j时刻被第m个测点所接收；

为第n个目标源侧等效源发出声波的时刻，且声波在t^j时刻被第m个测点所接收；

是第n个干扰源侧等效源在τ^k时刻的对时间t的时域源强积分；

是第n个目标源侧等效源在τ^k时刻的对时间t的时域源强积分；

是以

为参量的线性插值函数；

是以

为参量的线性插值函数；

基于以延迟时间τ为参量的线性插值函数

由式(9)所表征：

根据式(9)所示线性插值函数

的表达式分别获得线性插值函数

和

则所有M个测点在t^j时刻的声压P^j和质点振速V^j分别由式(10)和式(11)所表征：

式(10)和式(11)中：

P^j是由m＝1，2，...，M的p(r_m，t^j)组成的列向量；V^j是由m＝1，2，...，M的v(r_m，t^j)组成的列向量；

是由n＝1，2，...，N的

组成的列向量；

是由n＝1，2，...，N的

组成的列向量；

和

均为M×N的传递矩阵，各传递矩阵的元素分别为：

对于指定的m和n，R_Imn和R_Omn在任意时刻均为常数，则各线性插值函数

和

均具有式(16)所示的线性插值函数

的相同形式：

因此，各传递矩阵

和

均具有由式(17)所表达的传递矩阵G^jk的相同形式：

G^jk＝G^(j-1)(k-1)＝G^(j-2)(k-2)＝…＝G^(j-k+1)(1) (17)

为使G^(j-k+1)(1)中元素不为0，要求满足式(18)：

式(18)中：

τ²为第二个延迟步阶；

为第n个等效源发出声波的时刻，且声波在t^j-k+1时刻被第m个测点所接收；

由式(18)推导得式(19)：

式(19)中：Rmax是所有等效源与测点之间的最大距离；

令

则将式(10)变换为式(20)、将式(11)变换为式(21)：

式(20)和式(21)中：

步骤d、从初始延迟步阶开始依次获得各延迟步阶测量面两侧各等效源的时域源强积分：

首先，联立式(20)和式(21)，获得式(22)：

式(22)中：

和

分别表示j＝1、k＝1时的传递矩阵

和

和

分别是由n＝1，2，...，N的

和

所组成的列向量；

是第n个干扰源侧等效源在τ^j时刻的对时间t的时域源强积分；

是第n个目标源侧等效源在τ^j时刻的对时间t的时域源强积分；

τ^j为第j个离散延迟步阶；

由式(22)依次获得各延迟步阶的测量面两侧各等效源的时域源强积分，如式(23)所示：

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

由于求等效源时域源强积分的过程属于逆问题，故采用加权最小二乘法稳定求解过程，加权最小二乘解的形式如式(24)：

式(24)中，W是维度为2M的权重对角矩阵，并有W为：

W＝diag(1/||P_rms||，1/||V_rms||)

其中：

P_rms和V_rms分别表示对各测点的时域声压和时域质点振速求均方根值；||·||表示求向量2-范数；

diag(1/||P_rms||，1/||V_rms||)表示以1/||P_rms||，1/||V_rms||中各元素组成对角矩阵。

步骤e、根据获得的测量面两侧的各时刻等效源的时域源强积分以及等效源与任意测点之间的传递关系，计算获得目标源声场，从而实现实时声场分离，如式(25)和式(26)所示：

式(25)和式(26)中：

是由m＝1，2，...，M的p_O(r_m，t^j)组成的列向量；

是由m＝1，2，...，M的v_O(r_m，t^j)组成的列向量；

p_O(r_m，t^j)为目标源侧等效源在t^j时刻，在第m个测点处产生的声压；

v_O(r_m，t^j)为目标源侧等效源在t^j时刻，在第m个测点处产生的质点振速。

本实施例中各测点上的时域声压信号和时域质点振速信号是采用声压和质点振速传感器阵列在测量面上一次测量获得；干扰源为噪声源、反射源或散射源。

针对本发明方法按如下方式进行检验：

在空间放置两个单极子声源分别作为目标源和干扰源，以目标源为原点建立笛卡尔坐标系o(x，y，z)，干扰源的坐标为(-0.1m，-0.1m，0.4m)。给定声源的时域源强积分s(t)为：

其中：f₀为中心频率，目标源中心频率f₀＝600Hz，干扰源中心频率f₀＝1200Hz。

测量面形状分别按平面、半圆柱面和半球面进行检验。

图1所示平面测量阵列，孔径为0.4m×0.4m，相邻传声器在x和y方向间距均为0.05m；

图2所示半圆柱面测量阵列，圆柱半径为0.2m，由9个半圆环、每个半圆环中分布有13个传声器所组成，每个半圆环之间轴向间距为0.05cm，同一个半圆环中相邻传声器的周向角为15°；

图3所示半球面测量阵列，球面半径为0.2m，由6圈、每圈24个传声器以及顶部的1个传声器组成一个半球面，同一圈阵列中相邻传声器的极角相差15°，不同圈阵列相应位置传声器的方位角相差15°；

在图1、图2和图3中，b₁表示干扰源，b₂表示干扰源侧等效源，b₃表示测点，b₄表示目标源侧等效源，b₅表示目标源；

三种测量阵列各测点上单极子声源辐射的理论声压p(r_m，t)和法向质点振速v(r_m，t)由式(27)和式(28)给出：

式(27)和式(28)中：

R_ms为第m个测量点到单极子声源之间的距离；

s′(t-R_ms/c)表示声源源强在t-R_ms/c时刻的值；

s(t-R_ms/c)表示声源源强对时间t的积分在t-R_ms/c时刻的值；

声压和质点振速的采样频率为f_e＝12.8kHz，共采样128个采样点。为了模拟实际测量条件，在信号中加入信噪比为30dB的高斯白噪声。

对于平面测量阵列，虚源面与测量面的距离为0.08m；对于半圆柱面测量阵列，虚源面与测量面的距离为0.15m；对于半球面测量阵列，虚源面与测量面的距离为0.08m。

为了在时域内检验本发明方法分离声场的有效性，选取两个测点，即测点A₁和测点A₂，观察测点A₁和测点A₂位置处的时域信号；对于平面测量阵列，测点A₁的坐标为A₁(0，0，0.2)，测点A₂的坐标为A₂(-0.1，-0.1，0.2)；对于半圆柱面测量阵列，测点A₁的坐标为A₁(0，0，0.2)，测点A₂的坐标为A₂(-0.1，-0.1，0.173)；对于半球面测量阵列，测点A₁的坐标为A₁(0，0，0.2)，测点A₂坐标值为A₂(-0.1，-0.1，0.141)，各坐标值是以m为单位。

图4、图5和图6分别为本发明中采用平面阵列、半圆柱面阵列和半球面阵列测量时，测点A₁处的时域声压信号；其中虚线即曲线a₁₁为未分离的声压信号，空心点线即曲线a₁₂为应用本发明方法分离出的目标源辐射声压信号，实心点线即曲线a₁₃为理论目标源辐射声压信号；图7、图8和图9分别为本发明中采用平面阵列、半圆柱面阵列和半球面阵列测量时，测点A₂处的时域声压信号；其中虚线即曲线a₂₁为未分离的声压信号，空心点线即曲线a₂₂为应用本发明方法分离出的目标源辐射声压信号，实心点线即曲线a₂₃为理论目标源辐射声压信号；由图4至图9可见，本发明方法能够有效消除干扰源的影响，分离出的目标源辐射声场与目标源辐射声场理论值吻合良好。

图10为本发明方法中在t＝68Δt时刻的声场的空间分布，图10中第一行对应平面测量阵列，第二行对应半圆柱面测量阵列，第三行对应半球面测量阵列；图10中第一列代表目标源理论辐射声场，第二列代表未分离的混合声场，第三列代表应用本发明方法分离出的目标源辐射声场。图10可见干扰源对目标源声场的影响非常严重，应用本发明方法能够很大程度上消除干扰。

为了量化本发明方法分离声场的有效性，引入两个指标T₁和T₂，其中T₁表示各测点处分离出的目标源辐射声场与目标源理论辐射声场之间相位的吻合度，T₁的值越靠近1，表明相位吻合度越高；T₂表示各测点处分离出的目标源辐射声场与目标源理论辐射声场之间幅值的差异度，T₂的值越靠近0，表明幅值差异度越小。T₁和T₂分别由式(29)和式(30)表示：

式(29)和式(30)中：

上标″T″代表向量的转置，p_t代表目标源理论辐射声场，p_s代表分离出的目标源辐射声场。

图11a、图11b和图11c为本发明中衡量声场分离准确性的量化指标T₁和T₂的计算值，其中图11a、图11b和图11c一一对应为平面测量阵列、半圆柱面测量阵列和半球面测量阵列；由图11a、图11b和图11c可见，实心点线表示的T₁的值均大于0.8，空心点线表示T₂的值绝大部分均小于0.2，本发明方法分离声场结果有效。

Claims (4)

1.一种基于时域等效源法的实时声场分离方法，其特征是：在目标源和干扰源之间布置一个测量面，所述目标源是指目标声源，所述干扰源是指干扰声源，同步采集测量面上各测点处的时域声压和时域质点振速；在测量面两侧布置若干等效源，用于分别模拟目标源声场和干扰源声场；采用线性插值函数对各等效源的时域源强积分进行插值处理，建立任意时刻各测点的声压和质点振速与各等效源源强之间的传递关系；从初始时刻开始依次获得对应时刻测量面两侧各等效源的时域源强积分；根据获得的测量面两侧的各时刻等效源的时域源强积分以及等效源与任意测点之间的传递关系，计算获得目标源声场和干扰源声场，实现实时声场分离；所述测量面为任意形状，能够在混响环境中分离出任意形状目标源的辐射声场。

2.根据权利要求1所述的基于时域等效源法的实时声场分离方法，其特征是按如下步骤进行：

步骤a、在目标源和干扰源之间布置一个测量面；在所述测量面上分布M个测点；并有M个传声器一一对应设置在各测点上进行同步测量，获得测量面上各测点处的时域声压和时域质点振速；

步骤b、位于目标源和测量面之间，以及位于干扰源和测量面之间各布置一个虚源面，两个虚源面到测量面的距离相等，在两个虚源面上各分布N个等效源，N不大于M；各等效源的分布与各测点的分布共形；所述等效源为标准点源、面源或体源；测量面上的声场为分布在测量面两侧的N个等效源所产生声场的叠加，测量面上第m个测点处的时域声压p(r_m，t)和时域质点振速v(r_m，t)由式(1)和式(2)所表征：

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

符号“*”表示两个函数之间的卷积；t表示采样时间；

n表示N个等效源中任意一个等效源；m表示M个测点中任意一个测点；

I表示干扰源侧声场，O表示目标源侧声场；

p表示声压，v表示质点振速；r_m表示第m个测点的位置；

r_In表示第n个干扰源侧等效源位置，r_On表示第n个目标源侧等效源位置；

q_In(t)为第n个干扰源侧等效源的时域源强；

q_On(t)为第n个目标源侧等效源的时域源强；

g_p(r_m，r_In，t)是第m个测点与第n个干扰源侧等效源之间的时域声压格林函数；

g_p(r_m，r_On，t)是第m个测点与第n个目标源侧等效源之间的时域声压格林函数；

g_v(r_m，r_In，t)是第m个测点与第n个干扰源侧等效源之间的时域质点振速格林函数；

g_v(r_m，r_On，t)是第m个测点与第n个目标源侧等效源之间的时域质点振速格林函数；

基于时域声压格林函数g_p(r_m，r_n，t)和时域质点振速格林函数g_v(r_m，r_n，t)分别由式(3)和式(4)所示表征，根据按式(3)和式(4)的形式分别获得格林函数g_p(r_m，r_In，t)、g_p(r_m，r_On，t)、g_v(r_m，r_In，t)和g_v(r_m，r_On，t)；

在式(3)和式(4)中：

c是声速；R_mn是第m个测点和第n个等效源之间的距离；x表示测量面的单位法向向量；

ρ是介质密度；δ(t-R_mn/c)表示狄拉克函数；H(t-R_mn/c)表示赫维赛德函数；

通过对时间进行卷积计算，将式(1)变换为式(5)，将式(2)变换为式(6)：

式(5)和式(6)中：

R_Imn为第m个测点和第n个干扰源侧等效源之间的距离；

R_omn为第m个测点和第n个目标源侧等效源之间的距离；

τ_Imn和τ_Omn均为延迟时刻，并有：τ_Imn＝t-R_Imn/c，τ_Omn＝t-R_omn/c；

为第n个干扰源侧等效源在τ_Imn时刻的对时间t的时域源强积分；

为第n个干扰源侧等效源在τ_Imn时刻的时域源强；

为第n个目标源侧等效源在τ_Omn时刻的对时间t的时域源强积分；

为第n个目标源侧等效源在τ_Omn时刻的时域源强；

步骤c、采用线性插值函数对各等效源的时域源强积分进行插值处理，第m个测点处在t^j时刻的声压p(r_m，t^j)和质点振速v(r_m，t^j)分别表示为式(7)和式(8)：

式(7)和式(8)中：

τ表示延迟时间；

j表示采样时间离散点数，J为总采样步阶数，j＝1，2，...，J；

k表示延迟时间离散点数，K为总延迟步阶数，k＝1，2，...，K，且：J＝K；

t^j为第j个离散采样步阶，t^j＝t¹+(j-1)Δt，Δt为采样时间间隔，t¹为初始采样步阶；

τ^k为第k个离散延迟步阶，τ^k＝τ¹+(k-1)Δτ，τ¹＝t¹-R_min/c；

Δτ为延迟时间间隔，Δτ＝Δt，R_min为所有等效源和测点之间的最小距离；

为第n个干扰源侧等效源发出声波的时刻，且所述声波在t^j时刻被第m个测点所接收；

为第n个目标源侧等效源发出声波的时刻，且所述声波在t^j时刻被第m个测点所接收；

是第n个干扰源侧等效源在τ^k时刻的对时间t的时域源强积分；

是第n个目标源侧等效源在τ^k时刻的对时间t的时域源强积分；

是以

为参量的线性插值函数；

是以

为参量的线性插值函数；

基于以延迟时间τ为参量的线性插值函数

由式(9)所表征：

根据式(9)所示线性插值函数

的表达式分别获得线性插值函数

和

则所有M个测点在t^j时刻的声压P^j和质点振速V^j分别由式(10)和式(11)所表征：

式(10)和式(11)中：

P^j是由m＝1，2，...，M的p(r_m，t^j)组成的列向量；V^j是由m＝1，2，...，M的v(r_m，t^j)组成的列向量；

是由n＝1，2，...，N的

组成的列向量；

是由n＝1，2，...，N的

组成的列向量；

和

均为M×N的传递矩阵，各传递矩阵的元素分别为：

对于指定的m和n，R_Imn和Ro_mn在任意时刻均为常数，则各线性插值函数

和

均具有式(16)所示的线性插值函数

的相同形式：

因此，各传递矩阵

和

均具有由式(17)所表达的传递矩阵G^jk的相同形式：

G^jk＝G^(j-1)(k-1)＝G^(j-2)(k-2)＝…＝G^(j-k+1)(1) (17)

为使G^(j-k+1)(1)中元素不为0，要求满足式(18)：

式(18)中：

τ²为第二个延迟步阶；

为第n个等效源发出声波的时刻，且所述声波在t^j-k+1时刻被第m个测点所接收；

由式(18)推导得式(19)：

式(19)中：R_max是所有等效源与测点之间的最大距离；

令

则将式(10)变换为式(20)、将式(11)变换为式(21)：

式(20)和式(21)中：

步骤d、从初始延迟步阶开始依次获得各延迟步阶测量面两侧各等效源的时域源强积分：

首先，联立式(20)和式(21)，获得式(22)：

式(22)中：

和

分别表示j＝1、k＝1时的传递矩阵

和

和

分别是由n＝1，2，...，N的

和

所组成的列向量；

是第n个干扰源侧等效源在τ^j时刻的对时间t的时域源强积分；

是第n个目标源侧等效源在τ^j时刻的对时间t的时域源强积分；

τ^j为第j个离散延迟步阶；

由式(22)依次获得各延迟步阶的测量面两侧各等效源的时域源强积分，如式(23)所示：

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

由于求等效源时域源强积分的过程属于逆问题，故采用加权最小二乘法稳定求解过程，加权最小二乘解的形式如式(24)：

式(24)中，W为权重对角矩阵；

步骤e、根据获得的测量面两侧的各时刻等效源的时域源强积分以及等效源与任意测点之间的传递关系，计算获得目标源声场，实现实时声场分离如式(25)和式(26)所示：

式(25)和式(26)中：

是由m＝1，2，...，M的p_o(r_m，t^j)组成的列向量；

是由m＝1，2，...，M的v_O(r_m，t^j)组成的列向量；

p_O(r_m，t^j)为目标源侧等效源在t^j时刻，在第m个测点处产生的声压；

v_O(r_m，t^j)为目标源侧等效源在t^j时刻，在第m个测点处产生的质点振速。

3.根据权利要求1所述的基于时域等效源法的实时声场分离方法，其特征是：所述各测点上的时域声压信号和时域质点振速信号是采用声压和质点振速传感器阵列在测量面上一次测量获得。

4.根据权利要求1所述的基于时域等效源法的实时声场分离方法，其特征是：所述干扰源为噪声源、反射源或散射源。

Images