CN111812581A - 基于原子范数的球面阵列声源波达方向估计方法 - Google Patents

Abstract

本发明公开了基于原子范数的球面阵列声源波达方向估计方法，步骤为：1)搭建由Q个传声器(2)构成的球形传声器(2)阵列；2)声源(1)向球形传声器(2)阵列辐射声波；3)建立声源波达方向测量模型，并构建传声器(2)测量得到的声压信号矩阵P^★；4)建立协方差矩阵

5)利用球面ESPRIT算法对协方差矩阵

进行解算，确定声源波达方向。本发明能克服球面ESPRIT在高频、相干声源或少数据快拍工况下失效的缺陷，并显著提高低SNR工况下的声源DOA估计精度，即使在存在混响的普通测试环境中，本发明仍然有效。

Description

基于原子范数的球面阵列声源波达方向估计方法

技术领域

本发明涉及声源识别领域，具体是基于原子范数的球面阵列声源波达方向估计方法。

背景技术

声源波达方向(Direction-Of-Arrival,DOA)估计问题普遍存在于噪声源识别、目标探测、故障诊断等领域。基于球面传声器阵列测量的旋转不变信号参数估计(Estimationof Signal Parameters via Rotational Invariance Technique,ESPRIT)技术，简称球面ESPRIT，凭借能全景估计声源波达方向、计算复杂度低等优势而备受关注。球面ESPRIT以阵列传声器测量信号的协方差矩阵为输入，基于球谐函数的递归关系，将声源DOA估计问题转化为最小二乘求解和特征值分解问题。现有球面ESPRIT仅适用于传声器采样的球谐函数满足正交性的情形，传声器的离散性使满足该特性的球谐函数的阶都不高，这限制了可用的上限频率。其次，作为一种子空间方法，球面ESPRIT不可避免地承受子空间方法对相干声源、少数据快拍或低信噪比(Signal-to-Noise Ratio,SNR)工况失效的固有缺陷。这些问题和缺陷成为阻碍球面ESPRIT成功解决各种声源DOA估计问题的关键障碍。针对上限频率受限的缺陷，目前尚未见解决方法的相关报道。

针对相干声源、少数据快拍或低SNR工况失效的缺陷，现有技术采用正向反向平均方法和频率光滑技术来去除源相关性。但是，而且正向反向平均方法仅适用于相干声源数目为2的情形，频率光滑技术仅适用于声源辐射宽带信号的情形。

综上所述，球面ESPRIT对高频声源、相干声源、少数据快拍或低SNR工况失效的缺陷还悬而未决。

发明内容

本发明的目的是解决现有技术中存在的问题。

为实现本发明目的而采用的技术方案是这样的，基于原子范数的球面阵列声源波达方向估计方法，包括以下步骤：

1)搭建由Q个传声器构成的球形传声器阵列。球形传声器阵列中心记为坐标原点。其中，第q个传声器位置记为

a为阵列半径，q＝1,2,…,Q。Ω＝(θ,φ)表示球形传声器阵列所在三维空间内的任意方向。θ∈[0°,180°]为仰角，φ∈[0°,360°)为方位角。

2)声源向球形传声器阵列辐射声波。

3)建立声源波达方向测量模型，并构建传声器测量得到的声压信号矩阵P^★。

建立声源波达方向测量模型，包括以下步骤：

3.1)计算第i个声源到第q个传声器的传递函数

即：

式中，n和m分别为球谐函数的阶和次。b_n(ka)为模态强度。

为Ω方向的球谐函数。Ω_Si表示第i个声源波达方向。i＝1,2,…,I。I为声源总数。k为声源辐射声波的波数。上标*表示共轭。

模态强度b_n(ka)如下所示：

式中，j_n(ka)为n阶第一类球贝塞尔函数，

为n阶第二类球汉克尔函数。j'_n(ka)和

分别为n阶第一类球贝塞尔函数j_n(ka)和n阶第二类球汉克尔函数

的一阶导数。开口球表示传声器布置在开口球体表面。刚性球表示传声器布置在刚性球体表面。

3.2)计算Ω方向的球谐函数

即：

式中，

为连带勒让德函数。β_n,m,κ为连带勒让德函数系数。球谐函数第n阶m次项对应的球谐系数

3.3)建立Q个传声器所在方向的向量

和对应的球谐函数向量

建立I个声源所在方向的向量

和对应的球谐函数向量

3.4)建立每个声源到所有传声器的传递函数矩阵，记为

上标H表示转置共轭。

其中，传声器所在方向对应的球谐函数矩阵

如下所示：

式中，N₀＝∞表示球谐函数最高阶。

声源所在方向对应的球谐函数矩阵

如下所示：

阵列模态强度矩阵

如下所示：

3.5)建立各传声器测量得到的声压信号矩阵

即：

式中，

为噪声矩阵。信噪比SNR＝20lg(||P^★-N||_F/||N||_F)。声源强度矩阵

L为快拍总数。||||_F表示Frobenius范数。上标^★表示测量量。

3.6)对球谐函数矩阵

球谐函数矩阵

和模态强度矩阵

的阶n进行截断，即令球谐函数矩阵

球谐函数矩阵

和模态强度矩阵

的最高阶

表示将数值向正无穷方向圆整到次近的整数。

基于最高阶N，更新声压信号矩阵P^★如下：

4)建立协方差矩阵

建立协方差矩阵

包括以下步骤：

4.1)建立原子范数最小化模型，包括以下步骤：

4.1.1)建立连带勒让德函数表达式，即：

式中，x为函数输入。

其中，连带勒让德函数项(x²-1)ⁿ如下所示：

式中，多项式系数

4.1.2)将公式(11)代入公式(9)，得到连带勒让德函数次方m≥0时连带勒让德函数表达式，即：

4.1.3)确定仰角θ的正弦函数sinθ＝(e^jθ-e^-jθ)/(2j)和余弦函数cosθ＝(e^jθ+e^-jθ)/2。

4.1.4)基于步骤4.1.2)和步骤4.1.3)，更新连带勒让德函数表达式如下：

4.1.5)连带勒让德函数项(e^jθ-e^-jθ)^m、连带勒让德函数项(e^jθ+e^-jθ-2)^o-m分别如下所示：

4.1.6)基于公式(14)和公式(15)，更新连带勒让德函数表达式如下：

4.1.7)给定函数阶数n和级数m，令索引o从m增至n、索引u从0增至m、索引v从0增至o-m、索引w从0增至o-m-v。

在每组(o,u,v,w)下，根据公式(16)确定连带勒让德函数的三角多项式展开中的指数索引κ＝2u+v-m-w和连带勒让德函数的三角多项式中的系数。计算完所有组(o,u,v,w)后，同一κ值对应的系数相加即为β_n,m,κ。

