专利摘要

本发明公开的一种多球阵列多宽带声源快速定向方法，属于声源定向探测技术领域。本发明基于离散傅里叶变换，建立多球阵列声场的频率空间域模型；基于球傅里叶变换，构建多球阵列声场的球谐域模型；选择每一阶球谐阶数处模态强度最大值对应的球谐系数得到多球阵列融合后的球谐系数；去除球谐系数中和频率相关分量得到只包含角度相关分量的球傅里叶变换成分；构建球傅里叶变换成分的交叉谱，将球傅里叶变换成分或者其交叉谱进行分块，结合线性回归的方法得到传播算子的估计值；基于传播算子的估计值构建正交化的噪声子空间，利用信号和噪声相互正交得到波达方向的空间谱；扫描入射声源方向，空间谱中最大的峰值对应的声源方向即为声源的波达方向角。

权利要求

1.一种多球阵列多宽带声源快速定向方法，其特征在于：包括如下步骤，

步骤1：基于离散傅里叶变换，建立多球阵列声场的频率空间域模型；所述的多球阵列由D个同心的开放球组成，其中D≥2，定义D个球阵列的半径分别为r1,r2,…,rD；

步骤2：基于球傅里叶变换，构建多球阵列声场的球谐域模型；

步骤3：选择每一阶球谐阶数处的模态强度最大值对应的球谐系数得到多球阵列融合后的球谐系数；

步骤4：去除球谐系数中和频率相关分量得到只包含角度相关分量的球傅里叶变换成分；

步骤5：构建球傅里叶变换成分的交叉谱，将球傅里叶变换成分或者其交叉谱进行分块，结合线性回归的方法得到传播算子的估计值；

步骤6：基于传播算子的估计值构建正交化的噪声子空间，利用信号和噪声相互正交得到波达方向的空间谱；

步骤7：扫描入射声源方向Ωs，空间谱中最大的峰值对应的声源方向Ω即为声源的波达方向角；

步骤1具体实现方法为，

任一开放球的表面都均匀的分布着L个全指向的声压传感器，第l-th个声压传感器所在的空间方位为Ωl＝(θl,φl)(l＝1,2,…,L)，且在笛卡尔坐标系中第d-th个球上的第l-th个声压传感器的位置矢量表示为rdl＝(rdlsinθlcosφl,rdlsinθlsinφl,rdlcosθl)T，其中rdl表示从球心到第d-th个球上的第l-th个声压传感器的距离；其中，所有的方位角φ均在XOY平面沿X轴逆时针测得，所有的俯仰角θ均沿Z轴向下测得；

定义有S个波数为k的平面波入射到阵列，其中，k＝2πf/c，f表示声源频率，c表示声速，且第s-th个平面波的波达方向为Ωs＝(θs,φs)(s＝1,2…S)，那么第s-th个平面波的波矢被表示为ks＝-(kssinθscosφs,kssinθssinφs,kscosθs)T；

通过离散傅里叶变换获得声压传感器接收信号的频率空间域模型，写为：

xd(k)＝Ad(k)s(k)+nd(k)(1)

其中，且xd(k)＝[xd1(k),xd2(k),…xdL(k)]T表示第d-th个球阵列的观察向量，且s(k)＝[s1(k),s2(k),…sS(k)]T表示信号波形向量，且nd(k)＝[nd1(k),nd2(k),…,ndL(k)]T表示和信号向量s(k)无关的加性噪声向量,表示导向矩阵,写为

步骤2具体实现方法为，

在球谐域中式(2)中的元素被表示为

其中，n表示从0到N的球谐分解阶数，N表示最高的球谐分解阶数；m表示维度范围从-n到n；bn(krd)＝4πin×jn(krd)表示开放球的远场模态强度，其中，jn(krd)表示球贝塞尔函数；Ynm表示n阶m维度的球谐函数，且

其中，Pnm(cosθ)表示缔合勒让德函数，由式(3)导向矩阵Ad(k)写为：

Ad(k)＝Y(ΩL)B(krd)YH(ΩS) (5)

其中，表示阵列结构的球谐波矩阵，写为：

