CN101860779B - 用于球面阵的时域宽带谐波域波束形成器及波束形成方法 - Google Patents

Abstract

本发明涉及一种用于球面阵的时域宽带谐波域波束形成器及波束形成方法，该波束形成器包括：谐波变换模块、波束导向模块和波束图合成模块；所述的谐波变换模块，用于对传感器阵列接收的数据采样进行球谐波变换，得到谐波域数据；所述的波束导向模块，用于对谐波域数据进行加权与部分求和处理，输出导向后的谐波域数据；其特征在于，所述的波束图合成模块为基于FIR滤波器的波束图合成模块，波束导向后输出的各阶谐波域数据分别通过一个FIR滤波器，再对各FIR滤波器输出求和，获得波束输出时域数据。该方法对导向后的各阶谐波域数据进行FIR滤波求和处理，所述的各阶谐波域数据分别通过一个FIR滤波器，再对各FIR滤波器输出求和，获得时域波束输出信号。

Description

用于球面阵的时域宽带谐波域波束形成器及波束形成方法

技术领域

本发明涉及阵列信号处理领域，特别涉及一种针对球面阵的模式波束形成方法，尤其适用于宽带声信号的应用于球面麦克风阵列、球面水听器阵列的波束形成器及其波束形成方法。 

背景技术

波束形成处理广泛应用于麦克风阵列、声呐、雷达和无线通信等领域，一般的阵元域波束形成的处理过程为：采用空间分布的传感器阵列采集数据，然后对所采集的阵列数据进行线性加权组合处理得到一个标量波束输出，该处理器称为波束形成器。通过设计加权系数，可以使阵列接收系统响应具有方向性，因此，波束形成器可用于进行空域滤波，提高信噪比。 

近年来，文献中提出了针对球面阵的模式波束形成方法(J.Meyer and G.Elko，Ahighly scalable spherical microphone array based on an orthonormal decomposition of thesoundfield，in Proc.ICASSP，vol.2，May 2002，pp.1781-1784.)，其原理是对声场进行正交谐波分解，提取出各阶谐波，然后对各阶谐波进行加权综合，获得波束输出。与阵元域波束形成器相比，模式波束形成器处理通道减少，计算量也减少，且它可以将波束导向与加权分开，通过设计加权系数控制波束图，通过导向来控制波束主瓣方向，使得波束设计与扫描变得简便易行。 

文献中提出的模式波束形成方法通常采用窄带模型。借助于宽带信号傅里叶变换，利用设计出的窄带波束形成器加权系数，可以在频域实现宽带波束形成。但是，频域波束形成的主要缺点是：由于它是分块处理，使得波束输出与输入间不可避免具有延时，在某些要求实时收听的应用(例如：语音信号收听)中无法使用。 

在美国发明专利(G.W.Elko，R.A.Kubli，and J.Meyer，″Audio system based on atleast second-order eigenbeam，″United States Patent WO 03/061336 A1，July，2003.)中提出了一种基于谐波分解与滤波求和处理的时域宽带波束形成器框图。但，该专利中并未给出具体的实现方法，未指出采用何种滤波器，更没有给出滤波器的设计方法。而且，该方法在求取谐波时，引入了除法操作，这样，在某些频率除数特别小时，将造成波束形成器的稳健性较差。此外，它也没有考虑波束形成器的其他性能，如信噪比增益、旁瓣级等。 

发明内容

因此，本发明的目的是克服现有技术中基于窄带模型的模式波束形成方法在频域实现宽带波束形成的频域处理方法不适用于语音信号收听的缺点，从而提供了一种时域宽带谐波域波束形成器及其波束形成方法。 

为实现上述发明目的之一，本发明提供了一种用于球面阵的时域宽带谐波域波束形成器，该波束形成器包括：谐波变换模块、波束导向模块和波束图合成模块；所述的谐波变换模块，用于对传感器阵列接收的数据采样进行球谐波变换，得到谐波域数据；所述的波束导向模块，用于对谐波域数据进行加权与部分求和处理，输出导向后的谐波域数据；其特征在于，所述的波束图合成模块为基于FIR滤波器的波束图合成模块，波束导向后输出的各阶谐波域数据分别通过一个FIR滤波器，再对各FIR滤波器输出求和，获得波束输出时域数据。 

作为上述技术方案的一种改进，所述的球面阵为透明球或刚性球。 

作为上述技术方案的又一种改进，所述的波束导向模块，对于每一阶谐波域数据分为实部与虚部分别将各自由度(m＝0，…，n)谐波域数据进行单独加权求和处理，获得该阶谐波域数据导向后输出数据，其中，加权参数为： 

自由度m＝0时，此时谐波域该数据为实数，加权值为 

自由度0＜m≤n时，谐波域数据的实部加权值为： 

$2\underset{m=1}{\overset{n}{\mathrm{\Σ}}}\sqrt{\frac{(n-m)!}{(n+m)!}}{P}_n^m\left(\cos {\mathrm{\θ}}_0\right)\cos \left(m{\mathrm{\φ}}_0\right);$

谐波域数据的虚部加权值为： 

$2\underset{m=1}{\overset{n}{\mathrm{\Σ}}}\sqrt{\frac{(n-m)!}{(n+m)!}}{P}_n^m\left(\cos {\mathrm{\θ}}_0\right)\sin \left(m{\mathrm{\φ}}_0\right);$

式中，θ₀与φ₀分别是波束指向方向的垂直角与水平角， 表示伴随Legendre函数。 

作为上述技术方案的再一种改进，所述的FIR滤波器根据球面阵基阵模型的各阶模态构造，时域宽带谐波域波束形成器的各性能参数，包括：波束响应、指向性指数、白噪声增益、宽带主瓣响应随频率一致性和旁瓣大小； 

所述的基阵模型得到基阵的各阶模态b_n(ka)为： 

其中，a是球半径；k＝2πf/c表示波数，c表示声速，f表示频率； 

j_n与h_n分别是n阶球Bessel与Hankel函数，j_n′与h_n′分别是它们的导数； 

根据各阶模态b_n(ka)获得宽带谐波域阵列流形向量u(f，Θ)， 

$u(f,\mathrm{\Θ})=a(f,\mathrm{\Θ})\⊗e\left(f\right),$ 其中， 

