CN110554358A - 一种基于虚拟球阵列扩展技术的噪声源定位识别方法 - Google Patents

Abstract

本发明提供一种基于虚拟球阵列扩展技术的噪声源定位识别方法，具体地说是将基于波束形成正则化矩阵的虚拟阵列扩展技术与空心球阵相结合，得到虚拟双半径空心球阵列和通过‘声压‑振速’联合处理获得的具有心形指向性的虚拟阵元组成的空心球阵列。该方法可有效地克服了单层空心球阵在其径向函数零点(球贝塞尔函数极小值)对应的频率定位不稳定问题，从而提高了噪声源定位识别精度。

Description

一种基于虚拟球阵列扩展技术的噪声源定位识别方法

技术领域

本发明涉及一种基于虚拟球阵列扩展技术的噪声源定位识别方法。

背景技术

准确获取目标的主要噪声源空间分布特性，是有针对性开展减振降噪工作的重要环节。相比其它阵列形式，球形传声器阵列具有全空间指向性和三维对称性，可同时获得到三维空间的声场分布信息。另外，相较于刚性球阵，空心球阵对声场影响较小。因此，开展基于空心球阵的噪声源定位识别方法研究具有重大的工程应用价值。

基于平面波假设的球谐函数波束形成算法最早由Jens Meyer与Gary Elko于2002年提出(A highly scalable spherical microphone array based on an orthonormaldecomposition of the sound field[C]//IEEE International Conference onAcoustics,Speech and Signal Processing.Orlando,FL,United states,May 2002,2:1781-1784.)。在此基础上，T.D.Abhayapala和D.B.Ward给出了空心球阵的设计方法(Theory and design of high order sound field microphones using sphericalmicrophone array[C].IEEE International Conference on Acoustics,Speech,andSignal Processing,2002:II-1949-II-1952.)。单层空心球阵的径向函数在某些频率处存在零点，导致定位方法在这些频点处的稳健性差，此时较小的阵列位置误差、噪声便导致定位结果变差甚至失效。针对以上问题，2007年，I.Balmages和B.Rafaely等设计了双半径空心球阵和虚拟空心球阵列(阵元具有心形指向性)。通过双层阵联合处理方法克服了单层空心球阵的径向函数零点问题。(Open-sphere designs for spherical microphone arrays[J].IEEE Trans.Audio,Speech&Language Process.,2007,15(2):727-732.)，(High-resolution plane-wave decomposition in an auditorium using a dual-radiusscanning spherical microphone array[J].J.Acoust.Soc.Am.,2007,122(5):2661-2668.)。但是在实际工程应用中采用双层阵列形式，需要大量传感器，测量成本高，测量系统庞大复杂。

针对上述问题，本发明将基于波束形成正则化矩阵的虚拟阵列扩展技术与空心球阵相结合，构建虚拟双半径空心球阵列和虚拟空心球阵列(矢量传声器)列，采用球谐函数展开波束形成算法，提出了适用于球形阵列的三维空间噪声源定位识别方法。该方法显著提高了噪声源定位识别精度，同时降低测量系统复杂度，节约成本。

发明内容

本发明的目的是为了提供一种基于虚拟球阵列扩展技术的噪声源定位识别方法，将基于波束形成正则化矩阵的虚拟阵列扩展技术与空心球阵相结合，得到虚拟双半径空心球阵列和通过‘声压-振速’联合处理获得的具有心形指向性的虚拟阵元组成的空心球阵列。该方法可有效地克服了单层空心球阵在其径向函数零点(球贝塞尔函数极小值)对应的频率定位不稳定问题，从而提高了噪声源定位识别精度。

本发明的目的是这样实现的：步骤如下：

步骤一：获取半径为a阵元数为L的球阵列接收声压向量p，其中第l个阵元为：p(ka,θ_l,φ_l))；

在球阵外部和声源之间选取半径为r'的球面作为虚拟源面并在该面上设置K个虚拟源及位置，计算虚拟源和球阵列的传递矩阵G(a_l,r'_k)，利用p和G(a_l,r'_k)计算得到波束形成正则化矩阵L，选取正则化参数λ，得到虚拟源的源强q_BF；

在虚拟源面和球阵之间选取一个球面并设置M个虚拟阵元及位置，得到虚拟源和虚拟球阵间的声场传递矩阵G(r'_k,r_m)，最后通过q_BF和G(r'_k,r_m)计算得到虚拟球阵接收声压向量p'，其中第m个虚拟阵元为：p'(kr,θ_m,φ_m))；

