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

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

Info

Publication number
CN111812581B
CN111812581B CN202010545752.1A CN202010545752A CN111812581B CN 111812581 B CN111812581 B CN 111812581B CN 202010545752 A CN202010545752 A CN 202010545752A CN 111812581 B CN111812581 B CN 111812581B
Authority
CN
China
Prior art keywords
matrix
sound source
spherical
function
sound
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.)
Active
Application number
CN202010545752.1A
Other languages
English (en)
Other versions
CN111812581A (zh
Inventor
褚志刚
杨洋
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Chongqing University
Original Assignee
Chongqing University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Chongqing University filed Critical Chongqing University
Priority to CN202010545752.1A priority Critical patent/CN111812581B/zh
Publication of CN111812581A publication Critical patent/CN111812581A/zh
Application granted granted Critical
Publication of CN111812581B publication Critical patent/CN111812581B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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
    • G01S3/80—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 using ultrasonic, sonic or infrasonic waves
    • G01S3/802—Systems for determining direction or deviation from predetermined direction

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

本发明公开了基于原子范数的球面阵列声源波达方向估计方法，步骤为：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分别为球谐函数的阶和次。bn(ka)为模态强度。为Ω方向的球谐函数。ΩSi表示第i个声源波达方向。i＝1,2,…,I。I为声源总数。k为声源辐射声波的波数。上标*表示共轭。
模态强度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进行截断，即令球谐函数矩阵球谐函数矩阵和模态强度矩阵的最高阶 表示将数值向正无穷方向圆整到次近的整数。
基于最高阶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，则球谐函数矩阵的转置共轭各传声器测量得到的声压信号矩阵P★≈YMNBNGDS+N。
4.1.11)建立输入矩阵X，即：
式中，为S的第i行。||ψi||2＝1。
公式(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])；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)，初始化辅助矩阵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维单位矩阵。为Tb(·)的伴随算子。对任意给定矩阵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个特征向量写入矩阵US中。
5.2)基于球谐函数递归关系，建立与矩阵US相关的线性方程组，并采用最小二乘方法对与矩阵US相关的线性方程组进行求解，得到含有声源波达方向的转换矩阵。线性方程组如下：
其中，为按上标(u,v)从矩阵US中选取一部分行组合而成的矩阵，ψ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范数被定义为任意向量的l2范数被定义为tr(A)表示对矩阵A求迹，diag(a)表示形成以a中元素为对角元素的对角矩阵。
2)声源1向球形传声器阵列辐射声波。
3)建立声源波达方向测量模型，并构建传声器测量得到的声压信号矩阵P★。
建立声源波达方向测量模型，包括以下步骤：
3.1)计算第i个声源到第q个传声器的传递函数t((ka,ΩMq)|ΩSi)，即：
式中，n和m分别为球谐函数的阶和次。bn(ka)为模态强度。为Ω方向的球谐函数。ΩSi表示第i个声源波达方向。i＝1,2,…,I。I为声源总数。k为声源辐射声波的波数。上标*表示共轭。
模态强度bn(ka)如下所示：
式中，jn(ka)为n阶第一类球贝塞尔函数，为n阶第二类球汉克尔函数。j'n(ka)和分别为n阶第一类球贝塞尔函数jn(ka)和n阶第二类球汉克尔函数的一阶导数。开口球表示传声器布置在开口球体表面。刚性球表示传声器布置在刚性球体表面。jn为符号函数。
3.2)计算Ω方向的球谐函数即：
式中，为连带勒让德函数。βn,m,κ为连带勒让德函数系数。球谐函数第n阶m次项对应的球谐系数
3.3)建立Q个传声器所在方向的向量和对应的球谐函数向量建立I个声源所在方向的向量和对应的球谐函数向量ΩMQ为第Q个传声器所在方向的向量。表示I维复数集。
3.4)建立每个声源到所有传声器的传递函数矩阵，记为上标H表示转置共轭。
其中，传声器所在方向对应的球谐函数矩阵如下所示：
式中，N0＝∞表示球谐函数最高阶。
声源所在方向对应的球谐函数矩阵如下所示：
阵列模态强度矩阵如下所示：
3.5)建立各传声器测量得到的声压信号矩阵即：
式中，为噪声矩阵。信噪比SNR＝20lg(||P★-N||F/||N||F)。声源强度矩阵L为快拍总数。||||F表示Frobenius范数。上标★表示测量量。
3.6)各传声器测量的声压信号包含无穷多阶球谐函数。用符号表示将数值向正无穷方向圆整到次近的整数，时，bn(ka)的幅值很小，相应阶次的贡献可忽略。基于该事实，可对球谐函数的阶进行截断。对球谐函数矩阵球谐函数矩阵和阵列模态强度矩阵的阶n进行截断，即令球谐函数矩阵球谐函数矩阵和阵列模态强度矩阵的最高阶 表示将数值向正无穷方向圆整到次近的整数。
基于最高阶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，则球谐函数矩阵的转置共轭各传声器测量得到的声压信号矩阵P★≈YMNBNGDS+N。
4.1.11)建立输入矩阵X，即：
式中，为S的第i行。元素向量||ψi||2＝1。
公式(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])；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★-YMNBNGX||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)，初始化辅助矩阵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维单位矩阵。为Tb(·)的伴随算子。对任意给定矩阵
M为对角矩阵。矩阵 为第κ(m)个对角线的元素均为1其它元素均为0的基本Toeplitz矩阵。
4.3.6)基于步骤3.3)，更新公式(29)如下：
公式(36)表示将Hermitian矩阵投影到半正定锥上，即对该Hermitian矩阵进行特征值分解，并令所有负特征值为0。
4.4)基于求得的协方差矩阵和建立协方差矩阵
5)利用球面ESPRIT算法对协方差矩阵进行解算，确定声源波达方向，包括以下步骤：
5.1)特征值分解并以特征值大小对特征向量进行降序排列。将前s个特征向量写入矩阵US中。
5.2)基于球谐函数递归关系，建立与矩阵US相关的线性方程组，并采用最小二乘方法对与矩阵US相关的线性方程组进行求解，得到含有声源波达方向的转换矩阵。
线性方程组如下所示：
其中，为按上标(u,v)从矩阵US中选取一部分行组合而成的矩阵，ψ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个传声器的刚性球阵列被采用，对应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估计均方根误差为其中，为i号声源DOA的真实值ΩSi＝(θSi,φSi)与估计值间的角距离。该公式适用于估计的声源数目大于等于真实的声源数目的情形，此时，认为估计的前I个较强声源对应真实声源。当估计的声源数目小于真实的声源数目时，声源丢失，认为σd很大。定义D次计算的声源DOA估计均方根误差为图3(a)给出了100次蒙特卡罗计算时σ随频率的变化曲线，各频率对应的N也被标出。显然，球面ESPRIT在N≤ND的低频段具有较低误差，在N＞ND的高频段误差显著增大，而不论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≤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方法仍有效。