a(f，Θ)＝[a₀(f，Θ)，…，a_n(f，Θ)，…，a_N(f，Θ)]^T， 

$a_n(f,\mathrm{\Θ})=b_n\left(\mathrm{ka}\right)\sqrt{\frac{2n+1}{4\mathrm{\π}}}{P}_n^0\left(\cos \mathrm{\Θ}\right)e^{\mathrm{j\πf}(L-1)T_s},$

$e\left(f\right)={[1,e^{-j2\mathrm{\πf}{T}_s},...,e^{-j(L-1)2\mathrm{\πf}{T}_s}]}^T,$

上式中，Θ是波束指向方向Ω₀与球坐标角Ω之间的夹角，Ω＝(θ，φ)，其中，θ与φ分别为垂直角与水平角； 

表示Kronecker积， 

表示伴随Legendre函数，T_s是采样周期，(·)^T表示转置； 

构造基于FIR滤波器的谐波域波束形成器波束响应表达式B(f，Ω)： 

B(f，Ω)＝u^T(f，Θ)h＝h^T u(f，Θ)， 

其中， 是(N+1)L×1维列向量，h_n＝[h_n1，h_n2，…，h_nL]^T是第n阶谐波对应的FIR滤波器加权系数，L是滤波器长度； 

构造宽带各向同性噪声协方差矩阵 

再构造波束形成器宽带各向同性噪声波束输出功率 

构造宽带白噪声增益BWNG、构造主瓣空间响应差异向量γ_MSRV，以及旁瓣性能向量B_SL；其中， 

单位功率宽带各向同性噪声协方差矩阵 

为： 

这里积分可以用离散求和计算，f_L与f_U分别是信号下、上边界频率， 

(·)^H表示共轭转置， 

diag{·}表示构造一个对角元素为括号中的元素的对角矩阵； 

宽带各向同性噪声波束输出功率 为： 

宽带白噪声增益BWNG为： 

其中，M为麦克风个数，为了提 取N阶球谐波，麦克风个数M需要满足M≥(N+1)²； 

主瓣空间响应差异向量γ_MSRV是N_MLK×1列向量，为：[γ_MSRV]_k+(j-1)K＝γ_MSRV(f_k，Θ_j)； 

其中，将工作频带[f_L，f_U]，主瓣区域Θ_ML进行离散化为f_k∈[f_L，f_U](k＝1，2，…，K)，Θ_j∈Θ_ML(j＝1，…，N_ML)；γ_MSRV(f，Θ)＝|h^Tu(f，Θ)-h^Tu(f₀，Θ)|；f₀是选择的参考频率，它可以是工作频带[f_L，f_U]内的任意频率； 

旁瓣性能向量B_SL是N_SLK×1列向量为：[B_SL]_k+(i-1)K＝B(f_k，Θ_i)； 

其中，将旁瓣区域Θ_SL离散化为Θ_i∈Θ_SL(i＝1，…，N_SL)。 

为实现本发明的另一目的，本发明还提供了一种用于球面阵的时域宽带谐波域波束形成方法，包括如下步骤： 

1)对通过传感器阵列接收的数据采样，若第s(s＝1，…，M)个传感器接收到的宽带信号用x_s(t)表示，其中t表示时间，则对传感器阵列数据进行采样后数据为 

进行球谐波变换，得到谐波域数据；对该宽带数据进行实数域谐波变换，得到的谐波域数据实部与虚部分别表述为： 

${\tilde{x}}_{\mathrm{nm}}\left(l\right)=\underset{s=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\α}}_s{x}_s\left(l\right)\mathrm{Re}\left[Y_n^m\left({\mathrm{\Ω}}_s\right)\right],$

其中，Re(·)与Im(·)分别表示实部与虚部， 

是n阶自由度为m的球谐波： 

$Y_n^m\left(\mathrm{\Ω}\right)=Y_n^m(\mathrm{\θ},\mathrm{\φ})=\sqrt{\frac{(2n+1)}{4\mathrm{\π}}\frac{(n-m)!}{(n+m)!}}{P}_n^m\left(\cos \mathrm{\θ}\right)e^{\mathrm{im\φ}},$

α_s(s＝1，…，M)是一个实数，它使球谐波 

满足如下的离散正交条件： 

$\underset{s=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\α}}_s{Y}_{n^{\′}}^{m^{\′}}\left({\mathrm{\Ω}}_s\right){\left[Y_n^m\left({\mathrm{\Ω}}_s\right)\right]}^{*}={\mathrm{\δ}}_{n-n^{\′}}{\mathrm{\δ}}_{m-m^{\′}},$

这里δ_n-n′与δ_m-m′是克罗内克函数； 

2)对谐波域数据进行加权与部分求和处理，输出导向后的谐波域数据； 

3)对导向后的各阶谐波域数据进行FIR滤波求和处理，所述的各阶谐波域数据分别通过一个FIR滤波器，再对各FIR滤波器输出求和，获得时域波束输出信号y(l)。 

作为上述技术方案的一种改进，所述的传感器阵列为球面均匀分布的阵列，α_s≡4π/M；所述的传感器阵列为非均匀球面阵，α_s通过采用最小二乘法求解上式得到。 

作为上述技术方案的又一种改进，所述的步骤2)将谐波变换得到的谐波域数据分为实部与虚部分别将各自由度(m＝0，…，n)谐波域数据进行单独加权求和处理，进行波束导向得到导向后的各阶谐波域数据为： 

式中，Ω₀＝(θ₀，φ₀)是波束指向方向； 

其中，加权参数为： 

自由度m＝0时，此时谐波域该数据为实数，加权值为 

自由度0＜m≤n时，谐波域数据的实部加权值为： 

$2\underset{m=1}{\overset{n}{\mathrm{\Σ}}}\sqrt{\frac{(n-m)!}{(n+m)!}}{P}_n^m\left(\cos {\mathrm{\θ}}_0\right)\cos \left(m{\mathrm{\φ}}_0\right);$

谐波域数据的虚部加权值为： 

$2\underset{m=1}{\overset{n}{\mathrm{\Σ}}}\sqrt{\frac{(n-m)!}{(n+m)!}}{P}_n^m\left(\cos {\mathrm{\θ}}_0\right)\sin \left(m{\mathrm{\φ}}_0\right);$