步骤二：基于虚拟球阵列扩展技术的球谐波展开噪声源定位方法。

本发明还包括这样一些结构特征：

1.步骤二具体包括：

(1)将虚拟阵列和真实阵列接收数据进行球函数分解，按照‘分阶-选择’联合处理方法得到虚拟双半径空心球阵列的组合声压，采用球函数展开波束形成算法进行噪声源定位；

选取半径为r的球面作为虚拟球阵位置且虚拟阵元方位与实际阵列阵元方向相同，将其与实际的球阵视为整体构成一个虚拟的双半径空心球阵，利用选择参数β_n将p'(kr,θ_l,φ_l)以及p(ka,θ_l,φ_l)按照球谐函数分阶选择处理获得虚拟双半径空心球阵组合声压p_Z(ka,kr,θ_l,φ_l)；将p_Z(ka,kr,θ_l,φ_l)代入到球谐波展开波束形成算法中实现噪声源定位；

(2)将虚拟阵列和真实阵列相同方位阵元作有限差分处理求得两个阵元中间位置r_c处振速v(kr_c,θ_l,φ_l)，并采用‘p+v’联合处理方式得到虚拟空心球阵列的组合声压进而采用球函数展开波束形成算法进行噪声源定位；

在一个半径r大于a的虚拟球面上设置虚拟阵元，使虚拟空心球阵的阵元数目及方位均与真实阵列相同，通过两个球阵相同方位的阵元作有限差分求得两个阵元中间位置r_c处的振速v(kr_c,θ_l,φ_l)，同时通过步骤一得到振速位置处的声压p'(kr_c,θ_l,φ_l)，选取一阶参数γ，将v(kr_c,θ_l,φ_l)和p'(kr_c,θ_l,φ_l)通过‘p+v’的联合处理获得其声压的组合形式p_E(kr_c,θ_l,φ_l)，最后将p_E(kr_c,θ_l,φ_l)代入到球函数展开波束形成算法中实现噪声源定位。

2.步骤二中的(1)具体包括：

在真实空心球阵外部选择一个半径为r_z的虚拟球面，并在球面上设置虚拟阵元数目及位置，构成一个虚拟的大孔径空心球阵列，通过选择参数β_n在j_n(ka)和j_n(kr)之间选取最佳的径向函数来克服单个空心球阵径向函数的零点问题，选择参数表达式为：

声压的组合形式可以表示为：

p_Z(ka,kr,θ_l,φ_l)＝(1-β_n)p_l(ka,θ_l,φ_l)+β_np'_l(kr,θ_l,φ_l)

将上述组合声压代入到球谐函数展开波束形成算法中得到：

3.步骤二中的(2)具体包括：

在一个半径r_z大于a的虚拟球面上设置虚拟阵元，保证虚拟空心球阵的阵元数目和方位与真实阵列的阵元相同，通过两个球阵相同方位的阵元作有限差分求得两个阵元中间位置处的振速：

同时采用第二步过程计算振速点位置的声压p(kr_c,θ,φ)，最终得到新的虚拟球阵阵元声压的组合形式为：

p_E(r_c,θ,φ)＝γp'(r_c,θ,φ)-(1-γ)ρcv_r(r_c,θ,φ)

式中γ是一阶参数，一般取γ＝0.5；

将上述的组合声压代入到球谐函数展开波束形成算法中如下：

与现有技术相比，本发明的有益效果是：本发明以小孔径空心球形阵列作为测量系统的核心设备，将基于波束形成正则化矩阵的虚拟阵列扩展技术与空心球阵相结合，构建虚拟双半径空心球阵列和虚拟空心球阵列(阵元具有心形指向性)列，采用球谐函数展开波束形成算法实现噪声源定位与识别。一方面，本发明有效克服了空心球阵径向函数零点问题，提高了噪声源定位识别精度；另一方面，本发明所述方法采用基于波束形成正则化矩阵的虚拟阵列扩展技术，利用虚拟球阵列进行噪声源定位识别，相较于传统的双球层阵列噪声源定位识别方法，测量系统简单，成本低，具有良好的工程应用前景。