Claims (1)

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)建立协方差矩阵其中，矩阵矩阵是一组声源中，各个声源独自引起的信号的协方差矩阵的和，不包括不同声源引起的信号间的协方差成分；上标H表示转置共轭；
5)利用球面ESPRIT算法对协方差矩阵进行解算，确定声源波达方向；
建立声源波达方向测量模型，包括以下步骤：
3.1)计算第i个声源到第q个传声器(2)的传递函数t((ka,ΩMq)|ΩSi)，即：
式中，n和m分别为球谐函数的阶和次；bn(ka)为模态强度；为Ω方向的球谐函数；ΩSi表示第i个声源波达方向；i＝1,2,…,I；I为声源总数；k为声源辐射声波的波数；上标*表示共轭；
模态强度bn(ka)如下所示：
式中，jn(ka)为n阶第一类球贝塞尔函数，为n阶第二类球汉克尔函数；j'n(ka)和分别为n阶第一类球贝塞尔函数jn(ka)和n阶第二类球汉克尔函数的一阶导数；开口球表示传声器(2)布置在开口球体表面；刚性球表示传声器(2)布置在刚性球体表面；
3.2)计算Ω方向的球谐函数即：
式中，为连带勒让德函数；βn,m,κ为连带勒让德函数系数；球谐函数第n阶m次项对应的球谐系数
3.3)建立Q个传声器(2)所在方向的向量和对应的球谐函数向量建立I个声源所在方向的向量和对应的球谐函数向量
3.4)建立每个声源到所有传声器(2)的传递函数矩阵，记为上标H表示转置共轭；
其中，传声器所在方向对应的球谐函数矩阵如下所示：
式中，N0＝∞表示球谐函数最高阶；
声源所在方向对应的球谐函数矩阵如下所示：
阵列模态强度矩阵如下所示：
3.5)建立各传声器(2)测量得到的声压信号矩阵即：
式中，为噪声矩阵；信噪比SNR＝20lg(||P★-N||F/||N||F)；声源强度矩阵L为快拍总数；|| ||F表示Frobenius范数；上标★表示测量量；
3.6)对球谐函数矩阵球谐函数矩阵和模态强度矩阵的阶n进行截断，即令球谐函数矩阵球谐函数矩阵和模态强度矩阵的最高阶表示将数值向正无穷方向圆整到次近的整数；
基于最高阶N，更新声压信号矩阵P★如下：
建立协方差矩阵包括以下步骤：
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中元素
dθ(θSi)、dφ(φSi)为构成矩阵D的基础向量；
记向量的第(N+κ)(2N+1)+N-m+1个元素为An,mβn,m,κ，κ＝-n,…,0,…,n，并令其它元素为0，则球谐函数矩阵的转置共轭各传声器(2)测量得到的声压信号矩阵P★≈YMNBNGDS+N；
4.1.11)建立输入矩阵X，即：
式中，为S的第i行；||ψi||2＝1；
公式(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])；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)，初始化辅助矩阵Z0＝Λ0＝0，γ+1次迭代时的变量更新为：
4.3.4)对Hermitian拉格朗日乘子矩阵和辅助矩阵进行划分，得到：
4.3.5)基于步骤4.3.3)和步骤4.3.4)，更新公式(28)如下：
式中，I1和I2分别为L×L维和(2N+1)2×(2N+1)2维单位矩阵；为Tb(·)的伴随算子；对任意给定矩阵
M为对角矩阵；矩阵 为第κ(m)个对角线的元素均为1其它元素均为0的基本Toeplitz矩阵；
4.3.6)基于步骤4.3.3)，更新公式(29)如下：
公式(36)表示将Hermitian矩阵投影到半正定锥上，即对该Hermitian矩阵进行特征值分解，并令所有负特征值为0；
4.4)基于求得的协方差矩阵和建立协方差矩阵
确定声源波达方向，包括以下步骤：
5.1)特征值分解并以特征值大小对特征向量进行降序排列；将前s个特征向量写入矩阵US中；
5.2)基于球谐函数递归关系，建立与矩阵US相关的线性方程组，并采用最小二乘方法对与矩阵US相关的线性方程组进行求解，得到含有声源波达方向的转换矩阵；
线性方程组如下所示：
其中，为按上标(u,v)从矩阵US中选取一部分行组合而成的矩阵，ψxy+、ψxy-和ψz为含有声源波达方向的转换矩阵，和为系数矩阵；
其中，系数矩阵和系数矩阵分别如下所示：
5.3)采用广义特征值分解法对含有声源波达方向的矩阵进行特征值分解，确定声源波达方向。
CN202010545752.1A 2020-06-16 2020-06-16 基于原子范数的球面阵列声源波达方向估计方法 Active CN111812581B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010545752.1A CN111812581B (zh) 2020-06-16 2020-06-16 基于原子范数的球面阵列声源波达方向估计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010545752.1A CN111812581B (zh) 2020-06-16 2020-06-16 基于原子范数的球面阵列声源波达方向估计方法