表示远场模态强度的对角矩阵，写为：

B(krd)＝diag(b0(krd),b1(krd),b1(krd),…,bN(krd)) (7)

表示入射声源的球谐波矩阵，写为：

定义矩阵表示权重因子的对角矩阵，写为：

V＝diag(α1,α2,…,αL)(9)

其中，αl是一个将方程由积分转换为求和的近似系数，对于均匀分布的采样方式有αl＝4π/L；

给式(1)观察向量xd(k)左乘YH(ΩL)×V，得观察信号的球傅里叶变换系数xdnm(k)，写为：

xdnm(k)＝B(krd)YH(Ωs)s(k)+ndnm(k) (10)

步骤3具体实现方法为，

在多开放球阵列中，不同半径的小球对应的模态强度bn(kr)只是对自变量kr中的半径r进行了缩放，因此通过选取合适的半径比，并且在每一个阶数n和每一个波数k处选择各开放球模态强度bn(kr)的最大值，即能够避免开放球阵列模态强度中的零点；实际相当于对多球阵列的球傅里叶系数向量进行滤波，通过将各个阶数n和各个波数k处的模态强度bn(kr)的最大值挑选出来融合成为一个新的单球阵列的球傅里叶系数向量；新的融合后的单球阵列的模态强度bn(kr)被写为：

bn(kr)＝4πjn×max[|jn(kr1)|,|jn(kr2)|,…,|jn(krD)|](11)

其中，max[·]表示选择方括号中的最大值，模态强度bn(kr)的自变量kr中的r代表了多球阵列半径rd中的任意一个值，且有

r＝rd,|jn(krd)|＝max[|jn(kr1)|,|jn(kr2)|,…,|jn(krD)|] (12)

利用式(11)将式(10)融合可得多球阵列的球傅里叶变换系数向量为：

xnm(k)＝B(kr)YH(Ωs)s(k)+nnm(k)(13)

其中，B(kr)由式(11)中的bn(kr)作为其对角元素构成，表示多球阵列远场模态强度的对角矩阵；nnm(k)＝YH(ΩL)×V×n(k),且n(k)表示加性噪声向量；

步骤4具体实现方法为，

给式(13)左乘B-1(kr)即可将频率相关的分量去除获得只包含角度相关的分量，即球傅里叶变换成分，写为：

Pnm(k)＝YH(Ωs)s(k)+N(k)(14)

其中，且N(k)＝B-1(kr)×YH(ΩL)×V×n(k)；定义矩阵令

F(k)＝B-1(kr)×YH(ΩL)×V(15)

因此N(k)可被重写为

N(k)＝F(k)×n(k)(16)

步骤5具体实现方法为，

对于宽带声源信号，其频率覆盖了整个信号带宽；根据波数k＝2πf/c的范围由声信号的带宽决定，现将声信号分解成X个子频带，则基于式(14)中的球傅里叶变换成分Pnm(k)宽带信号的球傅里叶变换成分被写为：

Pnm＝[Pnm(k1),Pnm(k2),...Pnm(kX)](17)

因此X个子频带的交叉谱被写为：

其中，σ2表示噪声方差，和分别表示信号交叉谱和噪声交叉谱，定义如下：

注意到式(20)中的噪声交叉谱并非空间白噪声，因此对其白化后的交叉谱为：

传播算子的定义是基于对球傅里叶变换成分Pnm或者其交叉谱的分块获得，其分块方式如下：

Pnm＝[PA；PB] (22)

式(22)中表示由矩阵Pnm的前S行组成的子矩阵，其中U＝(N+1)2-S表示由矩阵Pnm的后U行组成的子矩阵；式(23)中表示由矩阵的前S列组成的子矩阵，表示由矩阵的后U列组成的子矩阵；此时，传播算子和的估计值通过最小二乘法获得，分别表示为：

2.如权利要求1所述的一种多球阵列多宽带声源快速定向方法，其特征在于：

步骤6具体实现方法为，