4.1.8)联合公式(9)和m≥0时的连带勒让德函数系数β_n,m,κ，确定m＜0时的连带勒让德函数系数β_n,m,κ。

4.1.9)确定连带勒让德函数系数β_n,m,κ后，构建Ω方向的球谐函数

即：

4.1.10)构建矩阵

和矩阵

式中，矩阵D中元素

记第(N+κ)(2N+1)+N-m+1个元素为A_n,mβ_n,m,κ，κ＝-n,…,0,…,n，并令其它元素为0，则球谐函数矩阵的转置共轭

各传声器测量得到的声压信号矩阵P^★≈Y_MNB_NGDS+N。

4.1.11)建立输入矩阵X，即：

式中，

为S的第i行。

||ψ_i||₂＝1。

公式(18)原子范数如下所示：

式中，“inf”表示下确界。

为原子集合。

原子集合

如下所示：

4.1.12)建立原子范数最小化模型，即：

其中，ε为噪声控制参数。

是连续域内声源稀疏性的一个度量。

表示原子范数最小化问题的最优解。

4.2)建立等价的半正定规划模型，包括以下步骤：

4.2.1)将式(21)转化为如下半正定规划模型：

式中，u和E为辅助量。N_u为辅助量u中元素的数目。T_b(·)为二重Toeplitz算子。对任意给定向量

是(2N,2N)的半空间，N_u＝8N²+4N+1。

4.2.2)利用二重Toeplitz算子T_b(u)将u映射为(2N+1)×(2N+1)维的块Toeplitz型Hermitian矩阵，即：

式中，每个矩阵分块T_κ均为(2N+1)×(2N+1)维的Toeplitz矩阵：0≤κ≤2N。

矩阵分块T_κ如下所示：

4.2.3)建立矩阵

的Vandermonde分解式，从而令公式(21)和公式(22)等价；

Vandermonde分解式如下所示：

式中，矩阵V＝[d(Ω_S1),d(Ω_S2),…,d(Ω_Sr)]；对角矩阵Σ＝diag([σ₁,σ₂,…,σ_r])；

i＝1,2,…,r；r为矩阵

的秩；r≤2N+1是公式(25)存在的充分条件。矩阵

视为是一组声源中，各个声源独自引起的信号的协方差矩阵的和，不包括不同声源引起的信号间的协方差成分。

4.3)使用交替方向乘子方法对半正定规划模型进行解算，包括以下步骤：

4.3.1)使用交替方向乘子方法更新半正定规划模型，得到：

式中，Z为辅助矩阵，τ为规则化参数。

4.3.2)建立公式(26)的增广拉格朗日函数表达式，即：

式中，

为Hermitian拉格朗日乘子矩阵。ρ＞0为惩罚参数。“<·,·>”表示内积。

4.3.3)利用交替方向乘子方法求解公式(26)，初始化辅助矩阵Z⁰＝Λ⁰＝0，γ+1次迭代时的变量更新为：

4.3.4)对Hermitian拉格朗日乘子矩阵和辅助矩阵进行划分，得到：

4.3.5)基于步骤3.3)和步骤3.4)，更新公式(28)如下：

式中，I₁和I₂分别为L×L维和(2N+1)²×(2N+1)²维单位矩阵。

为T_b(·)的伴随算子。对任意给定矩阵

M＝diag([(2N+1)×[2N+1,2N,…,1],

M为对角矩阵。矩阵

为第κ(m)个对角线的元素均为1其它元素均为0的基本Toeplitz矩阵。

4.3.6)基于步骤3.3)，更新公式(29)如下：

公式(36)表示将Hermitian矩阵

投影到半正定锥上，即对该Hermitian矩阵进行特征值分解，并令所有负特征值为。

4.4)基于求得的协方差矩阵和

建立协方差矩阵

5)利用球面ESPRIT算法对协方差矩阵

进行解算，确定声源波达方向，包括以下步骤：

5.1)特征值分解

并以特征值大小对特征向量进行降序排列。将前s个特征向量写入矩阵U_S中。

5.2)基于球谐函数递归关系，建立与矩阵U_S相关的线性方程组，并采用最小二乘方法对与矩阵U_S相关的线性方程组进行求解，得到含有声源波达方向的转换矩阵。线性方程组如下：

其中，

为按上标(u,v)从矩阵U_S中选取一部分行组合而成的矩阵，ψ_xy+、ψ_xy-和ψ_z为含有声源波达方向的转换矩阵，

和

为系数矩阵。

其中，系数矩阵

和系数矩阵

分别如下所示：

5.3)采用广义特征值分解法对含有声源波达方向的矩阵进行特征值分解，确定声源波达方向。

本发明的技术效果是毋庸置疑的，本发明建立了基于原子范数的新型球面ESPRIT技术，并基于仿真模拟和验证试验分析其性能，结果表明ANM+球面ESPRIT能够完美克服球面ESPRIT本发明能克服球面ESPRIT在高频、相干声源或少数据快拍工况下失效的缺陷，并显著提高低SNR工况下的声源DOA估计精度，即使在存在混响的普通测试环境中，本发明仍然有效。其次本发明中所推导的基于ADMM的求解算法相比基于IPM的SDPT3求解器更加高效。

附图说明

图1为球形传声器阵列和声源示意图；

图2为不同频率的单次蒙特卡罗计算的声源识别成像图；

图2(a)为1000Hz频率下球面ESPRIT的声源成像图；

图2(b)为1000Hz频率下ANM+球面ESPRIT的声源成像图，ANM由基于IPM的SDPT3求解器求解；

图2(c)为1000Hz频率下ANM+球面ESPRIT的声源成像图，ANM由基于ADMM的求解算法求解；

图2(d)为3000Hz频率下球面ESPRIT的声源成像图；

图2(e)为3000Hz频率下ANM+球面ESPRIT的声源成像图，ANM由基于IPM的SDPT3求解器求解；

图2(f)为3000Hz频率下ANM+球面ESPRIT的声源成像图，ANM由基于ADMM的求解算法求解；

图2(g)为1000Hz频率下球面ESPRIT的声源成像图；

图2(h)为1000Hz频率下ANM+球面ESPRIT的声源成像图，ANM由基于IPM的SDPT3求解器求解；

图2(i)为1000Hz频率下ANM+球面ESPRIT的声源成像图，ANM由基于ADMM的求解算法求解；

图3为声源DOA估计均方根误差和ANM求解耗时随频率的变化曲线；

图3(a)为100次蒙特卡罗计算时σ随频率的变化曲线

图3(b)为两种ANM求解方法耗时的对比曲线；

图4为不同声源相干性的单次蒙特卡罗计算的声源识别成像图；

图4(a)为声源互不相干时球面ESPRIT的声源成像图；

图4(b)为声源互不相干时ANM+球面ESPRIT的声源成像图，ANM由基于IPM的SDPT3求解器求解；

