基于原子范数的球面阵列声源波达方向估计方法
阅读说明：本技术 基于原子范数的球面阵列声源波达方向估计方法 (Spherical array sound source direction of arrival estimation method based on atomic norm ) 是由 褚志刚 杨洋 于 2020-06-16 设计创作，主要内容包括：本发明公开了基于原子范数的球面阵列声源波达方向估计方法,步骤为：1)搭建由Q个传声器(2)构成的球形传声器(2)阵列；2)声源(1)向球形传声器(2)阵列辐射声波；3)建立声源波达方向测量模型,并构建传声器(2)测量得到的声压信号矩阵P<Sup>★</Sup>；4)建立协方差矩阵<Image he="80" wi="260" file="DDA0002540654490000011.GIF" imgContent="drawing" imgFormat="GIF" orientation="portrait" inline="no"></Image>5)利用球面ESPRIT算法对协方差矩阵<Image he="81" wi="230" file="DDA0002540654490000012.GIF" imgContent="drawing" imgFormat="GIF" orientation="portrait" inline="no"></Image>进行解算,确定声源波达方向。本发明能克服球面ESPRIT在高频、相干声源或少数据快拍工况下失效的缺陷,并显著提高低SNR工况下的声源DOA估计精度,即使在存在混响的普通测试环境中,本发明仍然有效。(The invention discloses a spherical array sound source direction of arrival estimation method based on atomic norm, which comprises the following steps: 1) building a spherical microphone (2) array consisting of Q microphones (2); 2) a sound source (1) radiates sound waves to a spherical microphone (2) array; 3) establishing a sound source direction of arrival measurement model and constructing a sound pressure signal matrix P obtained by measuring a microphone (2) ★ (ii) a 4) Establishing a covariance matrix 5) Covariance matrix pair using spherical ESPRIT algorithm And resolving is carried out, and the arrival direction of the sound source is determined. The method can overcome the defect that the spherical ESPRIT fails under the working conditions of high frequency, coherent sound source or few data snapshots, obviously improve the DOA estimation precision of the sound source under the working condition of low SNR, and is still effective even in the common test environment with reverberation.)
技术领域
本发明涉及声源识别领域，具体是基于原子范数的球面阵列声源波达方向估计方法。
背景技术
声源波达方向(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分别为球谐函数的阶和次。bn(ka)为模态强度。
模态强度bn(ka)如下所示：
式中，jn(ka)为n阶第一类球贝塞尔函数，为n阶第二类球汉克尔函数。j'n(ka)和分别为n阶第一类球贝塞尔函数jn(ka)和n阶第二类球汉克尔函数的一阶导数。开口球表示传声器布置在开口球体表面。刚性球表示传声器布置在刚性球体表面。
3.2)计算Ω方向的球谐函数即：
式中，为连带勒让德函数。βn,m,κ为连带勒让德函数系数。球谐函数第n阶m次项对应的球谐系数
3.3)建立Q个传声器所在方向的向量和对应的球谐函数向量建立I个声源所在方向的向量
3.4)建立每个声源到所有传声器的传递函数矩阵，记为上标H表示转置共轭。
其中，传声器所在方向对应的球谐函数矩阵
式中，N0＝∞表示球谐函数最高阶。
声源所在方向对应的球谐函数矩阵如下所示：
阵列模态强度矩阵
3.5)建立各传声器测量得到的声压信号矩阵即：
式中，为噪声矩阵。信噪比SNR＝20lg(||P★-N||F/||N||F)。声源强度矩阵L为快拍总数。||||F表示Frobenius范数。上标★表示测量量。
3.6)对球谐函数矩阵
基于最高阶N，更新声压信号矩阵P★如下：
4)建立协方差矩阵
建立协方差矩阵包括以下步骤：
4.1)建立原子范数最小化模型，包括以下步骤：
4.1.1)建立连带勒让德函数表达式，即：
式中，x为函数输入。
其中，连带勒让德函数项(x2-1)n如下所示：
式中，多项式系数
4.1.2)将公式(11)代入公式(9)，得到连带勒让德函数次方m≥0时连带勒让德函数表达式，即：
4.1.3)确定仰角θ的正弦函数sinθ＝(ejθ-e-jθ)/(2j)和余弦函数cosθ＝(ejθ+e-jθ)/2。
4.1.4)基于步骤4.1.2)和步骤4.1.3)，更新连带勒让德函数表达式如下：
4.1.5)连带勒让德函数项(ejθ-e-jθ)m、连带勒让德函数项(ejθ+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个元素为An,mβn,m,κ，κ＝-n,…,0,…,n，并令其它元素为0，则球谐函数矩阵的转置共轭
4.1.11)建立输入矩阵X，即：
式中，
公式(18)原子范数如下所示：
式中，“inf”表示下确界。为原子集合。
原子集合如下所示：
4.1.12)建立原子范数最小化模型，即：
其中，ε为噪声控制参数。
4.2)建立等价的半正定规划模型，包括以下步骤：
4.2.1)将式(21)转化为如下半正定规划模型：
式中，u和E为辅助量。Nu为辅助量u中元素的数目。Tb(·)为二重Toeplitz算子。对任意给定向量 是(2N,2N)的半空间，Nu＝8N2+4N+1。
4.2.2)利用二重Toeplitz算子Tb(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([σ1,σ2,…,σr])；
4.3)使用交替方向乘子方法对半正定规划模型进行解算，包括以下步骤：
4.3.1)使用交替方向乘子方法更新半正定规划模型，得到：
式中，Z为辅助矩阵，τ为规则化参数。
4.3.2)建立公式(26)的增广拉格朗日函数表达式，即：
式中，
4.3.3)利用交替方向乘子方法求解公式(26)，初始化辅助矩阵Z0＝Λ0＝0，γ+1次迭代时的变量更新为：
4.3.4)对Hermitian拉格朗日乘子矩阵和辅助矩阵进行划分，得到：
4.3.5)基于步骤3.3)和步骤3.4)，更新公式(28)如下：
式中，I1和I2分别为L×L维和(2N+1)2×(2N+1)2维单位矩阵。
4.3.6)基于步骤3.3)，更新公式(29)如下：
公式(36)表示将Hermitian矩阵
4.4)基于求得的协方差矩阵和建立协方差矩阵
5)利用球面ESPRIT算法对协方差矩阵
5.1)特征值分解
5.2)基于球谐函数递归关系，建立与矩阵US相关的线性方程组，并采用最小二乘方法对与矩阵US相关的线性方程组进行求解，得到含有声源波达方向的转换矩阵。线性方程组如下：
其中，
其中，系数矩阵
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中“●”和“○”分别表示传声器和声源。
表示实数集，
2)声源1向球形传声器阵列辐射声波。
3)建立声源波达方向测量模型，并构建传声器测量得到的声压信号矩阵P★。
建立声源波达方向测量模型，包括以下步骤：
3.1)计算第i个声源到第q个传声器的传递函数t((ka,ΩMq)|ΩSi)，即：
式中，n和m分别为球谐函数的阶和次。bn(ka)为模态强度。
模态强度bn(ka)如下所示：
式中，jn(ka)为n阶第一类球贝塞尔函数，为n阶第二类球汉克尔函数。j'n(ka)和
3.2)计算Ω方向的球谐函数即：
式中，
3.3)建立Q个传声器所在方向的向量
3.4)建立每个声源到所有传声器的传递函数矩阵，记为上标H表示转置共轭。
其中，传声器所在方向对应的球谐函数矩阵如下所示：
式中，N0＝∞表示球谐函数最高阶。
声源所在方向对应的球谐函数矩阵如下所示：
阵列模态强度矩阵如下所示：
3.5)建立各传声器测量得到的声压信号矩阵
式中，
3.6)各传声器测量的声压信号包含无穷多阶球谐函数。用符号
基于最高阶N，更新声压信号矩阵P★如下：
4)建立协方差矩阵
建立协方差矩阵
4.1)建立原子范数最小化模型，包括以下步骤：
4.1.1)建立连带勒让德函数表达式，即：
式中，x为函数输入。
其中，连带勒让德函数项(x2-1)n如下所示：
式中，多项式系数
4.1.2)将公式(11)代入公式(9)，得到连带勒让德函数次方m≥0时连带勒让德函数表达式，即：
4.1.3)确定仰角θ的正弦函数sinθ＝(ejθ-e-jθ)/(2j)和余弦函数cosθ＝(ejθ+e-jθ)/2。
4.1.4)基于步骤4.1.2)和步骤4.1.3)，更新连带勒让德函数表达式如下：
4.1.5)连带勒让德函数项(ejθ-e-jθ)m、连带勒让德函数项(ejθ+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个元素为An,mβn,m,κ，κ＝-n,…,0,…,n，并令其它元素为0，则球谐函数矩阵的转置共轭
4.1.11)建立输入矩阵X，即：
式中，为S的第i行。元素向量
公式(18)原子范数如下所示：
式中，“inf”表示下确界。
原子集合
4.1.12)建立原子范数最小化模型，即：
其中，ε为噪声控制参数。是连续域内声源稀疏性的一个度量。最小化
4.2)建立等价的半正定规划模型，包括以下步骤：
4.2.1)将式(21)转化为如下半正定规划模型：
式中，u和E为辅助量。Nu为辅助量u中元素的数目。Tb(·)为二重Toeplitz算子。对任意给定向量
4.2.2)利用二重Toeplitz算子Tb(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([σ1,σ2,…,σr])；
是一组源中，各个源独自引起的信号的协方差矩阵的和，不包括不同源引起的信号间的协方差成分，即去除了源间的相关性，可看作是一组不相干源引起的信号的协方差矩阵，且该特性与采用的数据快拍数目无关。由于施加约束||P★-YMNBNGX||F≤ε，ANM具有噪声滤除功能。此外，ANM未涉及传声器采样的球谐函数的正交性。因此，本发明用
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)，初始化辅助矩阵Z0＝Λ0＝0，γ+1次迭代时的变量更新为：
4.3.4)对Hermitian拉格朗日乘子矩阵和辅助矩阵进行划分，得到：
4.3.5)基于步骤3.3)和步骤3.4)，更新公式(28)如下：
式中，I1和I2分别为L×L维和(2N+1)2×(2N+1)2维单位矩阵。
4.3.6)基于步骤3.3)，更新公式(29)如下：
公式(36)表示将Hermitian矩阵投影到半正定锥上，即对该Hermitian矩阵进行特征值分解，并令所有负特征值为0。
4.4)基于求得的协方差矩阵和建立协方差矩阵
5)利用球面ESPRIT算法对协方差矩阵
5.1)特征值分解并以特征值大小对特征向量进行降序排列。将前s个特征向量写入矩阵US中。
5.2)基于球谐函数递归关系，建立与矩阵US相关的线性方程组，并采用最小二乘方法对与矩阵US相关的线性方程组进行求解，得到含有声源波达方向的转换矩阵。
线性方程组如下所示：
其中，
其中，系数矩阵
式中，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个传声器的刚性球阵列被采用，对应ND为5。正向声场模拟时，由于无法计算无穷多项，N0取为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||2/||uγ-1||2，小于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≤ND，声源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估计均方根误差为
假设频率为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对少数据快拍或低信噪比工况失效的缺陷，显著提高其声源识别性能。
基于半消声室内的扬声器声源识别试验验证仿真结论的正确性并检验提出的方法在半消声测试环境中的有效性，然后基于普通房间内的扬声器声源识别试验检验提出的方法在普通测试环境中的有效性。试验采用
图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≤ND＝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＞ND＝5)。图7(d)证明球面ESPRIT在N＞ND的高频失效，图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＞ND的高频失效的缺陷。这些现象证明即使在普通测试环境中，提出的ANM+球面ESPRIT方法仍有效。