式中，θ₀与φ₀分别是波束指向方向的垂直角与水平角， 

表示伴随Legendre函数。 

作为上述技术方案的另一种改进，所述的步骤3)对各阶谐波进行滤波求和，时域波束输出信号y(l)为： 

$y\left(l\right)=\underset{n=0}{\overset{N}{\mathrm{\Σ}}}\{x_n(l,{\mathrm{\Ω}}_0)*h_n\};$

式中，*表示卷积； 

其中，所述的FIR滤波器根据球面阵基阵模型的各阶模态构造，时域宽带谐波域波束形成器的各性能参数，包括：波束响应、指向性指数、白噪声增益、宽带主瓣响应随频率一致性和旁瓣大小； 

所述的基阵模型得到基阵的各阶模态b_n(ka)为： 

其中，a是球半径；k＝2πf/c表示波数，c表示声速，f表示频率； 

j_n与h_n分别是n阶球Bessel与Hankel函数，j_n′与h_n′分别是它们的导数； 

根据各阶模态b_n(ka)获得宽带谐波域阵列流形向量u(f，Θ)， 

$u(f,\mathrm{\Θ})=a(f,\mathrm{\Θ})\⊗e\left(f\right),$

其中，a(f，Θ)＝[a₀(f，Θ)，…，a_n(f，Θ)，…，a_N(f，Θ)]^T， 

$a_n(f,\mathrm{\Θ})=b_n\left(\mathrm{ka}\right)\sqrt{\frac{2n+1}{4\mathrm{\π}}}{P}_n^0\left(\cos \mathrm{\Θ}\right)e^{\mathrm{j\πf}(L-1)T_s},$

$e\left(f\right)={[1,e^{-j2\mathrm{\πf}{T}_s},...,e^{-j(L-1)2\mathrm{\πf}{T}_s}]}^T,$

上式中，Θ是波束指向方向Ω₀与球坐标角Ω之间的夹角，Ω＝(θ，φ)，其中，θ与φ分别为垂直角与水平角； 

表示Kronecker积， 

表示伴随Legendre函数，T_s是采样周期，(·)^T表示转置； 

构造基于FIR滤波器的谐波域波束形成器波束响应表达式B(f，Ω)： 

B(f，Ω)＝u^T(f，Θ)h＝h^Tu(f，Θ)， 

其中， 是(N+1)L×1维列向量，h_n＝[h_n1，h_n2，…，h_nL]^T是第n阶谐波对应的FIR滤波器加权系数，L是滤波器长度； 

构造宽带各向同性噪声协方差矩阵 

再构造波束形成器宽带各向同性噪声波束输出功率 

构造宽带白噪声增益BWNG、构造主瓣空间响应差异向量γ_MSRV，以及旁瓣性能向量B_SL；其中， 

单位功率宽带各向同性噪声协方差矩阵 

为： 

这里积分可以用离散求和计算，f_L与f_U分别是信号下、上边界频率， 

(·)^H表示共轭转置， 

diag{·}表示构造一个对角元素为括号中的元素的对角矩阵； 

宽带各向同性噪声波束输出功率 

为： 

宽带白噪声增益BWNG为： 

其中，M为麦克风个数，为了提取N阶球谐波，麦克风个数M需要满足M≥(N+1)²； 

主瓣空间响应差异向量γ_MSRV是N_MLK×1列向量，为：[γ_MSRV]_k+(j-1)K＝γ_MSRV(f_k，Θ_j)； 

其中，将工作频带[f_L，f_U]，主瓣区域Θ_ML进行离散化为f_k∈[f_L，f_U](k＝1，2，…，K)，Θ_j∈Θ_ML(j＝1，…，N_ML)；γ_MSRV(f，Θ)＝|h^Tu(f，Θ)-h^Tu(f₀，Θ)|；f₀是选择的参考频率，它可以是工作频带[f_L，f_U]内的任意频率； 

旁瓣性能向量B_SL是N_SLK×1列向量为： 

[B_SL]_k+(i-1)K＝B(f_k，Θ_i)； 

其中，将旁瓣区域Θ_SL离散化为Θ_i∈Θ_SL(i＝1，…，N_SL)。 

作为上述技术方案的进一步的改进，所述的FIR滤波器设计还包括宽带谐波域波束优化设计表达式，根据实际工程应用需求所提出的设计指标，包括：旁瓣、指向性、稳健性或主瓣响应随频率一致性误差，将其中任意一个指标作为代价函数，其余指标作为约束值，求解获得优化的波束形成器参数，构造时域宽带谐波域波束形成器； 

所述的宽带谐波域波束优化设计表达式包括： 

$\underset{h}{\min }{\mathrm{\μ}}_l,$ l＝{1，2，3，4}，            (1a)和 

subject to B(f_k，Ω₀)＝4π/M，k＝1，2，…，K，(1b)； 

并进一步地包括： 

${\stackrel{\‾}{P}}_{\mathrm{isoout}}\≤{\mathrm{\μ}}_1,---\left(1c\right)$

${\left|\right|{\mathrm{\γ}}_{\mathrm{MSRV}}\left|\right|}_{q_1}\≤{\mathrm{\μ}}_2,---\left(1d\right)$

${\left|\right|B_{\mathrm{SL}}\left|\right|}_{q_2}\≤{\mathrm{\μ}}_3,---\left(1e\right)$ 或/和 

BWNG^-1≤μ₄，            (1f) 

式中，q₁，q₂∈{2，∞}，分别表示L₂与L_∞范数， 

包括1个代价函数与3个用户参数； 

式(1b)是无失真约束；式(1c)是为了约束各向同性噪声波束输出功率，用于提高波束指向性指数；式(1d)是约束波束主瓣空间响应差异，用于使波束在宽带内具有近似恒定的主瓣响应；式(1e)是用于控制波束旁瓣；式(1f)是约束波束白噪声增益，用于提高波束稳健性。上述的优化表达式中，式(1a)和式(1b)是必选的，式(1c)、式(1d)、式(1e)和式(1f)可以任选其中的一个或几个，进行优化。 

作为上述技术方案的再一种改进，采用二阶锥规划方法对宽带谐波域波束优化设计表达式进行约束优化问题求解。 

本发明的优点或有益效果在于， 