图4(c)为声源互不相干时ANM+球面ESPRIT的声源成像图，ANM由基于ADMM的求解算法求解；

图4(d)为声源部分相干时球面ESPRIT的声源成像图；

图4(e)为声源部分相干时ANM+球面ESPRIT的声源成像图，ANM由基于IPM的SDPT3求解器求解；

图4(f)为声源部分相干时ANM+球面ESPRIT的声源成像图，ANM由基于ADMM的求解算法求解；

图4(g)为声源完全相干时球面ESPRIT的声源成像图；

图4(h)为声源完全相干时ANM+球面ESPRIT的声源成像图，ANM由基于IPM的SDPT3求解器求解；

图4(i)为声源完全相干时ANM+球面ESPRIT的声源成像图，ANM由基于ADMM的求解算法求解；

图5为不同快拍总数不同SNR下声源DOA估计均方根误差的互补累积分布函数；

图5(a)为T＝2°、球面ESPRIT、100次蒙特卡罗计算时不同快拍总数不同SNR下声源DOA估计均方根误差的互补累积分布函数；

图5(b)为T＝2°、ANM+球面ESPRIT(ANM由基于IPM的SDPT3求解器求解)、100次蒙特卡罗计算时不同快拍总数不同SNR下声源DOA估计均方根误差的互补累积分布函数；

图5(c)为T＝2°、ANM+球面ESPRIT(ANM由基于ADMM的求解算法求解)、100次蒙特卡罗计算时不同快拍总数不同SNR下声源DOA估计均方根误差的互补累积分布函数；

图5(d)为T＝1°、球面ESPRIT、100次蒙特卡罗计算时不同快拍总数不同SNR下声源DOA估计均方根误差的互补累积分布函数；

图5(e)为T＝1°、ANM+球面ESPRIT(ANM由基于IPM的SDPT3求解器求解)、100次蒙特卡罗计算时不同快拍总数不同SNR下声源DOA估计均方根误差的互补累积分布函数；

图5(f)为T＝1°、ANM+球面ESPRIT(ANM由基于ADMM的求解算法求解)、100次蒙特卡罗计算时不同快拍总数不同SNR下声源DOA估计均方根误差的互补累积分布函数；

图6(a)为半消声室内的试验布局图；

图6(b)为半消声室内三维空间展开图；

图7为半消声室内的试验声源成像图；

图7(a)为扬声器由稳态白噪声信号激励、1500Hz频率、采用快拍总数为30时球面ESPRIT的声源成像图；

图7(b)为扬声器由稳态白噪声信号激励、1500Hz频率、采用快拍总数为30时ANM+球面ESPRIT的声源成像图，ANM由基于IPM的SDPT3求解器求解；

图7(c)为扬声器由稳态白噪声信号激励、1500Hz频率、采用快拍总数为30时ANM+球面ESPRIT的声源成像图，ANM由基于ADMM的求解算法求解；

图7(d)为扬声器由稳态白噪声信号激励、3000Hz频率、采用快拍总数为30时球面ESPRIT的声源成像图；

图7(e)为扬声器由稳态白噪声信号激励、3000Hz频率、采用快拍总数为30时ANM+球面ESPRIT的声源成像图，ANM由基于IPM的SDPT3求解器求解；

图7(f)为扬声器由稳态白噪声信号激励、3000Hz频率、采用快拍总数为30时ANM+球面ESPRIT的声源成像图，ANM由基于ADMM的求解算法求解；

图7(g)为扬声器由稳态白噪声信号激励、1500Hz频率、采用快拍总数为5时球面ESPRIT的声源成像图；

图7(h)为扬声器由稳态白噪声信号激励、1500Hz频率、采用快拍总数为5时ANM+球面ESPRIT的声源成像图，ANM由基于IPM的SDPT3求解器求解；

图7(i)为扬声器由稳态白噪声信号激励、1500Hz频率、采用快拍总数为5时ANM+球面ESPRIT的声源成像图，ANM由基于ADMM的求解算法求解；

图7(j)为扬声器由稳态白噪声信号激励、1500Hz频率、采用快拍总数为1时球面ESPRIT的声源成像图；

图7(k)为扬声器由稳态白噪声信号激励、1500Hz频率、采用快拍总数为1时ANM+球面ESPRIT的声源成像图，ANM由基于IPM的SDPT3求解器求解；

图7(l)为扬声器由稳态白噪声信号激励、1500Hz频率、采用快拍总数为1时ANM+球面ESPRIT的声源成像图，ANM由基于ADMM的求解算法求解；

图7(m)为扬声器由同一纯音信号激励、1500Hz频率、采用快拍总数为30时球面ESPRIT的声源成像图；

图7(n)为扬声器由同一纯音信号激励、1500Hz频率、采用快拍总数为30时ANM+球面ESPRIT的声源成像图，ANM由基于IPM的SDPT3求解器求解；

图7(o)为扬声器由同一纯音信号激励、1500Hz频率、采用快拍总数为30时ANM+球面ESPRIT的声源成像图，ANM由基于ADMM的求解算法求解；

图8(a)为普通房间内的试验布局图；

图8(b)为普通房间三维空间展开图；

图9为普通房间内的试验声源成像图；

图9(a)为扬声器由稳态白噪声信号激励、1000Hz频率、采用快拍总数为30时球面ESPRIT的声源成像图；

图9(b)为扬声器由稳态白噪声信号激励、1000Hz频率、采用快拍总数为30时ANM+球面ESPRIT的声源成像图，ANM由基于IPM的SDPT3求解器求解；

图9(c)为扬声器由稳态白噪声信号激励、1000Hz频率、采用快拍总数为30时ANM+球面ESPRIT的声源成像图，ANM由基于ADMM的求解算法求解；

图9(d)为扬声器由稳态白噪声信号激励、2000Hz频率、采用快拍总数为30时球面ESPRIT的声源成像图；

图9(e)为扬声器由稳态白噪声信号激励、2000Hz频率、采用快拍总数为30时ANM+球面ESPRIT的声源成像图，ANM由基于IPM的SDPT3求解器求解；

图9(f)为扬声器由稳态白噪声信号激励、2000Hz频率、采用快拍总数为30时ANM+球面ESPRIT的声源成像图，ANM由基于ADMM的求解算法求解；

图9(g)为扬声器由稳态白噪声信号激励、3000Hz频率、采用快拍总数为30时球面ESPRIT的声源成像图；

图9(h)为扬声器由稳态白噪声信号激励、3000Hz频率、采用快拍总数为30时ANM+球面ESPRIT的声源成像图，ANM由基于IPM的SDPT3求解器求解；

图9(i)为扬声器由稳态白噪声信号激励、3000Hz频率、采用快拍总数为30时ANM+球面ESPRIT的声源成像图，ANM由基于ADMM的求解算法求解；