构建矩阵其中或则有QHPnm＝0或者成立；因此矩阵Q的向量空间就会正交于方向向量a(Ωs)；其中，a(Ωs)表示导向向量，由YH(Ωs)或者的列向量组成；

为了引入噪声空间的投影算子，对矩阵Q正交化得其正交化矩阵Q0，即:

Q0＝Q(QHQ)-1/2(26)

则波达方向估计的空间谱为：

3.如权利要求1所述的一种多球阵列多宽带声源快速定向方法，其特征在于：还包括步骤8：根据步骤1建立的多球阵列声场的频率空间域模型，步骤2到步骤4得到的声场的球谐域模型，以及步骤5到步骤7的传播算子方法对声源的入射方向进行估计，进而提高多个宽带声源信号波达方向估计精度，减少算法运算量，解决实际工程问题。

4.如权利要求3所述的一种多球阵列多宽带声源快速定向方法，其特征在于：步骤8所述的解决实际工程问题包括：(1)球形阵列的尺寸设计；(2)提高多目标定向的空间分辨率和定向精度；(3)重建声场被映射到球谐域后其频率分量和角度分量分离，便于采用频率分析方法解决宽带问题。

说明书

技术领域

本发明涉及一种多宽带声源定向方法，尤其涉及一种基于多球阵列的多宽带源快速定向方法，属于声源定向探测技术领域。

背景技术

被动声探测技术是阵列信号处理应用中的重要发展方向之一，其主要任务是发现声源目标并确定声源波达方向(Direction of Arrival,DOA)，是目标探测的一个重要手段和技术分支。在很多实际的应用场合中多目标共存的情况较为普遍，同时声源信号频带较宽，不同信号所处频段可能不同。因此，要求被动声探测系统能够同时对多个宽带目标的方位角、俯仰角进行准确定向，进而估计目标位置。

通常的解决宽带信号定向估计的方法是把宽带信号转换为窄带信号再应用窄带信号的相关理论进行处理。该方法不但运算复杂，而且当各窄带频率范围内信噪比不同时，定向性能将会大幅下降。球形声压传感器阵列(Spherical Microphone Arrays,SMA)已被广泛应用于声场再现，波束形成和声源波达方向估计，其中最著名的高分辨率方法是基于球谐分解的多重信号分类算法(Spherical Harmonic Multiple Signal Classification,SH-MUSIC)。与传统的多重信号分类(Multiple Signal Classification,MUSIC)算法相比，SH-MUSIC对三维空间中的DOA估计性能具有显著的提升。然而，当声压传感器数目很多时进行特征值分解获得特征向量需要消耗巨大的计算量。

发明内容

针对目前对多宽带声源信号定向中存在的运算复杂、定位精度不足的问题。本发明公开的一种多球阵列多宽带声源快速定向方法要解决的技术问题是，解决现有多宽带声源定向方法中存在的上述问题，进而提高声源定向方法的运算速度和定向精度；同时相比于现有的单球阵列定向技术，所设计的多球阵列可增加阵列可处理的声源信号的频带宽度，更有利于对宽带声源信号进行波达方向估计。

本发明的目的是通过以下技术方案实现的：

本发明公开的一种多球阵列多宽带声源快速定向方法，首先通过对多球阵列各阵元接收到的声场声压信号进行离散傅里叶变换获得声场的频率空间域模型；其次通过球傅里叶变换(Spherical Fourier Transform，SFT)将频率相关的分量和角度相关的分量相互解耦得到球谐系数，通过选择每一阶球谐阶数处的模态强度最大值对应的球谐系数获得多球阵列融合后的球谐系数；然后将球谐系数中和频率相关分量去除得到只包含角度相关分量的球傅里叶变换成分，构建球傅里叶变换成分的交叉谱，将球傅里叶变换成分或者其交叉谱进行分块，结合线性回归的方法得到传播算子的估计值；接着基于传播算子的估计值构建正交化的噪声子空间，利用信号和噪声相互正交得到波达方向的空间谱；最后扫描入射声源方向，空间谱中最大的峰值对应的声源方向即为声源的波达方向角。