本发明设计了宽带谐波域波束形成器时域实现结构及其参数设计统一优化设计方法，该参数统一优化设计法设计非常灵活，能够兼顾波束形成器多个性能指标，如指向性、旁瓣级、稳健性、宽带主瓣响应偏差等，可以选择不同约束组合构成不同的波束优化准则，具有很强的实用性和广泛应用的适应性。 

相对于现有的频域方法的分块处理，本发明的时域方法是连续处理，其波束输出更适合于语音信号收听。 

本发明将经典阵元域波束形成中的通道数M降为谐波域处理的N+1(其中M≥(N+1)²)，降低了计算量。 

本发明构造的波束合成单元FIR滤波器设计优化问题，能够兼顾波束形成器多个性能指标，如指向性、旁瓣级、稳健性、宽带主瓣响应偏差等，是一种统一优化法。 

本发明可以适当减少需要兼顾的性能指标数目，根据实际需要对代价函数与约束进行适当选择，不同的约束组合构成不同的波束优化准则，获得的波束形成器都是该统一优化法的特例，导致本发明可以应对不同的设计要求。 

本发明统一优化法设计非常灵活，具有很强的实用性和广泛应用的适应性。 

附图说明

以下，结合附图来详细说明本发明的实施例，其中： 

图1是宽带谐波域波束形成器时域实现框图。 

图2是示例32元球面麦克风阵列阵元位置图。 

图3(a)是典型实施例计算得到的波束形成器中FIR滤波器系数显示图。 

图3(b)是典型实施例计算得到的波束形成器中FIR滤波器频率响应幅度。 

图3(c)是典型实施例波束响应图。 

图3(d)是波束形成器在各频率的指向性指数与白噪声增益显示图。 

图4(a)是线性调频声源信号波形。 

图4(b)是时域宽带谐波域波束形成器输出波形。 

图5是设计的宽带恒定主瓣响应波束图。 

具体实施方式

下面结合附图和具体实施例对本发明作进一步地描述。 

本实施例的时域宽带谐波域波束形成器，该波束形成器包括如下步骤： 

步骤1)：构造时域宽带谐波域波束形成器各性能参数，包括：波束响应、指向性指数、白噪声增益、宽带主瓣响应随频率一致性和旁瓣大小等)，具体为： 

首先，根据基阵模型得到基阵的各阶模态b_n(ka)，其中a是球半径，k＝2πf/c表示波数，c表示声速，f表示频率； 

然后，根据各阶模态b_n(ka)计算宽带谐波域阵列流形向量u(f，Θ)，继而构造波束响应表达式B(f，Ω)，其中，Θ是波束指向方向Ω₀与球坐标角Ω之间的夹角，Ω＝(θ，φ)，其中，θ与φ分别为垂直角与水平角； 

再，计算宽带各向同性噪声协方差矩阵 

继而构造宽带各向同性噪声波束输出功率 

构造宽带白噪声增益BWNG、构造主瓣空间响应差异向量γ_MSRV，以及旁瓣性能向量B_SL。 

步骤2)：根据实际工程应用需求提出设计指标，包括：旁瓣、指向性、稳健性、主瓣响应随频率一致性误差等，将其中任意一个指标作为代价函数，其余指标作为约束值，构造波束形成器优化设计表达式。 

步骤3)：最后，对构造的约束优化问题求解，获得优化的波束形成器参数，即得到时域宽带谐波域波束形成器。 

步骤4)：首先通过传感器(以下不防假设为麦克风)阵列接收数据采样，进行球谐波变换，得到谐波域数据； 

步骤5)：对谐波域数据进行加权与部分求和处理，输出导向后的谐波域数据； 

步骤6)：根据设计出的波束形成各参数，对导向后的谐波域数据进行加权求和处理，获得波束输出数据； 

上述技术方案中，所述步骤1)中，以透明球与刚性球为例，其各阶模态表示为： 

其中 

j_n与h_n分别是n阶球Bessel与Hankel函数，j_n′与h_n′分别是它们的导数。 

上述技术方案中，所述步骤1)中，宽带谐波域阵列流形向量u(f，Θ)表示为： 

$u(f,\mathrm{\Θ})=a(f,\mathrm{\Θ})\⊗e\left(f\right),$

其中 表示Kronecker积，a(f，Θ)＝[a₀(f，Θ)，…，a_n(f，Θ)，…，a_N(f，Θ)]^T， 

表示伴随Legendre函数，T_s是采样周期，(·)^T表示转置。 

上述技术方案中，所述步骤1)中，波束响应表示为： 

B(f，Ω)＝u^T(f，Θ)h＝h^Tu(f，Θ)， 

其中 

是(N+1)L×1维列向量，h_n＝[h_n1，h_n2，…，h_nL]^T是第n阶谐波对应的FIR滤波器加权系数，L是滤波器长度。 

上述技术方案中，所述步骤1)中，单位功率宽带各向同性噪声协方差矩阵 

表示为： 

${\stackrel{\‾}{Q}}_{\mathrm{hiso}}={\∫}_{f_L}^{f_U}{Q}_{\mathrm{hiso}}\left(f\right).$

这里积分可以用离散求和计算，f_L与f_U分别是信号下、上边界频率， 

(·)^H表示共轭转置， diag{·}表示构造一个对角元素为括号中的元素的对角矩阵。 

上述技术方案中，所述步骤1)中，宽带各向同性噪声波束输出功率 表示为： 

${\stackrel{\‾}{P}}_{\mathrm{isoout}}=h^T{\stackrel{\‾}{Q}}_{\mathrm{hiso}}h.$

上述技术方案中，所述步骤1)中，宽带白噪声增益BWNG表示为： 

$\mathrm{BWNG}=\frac{4\mathrm{\π}/M}{h^{T}h}.$

其中M为麦克风个数，为了提取N阶球谐波，麦克风个数M需要满足 M≥(N+1)²。 

上述技术方案中，所述步骤1)中，将工作频带[f_L，f_U]，主瓣区域Θ_ML进行离散化为f_k∈[f_L，f_U](k＝1，2，…，K)，Θ_j∈Θ_ML(j＝1，…，N_ML)。主瓣空间响应差异向量γ_MSRV是N_MLK×1列向量，定义为： 

[γ_MSRV]_k+(j-1)K＝γ_MSRV(f_k，Θ_j)； 