图9(j)为扬声器由稳态白噪声信号激励、4000Hz频率、采用快拍总数为30时球面ESPRIT的声源成像图；

图9(k)为扬声器由稳态白噪声信号激励、4000Hz频率、采用快拍总数为30时ANM+球面ESPRIT的声源成像图，ANM由基于IPM的SDPT3求解器求解；

图9(l)为扬声器由稳态白噪声信号激励、4000Hz频率、采用快拍总数为30时ANM+球面ESPRIT的声源成像图，ANM由基于ADMM的求解算法求解；

图9(m)为扬声器由稳态白噪声信号激励、5000Hz频率、采用快拍总数为30时球面ESPRIT的声源成像图；

图9(n)为扬声器由稳态白噪声信号激励、5000Hz频率、采用快拍总数为30时ANM+球面ESPRIT的声源成像图，ANM由基于IPM的SDPT3求解器求解；

图9(o)为扬声器由稳态白噪声信号激励、5000Hz频率、采用快拍总数为30时ANM+球面ESPRIT的声源成像图，ANM由基于ADMM的求解算法求解；

图中，声源1、传声器2。

具体实施方式

下面结合实施例对本发明作进一步说明，但不应该理解为本发明上述主题范围仅限于下述实施例。在不脱离本发明上述技术思想的情况下，根据本领域普通技术知识和惯用手段，做出各种替换和变更，均应包括在本发明的保护范围内。

实施例1：

参见图1至图9，基于原子范数的球面阵列声源波达方向估计方法，包括以下步骤：

1)搭建由Q个传声器2构成的球形传声器阵列。球形传声器阵列中心记为坐标原点。其中，第q个传声器位置记为(a,Ω_Mq)。a为阵列半径，q＝1,2,…,Q。Ω＝(θ,φ)表示球形传声器阵列所在三维空间内的任意方向。θ∈[0°,180°]为仰角，φ∈[0°,360°)为方位角。附图1中“●”和“○”分别表示传声器和声源。

表示实数集，

表示正实数集，

表示复数集，向量用黑体小写字母表示，矩阵用黑体大写字母表示，上标“^★”表示测量量，上标“^T”、“*”和“^H”分别表示转置、共轭和转置共轭，符号

表示Kronecker乘积，符号“|·|”表示求标量的模或集合的势，广义不等式A≥0表示矩阵A半正定，任意矩阵

的Frobenius范数被定义为

任意向量

的l₂范数被定义为

tr(A)表示对矩阵A求迹，diag(a)表示形成以a中元素为对角元素的对角矩阵。

2)声源1向球形传声器阵列辐射声波。

3)建立声源波达方向测量模型，并构建传声器测量得到的声压信号矩阵P^★。

建立声源波达方向测量模型，包括以下步骤：

3.1)计算第i个声源到第q个传声器的传递函数t((ka,Ω_Mq)|Ω_Si)，即：

式中，n和m分别为球谐函数的阶和次。b_n(ka)为模态强度。

为Ω方向的球谐函数。Ω_Si表示第i个声源波达方向。i＝1,2,…,I。I为声源总数。k为声源辐射声波的波数。上标*表示共轭。

模态强度b_n(ka)如下所示：

式中，j_n(ka)为n阶第一类球贝塞尔函数，

为n阶第二类球汉克尔函数。j'_n(ka)和

分别为n阶第一类球贝塞尔函数j_n(ka)和n阶第二类球汉克尔函数

的一阶导数。开口球表示传声器布置在开口球体表面。刚性球表示传声器布置在刚性球体表面。jⁿ为符号函数。

3.2)计算Ω方向的球谐函数

即：

式中，

为连带勒让德函数。β_n,m,κ为连带勒让德函数系数。球谐函数第n阶m次项对应的球谐系数

3.3)建立Q个传声器所在方向的向量

和对应的球谐函数向量

建立I个声源所在方向的向量

和对应的球谐函数向量

Ω_MQ为第Q个传声器所在方向的向量。

表示I维复数集。

3.4)建立每个声源到所有传声器的传递函数矩阵，记为

上标H表示转置共轭。

其中，传声器所在方向对应的球谐函数矩阵

如下所示：

式中，N₀＝∞表示球谐函数最高阶。

声源所在方向对应的球谐函数矩阵

如下所示：

阵列模态强度矩阵

如下所示：

3.5)建立各传声器测量得到的声压信号矩阵

即：

式中，

为噪声矩阵。信噪比SNR＝20lg(||P^★-N||_F/||N||_F)。声源强度矩阵

L为快拍总数。||||_F表示Frobenius范数。上标^★表示测量量。

3.6)各传声器测量的声压信号包含无穷多阶球谐函数。用符号

表示将数值向正无穷方向圆整到次近的整数，

时，b_n(ka)的幅值很小，相应阶次的贡献可忽略。基于该事实，可对球谐函数的阶进行截断。对球谐函数矩阵

球谐函数矩阵

和阵列模态强度矩阵

的阶n进行截断，即令球谐函数矩阵

球谐函数矩阵

和阵列模态强度矩阵

的最高阶

表示将数值向正无穷方向圆整到次近的整数。

基于最高阶N，更新声压信号矩阵P^★如下：

4)建立协方差矩阵

建立协方差矩阵

包括以下步骤：

4.1)建立原子范数最小化模型，包括以下步骤：

4.1.1)建立连带勒让德函数表达式，即：

式中，x为函数输入。

其中，连带勒让德函数项(x²-1)ⁿ如下所示：

式中，多项式系数

4.1.2)将公式(11)代入公式(9)，得到连带勒让德函数次方m≥0时连带勒让德函数表达式，即：

4.1.3)确定仰角θ的正弦函数sinθ＝(e^jθ-e^-jθ)/(2j)和余弦函数cosθ＝(e^jθ+e^-jθ)/2。

4.1.4)基于步骤4.1.2)和步骤4.1.3)，更新连带勒让德函数表达式如下：

4.1.5)连带勒让德函数项(e^jθ-e^-jθ)^m、连带勒让德函数项(e^jθ+e^-jθ-2)^o-m分别如下所示：

式中，j为虚数。

4.1.6)基于公式(14)和公式(15)，更新连带勒让德函数表达式如下：

4.1.7)给定函数阶数n和级数m，令索引o从m增至n、索引u从0增至m、索引v从0增至o-m、索引w从0增至o-m-v。

在每组(o,u,v,w)下，根据公式(16)确定连带勒让德函数的三角多项式展开中的指数索引κ＝2u+v-m-w和连带勒让德函数的三角多项式中的系数。计算完所有组(o,u,v,w)后，同一κ值对应的系数相加即为β_n,m,κ。

4.1.8)联合公式(9)和m≥0时的连带勒让德函数系数β_n,m,κ，确定m＜0时的连带勒让德函数系数β_n,m,κ。

