CN107884741A - 一种多球阵列多宽带声源快速定向方法 - Google Patents
一种多球阵列多宽带声源快速定向方法 Download PDFInfo
- Publication number
- CN107884741A CN107884741A CN201711031828.3A CN201711031828A CN107884741A CN 107884741 A CN107884741 A CN 107884741A CN 201711031828 A CN201711031828 A CN 201711031828A CN 107884741 A CN107884741 A CN 107884741A
- Authority
- CN
- China
- Prior art keywords
- array
- spherical
- sphere
- matrix
- sound source
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 47
- 238000003491 array Methods 0.000 title abstract description 10
- 238000001228 spectrum Methods 0.000 claims abstract description 51
- 238000012417 linear regression Methods 0.000 claims abstract description 6
- 239000011159 matrix material Substances 0.000 claims description 54
- 238000000354 decomposition reaction Methods 0.000 claims description 11
- 230000001419 dependent effect Effects 0.000 claims description 10
- 239000000654 additive Substances 0.000 claims description 6
- 230000000996 additive effect Effects 0.000 claims description 6
- 238000005070 sampling Methods 0.000 claims description 4
- 238000000638 solvent extraction Methods 0.000 claims description 4
- 238000004458 analytical method Methods 0.000 claims description 3
- 230000000903 blocking effect Effects 0.000 claims description 3
- 239000004576 sand Substances 0.000 claims description 3
- 238000001914 filtration Methods 0.000 claims description 2
- 238000001514 detection method Methods 0.000 abstract description 6
- 238000005516 engineering process Methods 0.000 abstract description 4
- 230000004927 fusion Effects 0.000 abstract description 3
- 230000009466 transformation Effects 0.000 abstract 4
- 230000000875 corresponding effect Effects 0.000 abstract 2
- 230000002596 correlated effect Effects 0.000 abstract 1
- 238000004422 calculation algorithm Methods 0.000 description 31
- 238000004364 calculation method Methods 0.000 description 4
- 238000012545 processing Methods 0.000 description 4
- 238000007635 classification algorithm Methods 0.000 description 3
- 230000009977 dual effect Effects 0.000 description 3
- 239000011324 bead Substances 0.000 description 2
- 230000002969 morbid Effects 0.000 description 2
- 229940037201 oris Drugs 0.000 description 2
- 238000013459 approach Methods 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 230000002349 favourable effect Effects 0.000 description 1
- 238000009434 installation Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 239000008188 pellet Substances 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 238000004513 sizing Methods 0.000 description 1
- 230000003595 spectral effect Effects 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000002087 whitening effect Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S3/00—Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received
Landscapes
- Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Circuit For Audible Band Transducer (AREA)
Abstract
本发明公开的一种多球阵列多宽带声源快速定向方法，属于声源定向探测技术领域。本发明基于离散傅里叶变换，建立多球阵列声场的频率空间域模型；基于球傅里叶变换，构建多球阵列声场的球谐域模型；选择每一阶球谐阶数处模态强度最大值对应的球谐系数得到多球阵列融合后的球谐系数；去除球谐系数中和频率相关分量得到只包含角度相关分量的球傅里叶变换成分；构建球傅里叶变换成分的交叉谱，将球傅里叶变换成分或者其交叉谱进行分块，结合线性回归的方法得到传播算子的估计值；基于传播算子的估计值构建正交化的噪声子空间，利用信号和噪声相互正交得到波达方向的空间谱；扫描入射声源方向，空间谱中最大的峰值对应的声源方向即为声源的波达方向角。
Description
技术领域
本发明涉及一种多宽带声源定向方法，尤其涉及一种基于多球阵列的多宽带源快速定向方法，属于声源定向探测技术领域。
背景技术
被动声探测技术是阵列信号处理应用中的重要发展方向之一，其主要任务是发现声源目标并确定声源波达方向(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频率变化率线性扫频。所述的三个声信号的入射方向分别为(120°，90°),(90°,180°)和(60°,270°)。麦克风接收信号的信噪比设置为0dB。快拍数为2048个，采样率为20.48kHz，即将声源信号分为X＝200个10Hz带宽的子频带。
通过离散傅里叶变换获得声压传感器接收信号的频率空间域模型，写为：
x1(k)＝A1(k)s(k)+n1(k) (28)
x2(k)＝A2(k)s(k)+n2(k) (29)
其中，且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：基于球傅里叶变换，构建多球阵列声场的球谐域模型。
给式(28)和式(29)观察向量xd(k)左乘YH(ΩL)×V，得观察信号的球傅里叶变换系数xdnm(k)，写为：
x1nm(k)＝B(kr1)YH(Ωs)s(k)+n1nm(k) (30)
x2nm(k)＝B(kr2)YH(Ωs)s(k)+n2nm(k) (31)
步骤3：选择每一阶球谐阶数处的模态强度最大值对应的球谐系数得到多球阵列融合后的球谐系数。
新的融合后的单球阵列的模态强度bn(kr)被写为：
bn(kr)＝4πjn×max[|jn(kr1)|,|jn(kr2)|] (32)
其中，max[·]表示选择方括号中的最大值，模态强度bn(kr)的自变量kr中的r代表了双球阵列半径r1，r2中的任意一个值，且有
r＝rd,|jn(krd)|＝max[|jn(kr1)|,|jn(kr2)|] (33)
利用式(32)将式(30)和式(31)融合可得双球阵列的球傅里叶变换系数向量为：
xnm(k)＝B(kr)YH(Ωs)s(k)+nnm(k) (34)
其中，B(kr)由式(32)中的bn(kr)作为其对角元素构成，表示双球阵列远场模态强度的对角矩阵；nnm(k)＝YH(ΩL)×V×n(k),且n(k)表示加性噪声向量。
步骤4：去除球谐系数中和频率相关分量得到只包含角度相关分量的球傅里叶变换成分。
给式(34)左乘B-1(kr)即可将频率相关的分量去除获得只包含角度相关的分量，即球傅里叶变换成分，写为：
Pnm(k)＝YH(Ωs)s(k)+N(k) (35)
步骤5：构建球傅里叶变换成分的交叉谱，将球傅里叶变换成分或者其交叉谱进行分块，结合线性回归的方法得到传播算子的估计值。
基于式(35)中的球傅里叶变换成分Pnm(k)宽带信号的球傅里叶变换成分被写为：
Pnm＝[Pnm(k1),Pnm(k2),...Pnm(k200)] (36)
因此X＝200个子频带的交叉谱被写为：
其中，σ2表示噪声方差，和分别表示信号交叉谱和噪声交叉谱，定义如下：
注意到式(39)中的噪声交叉谱并非空间白噪声，因此对其白化后的交叉谱为：
传播算子的定义是基于对球傅里叶变换成分Pnm或者其交叉谱的分块获得，其分块方式如下：
Pnm＝[PA；PB] (41)
式(41)中表示由矩阵Pnm的前S行组成的子矩阵，其中U＝(N+1)2-S表示由矩阵Pnm的后U行组成的子矩阵；式(42)中表示由矩阵的前S列组成的子矩阵，表示由矩阵的后U列组成的子矩阵；此时，传播算子和的估计值通过最小二乘法获得，分别表示为：
步骤6：基于传播算子的估计值构建正交化的噪声子空间，利用信号和噪声相互正交得到波达方向的空间谱。
构建矩阵其中或则有QHPnm＝0或者成立。因此矩阵Q的向量空间就会正交于方向向量a(Ωs)。其中，a(Ωs)表示导向向量，由YH(Ωs)或者的列向量组成。
为了引入噪声空间的投影算子，对矩阵Q正交化得其正交化矩阵Q0，即:
Q0＝Q(QHQ)-1/2 (45)
则波达方向估计的空间谱为：
步骤7：扫描入射声源方向Ωs，空间谱中最大的峰值对应的声源方向Ω即为声源的波达方向角。图3给出了SH-SFT-OPM算法SH-CSM-OPM算法和SH-MUSIC算法三种波达方向估计算法的空间谱，其中俯仰角和方位角的扫描间隔为(1°,2°)。通过比较可以得出：图3(a)中使用SH-SFT-OPM算法的波达方向估计准确性最差；图3(b)和图3(c)中分别使用SH-CSM-OPM算法和SH-MUSIC算法的波达方向估计结果非常接近；和图3(a)中使用SH-SFT-OPM算法得到的空间谱相比，所述两种算法的空间谱的谱峰均很尖锐。
基于MATLAB软件测试使用双球阵列的SH-SFT-OPM,SH-CSM-OPM和SH-MUSIC三种波达方向估计算法的计算时间。其方位角和俯仰角扫描间隔分别设置为(5°,10°),(2°,4°)和(1°,2°)。下表为经过十次实验的平均计算时间。
由上表结果可以得出，随着扫描间隔变小三种算法波达方向估计的计算时间快速的增大。表明所述三种算法的主要的计算量都集中在步骤7扫描入射声源方向的过程中。其中SH-SFT-OPM算法的计算时间最少，因为该算法不需要计算交叉谱也不需要噪声白化。而SH-CSM-OPM算法和SH-MUSIC算法的计算时间非常相近。
实施例2
定义球傅里叶变换成分的信噪比(Signal-to-Noise Ratio of the SFTcomponents，SNR-SFT)，且SNR-SFT被写为：
其中，||·||F表示Frobenius范数，Pnm(k)和分别表示球傅里叶变换成分的理论值和估计值。根据实施例1中的定向算法程序执行的结果，所述的SNR-SFT参数可对球麦克风阵列的DOA估计性能进行评价。
使用三个非相干的远场正弦线性扫频信号，三个声信号频率范围从0.1kHz到9.1kHz并以10kHz/s频率变化率线性扫频。所述的三个声信号的入射方向分别为(120°,90°),(90°,180°)和(60°,270°)。麦克风接收信号的信噪比设置为0dB。快拍数为2048个，采样率为20.48kHz，即将声源信号分为900个10Hz带宽的子频带。
使用的三个球形阵列的配置：阵列1使用了半径为r＝0.1m的单开放球阵列，阵列2使用了半径为r＝0.14m的单开放球，阵列3使用了半径为r＝0.1m&0.14m的双开放球阵列，且在每个球的表面都以相同的方式均匀分布了32个全指向麦克风。最高球谐分解阶数此处N＝4。
将所述的参数代入步骤1到步骤4，得到球傅里叶变换成分的估计值然后将球傅里叶变换成分的估计值代入式(47)，所得结果如图4所示。图4中给出了三个不同配置的球形阵列对应的以频率f为自变量的球傅里叶变换成分的信噪比幅度。该图结果表明所述三个阵列的性能具有明显的差异。单开放球阵列1和阵列2中，在一些固定的频点处对应的SNR-SFT的值出现了病态的偏小值，所述的现象是由于在这些固定的频点处对应的模态强度为零引起的。同时由于麦克风的安装误差以及传感器引入的噪声等因素的存在将会导致阵列性能的鲁棒性非常差。而阵列3所使用的双开放球阵列，通过选择合适的半径比避免了模态强度bn(kr)中的零值，因此SNR-SFT的值中不再有病态的偏小值。此外，双球阵列中的大半径球阵列降低了可处理频率范围的下边界，小半径球阵列增加了可处理频率范围的上边界，即和单开放球阵列相比双开放球阵列具有更宽的可处理频率范围。
将SFT-SNR大于-5dB所对应的频率值定义为阵列可处理频率范围。图5中给出了不同半径的双球阵列对应的可处理频率范围。所述的分析使用的双球阵列的配置和阵列3的配置均相同只是改变了球阵列的半径大小，双球阵列的半径比ρ仍然等于1.4。考虑到实际使用的环境将小球的半径限定为0.05米到0.5米的范围内。其它所有仿真条件均和前面图4中所使用的相同。从图5中可以看出，小球半径越小可处理的频率范围越高，小球半径越大可处理的频率范围越低。与所述的结论相一致。此外，针对不同频率范围的目标声源，可以用所述的结果来指导球形阵列的尺寸设计。
以上所述的具体描述，对发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
Claims (9)
1.一种多球阵列多宽带声源快速定向方法，其特征在于：包括如下步骤，
步骤1：基于离散傅里叶变换，建立多球阵列声场的频率空间域模型；
步骤2：基于球傅里叶变换，构建多球阵列声场的球谐域模型；
步骤3：选择每一阶球谐阶数处的模态强度最大值对应的球谐系数得到多球阵列融合后的球谐系数；
步骤4：去除球谐系数中和频率相关分量得到只包含角度相关分量的球傅里叶变换成分；
步骤5：构建球傅里叶变换成分的交叉谱，将球傅里叶变换成分或者其交叉谱进行分块，结合线性回归的方法得到传播算子的估计值；
步骤6：基于传播算子的估计值构建正交化的噪声子空间，利用信号和噪声相互正交得到波达方向的空间谱；
步骤7：扫描入射声源方向Ωs，空间谱中最大的峰值对应的声源方向Ω即为声源的波达方向角。
2.如权利要求1所述的一种多球阵列多宽带声源快速定向方法，其特征在于：还包括步骤8：根据步骤1建立的多球阵列声场的频率空间域模型，步骤2到步骤4得到的声场的球谐域模型，以及步骤5到步骤7的传播算子方法对声源的入射方向进行估计，进而提高多个宽带声源信号波达方向估计精度，减少算法运算量，解决实际工程问题。
3.如权利要求2所述的一种多球阵列多宽带声源快速定向方法，其特征在于：步骤8所述的解决实际工程问题包括：(1)球形阵列的尺寸设计；(2)提高多目标定向的空间分辨率和定向精度；(3)重建声场被映射到球谐域后其频率分量和角度分量分离，便于采用频率分析方法解决宽带问题。
4.如权利要求1、2或3所述的一种多球阵列多宽带声源快速定向方法，其特征在于：
步骤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)无关的加性噪声向量, 表示导向矩阵,写为
5.如权利要求4所述的一种多球阵列多宽带声源快速定向方法，其特征在于：
步骤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)
6.如权利要求5所述的一种多球阵列多宽带声源快速定向方法，其特征在于：
步骤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)表示加性噪声向量。
7.如权利要求6所述的一种多球阵列多宽带声源快速定向方法，其特征在于：
步骤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)
8.如权利要求7所述的一种多球阵列多宽带声源快速定向方法，其特征在于：
步骤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列组成的子矩阵；此时，传播算子和的估计值通过最小二乘法获得，分别表示为：
9.如权利要求8所述的一种多球阵列多宽带声源快速定向方法，其特征在于：
步骤6具体实现方法为，
构建矩阵其中或则有QHPnm＝0或者成立；因此矩阵Q的向量空间就会正交于方向向量a(Ωs)；其中，a(Ωs)表示导向向量，由YH(Ωs)或者的列向量组成；
为了引入噪声空间的投影算子，对矩阵Q正交化得其正交化矩阵Q0，即:
Q0＝Q(QHQ)-1/2 (26)
则波达方向估计的空间谱为：
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201711031828.3A CN107884741B (zh) | 2017-10-30 | 2017-10-30 | 一种多球阵列多宽带声源快速定向方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201711031828.3A CN107884741B (zh) | 2017-10-30 | 2017-10-30 | 一种多球阵列多宽带声源快速定向方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107884741A true CN107884741A (zh) | 2018-04-06 |
| CN107884741B CN107884741B (zh) | 2021-01-19 |
Family
ID=61782840
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201711031828.3A Expired - Fee Related CN107884741B (zh) | 2017-10-30 | 2017-10-30 | 一种多球阵列多宽带声源快速定向方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107884741B (zh) |
Cited By (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108594166A (zh) * | 2018-04-19 | 2018-09-28 | 广东工业大学 | 一种二维波达方向估计方法及装置 |
| CN109671439A (zh) * | 2018-12-19 | 2019-04-23 | 成都大学 | 一种智能化果林鸟害防治设备及其鸟类定位方法 |
| CN110006516A (zh) * | 2019-03-25 | 2019-07-12 | 中国船舶重工集团公司第七一五研究所 | 一种光纤水听器阵列的灵敏度快速标定装置及方法 |
| CN110133579A (zh) * | 2019-04-11 | 2019-08-16 | 南京航空航天大学 | 适用于球面麦克风阵列声源定向的球谐波阶数自适应选择方法 |
| CN110907892A (zh) * | 2019-12-05 | 2020-03-24 | 扬州大学 | 一种球麦克风阵列语音信号到达角估计方法 |
| CN111812581A (zh) * | 2020-06-16 | 2020-10-23 | 重庆大学 | 基于原子范数的球面阵列声源波达方向估计方法 |
| CN112698263A (zh) * | 2020-11-10 | 2021-04-23 | 重庆邮电大学 | 一种基于正交传播算子的单基地互质mimo阵列doa估计算法 |
| CN113138367A (zh) * | 2020-01-20 | 2021-07-20 | 中国科学院上海微系统与信息技术研究所 | 一种目标定位方法、装置、电子设备及存储介质 |
| CN114527427A (zh) * | 2022-01-27 | 2022-05-24 | 华南理工大学 | 一种基于球形麦克风阵列的低频波束形成声源定位方法 |
| CN115049813A (zh) * | 2022-08-17 | 2022-09-13 | 南京航空航天大学 | 一种基于一阶球谐的粗配准方法、装置及系统 |
| CN115061089A (zh) * | 2022-05-12 | 2022-09-16 | 苏州清听声学科技有限公司 | 一种声源定位方法、系统、介质、设备及装置 |
Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102855409A (zh) * | 2012-09-20 | 2013-01-02 | 辽宁工业大学 | 近场径向干扰源抑制算法及其球麦克风阵列语音采集装置 |
| CN102866385A (zh) * | 2012-09-10 | 2013-01-09 | 上海大学 | 一种基于球麦克风阵列的多声源定位方法 |
| CN103546221A (zh) * | 2013-10-25 | 2014-01-29 | 东南大学 | 一种宽带相干信号波达角估计方法 |
| CN103592628A (zh) * | 2013-11-12 | 2014-02-19 | 上海大学 | 一种基于球谐域实值权重波束形成的多声源定位方法 |
| CN103777197A (zh) * | 2013-12-24 | 2014-05-07 | 南京航空航天大学 | 单基地mimo雷达中降维传播算子的方位估计方法 |
| CN104407318A (zh) * | 2014-11-20 | 2015-03-11 | 中国船舶重工集团公司第七〇五研究所 | 一种48元均匀圆柱阵超短基线精确定向方法 |
| CN105301563A (zh) * | 2015-11-10 | 2016-02-03 | 南京信息工程大学 | 一种基于一致聚焦变换最小二乘法的双声源定位方法 |
| CN106596088A (zh) * | 2016-12-13 | 2017-04-26 | 东南大学 | 基于近场声源聚焦定位的碰摩声发射故障位置识别方法 |
-
2017
- 2017-10-30 CN CN201711031828.3A patent/CN107884741B/zh not_active Expired - Fee Related
Patent Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102866385A (zh) * | 2012-09-10 | 2013-01-09 | 上海大学 | 一种基于球麦克风阵列的多声源定位方法 |
| CN102855409A (zh) * | 2012-09-20 | 2013-01-02 | 辽宁工业大学 | 近场径向干扰源抑制算法及其球麦克风阵列语音采集装置 |
| CN103546221A (zh) * | 2013-10-25 | 2014-01-29 | 东南大学 | 一种宽带相干信号波达角估计方法 |
| CN103592628A (zh) * | 2013-11-12 | 2014-02-19 | 上海大学 | 一种基于球谐域实值权重波束形成的多声源定位方法 |
| CN103777197A (zh) * | 2013-12-24 | 2014-05-07 | 南京航空航天大学 | 单基地mimo雷达中降维传播算子的方位估计方法 |
| CN104407318A (zh) * | 2014-11-20 | 2015-03-11 | 中国船舶重工集团公司第七〇五研究所 | 一种48元均匀圆柱阵超短基线精确定向方法 |
| CN105301563A (zh) * | 2015-11-10 | 2016-02-03 | 南京信息工程大学 | 一种基于一致聚焦变换最小二乘法的双声源定位方法 |
| CN106596088A (zh) * | 2016-12-13 | 2017-04-26 | 东南大学 | 基于近场声源聚焦定位的碰摩声发射故障位置识别方法 |
Non-Patent Citations (3)
| Title |
|---|
| XI PAN等: "Multiple Spherical Arrays Design for Acoustic Source localization", 《2016 SENSOR SIGNAL PROCESSING FOR DEFENCE (SSPD)》 * |
| 汤辰: "基于球阵列测量的近场波束形成技术研究", 《中国优秀硕士学位论文全文数据库工程科技Ⅱ辑》 * |
| 黄青华等: "基于球谐域MVDR波束形成技术的远场多", 《新型工业化》 * |
Cited By (18)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108594166B (zh) * | 2018-04-19 | 2022-03-25 | 广东工业大学 | 一种二维波达方向估计方法及装置 |
| CN108594166A (zh) * | 2018-04-19 | 2018-09-28 | 广东工业大学 | 一种二维波达方向估计方法及装置 |
| CN109671439A (zh) * | 2018-12-19 | 2019-04-23 | 成都大学 | 一种智能化果林鸟害防治设备及其鸟类定位方法 |
| CN109671439B (zh) * | 2018-12-19 | 2024-01-19 | 成都大学 | 一种智能化果林鸟害防治设备及其鸟类定位方法 |
| CN110006516A (zh) * | 2019-03-25 | 2019-07-12 | 中国船舶重工集团公司第七一五研究所 | 一种光纤水听器阵列的灵敏度快速标定装置及方法 |
| CN110133579A (zh) * | 2019-04-11 | 2019-08-16 | 南京航空航天大学 | 适用于球面麦克风阵列声源定向的球谐波阶数自适应选择方法 |
| CN110907892A (zh) * | 2019-12-05 | 2020-03-24 | 扬州大学 | 一种球麦克风阵列语音信号到达角估计方法 |
| CN113138367A (zh) * | 2020-01-20 | 2021-07-20 | 中国科学院上海微系统与信息技术研究所 | 一种目标定位方法、装置、电子设备及存储介质 |
| CN111812581B (zh) * | 2020-06-16 | 2023-11-14 | 重庆大学 | 基于原子范数的球面阵列声源波达方向估计方法 |
| CN111812581A (zh) * | 2020-06-16 | 2020-10-23 | 重庆大学 | 基于原子范数的球面阵列声源波达方向估计方法 |
| CN112698263A (zh) * | 2020-11-10 | 2021-04-23 | 重庆邮电大学 | 一种基于正交传播算子的单基地互质mimo阵列doa估计算法 |
| CN114527427A (zh) * | 2022-01-27 | 2022-05-24 | 华南理工大学 | 一种基于球形麦克风阵列的低频波束形成声源定位方法 |
| CN114527427B (zh) * | 2022-01-27 | 2024-03-29 | 华南理工大学 | 一种基于球形麦克风阵列的低频波束形成声源定位方法 |
| CN115061089A (zh) * | 2022-05-12 | 2022-09-16 | 苏州清听声学科技有限公司 | 一种声源定位方法、系统、介质、设备及装置 |
| WO2023217082A1 (zh) * | 2022-05-12 | 2023-11-16 | 苏州清听声学科技有限公司 | 一种声源定位方法、系统、介质、设备及装置 |
| CN115061089B (zh) * | 2022-05-12 | 2024-02-23 | 苏州清听声学科技有限公司 | 一种声源定位方法、系统、介质、设备及装置 |
| CN115049813A (zh) * | 2022-08-17 | 2022-09-13 | 南京航空航天大学 | 一种基于一阶球谐的粗配准方法、装置及系统 |
| CN115049813B (zh) * | 2022-08-17 | 2022-11-01 | 南京航空航天大学 | 一种基于一阶球谐的粗配准方法、装置及系统 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107884741B (zh) | 2021-01-19 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107884741B (zh) | 一种多球阵列多宽带声源快速定向方法 | |
| CN112180329B (zh) | 一种基于阵元随机均匀分布球阵反卷积波束形成的汽车噪声源声成像方法 | |
| CN106443587B (zh) | 一种高分辨率的快速反卷积声源成像算法 | |
| CN106125041B (zh) | 基于子空间加权稀疏恢复的宽带信号源定位方法 | |
| AU2015292238B2 (en) | Planar sensor array | |
| CN106646376A (zh) | 基于加权修正参数的p范数噪声源定位识别方法 | |
| CN104166120B (zh) | 一种声矢量圆阵稳健宽带mvdr方位估计方法 | |
| CN105334508B (zh) | 一种稀疏阵列宽带波束形成的栅瓣抑制方法 | |
| CN110133579B (zh) | 适用于球面麦克风阵列声源定向的球谐波阶数自适应选择方法 | |
| CN106199505B (zh) | 一种声矢量圆阵模态域稳健方位估计方法 | |
| Yang et al. | Functional delay and sum beamforming for three-dimensional acoustic source identification with solid spherical arrays | |
| CN109459744B (zh) | 一种实现多干扰抑制的稳健自适应波束形成方法 | |
| CN110346752B (zh) | 基于互质稀疏阵的无模糊测向方法 | |
| CN106680762A (zh) | 一种基于互协方差稀疏重构的声矢量阵方位估计方法 | |
| CN114527427B (zh) | 一种基于球形麦克风阵列的低频波束形成声源定位方法 | |
| CN112285647B (zh) | 一种基于稀疏表示与重构的信号方位高分辨估计方法 | |
| CN108447499A (zh) | 一种双层圆环麦克风阵列语音增强方法 | |
| CN110837076A (zh) | 一种基于张量分解的矢量水听器阵列方位估计方法 | |
| CN106997037A (zh) | 声矢量传感器阵列空间旋转解相干到达角估计方法 | |
| CN104035064B (zh) | 适用于任意阵型的稳健宽带导向最小方差波束形成方法 | |
| Hafezi et al. | Multiple source localization in the spherical harmonic domain using augmented intensity vectors based on grid search | |
| CN106526563A (zh) | 一种基于互相关虚拟阵的五元体积阵多目标方位估计方法 | |
| CN110501669B (zh) | 一种中心对称声矢量圆阵快速空间谱压缩超分辨方位估计方法 | |
| CN106908754B (zh) | L型声矢量传感器阵列esprit解相干参数估计方法 | |
| CN109541526B (zh) | 一种利用矩阵变换的圆环阵方位估计方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20210119 |