其中， 

γ_MSRV(f，Θ)＝|h^Tu(f，Θ)-h^Tu(f₀，Θ)|； 

其中f₀是选择的参考频率，它可以是工作频带[f_L，f_U]内的任意频率。 

上述技术方案中，所述步骤1)中，将旁瓣区域Θ_SL离散化为Θ_i∈Θ_SL(i＝1，…，N_SL)，旁瓣性能向量B_SL是N_SLK×1列向量，表示为 

[B_SL]_k+(i-1)K＝B(f_k，Θ_i)。 

上述技术方案中，所述步骤2)中，构造宽带谐波域波束优化设计问题表达式为： 

$\underset{h}{\min }{\mathrm{\μ}}_l,$ l＝{1，2，3，4}，(1a) 

subject toB(f_k，Ω₀)＝4π/M，k＝1，2，…，K，(1b) 

${\stackrel{\‾}{P}}_{\mathrm{isoout}}\≤{\mathrm{\μ}}_1,---\left(1c\right)$

${\left|\right|{\mathrm{\γ}}_{\mathrm{MSRV}}\left|\right|}_{q_1}\≤{\mathrm{\μ}}_2,---\left(1d\right)$

${\left|\right|B_{\mathrm{SL}}\left|\right|}_{q_2}\≤{\mathrm{\μ}}_3,---\left(1e\right)$ 和 

BWNG^-1≤μ₄，                                (1f) 

式中q₁，q₂∈{2，∞}，分别表示L₂与L_∞范数， 

包括1个代价函数与3个用户参数。式(1b)是无失真约束；式(1c)是为了约束各向同性噪声波束输出功率，用于提高波束指向性指数；式(1d)是约束波束主瓣空间响应差异，用于使波束在宽带内具有近似恒定的主瓣响应；式(1e)是用于控制波束旁瓣；式(1f)是约束波束白噪声增益，用于提高波束稳健性。本例中，同时采用式(1c)、式(1d)、式(1e)和式(1f)分别进行优化。 

上述技术方案中，所述步骤3)中，可以采用某些数学方法对式(1a～f)所示多约束优化问题求解，例如采用二阶锥规划方法求解。 

上述技术方案中，所述步骤4)中，麦克风阵列接收数据，若第s(s＝1，…，M)个麦克风接收到的宽带信号用x_s(t)表示，其中t表示时间，对麦克风阵列数据进行采样后数据为 对该宽带数据进行实数域谐波变换，得到的谐波域数据实部与虚部分别表述为： 

${\tilde{x}}_{\mathrm{nm}}\left(l\right)=\underset{s=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\α}}_s{x}_s\left(l\right)\mathrm{Re}\left[Y_n^m\left({\mathrm{\Ω}}_s\right)\right];$

其中，Re(·)与Im(·)分别表示实部与虚部， 

是n阶自由度为m的球谐波： 

$Y_n^m\left(\mathrm{\Ω}\right)=Y_n^m(\mathrm{\θ},\mathrm{\φ})=\sqrt{\frac{(2n+1)}{4\mathrm{\π}}\frac{(n-m)!}{(n+m)!}}{P}_n^m\left(\cos \mathrm{\θ}\right)e^{\mathrm{im\φ}},$

α_s(s＝1，…，M)是一个实数，它使球谐波 

满足如下的离散正交条件： 

$\underset{s=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\α}}_s{Y}_{n^{\′}}^{m^{\′}}\left({\mathrm{\Ω}}_s\right){\left[Y_n^m\left({\mathrm{\Ω}}_s\right)\right]}^{*}={\mathrm{\δ}}_{n-n^{\′}}{\mathrm{\δ}}_{m-m^{\′}},$

这里δ_n-n′与δ_m-m′是克罗内克函数。对于麦克风在球面均匀分布的阵列α_s≡4π/M。对于非均匀球面阵，α_s可以通过采用最小二乘法求解上式得到。 

上述技术方案中，所述步骤5)中，对谐波域数据进行加权与部分求和处理，输出导向后的谐波域数据 

上述技术方案中，所述步骤6)中，对导向后谐波域数据进行波束形成，即对各阶导向后谐波进行滤波求和，时域波束输出y(l)可以写成： 

$y\left(l\right)=\underset{n=0}{\overset{N}{\mathrm{\Σ}}}\{x_n(l,{\mathrm{\Ω}}_0)*h_n\},$

式中*表示卷积。 

本发明设计的时域宽带谐波域波束形成方法可用球面麦克风阵列系统来实现，图1给出了宽带谐波域波束形成器时域实现框图。 

如图1所示，该波束形成器框图包括三部分：谐波变换模块、波束导向模块与波束图合成模块。本实施例中，前两者(变换模块与波束导向模块)与Elko等人发明(WO 03/061336 A1)的方法在结构上类似，但波束导向模块的参数不同；第三部分(基于FIR滤波器的波束合成单元)是本发明独创。 

考虑一个半径为4.2厘米的球面麦克风阵，M＝32个麦克风位于球面截断20面体中心。如图2所示，示例性给出了基阵阵元位置。 

采用该麦克风阵列采集声场数据，采样频率为f_s＝14700Hz。该球面阵麦克风近似均匀分布，可以近似取α_s＝4π/M。假设进行谐波分解时取N＝4，波束观察方向为Ω₀＝0°，FIR滤波器长度为L＝65。 

典型方案1：假设该球面阵为透明球阵，我们要设计一个稳健最大指向性指数波束形成器。一个工作频带为[f_L，f_U]＝[500，5000]Hz的线性调频信号从0°方向入射到该麦克风阵，用2°间隔对主瓣与旁瓣区域进行离散化：Θ_ML＝[0°:2°:40°]，Θ_SL＝[48°:2°:180°]。 

在步骤1)中，对于该刚性球，取 用K＝51个频率网格离散化工作频带[f_L，f_U]： 

k＝1，2，…，K。构造B(f_k，Ω₀)(k＝1，2，…，K)、 BWNG、γ_MSRV与B_SL。 

在步骤2)中，取l＝1，μ₂＝∞，μ₃＝∞，μ₄＝4π/M。按照式(1a～f)构造优化问题。本例中，同时采用式(1c)和式(1f)分别进行优化。 