4.1.9)确定连带勒让德函数系数β_n,m,κ后，构建Ω方向的球谐函数

即：

4.1.10)构建矩阵

和矩阵

式中，矩阵D中元素

记第(N+κ)(2N+1)+N-m+1个元素为A_n,mβ_n,m,κ，κ＝-n,…,0,…,n，并令其它元素为0，则球谐函数矩阵的转置共轭

各传声器测量得到的声压信号矩阵P^★≈Y_MNB_NGDS+N。

4.1.11)建立输入矩阵X，即：

式中，

为S的第i行。元素

向量

||ψ_i||₂＝1。

公式(18)原子范数如下所示：

式中，“inf”表示下确界。

为原子集合。

原子集合

如下所示：

4.1.12)建立原子范数最小化模型，即：

其中，ε为噪声控制参数。

是连续域内声源稀疏性的一个度量。最小化

即对声源分布施加稀疏约束。

表示原子范数最小化问题的最优解。

4.2)建立等价的半正定规划模型，包括以下步骤：

4.2.1)将式(21)转化为如下半正定规划模型：

式中，u和E为辅助量。N_u为辅助量u中元素的数目。T_b(·)为二重Toeplitz算子。对任意给定向量

是(2N,2N)的半空间，N_u＝8N²+4N+1。

4.2.2)利用二重Toeplitz算子T_b(u)将u映射为(2N+1)×(2N+1)维的块Toeplitz型Hermitian矩阵，即：

式中，每个矩阵分块T_κ均为(2N+1)×(2N+1)维的Toeplitz矩阵：0≤κ≤2N。

矩阵分块T_κ如下所示：

4.2.3)建立矩阵

的Vandermonde分解式，从而令公式(21)和公式(22)等价；

Vandermonde分解式如下所示：

式中，矩阵V＝[d(Ω_S1),d(Ω_S2),…,d(Ω_Sr)]；对角矩阵Σ＝diag([σ₁,σ₂,…,σ_r])；

i＝1,2,…,r；r为矩阵

的秩。文献“Compressive two-dimensionalharmonic retrieval via atomic norm minimization”,Y.Chi,Y.Chen.IEEETransactions on Signal Processing,volume 63,issue No.4,pages 1030-1042,February 2015(基于原子范数最小化的二维压缩谐波检测,Y.Chi,Y.Chen.IEEETransactions on Signal Processing,第63卷,第4期,第1030-1042页,2015年2月)已证明：r≤2N+1时，式(25)成立。

是一组源中，各个源独自引起的信号的协方差矩阵的和，不包括不同源引起的信号间的协方差成分，即

去除了源间的相关性，可看作是一组不相干源引起的信号的协方差矩阵，且该特性与采用的数据快拍数目无关。由于施加约束||P^★-Y_MNB_NGX||_F≤ε，ANM具有噪声滤除功能。此外，ANM未涉及传声器采样的球谐函数的正交性。因此，本发明用

作为球面ESPRIT的输入矩阵能够克服其对高频声源、相干声源、少数据快拍或低SNR工况失效的缺陷，这里，用

而不直接用

是因为前者的前I个较大特征值对应的特征向量构成的矩阵才与

跨越相同子空间。

4.3)半正定规划为标准凸优化问题，可用基于内点方法(Interior PointMethod,IPM)的CVX工具箱中现成的SDPT3求解器求解，在本发明中给出了更加高效的基于交替方向乘子方法(Alternating Direction Method of Multipliers,ADMM)的求解算法。使用交替方向乘子方法对半正定规划模型进行解算，包括以下步骤：

4.3.1)使用交替方向乘子方法更新半正定规划模型，得到：

式中，Z为辅助矩阵，τ为规则化参数。

4.3.2)建立公式(26)的增广拉格朗日函数表达式，即：

式中，

为Hermitian拉格朗日乘子矩阵。ρ＞0为惩罚参数。“<·,·>”表示内积。

4.3.3)利用交替方向乘子方法求解公式(26)，初始化辅助矩阵Z⁰＝Λ⁰＝0，γ+1次迭代时的变量更新为：

4.3.4)对Hermitian拉格朗日乘子矩阵和辅助矩阵进行划分，得到：

4.3.5)基于步骤3.3)和步骤3.4)，更新公式(28)如下：

式中，I₁和I₂分别为L×L维和(2N+1)²×(2N+1)²维单位矩阵。

为T_b(·)的伴随算子。对任意给定矩阵

M为对角矩阵。矩阵

为第κ(m)个对角线的元素均为1其它元素均为0的基本Toeplitz矩阵。

4.3.6)基于步骤3.3)，更新公式(29)如下：

公式(36)表示将Hermitian矩阵

投影到半正定锥上，即对该Hermitian矩阵进行特征值分解，并令所有负特征值为0。

4.4)基于求得的协方差矩阵和

建立协方差矩阵

5)利用球面ESPRIT算法对协方差矩阵

进行解算，确定声源波达方向，包括以下步骤：

5.1)特征值分解

并以特征值大小对特征向量进行降序排列。将前s个特征向量写入矩阵U_S中。

5.2)基于球谐函数递归关系，建立与矩阵U_S相关的线性方程组，并采用最小二乘方法对与矩阵U_S相关的线性方程组进行求解，得到含有声源波达方向的转换矩阵。

线性方程组如下所示：

其中，

为按上标(u,v)从矩阵U_S中选取一部分行组合而成的矩阵，ψ_xy+、ψ_xy-和ψ_z为含有声源波达方向的转换矩阵，

和

为系数矩阵。

其中，系数矩阵

和系数矩阵

分别如下所示：

式中，w、v为系数矩阵

和系数矩阵

中的元素。

5.3)采用广义特征值分解法对含有声源波达方向的矩阵进行特征值分解，确定声源波达方向。

详细步骤可参考文献“Parametric direction-of-arrival estimation withthree recurrence relations of spherical harmonics”,B.Jo,J.W.Choi.Journal ofAcoustical Society of America,volume 145,issue No.1,pages 480-488,January2019(基于三个球谐函数递归关系的参数化波达方向估计,B.Jo,J.W.Choi.Journal ofAcoustical Society of America,第145卷,第1期,第1030-1042页,2019年1月)”。

实施例2：

基于原子范数的球面阵列声源波达方向估计方法的验证试验：

进行声源识别仿真模拟。5个声源被假设，DOA依次为(30°,90°)、(150°,270°)、(120°,80°)、(60°,290°)和(90°,180°)，强度

依次为100dB、97.5dB、95dB、92.5dB和90dB(参考2×10-5Pa)。