Publications (2)

Publication Number Publication Date
CN111812581A CN111812581A (zh) 2020-10-23
CN111812581B true CN111812581B (zh) 2023-11-14

Family

ID=72845111

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010545752.1A Active CN111812581B (zh) 2020-06-16 2020-06-16 基于原子范数的球面阵列声源波达方向估计方法

Country Status (1)

Country Link
CN (1) CN111812581B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112363110B (zh) * 2020-11-30 2021-12-21 海南大学 一种基于嵌套交叉偶极子阵列的无网格单比特doa估计方法
CN113673317B (zh) * 2021-07-12 2023-04-07 电子科技大学 基于原子范数最小化可降维的二维离格doa估计方法
CN115407270B (zh) * 2022-08-19 2023-11-17 苏州清听声学科技有限公司 一种分布式阵列的声源定位方法

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7872583B1 (en) * 2005-12-15 2011-01-18 Invisitrack, Inc. Methods and system for multi-path mitigation in tracking objects using reduced attenuation RF technology
CN103888889A (zh) * 2014-04-07 2014-06-25 北京工业大学 一种基于球谐展开的多声道转换方法
CN107817465A (zh) * 2017-10-12 2018-03-20 中国人民解放军陆军工程大学 超高斯噪声背景下的基于无网格压缩感知的doa估计方法
CN107884741A (zh) * 2017-10-30 2018-04-06 北京理工大学 一种多球阵列多宽带声源快速定向方法
CN109932682A (zh) * 2019-02-19 2019-06-25 重庆工业职业技术学院 二维多快拍无网格压缩波束形成声源识别方法
CN110727914A (zh) * 2019-09-30 2020-01-24 中国人民解放军战略支援部队信息工程大学 一种基于向量运算的垂线偏差单点计算方法
KR20200020233A (ko) * 2018-08-16 2020-02-26 국방과학연구소 구형 마이크로폰 어레이를 이용한 음원의 입사 방향 추정방법
CN110907892A (zh) * 2019-12-05 2020-03-24 扬州大学 一种球麦克风阵列语音信号到达角估计方法

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7872583B1 (en) * 2005-12-15 2011-01-18 Invisitrack, Inc. Methods and system for multi-path mitigation in tracking objects using reduced attenuation RF technology
CN103888889A (zh) * 2014-04-07 2014-06-25 北京工业大学 一种基于球谐展开的多声道转换方法
CN107817465A (zh) * 2017-10-12 2018-03-20 中国人民解放军陆军工程大学 超高斯噪声背景下的基于无网格压缩感知的doa估计方法
CN107884741A (zh) * 2017-10-30 2018-04-06 北京理工大学 一种多球阵列多宽带声源快速定向方法
KR20200020233A (ko) * 2018-08-16 2020-02-26 국방과학연구소 구형 마이크로폰 어레이를 이용한 음원의 입사 방향 추정방법
CN109932682A (zh) * 2019-02-19 2019-06-25 重庆工业职业技术学院 二维多快拍无网格压缩波束形成声源识别方法
CN110727914A (zh) * 2019-09-30 2020-01-24 中国人民解放军战略支援部队信息工程大学 一种基于向量运算的垂线偏差单点计算方法
CN110907892A (zh) * 2019-12-05 2020-03-24 扬州大学 一种球麦克风阵列语音信号到达角估计方法