在步骤3)中，采用二阶锥规划方法求解FIR滤波器h。得到的FIR滤波器系数h₀，h₁，…，h_N显示于图3(a)中。这些FIR对应的频率响应幅度显示于图3(b)中。由这些FIR滤波器构成的宽带波束形成器产生的波束图显示于图3(c)中。该波束形成器在各频率的指向性指数DI与白噪声增益WNG显示于图3(d)中。从图3(d)中可以看出，波束白噪声增益都高于-3dB，表明该波束具有较好的稳健性，其指向性指数也维持在8～14dB。 

在步骤4)中，对入射信号进行波束形成，获得的波束输出时间序列显示于图4(b)中。作为比较，入射信号波形显示于图4(a)中，比较图4(a)与图4(b)可见，波束形成器保证了输入信号无失真输出。 

典型方案2：假设该球面阵为刚性球阵，假设我们要设计一个频率不变波束形成器，信号频带为[f_L，f_U]＝[1250，5000]Hz，取Θ_ML＝[0°:2°:180°]，l＝1， μ₂＝10^-1.5·4π/M，q₁＝2，μ₃＝∞，μ₄＝2π/M。采用同样的步骤，获得的宽带波束图显示于图5中。本例中，同时采用式(1c)、式(1d)和式(1f)分别进行优化。 

最后所应说明的是，以上实施例仅用以说明本发明的技术方案而非限制。尽管参照实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，对本发明的技术方案进行修改或者等同替换，都不脱离本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围当中。 

Claims (9)

1.一种用于球面阵的时域宽带谐波域波束形成器，该波束形成器包括：谐波变换模块、波束导向模块和波束图合成模块；所述的谐波变换模块，用于对传感器阵列接收的数据采样进行球谐波变换，得到谐波域数据；所述的波束导向模块，用于对谐波域数据进行加权与部分求和处理，输出导向后的谐波域数据；其特征在于，所述的波束图合成模块为基于FIR滤波器的波束图合成模块，波束导向后输出的各阶谐波域数据分别通过一个FIR滤波器，再对各FIR滤波器输出求和，获得波束输出时域数据；

所述的FIR滤波器根据球面阵基阵模型的各阶模态构造，时域宽带谐波域波束形成器的各性能参数，包括：波束响应、指向性指数、白噪声增益、宽带主瓣响应随频率一致性和旁瓣大小；

所述的基阵模型得到基阵的各阶模态b_n(ka)为：

其中，a是球半径；k=2πf/c表示波数，c表示声速，f表示频率；

j_n与h_n分别是n阶球Bessel与Hankel函数，j_n′与h_n′分别是它们的导数；

根据各阶模态b_n(ka)获得宽带谐波域阵列流形向量u(f,Θ)， $u(f,\mathrm{\Θ})=a(f,\mathrm{\Θ})\⊗e\left(f\right),$ 其中，

a(f,Θ)=[a₀(f,Θ),…,a_n(f,Θ),…,a_N(f,Θ)]^T，

$a_n(f,\mathrm{\Θ})=b_n\left(\mathrm{ka}\right)\sqrt{\frac{2n+1}{4\mathrm{\π}}}{P}_n^0\left(\cos \mathrm{\Θ}\right)e^{\mathrm{j\πf}(L-1)T_s},$

$e\left(f\right)={[1,e^{-j2\mathrm{\πf}{T}_s},\·\·\·,e^{-j(L-1)2\mathrm{\πf}{T}_s}]}^T,$

上式中，Θ是波束指向方向Ω₀与球坐标角Ω之间的夹角，Ω=(θ，φ)，其中，θ与φ分别为垂直角与水平角；

表示Kronecker积，表示伴随Legendre函数，T_s是采样周期，(·)^T表示转置；

构造FIR滤波器的波束响应表达式B(f,Ω)：

B(f,Ω)=u^T(f,Θ)h=h^Tu(f,Θ)，

其中，

是(N+1)L×1维列向量，h_n=[h_n1，h_n2，…,h_nL]^T是第n阶谐波对应的FIR滤波器加权系数，L是滤波器长度；

构造宽带各向同性噪声协方差矩阵

再构造FIR滤波器的宽带各向同性噪声波束输出功率

构造宽带白噪声增益BWNG、构造主瓣空间响应差异向量γ_MSRV，以及旁瓣性能向量B_SL；其中，

单位功率宽带各向同性噪声协方差矩阵

为：

这里积分用离散求和计算，f_L与f_U分别是信号下、上边界频率， $Q_{\mathrm{hiso}}\left(f\right)=[I_{(N+1)\×(N+1)}\⊗e\left(f\right)]Q_{\mathrm{ciso}}\left(f\right){[I_{(N+1)\×(N+1)}\⊗e\left(f\right)]}^H,$ (·)^H表示共轭转置， $Q_{\mathrm{ciso}}\left(f\right)=\frac{1}{4\mathrm{\π}}\mathrm{diag}\{{\left|b_0\left(\mathrm{ka}\right)\right|}^2,{\left|b_1\left(\mathrm{ka}\right)\right|}^2,{\left|b_2\left(\mathrm{ka}\right)\right|}^2,\·\·\·,{\left|b_N\left(\mathrm{ka}\right)\right|}^2\},$ diag{·}表示构造一个对角元素为括号中的元素的对角矩阵；宽带各向同性噪声波束输出功率

为：

宽带白噪声增益BWNG为：

其中，M为麦克风个数，为了提取N阶球谐波，麦克风个数M需要满足M≥(N+1)²；

主瓣空间响应差异向量γ_MSRV是N_MLK×1列向量，为：[γ_MSRV]_k+(j-1)K=γ_MSRV(f_k,Θ_j)；

其中，将工作频带[f_L,f_U]，主瓣区域Θ_ML进行离散化为f_k∈[f_L,f_U](k=1,2,…,K)，Θ_j∈Θ_ML(j=1,…,N_ML)；γ_MSRV(f,Θ)=|h^Tu(f,Θ)-h^Tu(f₀,Θ)|；f₀是选择的参考频率，它是工作频带[f_L,f_U]内的任意频率；

旁瓣性能向量B_SL是N_SLK×1列向量为：[B_SL]_k+(i-1)K=B(f_k,Θ_i)；

其中，将旁瓣区域Θ_SL离散化为Θ_i∈Θ_SL(i=1,…,N_SL)。