公司的半径为0.0975m的包含36个传声器的刚性球阵列被采用，对应N_D为5。正向声场模拟时，由于无法计算无穷多项，N₀取为20。用基于ADMM的求解算法求解式(22)所示ANM的等价半正定规划时，相关参数设置如下：ρ取1，τ按文献“Off-the-grid line spectrum denoising and estimation with multiplemeasurement vectors”,Y.Li,Y.Chi.IEEE Transactions on Signal Processing,volume64,issue No.5,pages 1157-1269,March 2016(基于多重测量矢量的一维无网格线谱去噪和估计,Y.Li,Y.Chi.IEEE Transactions on Signal Processing,第64卷,第5期,第1257-1269页,2016年3月)中给出的原则取值，最大迭代次数取1000，当连续两次迭代的u间的相对变化量，即||u^γ-u^γ-1||₂/||u^γ-1||₂，小于10-3或最大迭代次数被完成时，迭代终止。所有仿真均在3.70GHz Intel(R)Core(TM)i7-8700K的CPU上用MATLAB R2014a执行。

假设声源互不相干，快拍总数为30，无噪声干扰，改变声源辐射声波的频率进行仿真，每个频率下进行多次蒙特卡罗计算，每次计算随机生成S。图2为1000Hz、3000Hz和5000Hz时具有代表性的单次蒙特卡罗计算的声源成像图，每幅子图中，“*”和“○”分别指示重构的和真实的声源分布，重构结果和真实结果分别参考各自的最大值进行dB缩放，上方标注为重构的参考2×10^-5Pa的最大值，后续成像图亦如此。由图2可见，三个频率下，球谐函数阶截断后的最高阶N依次为3、7和10。第1列(图2(a)、(d)和(g))对应球面ESPRIT，仅1000Hz时(图2(a))，N≤N_D，声源DOA被准确估计。第2和3列(图2(b)、(c)、(e)、(f)、(h)和(i))对应ANM+球面ESPRIT，不论ANM由基于IPM的SDPT3求解器求解(图2(b)、(e)和(h))，还是由基于ADMM的求解算法求解(图2(c)、(f)和(i))，各频率下，声源DOA均被准确估计。定义第d次计算的声源DOA估计均方根误差为

其中，

为i号声源DOA的真实值Ω_Si＝(θ_Si,φ_Si)与估计值

间的角距离。该公式适用于估计的声源数目大于等于真实的声源数目的情形，此时，认为估计的前I个较强声源对应真实声源。当估计的声源数目小于真实的声源数目时，声源丢失，认为σ_d很大。定义D次计算的声源DOA估计均方根误差为

图3(a)给出了100次蒙特卡罗计算时σ随频率的变化曲线，各频率对应的N也被标出。显然，球面ESPRIT在N≤N_D的低频段具有较低误差，在N＞N_D的高频段误差显著增大，而不论ANM由基于IPM的SDPT3求解器求解，还是由基于ADMM的求解算法求解，各频率下，ANM+球面ESPRIT的误差均极低。图2和图3(a)均证明：运用ANM能完美克服球面ESPRIT高频失效的缺陷。图3(b)对比了两种ANM求解方法的耗时。显然，基于ADMM的求解算法比基于IPM的SDPT3求解器更高效；随频率增高，N变大，ANM问题涉及矩阵的维度变大，两种方法的耗时均呈增长趋势，在更高频率(更大N)，例如5200Hz(N＝11)，基于IPM的SDPT3求解器受仅适用于小维度矩阵问题的固有特性的限制将无法工作，而基于ADMM的求解算法仍然能够应用。

假设频率为1500Hz，快拍总数为30，无噪声干扰，改变声源间的相干性进行仿真。图4为具有代表性的单次蒙特卡罗计算的声源成像图。第1列(图4(a)、(d)和(g))对应球面ESPRIT，仅声源互不相干时(图4(a))，各声源的DOA被准确估计，声源部分相干(图4(d))或完全相干(图4(g))时，球面ESPRIT均失效。第2和3列(图4(b)、(c)、(e)、(f)、(h)和(i))对应ANM+球面ESPRIT，不论声源互不相干(图4(b)和(c))，还是部分相干(图4(e)和(f))，或是完全相干(图4(h)和(i))，不论ANM由基于IPM的SDPT3求解器求解(图4(b)、(e)和(h))，还是由基于ADMM的求解算法求解(图4(c)、(f)和(i))，ANM+球面ESPRIT均能准确估计各声源的DOA。这证明：运用ANM能完美克服球面ESPRIT对相干声源失效的缺陷。获得图4(b)、(e)和(h)时，由基于IPM的SDPT3求解器求解ANM的耗时分别约为32s、32s和35s，获得图4(c)、(f)和(i)时，由基于ADMM的求解算法求解ANM的耗时分别约为6s、8s和9s，再次证明了后者相比于前者的高效性。

假设声源互不相干，频率为1500Hz，改变快拍总数和SNR进行仿真，每对快拍总数和SNR下进行多次蒙特卡罗计算，每次计算随机生成S和N。声源DOA估计均方根误差的互补累积分布函数为F(T)＝P(σ_d＞T)＝|{σ_d|σ_d＞T}|/D，其中，T为自变量，P(·)表示括号内事件发生的概率。T取小值时，σ_d≤T意味着声源DOA估计成功，F(T)表示声源DOA估计失败的概率，F(T)≈0意味着声源DOA被高概率成功估计。图5给出了100次蒙特卡罗计算时不同快拍总数不同SNR下声源DOA估计均方根误差的互补累积分布函数。图5(a)-(c)对应T＝2°，相比图5(a)所示的球面ESPRIT，图5(b)和(c)所示的ANM+球面ESPRIT能高概率成功估计声源DOA的区域明显更大。当SNR足够高时，ANM+球面ESPRIT在少快拍甚至单快拍下可高概率成功估计各声源的DOA，而球面ESPRIT不能；在低SNR下，前者声源DOA估计失败的概率整体也低于后者。图5(d)-(f)对应T＝1°，即认为声源DOA估计成功的标准更严格。对比图5(a)和(d)显见，球面ESPRIT的F(T)≈0的区域从有变无，对比图5(b)和(e)、图5(c)和(f)显见，ANM+球面ESPRIT的F(T)≈0的区域仅轻微变小，说明ANM+球面ESPRIT的声源DOA估计精度更高。仿真结果表明ANM+球面ESPRIT能够解决传统球面ESPRIT对少数据快拍或低信噪比工况失效的缺陷，显著提高其声源识别性能。

基于半消声室内的扬声器声源识别试验验证仿真结论的正确性并检验提出的方法在半消声测试环境中的有效性，然后基于普通房间内的扬声器声源识别试验检验提出的方法在普通测试环境中的有效性。试验采用

公司的半径为0.0975m的包含36个4958型传声器和12个uEye UI-122xLE型摄像头的刚性球阵列进行测量。各传声器测量的声信号经PULSE 3560D型数据采集系统同时采集并传输到BKCONNECT中进行频谱分析。测量时长为5s，采样频率为16384Hz，信号添加汉宁窗，每个快拍长1s、对应频率分辨率为1Hz，重叠率为90％，共40个快拍被获得。基于ADMM的求解算法中相关参数的设置与仿真中一致。