本发明公开的一种多球阵列多宽带声源快速定向方法，通过将具有更宽的可处理频率范围的多球阵列和具有更快速的波达方向估计能力的正交传播算子方法(Orthonormal Propagator Method,OPM)相结合的方法，有效减小声源定向算法的运算量，提高定向精度，同时增加阵列可处理的声源信号的频带宽度。

所述的线性回归的方法使用最小二乘法。

所述的传播算子的估计值可以通过球傅里叶变换成分及其交叉谱两种方法获得，分别称为球傅立叶变换成分传播算子算法(SH-SFT-OPM)和交叉谱传播算子算法(SH-CSM-OPM)。

本发明公开的一种多球阵列多宽带声源快速定向方法，包括如下步骤：

步骤1：基于离散傅里叶变换，建立多球阵列声场的频率空间域模型。

所述的多球列由D个同心的开放球组成(D≥2)，定义D个球阵列的半径分别为r1,r2,…,rD。任一开放球的表面都均匀的分布着L个全指向的声压传感器，第l-th个声压传感器所在的空间方位为Ωl＝(θl,φl)(l＝1,2,…,L)，且在笛卡尔坐标系中第d-th个球上的第l-th个声压传感器的位置矢量表示为rdl＝(rdlsinθlcosφl,rdlsinθlsinφl,rdlcosθl)T，其中rdl表示从球心到第d-th个球上的第l-th个声压传感器的距离。其中，所有的方位角φ均在XOY平面沿X轴逆时针测得，所有的俯仰角θ均沿Z轴向下测得。

定义有S个波数为k的平面波入射到阵列，其中，k＝2πf/c，f表示声源频率，c表示声速，且第s-th个平面波的波达方向为Ωs＝(θs,φs)(s＝1,2…S)，那么第s-th个平面波的波矢被表示为ks＝-(kssinθscosφs,kssinθssinφs,kscosθs)T。

通过离散傅里叶变换获得声压传感器接收信号的频率空间域模型，写为：

xd(k)＝Ad(k)s(k)+nd(k)(1)其中， 且xd(k)＝[xd1(k),xd2(k),…xdL(k)]T表示第d-th个球阵列的观察向量， 且s(k)＝[s1(k),s2(k),…sS(k)]T表示信号波形向量， 且nd(k)＝[nd1(k),nd2(k),…,ndL(k)]T表示和信号向量s(k)无关的加性噪声向量, 表示导向矩阵,写为

步骤2：基于球傅里叶变换，ad1(构)建ad多2(球)阵列ad声S(场)的球谐域模型。

在球谐域中式(2)中的元素被表示为

其中， n表示从0到N的球谐分解阶数，N表示最高的球谐分解阶数。m表示维度范围从-n到n。bn(krd)＝4πin×jn(krd)表示开放球的远场模态强度，其中，jn(krd)表示球贝塞尔函数。Ynm表示n阶m维度的球谐函数，且

其中，Pnm(cosθ)表示缔合勒让德函数，由式(3)导向矩阵Ad(k)写为：

Ad(k)＝Y(ΩL)B(krd)YH(ΩS)(5)其中， 表示阵列结构的球谐波矩阵，写为：

表示远场模态强度的对角矩阵，写为：

B(krd)＝diag(b0(krd),b1(krd),b1(krd),…,bN(krd))(7)

表示入射声源的球谐波矩阵，写为：

定义矩阵 表示权重因子的对角矩阵，写为：

V＝diag(α1,α2,…,αL)(9)

其中，αl是一个将方程由积分转换为求和的近似系数，对于均匀分布的采样方式有αl＝4π/L。

给式(1)观察向量xd(k)左乘YH(ΩL)×V，得观察信号的球傅里叶变换系数xdnm(k)，写为：

xdnm(k)＝B(krd)YH(Ωs)s(k)+ndnm(k)(10)

步骤3：选择每一阶球谐阶数处的模态强度最大值对应的球谐系数得到多球阵列融合后的球谐系数。