2.根据权利要求1所述的时域宽带谐波域波束形成器，其特征在于，所述的球面阵为透明球或刚性球。

3.根据权利要求1所述的时域宽带谐波域波束形成器，其特征在于，所述的波束导向模块，对于每一阶谐波域数据分为实部与虚部分别将各自由度m谐波域数据进行单独加权求和处理，获得该阶谐波域数据导向后输出数据，其中，加权参数为：自由度m=0时，此时谐波域该数据为实数，加权值为

自由度0<m≤n时，谐波域数据的实部加权值为：

$2\stackrel{n}{\underset{m=1}{\mathrm{\&Sigma;}}}\sqrt{\frac{(n-m)!}{(n+m)!}}{P}_n^m\left(\cos {\mathrm{\&theta;}}_0\right)\cos \left(m{\mathrm{\&phi;}}_0\right);$

谐波域数据的虚部加权值为：

$2\underset{m=1}{\overset{n}{\mathrm{\&Sigma;}}}\sqrt{\frac{(n-m)!}{(n+m)!}}{P}_n^m\left(\cos {\mathrm{\&theta;}}_0\right)\sin \left(m{\mathrm{\&phi;}}_0\right);$

式中，θ₀与φ₀分别是波束指向方向的垂直角与水平角，

表示伴随Legendre函数。

4.一种用于球面阵的时域宽带谐波域波束形成方法，包括如下步骤：

1）对通过传感器阵列接收的数据采样，若第s个传感器接收到的宽带信号用x_s(t)表示，其中s=1,…,M，t表示时间，则对传感器阵列数据进行采样后数据为

进行球谐波变换，得到谐波域数据；对该宽带数据进行实数域谐波变换，得到的谐波域数据实部与虚部分别表述为：