Non-Patent Citations (9)

* Cited by examiner, † Cited by third party
Title
Alternating direction method of multipliers for weighted atomic norm minimization in two-dimensional grid-free compressive beamforming;Yang 等;Journal of the acoustical society of America;第144卷(第5期);361-366 *
Direction of Arrive Estimation in Spherical Harmonic Domain Using Super Resolution Approach;Pan 等;Machine learning and intelligence communications second international conference;179-187 *
parametric direction-of-arrival estimation with three recurrence relations of spherical harmonics;Byeongho 等;The journal of the acoustical society of America;第145卷(第1期);480-488 *
Spherical harmonic atomic norm and its application to DOA estimation;Pan 等;IEEE ACCESS;第7卷;156555-156568 *
基于Karhunen-Loeve变换的正交椭圆球面波脉冲波形设计;王红星 等;电子与信息学报;第34卷(第10期);2415-2420 *
基于声压球谐函数分解的球面波束形成噪声源识别;褚志刚 等;农业工程学报;第28卷(第增1期);146-151 *
朱亚林.球谐波域连续压缩感知DOA估计.中国优秀硕士学位论文全文数据库信息科技辑.2020,(2020年第01期),I135-564. *
球谐域传播算子快速声定向算法;潘曦 等;兵工学报;第39卷(第10期);1936-1944 *
球谐波域连续压缩感知DOA估计;朱亚林;中国优秀硕士学位论文全文数据库信息科技辑(2020年第01期);I135-564 *

Also Published As

Publication number Publication date
CN111812581A (zh) 2020-10-23

Similar Documents

Publication Publication Date Title
CN111812581B (zh) 基于原子范数的球面阵列声源波达方向估计方法
CN111190136B (zh) 一种基于特定频率组合信号的一维doa估计方法
CN107247251B (zh) 基于压缩感知的三维声源定位方法
CN109655799B (zh) 基于iaa的协方差矩阵向量化的非均匀稀疏阵列测向方法
CN107544051A (zh) 嵌套阵列基于k‑r子空间的波达方向估计方法
CN106125041B (zh) 基于子空间加权稀疏恢复的宽带信号源定位方法
CN110082708A (zh) 非均匀阵列设计和波达方向估计方法
CN107290709A (zh) 基于范德蒙分解的互质阵列波达方向估计方法
CN106443587A (zh) 一种高分辨率的快速反卷积声源成像算法
CN112285647B (zh) 一种基于稀疏表示与重构的信号方位高分辨估计方法
CN106483503A (zh) 实心球阵列三维声源识别的快速反卷积方法
CN107576931A (zh) 一种基于协方差低维度迭代稀疏重构的相关/相干信号波达方向估计方法
CN109375152A (zh) 电磁矢量嵌套l阵下低复杂度的doa与极化联合估计方法
CN109343003B (zh) 一种快速迭代收缩波束形成声源识别方法
CN110837076A (zh) 一种基于张量分解的矢量水听器阵列方位估计方法
CN109870669A (zh) 一种二维多快拍无网格压缩波束形成声源识别方法
Yang et al. Enhancement of direction-of-arrival estimation performance of spherical ESPRIT via atomic norm minimisation
CN106908754B (zh) L型声矢量传感器阵列esprit解相干参数估计方法
CN115291169A (zh) 一种声源成像方法、系统、设备及存储介质
Yu et al. Adaptive imaging of sound source based on total variation prior and a subspace iteration integrated variational bayesian method
CN113866718A (zh) 一种基于互质阵的匹配场被动定位方法
CN109932682B (zh) 二维多快拍无网格压缩波束形成声源识别方法
Yin et al. Super-resolution compressive spherical beamforming based on off-grid sparse Bayesian inference
CN109061551B (zh) 一种基于多项式求根的无网格稀疏谱估计方法
CN113381793B (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