图6(a)为试验布局，五个扬声器围绕球阵列布置，图6(b)为12个摄像头拍摄的照片组合而成的三维空间展开图，与θ和φ的对应关系被标出，两幅图中均用圆圈标记扬声器位置并对应编号。图7为声源成像图，第1列(图7(a)、(d)、(g)、(j)和(m))对应球面ESPRIT，第2列(图7(b)、(e)、(h)、(k)和(n))对应ANM+球面ESPRIT且ANM由基于IPM的SDPT3求解器求解，第3列(图7(c)、(f)、(i)、(l)和(o))对应ANM+球面ESPRIT且ANM由基于ADMM的求解算法求解。图7(a)-(l)对应五个扬声器均由稳态白噪声信号激励的工况，声源互不相干。计算图7(a)-(c)时，30个快拍被采用，频率为1500Hz(4＝N≤N_D＝5)。如图7(a)示，球面ESPRIT估计的DOA中有四个接近真实声源(1、3、4和5号声源)的DOA，2号声源的DOA未被估计出，这可能是受地面反射干扰的缘故。如图7(b)和(c)示，ANM+球面ESPRIT准确估计了每个声源的DOA。计算图7(d)-(f)时，30个快拍被采用，频率为3000Hz(7＝N＞N_D＝5)。图7(d)证明球面ESPRIT在N＞N_D的高频失效，图7(e)和(f)证明ANM+球面ESPRIT在高频时仍能准确估计声源DOA。计算图7(g)-(i)时，仅5个快拍被采用，计算图7(j)-(l)时，仅单个快拍被采用，频率均为1500Hz。图7(g)和(j)证明球面ESPRIT在少数据快拍时失效，图7(h)、(i)、(k)和(l)证明ANM+球面ESPRIT在少数据快拍时仍能准确估计声源DOA。图7(m)-(o)对应五个扬声器均由同一纯音信号激励的工况，声源彼此相干，计算时采用30个快拍和1500Hz频率。图7(m)证明球面ESPRIT对相干声源失效，图7(n)和(o)证明ANM+球面ESPRIT仍能准确估计相干声源的DOA。获得图7(b)、(e)、(h)、(k)和(n)时，求解ANM的耗时分别约为34s、392s、5s、3s和46s，获得图7(c)、(f)、(i)、(l)和(o)时，求解ANM的耗时分别约为2s、22s、1s、1s和6s，证明基于ADMM的求解算法比基于IPM的SDPT3求解器更高效。这些试验结果与仿真结果一致，证明仿真结论正确，同时表明提出的ANM+球面ESPRIT方法在半消声测试环境中有效。

相比半消声室，普通房间存在严重混响干扰，且频率越低越严重。图8给出了普通房间内的试验布局和三维空间展开图。五个扬声器均由稳态白噪声信号激励，计算时30个快拍被采用。图9为声源成像图，第一行(图9(a)-(c))对应1000Hz，第二行(图9(d)-(f))对应2000Hz，第三行(图9(g)-(i))对应3000Hz，第四行(图9(j)-(l))对应4000Hz，第五行(图9(m)-(o))对应5000Hz，第1列(图9(a)、(d)、(g)、(j)和(m))对应球面ESPRIT，第2列(图9(b)、(e)、(h)、(k)和(n))对应ANM+球面ESPRIT且ANM由基于IPM的SDPT3求解器求解，第3列(图9(c)、(f)、(i)、(l)和(o))对应ANM+球面ESPRIT且ANM由基于ADMM的求解算法求解。显见，五个频率下，球面ESPRIT均不能有效估计声源DOA，除1000Hz外的其它频率下，不论ANM由哪种方法求解，ANM+球面ESPRIT均能准确估计每个声源的DOA，尽管估计的声源数目多于真实的声源数目，估计结果中多出一些弱的虚假源。1000Hz和2000Hz时球面ESPRIT失效及1000Hz时ANM+球面ESPRIT的低声源DOA估计准确度主要归咎于低频时比较严重的房间混响干扰，3000Hz、4000Hz和5000Hz时球面ESPRIT失效主要归咎于其对N＞N_D的高频失效的缺陷。这些现象证明即使在普通测试环境中，提出的ANM+球面ESPRIT方法仍有效。

Claims (4)

1.基于原子范数的球面阵列声源波达方向估计方法，其特征在于，包括以下步骤：

1)搭建由所述Q个传声器(2)构成的球形传声器(2)阵列；球形传声器(2)阵列中心记为坐标原点；其中，第q个传声器(2)位置记为(a,Ω_Mq)；a为阵列半径，q＝1,2,…,Q；Ω＝(θ,φ)表示球形传声器(2)阵列所在三维空间内的任意方向；θ∈[0°,180°]为仰角，φ∈[0°,360°)为方位角；

2)声源(1)向球形传声器(2)阵列辐射声波。

3)建立声源波达方向测量模型，并构建传声器(2)测量得到的声压信号矩阵P^★；

4)建立协方差矩阵

5)利用球面ESPRIT算法对协方差矩阵

进行解算，确定声源波达方向。

2.根据权利要求1或2所述的基于原子范数的球面阵列声源波达方向估计方法，其特征在于：建立声源波达方向测量模型，包括以下步骤：

1)计算第i个声源到第q个传声器(2)的传递函数t((ka,Ω_Mq)|Ω_Si)，即：

式中，n和m分别为球谐函数的阶和次；b_n(ka)为模态强度；

为Ω方向的球谐函数；Ω_Si表示第i个声源波达方向；i＝1,2,…,I；I为声源总数；k为声源辐射声波的波数；上标*表示共轭；

模态强度b_n(ka)如下所示：

式中，j_n(ka)为n阶第一类球贝塞尔函数，

为n阶第二类球汉克尔函数；j'_n(ka)和

分别为n阶第一类球贝塞尔函数j_n(ka)和n阶第二类球汉克尔函数

的一阶导数；开口球表示传声器(2)布置在开口球体表面；刚性球表示传声器(2)布置在刚性球体表面；

2)计算Ω方向的球谐函数

即：

式中，

为连带勒让德函数；β_n,m,κ为连带勒让德函数系数；球谐函数第n阶m次项对应的球谐系数

3)建立Q个传声器(2)所在方向的向量

和对应的球谐函数向量

建立I个声源所在方向的向量

和对应的球谐函数向量

4)建立每个声源到所有传声器(2)的传递函数矩阵，记为

上标H表示转置共轭；

其中，传声器所在方向对应的球谐函数矩阵

如下所示：

式中，N₀＝∞表示球谐函数最高阶；