${\tilde{x}}_{\mathrm{nm}}\left(l\right)=\underset{s=1}{\overset{M}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_s{x}_s\left(l\right)\mathrm{Re}\left[Y_n^m\left({\mathrm{\&Omega;}}_s\right)\right],$

其中，Re(·)与Im(·)分别表示实部与虚部，

是n阶自由度为m的球谐波：

$Y_n^m\left(\mathrm{\&Omega;}\right)=Y_n^m(\mathrm{\&theta;},\mathrm{\&phi;})=\sqrt{\frac{(2n+1)}{4\mathrm{\&pi;}}\frac{(n-m)!}{(n+m)!}}{P}_n^m\left(\cos \mathrm{\&theta;}\right)e^{\mathrm{im\&phi;}},$

其中，θ与φ分别为垂直角与水平角；表示伴随Legendre函数；α_s是一个实数，s=1,…,M，它使球谐波

满足如下的离散正交条件：

$\underset{s=1}{\overset{M}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_s{Y}_{n^{\&prime;}}^{m^{\&prime;}}\left({\mathrm{\&Omega;}}_s\right){\left[Y_n^m\left({\mathrm{\&Omega;}}_s\right)\right]}^{*}={\mathrm{\&delta;}}_{n-n^{\&prime;}}{\mathrm{\&delta;}}_{m-m^{\&prime;}},$

这里

与

是克罗内克函数；

2）对谐波域数据进行加权与部分求和处理，输出导向后的谐波域数据；

3）对导向后的各阶谐波域数据进行FIR滤波求和处理，即所述的各阶谐波域数据分别通过一个FIR滤波器，再对各FIR滤波器输出求和，获得时域波束输出信号y(l)；

所述的步骤3）对各阶导向后谐波进行滤波求和，时域波束输出信号y(l)为：

$y\left(l\right)=\underset{n=0}{\overset{N}{\mathrm{\&Sigma;}}}\{x_n(l,{\mathrm{\&Omega;}}_0)*h_n\};$

式中，*表示卷积；Ω₀是指波束指向方向；x_n(l,Ω₀)是指导向后的各阶谐波域数据；

其中，所述的FIR滤波器根据球面阵基阵模型的各阶模态构造，时域宽带谐波域波束形成器的各性能参数，包括：波束响应、指向性指数、白噪声增益、宽带主瓣响应随频率一致性和旁瓣大小；

所述的基阵模型得到基阵的各阶模态b_n(ka)为：

其中，a是球半径；k=2πf/c表示波数，c表示声速，f表示频率；

j_n与h_n分别是n阶球Bessel与Hankel函数，

与

分别是它们的导数；根据各阶模态b_n(ka)获得宽带谐波域阵列流形向量u(f,Θ)， $u(f,\mathrm{\&Theta;})=a(f,\mathrm{\&Theta;})\&CircleTimes;e\left(f\right),$

其中，a(f,Θ)=[a₀(f,Θ),…,a_n(f,Θ),…,a_N(f,Θ)]^T，

$a_n(f,\mathrm{\&Theta;})=b_n\left(\mathrm{ka}\right)\sqrt{\frac{2n+1}{4\mathrm{\&pi;}}}{P}_n^0\left(\cos \mathrm{\&Theta;}\right)e^{\mathrm{j\&pi;f}(L-1)T_s},$

$e\left(f\right)={[1,e^{-j2\mathrm{\&pi;f}{T}_s},\&CenterDot;\&CenterDot;\&CenterDot;,e^{-j(L-1)2\mathrm{\&pi;f}{T}_s}]}^T,$

上式中，Θ是波束指向方向Ω₀与球坐标角Ω之间的夹角，Ω=(θ，φ)，其中，θ与φ分别为垂直角与水平角；

表示Kronecker积，表示伴随Legendre函数，T_s是采样周期，(·)^T表示转置；

构造FIR滤波器的波束响应表达式B(f,Ω)：

B(f,Ω)=u^T(f,Θ)h=h^Tu(f,Θ)，

其中，

是(N+1)L×1维列向量，h_n=[h_n1，h_n2，…,h_nL]^T是第n阶谐波对应的FIR滤波器加权系数，L是滤波器长度；构造宽带各向同性噪声协方差矩阵

再构造FIR滤波器的宽带各向同性噪声波束输出功率构造宽带白噪声增益BWNG、构造主瓣空间响应差异向量γ_MSRV，以及旁瓣性能向量B_SL；其中，

单位功率宽带各向同性噪声协方差矩阵为：这里积分用离散求和计算，f_L与f_U分别是信号下、上边界频率，

$Q_{\mathrm{hiso}}\left(f\right)=[I_{(N+1)\&times;(N+1)}\&CircleTimes;e\left(f\right)]Q_{\mathrm{ciso}}\left(f\right){[I_{(N+1)\&times;(N+1)}\&CircleTimes;e\left(f\right)]}^H,$ (·)^H表示共轭转置，

$Q_{\mathrm{ciso}}\left(f\right)=\frac{1}{4\mathrm{\&pi;}}\mathrm{diag}\{{\left|b_0\left(\mathrm{ka}\right)\right|}^2,{\left|b_1\left(\mathrm{ka}\right)\right|}^2,{\left|b_2\left(\mathrm{ka}\right)\right|}^2,\&CenterDot;\&CenterDot;\&CenterDot;,{\left|b_N\left(\mathrm{ka}\right)\right|}^2\},$ diag{·}表示构造一个对角元素为括号中的元素的对角矩阵；

宽带各向同性噪声波束输出功率

为：

宽带白噪声增益BWNG为：

其中，M为麦克风个数，为了提取N阶球谐波，麦克风个数M需要满足M≥(N+1)²；

主瓣空间响应差异向量γ_MSRV是N_MLK×1列向量，为：[γ_MSRV]_k+(j-1)K=γ_MSRV(f_k,Θ_j)；

其中，将工作频带[f_L,f_U]，主瓣区域Θ_ML进行离散化为f_k∈[f_L,f_U](k=1,2,…,K)，Θ_j∈Θ_ML(j=1,…,N_ML)；γ_MSRV(f,Θ)=|h^Tu(f,Θ)-h^Tu(f₀,Θ)|；f₀是选择的参考频率，它是工作频带[f_L,f_U]内的任意频率；

旁瓣性能向量B_SL是N_SLK×1列向量为：

[B_SL]_k+(i-1)K=B(f_k,Θ_i)；

其中，将旁瓣区域Θ_SL离散化为Θ_i∈Θ_SL(i=1,…,N_SL)。

5.根据权利要求4所述的时域宽带谐波域波束形成方法，其特征在于，所述的传感器阵列为球面均匀分布的阵列，α_s≡4π/M。

6.根据权利要求4所述的时域宽带谐波域波束形成方法，其特征在于，所述的传感器阵列为非均匀球面阵，α_s通过采用最小二乘法求解离散正交条件 $\underset{s=1}{\overset{M}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_s{Y}_{n^{\&prime;}}^{m^{\&prime;}}\left({\mathrm{\&Omega;}}_s\right){\left[Y_n^m\left({\mathrm{\&Omega;}}_s\right)\right]}^{*}={\mathrm{\&delta;}}_{n-n^{\&prime;}}{\mathrm{\&delta;}}_{m-m^{\&prime;}}$ 得到。

7.根据权利要求4所述的时域宽带谐波域波束形成方法，其特征在于，所述的步骤2）将谐波变换得到的谐波域数据分为实部与虚部分别将各自由度m谐波域数据进行单独加权求和处理，进行波束导向得到导向后的各阶谐波域数据为：

式中，Ω₀＝(θ₀，φ₀)是波束指向方向；

其中，加权参数为：

自由度m＝0时，此时谐波域该数据为实数，加权值为

自由度0＜m≤n时，谐波域数据的实部加权值为：

$2\underset{m=1}{\overset{n}{\mathrm{\&Sigma;}}}\sqrt{\frac{(n-m)!}{(n+m)!}}{P}_n^m\left(\cos {\mathrm{\&theta;}}_0\right)\cos \left(m{\mathrm{\&phi;}}_0\right);$

谐波域数据的虚部加权值为：

$2\underset{m=1}{\overset{n}{\mathrm{\&Sigma;}}}\sqrt{\frac{(n-m)!}{(n+m)!}}{P}_n^m\left(\cos {\mathrm{\&theta;}}_0\right)\sin \left(m{\mathrm{\&phi;}}_0\right);$

式中，θ₀与φ₀分别是波束指向方向的垂直角与水平角，

表示伴随Legendre函数。

8.根据权利要求4所述的时域宽带谐波域波束形成方法，其特征在于，所述的FlR滤波器还包括宽带谐波域波束优化设计表达式，根据实际工程应用需求所提出的设计指标，包括：旁瓣、指向性、稳健性或主瓣响应随频率一致性误差，将其中任意一个指标作为代价函数，其余指标作为约束值，求解获得优化的波束形成器参数，构造时域宽带谐波域波束形成器；

所述的宽带谐波域波束优化设计表达式包括：

$\underset{h}{\min }{\mathrm{\&mu;}}_l,l=\left\{\mathrm{1,2,3,4}\right\},$ (1a)和

subject to B(f_k，Ω₀)＝4π/M，k＝1，2，…，K，    (1b)；

并进一步地包括：

${\stackrel{\&OverBar;}{P}}_{\mathrm{isoout}}\&le;{\mathrm{\&mu;}}_1,---\left(1c\right)$

${\left|\right|{\mathrm{\&gamma;}}_{\mathrm{MSRV}}\left|\right|}_{q_1}\&le;{\mathrm{\&mu;}}_2,---\left(1d\right)$

${\left|\right|B_{\mathrm{SL}}\left|\right|}_{q_2}\&le;{\mathrm{\&mu;}}_3,$                     (1e)或/和

BWNG^-1≤μ₄，    (1f)

式中，q₁，q₂∈{2，∞}，分别表示L₂与L_∞范数，

包括1个代价函数与3个用户参数；

式(1b)是无失真约束；式(1c)是为了约束各向同性噪声波束输出功率，用于提高波束指向性指数；式(1d)是约束波束主瓣空间响应差异，用于使波束在宽带内具有近似恒定的主瓣响应；式(1e)是用于控制波束旁瓣；式(1f)是约束波束白噪声增益，用于提高波束稳健性；式(1a)和式(1b)是必选的，式(1c)、式(1d)、式(1e)和式(1f)任选其中的一个或几个进行优化。

9.根据权利要求8所述的时域宽带谐波域波束形成方法，其特征在于，采用二阶锥规划方法对宽带谐波域波束优化设计表达式进行约束优化问题求解。

Images