在多开放球阵列中，不同半径的小球对应的模态强度bn(kr)只是对自变量kr中的半径r进行了缩放，因此通过选取合适的半径比，并且在每一个阶数n和每一个波数k处选择各开放球模态强度bn(kr)的最大值，即能够避免开放球阵列模态强度中的零点。该过程实际相当于对多球阵列的球傅里叶系数向量进行滤波，通过将各个阶数n和各个波数k处的模态强度bn(kr)的最大值挑选出来融合成为一个新的单球阵列的球傅里叶系数向量。新的融合后的单球阵列的模态强度bn(kr)被写为：

bn(kr)＝4πjn×max[|jn(kr1)|,|jn(kr2)|,…,|jn(krD)|](11)

其中，max[·]表示选择方括号中的最大值，模态强度bn(kr)的自变量kr中的r代表了多球阵列半径rd中的任意一个值，且有

r＝rd,|jn(krd)|＝max[|jn(kr1)|,|jn(kr2)|,…,|jn(krD)|](12)

利用式(11)将式(10)融合可得多球阵列的球傅里叶变换系数向量为：

xnm(k)＝B(kr)YH(Ωs)s(k)+nnm(k)(13)

其中，B(kr)由式(11)中的bn(kr)作为其对角元素构成，表示多球阵列远场模态强度的对角矩阵；nnm(k)＝YH(ΩL)×V×n(k),且n(k)表示加性噪声向量。

步骤4：去除球谐系数中和频率相关分量得到只包含角度相关分量的球傅里叶变换成分。

给式(13)左乘B-1(kr)即可将频率相关的分量去除获得只包含角度相关的分量，即球傅里叶变换成分，写为：

Pnm(k)＝YH(Ωs)s(k)+N(k)(14)

其中， 且N(k)＝B-1(kr)×YH(ΩL)×V×n(k)。定义矩阵 令

F(k)＝B-1(kr)×YH(ΩL)×V(15)

因此N(k)可被重写为

N(k)＝F(k)×n(k)(16)

步骤5：构建球傅里叶变换成分的交叉谱，将球傅里叶变换成分或者其交叉谱进行分块，结合线性回归的方法得到传播算子的估计值。

对于宽带声源信号，其频率覆盖了整个信号带宽。根据波数k＝2πf/c的范围由声信号的带宽决定，现将声信号分解成X个子频带，则基于式(14)中的球傅里叶变换成分Pnm(k)宽带信号的球傅里叶变换成分 被写为：

Pnm＝[Pnm(k1),Pnm(k2),...Pnm(kX)](17)

因此X个子频带的交叉谱 被写为：

其中，σ2表示噪声方差， 和 分别表示信号交叉谱和噪声交叉谱，定义如下：

注意到式(20)中的噪声交叉谱 并非空间白噪声，因此对其白化后的交叉谱 为：

传播算子的定义是基于对球傅里叶变换成分Pnm或者其交叉谱 的分块获得，其分块方式如下：

Pnm＝[PA；PB](22)

式(22)中 表示由矩阵Pnm的前S行组成的子矩阵， 其中U＝(N+1)2-S表示由矩阵Pnm的后U行组成的子矩阵；式(23)中 表示由矩阵 的前S列组成的子矩阵， 表示由矩阵 的后U列组成的子矩阵；此时，传播算子 和 的估计值通过最小二乘法获得，分别表示为：

步骤6：基于传播算子的估计值构建正交化的噪声子空间，利用信号和噪声相互正交得到波达方向的空间谱。

构建矩阵 其中 或 则有QHPnm＝0或者 成立。因此矩阵Q的向量空间就会正交于方向向量a(Ωs)。其中，a(Ωs)表示导向向量，由YH(Ωs)或者 的列向量组成。

为了引入噪声空间的投影算子，对矩阵Q正交化得其正交化矩阵Q0，即:

Q0＝Q(QHQ)-1/2(26)

则波达方向估计的空间谱为：

步骤7：扫描入射声源方向Ωs，空间谱中最大的峰值对应的声源方向Ω即为声源的波达方向角。

还包括步骤8：根据步骤1建立的多球阵列声场的频率空间域模型，步骤2到步骤4得到的声场的球谐域模型，以及步骤5到步骤7的传播算子方法对声源的入射方向进行估计，进而提高多个宽带声源信号波达方向估计精度，减少算法运算量，解决实际工程问题。