声源所在方向对应的球谐函数矩阵

如下所示：

阵列模态强度矩阵

如下所示：

5)建立各传声器(2)测量得到的声压信号矩阵

即：

式中，

为噪声矩阵；信噪比SNR＝20lg(||P^★-N||_F/||N||_F)；声源强度矩阵

L为快拍总数；|| ||_F表示Frobenius范数；上标★表示测量量；

6)对球谐函数矩阵

球谐函数矩阵

和模态强度矩阵

的阶n进行截断，即令球谐函数矩阵

球谐函数矩阵

和模态强度矩阵

的最高阶

表示将数值向正无穷方向圆整到次近的整数；

基于最高阶N，更新声压信号矩阵P^★如下：

3.根据权利要求1或2所述的基于原子范数的球面阵列声源波达方向估计方法，其特征在于，建立协方差矩阵

包括以下步骤：

1)建立原子范数最小化模型，包括以下步骤：

1.1)建立连带勒让德函数表达式，即：

式中，x为函数输入；

其中，连带勒让德函数项(x²-1)ⁿ如下所示：

式中，多项式系数

1.2)将公式(11)代入公式(9)，得到连带勒让德函数次方m≥0时连带勒让德函数表达式，即：

1.3)确定仰角θ的正弦函数sinθ＝(e^jθ-e^-jθ)/(2j)和余弦函数cosθ＝(e^jθ+e^-jθ)/2；

1.4)基于步骤1.2)和步骤1.3)，更新连带勒让德函数表达式如下：

1.5)连带勒让德函数项(e^jθ-e^-jθ)^m、连带勒让德函数项(e^jθ+e^-jθ-2)^o-m分别如下所示：

1.6)基于公式(14)和公式(15)，更新连带勒让德函数表达式如下：

1.7)给定函数阶数n和级数m，令索引o从m增至n、索引u从0增至m、索引v从0增至o-m、索引w从0增至o-m-v；

在每组(o,u,v,w)下，根据公式(16)确定连带勒让德函数的三角多项式展开中的指数索引κ＝2u+v-m-w和连带勒让德函数的三角多项式中的系数；计算完所有组(o,u,v,w)后，同一κ值对应的系数相加即为β_n,m,κ；

1.8)联合公式(9)和m≥0时的连带勒让德函数系数β_n,m,κ，确定m＜0时的连带勒让德函数系数β_n,m,κ；

1.9)确定连带勒让德函数系数β_n,m,κ后，构建Ω方向的球谐函数

即：

1.10)构建矩阵

和矩阵

式中，矩阵D中元素

d_θ(θ_Si)、d_φ(φ_Si)为构成矩阵D的基础向量；

记向量

的第(N+κ)(2N+1)+N-m+1个元素为A_n,mβ_n,m,κ，κ＝-n,…,0,…,n，并令其它元素为0，则球谐函数矩阵的转置共轭

各传声器(2)测量得到的声压信号矩阵P^★≈Y_MNB_NGDS+N；

1.11)建立输入矩阵X，即：

式中，

为S的第i行；

||ψ_i||₂＝1；

公式(18)原子范数如下所示：

式中，“inf”表示下确界；

为原子集合；

原子集合

如下所示：

1.12)建立原子范数最小化模型，即：

其中，ε为噪声控制参数；

是连续域内声源稀疏性的一个度量；

表示原子范数最小化问题的最优解；

2)建立等价的半正定规划模型，包括以下步骤：

2.1)将式(21)转化为如下半正定规划模型：

式中，u和E为辅助量；N_u为辅助量u中元素的数目；T_b(·)为二重Toeplitz算子；对任意给定向量

是(2N,2N)的半空间，N_u＝8N²+4N+1；

2.2)利用二重Toeplitz算子T_b(u)将u映射为(2N+1)×(2N+1)维的块Toeplitz型Hermitian矩阵，即：

式中，每个矩阵分块T_κ均为(2N+1)×(2N+1)维的Toeplitz矩阵：0≤κ≤2N；

矩阵分块T_κ如下所示：

23)建立矩阵

的Vandermonde分解式，从而令公式(21)和公式(22)等价；

Vandermonde分解式如下所示：

式中，矩阵V＝[d(Ω_S1),d(Ω_S2),…,d(Ω_Sr)]；对角矩阵Σ＝diag([σ₁,σ₂,…,σ_r])；

i＝1,2,…,r；r为矩阵

的秩；r≤2N+1是公式(25)存在的充分条件；矩阵

视为是一组声源中，各个声源独自引起的信号的协方差矩阵的和，不包括不同声源引起的信号间的协方差成分；

3)使用交替方向乘子方法对半正定规划模型进行解算，包括以下步骤：

3.1)使用交替方向乘子方法更新半正定规划模型，得到：

式中，Z为辅助矩阵，τ为规则化参数；

3.2)建立公式(26)的增广拉格朗日函数表达式，即：

式中，

为Hermitian拉格朗日乘子矩阵；ρ＞0为惩罚参数；“<·,·>”表示内积；

3.3)利用交替方向乘子方法求解公式(26)，初始化辅助矩阵Z⁰＝Λ⁰＝0，γ+1次迭代时的变量更新为：

3.4)对Hermitian拉格朗日乘子矩阵和辅助矩阵进行划分，得到：

3.5)基于步骤3.3)和步骤3.4)，更新公式(28)如下：

式中，I₁和I₂分别为L×L维和(2N+1)²×(2N+1)²维单位矩阵；

为T_b(·)的伴随算子。对任意给定矩阵

M为对角矩阵；矩阵

为第κ(m)个对角线的元素均为1其它元素均为0的基本Toeplitz矩阵；

3.6)基于步骤3.3)，更新公式(29)如下：

公式(36)表示将Hermitian矩阵

投影到半正定锥上，即对该Hermitian矩阵进行特征值分解，并令所有负特征值为0；

4)基于求得的协方差矩阵和

建立协方差矩阵

4.根据权利要求1所述的基于原子范数的球面阵列声源波达方向估计方法，其特征在于，确定声源波达方向，包括以下步骤：

1)特征值分解

并以特征值大小对特征向量进行降序排列；将前s个特征向量写入矩阵U_S中；

2)基于球谐函数递归关系，建立与矩阵U_S相关的线性方程组，并采用最小二乘方法对与矩阵U_S相关的线性方程组进行求解，得到含有声源波达方向的转换矩阵；

线性方程组如下所示：

其中，

为按上标(u,v)从矩阵U_S中选取一部分行组合而成的矩阵，ψ_xy+、ψ_xy-和ψ_z为含有声源波达方向的转换矩阵，

和

为系数矩阵；

其中，系数矩阵

和系数矩阵

分别如下所示：

3)采用广义特征值分解法对含有声源波达方向的矩阵进行特征值分解，确定声源波达方向。

Images