[发明专利]一种多球阵列多宽带声源快速定向方法

权利要求书

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

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

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

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

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

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

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

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

步骤1具体实现方法为，

任一开放球的表面都均匀的分布着L个全指向的声压传感器，第l-th个声压传感器所在的空间方位为Ω_l＝(θ_l,φ_l)(l＝1,2,…,L)，且在笛卡尔坐标系中第d-th个球上的第l-th个声压传感器的位置矢量表示为r_dl＝(r_dlsinθ_lcosφ_l,r_dlsinθ_lsinφ_l,r_dlcosθ_l)^T，其中r_dl表示从球心到第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个平面波的波矢被表示为k_s＝-(k_ssinθ_scosφ_s,k_ssinθ_ssinφ_s,k_scosθ_s)^T；

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

x_d(k)＝A_d(k)s(k)+n_d(k) (1)

其中，且x_d(k)＝[x_d1(k),x_d2(k),…x_dL(k)]^T表示第d-th个球阵列的观察向量，且s(k)＝[s₁(k),s₂(k),…s_S(k)]^T表示信号波形向量，且n_d(k)＝[n_d1(k),n_d2(k),…,n_dL(k)]^T表示和信号向量s(k)无关的加性噪声向量,表示导向矩阵,写为

步骤2具体实现方法为，

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

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

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

A_d(k)＝Y(Ω_L)B(kr_d)Y^H(Ω_S) (5)

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

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

B(kr_d)＝diag(b₀(kr_d),b₁(kr_d),b₁(kr_d),…,b_N(kr_d)) (7)

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

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

V＝diag(α₁,α₂,…,α_L) (9)

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

给式(1)观察向量x_d(k)左乘Y^H(Ω_L)×V，得观察信号的球傅里叶变换系数x_dnm(k)，写为：

x_dnm(k)＝B(kr_d)Y^H(Ω_s)s(k)+n_dnm(k) (10)

步骤3具体实现方法为，

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

b_n(kr)＝4πjⁿ×max[|j_n(kr₁)|,|j_n(kr₂)|,…,|j_n(kr_D)|] (11)

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

r＝r_d,|j_n(kr_d)|＝max[|j_n(kr₁)|,|j_n(kr₂)|,…,|j_n(kr_D)|] (12)

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

x_nm(k)＝B(kr)Y^H(Ω_s)s(k)+n_nm(k) (13)

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

步骤4具体实现方法为，

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

P_nm(k)＝Y^H(Ω_s)s(k)+N(k) (14)

其中，且N(k)＝B^-1(kr)×Y^H(Ω_L)×V×n(k)；定义矩阵令

F(k)＝B^-1(kr)×Y^H(Ω_L)×V (15)

因此N(k)可被重写为

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

步骤5具体实现方法为，

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

P_nm＝[P_nm(k₁),P_nm(k₂),...P_nm(k_X)] (17)

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

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

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

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

P_nm＝[P_A；P_B] (22)

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