步骤8所述的解决实际工程问题包括：(1)球形阵列的尺寸设计；(2)提高多目标定向的空间分辨率和定向精度；(3)重建声场被映射到球谐域后其频率分量和角度分量分离，便于采用频率分析方法解决宽带问题。

有益效果：

1、本发明公开的一种多球阵列多宽带声源快速定向方法，采用多球阵列结构。该结构中小半径的球阵列适用于处理高频率的声源信号，大半径的球阵列适用于处理低频率的声源信号；通过将多个半径不同的开放球以同心的方式置于同一个声场中，能够增加阵列可处理的声源信号频率范围；同时采用信号融合的方式提取出最佳球谐系数，还可以提高声源定向精度。

2、本发明公开的一种多球阵列多宽带声源快速定向方法，采用球傅里叶变换将声场映射到球谐域后其频率分量和角度分量相互解耦，然后将频率分量去除获得球傅里叶变换成分，有利于解决宽带声源定向的问题。

3、本发明公开的一种多球阵列多宽带声源快速定向方法，采用的球谐正交传播算子算法不需要对传感器接收信号的交叉谱矩阵进行特征值分解，避免现有的基于球谐分解的多重信号分类算法在多声压传感器时对特征值分解获得特征向量所消耗的巨大计算。其中，交叉谱传播算子算法(SH-CSM-OPM)和现有的基于球谐分解的多重信号分类算法(SH-MUSIC)的计算增益的阶数为O(S/(N+1)2)。而球傅立叶变换成分传播算子算法(SH-SFT-OPM)不需要进行步骤5中的构建球傅里叶变换成分的交叉谱，因此和SH-CSM-OPM，SH-MUSIC两种算法相比，该算法计算复杂度最小。

附图说明

图1多球阵列多宽带声源快速定向方法流程图；

图2为本发明的双球阵列几何模型示意图；

图3是具体实施方式中三种波达方向估计算法的空间谱，其中俯仰角和方位角的扫描间隔为(1°,2°)，(a)基于SH-SFT-OPM算法，(b)基于SH-CSM-OPM算法，(c)基于SH-MUSIC算法；

图4是具体实施方式中三个不同配置的球形阵列对应的以频率f为自变量的球傅里叶变换成分的信噪比幅度；

图5不同半径的双球阵列对应的可处理频率范围。

具体实施方式

为了更好的说明本发明的目的和优点，下面结合附图和实施例对本发明内容做进一步说明。

实施例1

步骤1：基于离散傅里叶变换，建立多球阵列声场的频率空间域模型。

所述的多球阵列由两个同心的开放球组成，定义两个球阵列的半径分别为r1＝0.1m，r2＝0.14m。任一开放球的表面都均匀的分布着L＝32个全指向的声压传感器，所述双球阵列的几何模型示意图如图2所示。最高球谐分解阶数 此处N＝4。

定义三个非相干的远场正弦线性扫频信号，三个声信号频率范围从1.2kHz到3.2kHz并以8kHz/s频率变化率线性扫频。所述的三

一种多球阵列多宽带声源快速定向方法专利购买费用说明

Q：办理专利转让的流程及所需资料

A：专利权人变更需要办理著录项目变更手续，有代理机构的，变更手续应当由代理机构办理。

1：专利变更应当使用专利局统一制作的“著录项目变更申报书”提出。

2：按规定缴纳著录项目变更手续费。

3：同时提交相关证明文件原件。

4：专利权转移的，变更后的专利权人委托新专利代理机构的，应当提交变更后的全体专利申请人签字或者盖章的委托书。

Q:专利著录项目变更费用如何缴交

A：（1）直接到国家知识产权局受理大厅收费窗口缴纳,（2）通过代办处缴纳，（3）通过邮局或者银行汇款,更多缴纳方式

Q：专利转让变更，多久能出结果

A：著录项目变更请求书递交后，一般1-2个月左右就会收到通知，国家知识产权局会下达《转让手续合格通知书》。