附图说明

图1基于虚拟球阵列扩展技术的噪声源定位识别方法总体框图；

图2阵元随机均匀分布空心球阵及球坐标系；

图3虚拟阵列扩展示意图；

图4a-c是不同球阵形式径向函数(阶数n＝0,1,2,3)对比图；

图5a-c是在ka＝2时，分别三种球阵结合球谐函数展开波束形成算法仿真结果图；

图6a-c是在ka＝6时，分别三种球阵结合球谐函数展开波束形成算法仿真结果图；

图7a-c是在空心球径向函数第一个零点(ka＝3.14)处，阵列扩展前后声源定位识别仿真结果对比图；

图8a-c是在空心球径向函数第二个零点(ka＝4.49)处，阵列扩展前后声源定位识别仿真结果对比图。

具体实施方式

下面结合附图与具体实施方式对本发明作进一步详细描述。

结合图1至图8c，本发明可有效地克服了单层空心球阵在其径向函数零点(球贝塞尔函数极小值)对应的频率定位不稳定问题，提高了噪声源定位精度，包括如下步骤：

第一步：获取半径为a阵元数为L的球阵列接收声压向量p(其中第l个阵元为p(ka,θ_l,φ_l))；在球阵外部和声源之间选取半径为r'的球面(与球阵同球心)作为虚拟源面并在该面上设置K个虚拟源及位置，计算虚拟源和球阵列的传递矩阵G(a_l,r'_k)，利用p和G(a_l,r'_k)计算得到波束形成正则化矩阵L，选取正则化参数λ，计算虚拟源的源强q_BF；在虚拟源面和球阵之间选取一个球面并设置M个虚拟阵元及位置，计算虚拟源和虚拟球阵间的声场传递矩阵G(r'_k,r_m)，最后通过q_BF和G(r'_k,r_m)计算得到虚拟球阵接收声压向量p'(其中第m个虚拟阵元写成p'(kr,θ_m,φ_m))。

如图2所示，以球阵中心为坐标原点建立坐标系，假设一点声源位于球坐标系下(r_s,θ_s,φ_s)处，该点声源在半径为a的空心球阵表面产生的声压：

式中j_n(ka)为n阶第一类球贝塞尔函数；为n阶第一类球汉克尔函数。阶次N的选取和球面阵元个数有关，令N_D为最大截断系数，通常(N_D+1)²对应的值不大于阵列传声器的数量。实际经验表明，按照以下方式确定N的取值，能够平衡阵列的分辨率与旁瓣抑制的关系：

其中[*]为取整函数。则空心球阵接收声压向量可以表示为p＝[p₁,p₂,…,p_L]。

如图3所示，在声源和球阵之间选取一个球面作为包络面，假设包络面S上存在外部实际声源辐射到该包络面的虚拟声源，则r′处的虚拟声源在S任意一点r的声场p_r′(r)满足Helmholtz方程：

那么可得包络面S上虚拟声源在S内的声场为

其中G(a,r')表示自由场格林函数。

根据以上理论分析，只要得到包络面S上的虚拟声源强度就可以获得包络面内任意场点声压值，也就是得到了S内任意位置的数据而无需增大阵列孔径，从而实现了虚拟球阵列扩展。

实际计算过程如图1所示，假设K个虚拟声源均匀分布在包络面上，将上述积分公式离散化为：

p(a_l)＝G(a_l,r'_k)q(r'_k) (5)

其中p(a_l)表示L个空心球阵阵元的接收声压向量，G(a_l,r'_k)表示虚拟声源和阵元间的传递矩阵。通过上式可以获得虚拟声源的源强为:

q＝G^-1p (6)

在实际计算过程中，由于声压传感器个数与包络面上次级声源个数不匹配，G是非方阵矩阵，直接求G^-1是病态的，得到的结果将与真实值有严重偏差，为了获得上式的稳定解，需要在求解过程中采取正则化方法。采用Tikhonov正则化方法。该方法是基于上式解的残余范数和解范数之间的联合加权达到最小的思想，即：

其中，Ω(q)＝||Lq||²称为离散平滑范数，L为正则化矩阵，约束了待定解的光滑性，当L＝I_N时该方法称为标准的Tikhonov正则化方法，λ为正则化参数。人们发现传统的正则化方法趋向于产生一个过光滑的解。对于一个实际问题，其精确解很可能具有不连续性或尖角等不光滑的性质，就会导致部分信息的丢失，从而无法获得精确的正则化解，为了克服传统正则化解过分光滑的缺点，在稳定求解的同时，可通过修改Tikhonov正则化过程中的罚项来提高正则化解的精度。当L≠I_N时，源强度的解为：

式中,

矩阵L为波束形成的归一化输出结果的倒数，称为波束形成正则化矩阵，它涵盖了次级声源的方位信息及相对源强信息，它在源强较大处施加小的惩罚项，在源强较弱处的惩罚项则较大，避免了由均一罚项导致的信息缺失。

在真实空心球阵外选取一个半径为r(r＞a)的虚拟球面，在球面上设置虚拟空心球阵的阵元方位和数目。取虚源和虚拟球阵阵元之间的声场传递函数为自由场格林函数。则虚拟球阵接收声压近似值向量为：

p_x(r_m)＝G(r_m,r'_k)q_BF(r'_k) (10)

第二步：基于虚拟球阵列扩展技术的球谐波展开噪声源定位方法。

由第一步可知，在得到包络面上虚拟源的源强度以后，就可以近似计算包络面内部任意一点声压的近似值，因此在包络面内选择不同的球面半径r即可获得不同阵孔径的虚拟空心球阵列。在本专利中采用两种不同的阵孔径的虚拟空心球阵实现两种不同的联合处理方式。

①在真实空心球阵外部选择一个半径为r_z的虚拟球面，并在球面上设置虚拟阵元数目及位置，构成一个虚拟的大孔径空心球阵列。将此虚拟空心球阵和真实的空心球阵联合视为一个组合空心球阵列，即双半径空心球阵。通过选择参数β_n在j_n(ka)和j_n(kr)之间选取最佳的径向函数来克服单个空心球阵径向函数的零点问题，选择参数表达式为：

声压的组合形式可以表示为：

p_Z(ka,kr,θ_l,φ_l)＝(1-β_n)p_l(ka,θ_l,φ_l)+β_np'_l(kr,θ_l,φ_l) (12)

将上述的组合声压代入到球谐函数展开波束形成算法中如下：

②在一个半径r_z稍大于a的虚拟球面上设置虚拟阵元，保证虚拟空心球阵的阵元数目和方位与真实阵列的阵元相同。通过两个球阵相同方位的阵元作有限差分求得两个阵元中间位置处的振速：

同时采用第二步过程计算振速点位置的声压p(kr_c,θ,φ)，最终得到新的虚拟球阵阵元声压的组合形式为：

p_E(r_c,θ,φ)＝γp'(r_c,θ,φ)-(1-γ)ρcv_r(r_c,θ,φ) (15)

式中γ是一阶参数，一般取γ＝0.5。

将上述的组合声压代入到球谐函数展开波束形成算法中如下：

式中α_l表示第l个阵元的权重。

下面结合具体数值给出本发明的仿真算例：

仿真中采用了半径a＝0.3m的64阵元随机均匀分布空心球阵列(α_l＝4π/L)，以球阵几何中心为坐标原点建立球坐标系。声源距离球阵中心1m，球坐标系下的坐标为(1m,89°,160°)，包络面半径为0.9m的球面，200个单极子点声源均匀分布在包络面上。空气中声速c＝340m/s，方位角为[0,360°]和俯仰角为[0,180°]扫面间隔均为1°。双半径空心球阵中的参数设置为r_z＝1.2a，虚拟空心球阵列(阵元具有心形指向性)中的外层虚拟球阵半径设置为r_z＝0.301m(两个球阵相同方位阵元间距0.001m)。

图4a-c给出了随着声源频率变化时，几种球阵形式不同阶径向函数(n＝0,1,2,3)的幅值的对比图，给出了真实空心球阵径向函数的前两个零点位置。从图4(a)中可以看出，真实球阵的径向函数存在零点。从图4(b)和4(c)可以看出，采用双半径空心球阵和虚拟空心球阵列(阵元具有心形的指向性)的径向函数在此频段范围内均不存在零点。

图5a-c和图6a-c分别给出了在ka＝2和ka＝6时，分别采用真实空心球阵、虚拟双半径空心球阵联合处理以及虚拟心形传声器球形阵列联合处理方法的噪声源定位结果。由图中可以看出采用三种不同的球阵形式均能达到良好的定位效果。

图7a-c给出了在ka＝3.14(真实空心球阵径向函数的第一个零点对应的频率)时，分别采用真实空心球阵、虚拟双半径空心球阵联合处理以及虚拟心形传声器球形阵列联合处理方法的噪声源定位结果。由图7(a)可以看出当ka＝3.14时，由于真实空心球阵的径向函数出现零点，导致噪声源定位方法失效，无法正确给出噪声源的位置。相比而言，基于阵列扩展技术的其他两种方法均可以有效给出噪声源所在位置，有效地克服了真实球阵零点处噪声源定位结果不稳定的问题。而且对比图7(b)、5(c)发现，采用虚拟双半径空心球阵联合处理方法的旁瓣较高，相比之下，采用(声压+振速)联合处理方式得到虚拟空心球阵列(阵元具有心形指向性)具有较低的旁瓣。

图8a-c给出了ka＝4.49(真实空心球阵径向函数的第二个零点对应的频率)条件下分别采用真实空心球阵、虚拟双半径空心球阵联合处理以及虚拟空心球阵列(阵元具有心形的指向性)的噪声源定位结果。由图8(a)所示，真实空心球阵在其径向函数的零点处噪声源定位结果失效，无法正确给出噪声源的真实位置。相比之下，采用两种虚拟球阵可以正确的给出噪声源的位置，克服了真实空心球阵在该频率下无法识别噪声源的问题。对比图8(b)、(c)可知，采用虚拟空心球阵列(阵元具有心形的指向性)可以有效抑制旁瓣，相比虚拟双半径球阵具有更低的旁瓣，更适合应用于实际环境中。

本发明通过虚拟阵列扩展技术，可以灵活采用不同的虚拟阵元处理方式即虚拟双半径空心球阵和虚拟空心球阵列(阵元具有心形指向性)，结合球谐波展开波束形成方法，有效地解决了真实空心球阵中径向函数零点频率出噪声源定位结果不稳定的问题。同时我们发现采用虚拟空心球阵列(阵元具有心形的指向性)对声源定位时的效果更好，具有良好的应用前景。

综上，本发明将基于波束形成正则化矩阵的虚拟阵列扩展技术与空心球阵相结合，构建虚拟大孔径空心球阵列、虚拟双半径空心球阵列和虚拟心形传声器球阵列，采用球谐函数展开波束形成算法，提出了一种适用于球形阵列的三维空间噪声源定位方法。该方法首先获取球形阵列的接收信号p。其次，选取包络面的位置，并在包络面上配置等效虚源，计算虚源和球阵间的声场传递矩阵G(r_l,r'_k)。根据波束形成正则化矩阵L＝diag(G^Hp/||G^Hp||_∞)^-1求解虚源强度q_BF。根据虚源和虚拟阵元间的传递矩阵G(r'_k,r_m)和虚源强度q_BF即可获得包络面内半径为r的球面上的声压近似值,即虚拟空心球阵列的近似接收声压。并通过选择不同的虚拟球面半径和处理方式获得：①获得虚拟双半径空心球阵声压的组合形式pZ(ka,kr,θ_l,φ_l)；②获得虚拟空心球阵列(阵元具有心形指向性)的组合声压p_E(kr_c,θ,φ)。最后采用球谐函数展开波束形成算法实现噪声源定位与识别。仿真分析结果表明，本发明可有效地获得虚拟大孔径空心球阵列、虚拟双半径空心球阵列及虚拟空心球阵列(阵元具有心形指向性)，从而克服了单层空心球阵列径向函数零点问题，提高了噪声源定位识别精度；同时，相较于传统的双球层阵列噪声源定位识别方法，本发明测量系统简单，成本低，具有良好的工程应用前景。

Claims (4)

1.一种基于虚拟球阵列扩展技术的噪声源定位识别方法，其特征在于：步骤如下：

步骤一：获取半径为a阵元数为L的球阵列接收声压向量p，其中第l个阵元为：p(ka,θ_l,φ_l))；

在球阵外部和声源之间选取半径为r'的球面作为虚拟源面并在该面上设置K个虚拟源及位置，计算虚拟源和球阵列的传递矩阵G(a_l,r'_k)，利用p和G(a_l,r'_k)计算得到波束形成正则化矩阵L，选取正则化参数λ，得到虚拟源的源强q_BF；

在虚拟源面和球阵之间选取一个球面并设置M个虚拟阵元及位置，得到虚拟源和虚拟球阵间的声场传递矩阵G(r'_k,r_m)，最后通过q_BF和G(r'_k,r_m)计算得到虚拟球阵接收声压向量p'，其中第m个虚拟阵元为：p'(kr,θ_m,φ_m))；

步骤二：基于虚拟球阵列扩展技术的球谐波展开噪声源定位方法。

2.根据权利要求1所述的一种基于虚拟球阵列扩展技术的噪声源定位识别方法，其特征在于：步骤二具体包括：

(1)将虚拟阵列和真实阵列接收数据进行球函数分解，按照‘分阶-选择’联合处理方法得到虚拟双半径空心球阵列的组合声压，采用球函数展开波束形成算法进行噪声源定位；

选取半径为r的球面作为虚拟球阵位置且虚拟阵元方位与实际阵列阵元方向相同，将其与实际的球阵视为整体构成一个虚拟的双半径空心球阵，利用选择参数β_n将p'(kr,θ_l,φ_l)以及p(ka,θ_l,φ_l)按照球谐函数分阶选择处理获得虚拟双半径空心球阵组合声压p_Z(ka,kr,θ_l,φ_l)；将p_Z(ka,kr,θ_l,φ_l)代入到球谐波展开波束形成算法中实现噪声源定位；

(2)将虚拟阵列和真实阵列相同方位阵元作有限差分处理求得两个阵元中间位置r_c处振速v(kr_c,θ_l,φ_l)，并采用‘p+v’联合处理方式得到虚拟空心球阵列的组合声压进而采用球函数展开波束形成算法进行噪声源定位；

在一个半径r的虚拟球面上设置虚拟阵元，使虚拟空心球阵的阵元数目及方位均与真实阵列相同，通过两个球阵相同方位的阵元作有限差分求得两个阵元中间位置r_c处的振速v(kr_c,θ_l,φ_l)，同时通过步骤一得到振速位置处的声压p'(kr_c,θ_l,φ_l)，选取一阶参数γ，将v(kr_c,θ_l,φ_l)和p'(kr_c,θ_l,φ_l)通过‘p+v’的联合处理获得其声压的组合形式p_E(kr_c,θ_l,φ_l)，最后将p_E(kr_c,θ_l,φ_l)代入到球函数展开波束形成算法中实现噪声源定位。

3.根据权利要求2所述的一种基于虚拟球阵列扩展技术的噪声源定位识别方法，其特征在于：步骤二中的(1)具体包括：

在真实空心球阵外部选择一个半径为r_z的虚拟球面，并在球面上设置虚拟阵元数目及位置，构成一个虚拟的大孔径空心球阵列，通过选择参数β_n在j_n(ka)和j_n(kr)之间选取最佳的径向函数来克服单个空心球阵径向函数的零点问题，选择参数表达式为：

声压的组合形式可以表示为：

p_Z(ka,kr,θ_l,φ_l)＝(1-β_n)p_l(ka,θ_l,φ_l)+β_np'_l(kr,θ_l,φ_l)

将上述组合声压代入到球谐函数展开波束形成算法中得到：

4.根据权利要求3所述的一种基于虚拟球阵列扩展技术的噪声源定位识别方法，其特征在于：步骤二中的(2)具体包括：

在一个半径r_z大于a的虚拟球面上设置虚拟阵元，保证虚拟空心球阵的阵元数目和方位与真实阵列的阵元相同，通过两个球阵相同方位的阵元作有限差分求得两个阵元中间位置处的振速：

同时采用第二步过程计算振速点位置的声压p(kr_c,θ,φ)，最终得到新的虚拟球阵阵元声压的组合形式为：

p_E(r_c,θ,φ)＝γp'(r_c,θ,φ)-(1-γ)ρcv_r(r_c,θ,φ)

式中γ是一阶参数，一般取γ＝0.5；

将上述的组合声压代入到球谐函数展开波束形成算法中如下：