阵列互耦情况下基于稀疏贝叶斯学习的离网格DOA估计

王绪虎，白浩东，张群飞，田 雨

(1.青岛理工大学 信息与控制工程学院，青岛 266520；2.西北工业大学 航海学院，西安 710072)

波达方向(direction of arrival,DOA)估计问题(在很多领域受到极大的关注，例如，雷达[1]、声呐[2-4]，生物医学[5]等。在过去几十年已经提出了许多高分辨的DOA估计算法，例如，多重信号分类算法[6](multiple signal classification,MUSIC)、求根MUSIC算法[7](Root-MUSIC)、旋转不变子空间算法[8](estimating signal parameters via rotational invariance techniques,ESPRIT)等，但上述高分辨算法在低信噪比、少快拍数的情况下，估计性能会严重下降，甚至失效。

随着压缩感知(compress sensing,CS)理论和稀疏重构技术的不断发展和成熟，许多学者将其与DOA估计联系起来，使DOA估计技术进入了新的发展阶段。相较于传统的高分辨算法，基于CS和稀疏重构技术的DOA估计方法，在低信噪比、少快拍的情况下表现出良好的估计性能。该种方法大致分为三类：凸优化方法、贪婪方法[9]、稀疏贝叶斯学习方法[10-11]。例如，文献[12]研究稀疏信号的表示以及DOA估计问题，提出L1范数奇异值分解(L1norm-singular value decomposition,L1-SVD)方法，把测向问题转化为求解L1范数的问题，但是其优化问题受到模型正则化参数的影响，影响估计的精度。王伟东等提出的基于稀疏信号功率迭代补偿的DOA估计算法，该算法基于补偿原理，使信号的稀疏表示近似于L0范数，把DOA估计问题转化为求解L0范数问题。相比于L1-SVD算法，该算法不需要设置正则化参数，对非相干源具有较高的估计精度。然而，凸优化方法的计算效率限制了其进一步发展。

随着CS理论研究的不断深入，稀疏贝叶斯学习方法(sparse Bayesian learning,SBL)被认为与凸优化方法具有相同的全局收敛性，而它计算效率又远远优于后者。Ji等将SBL引入CS领域中，提出贝叶斯压缩感知算法。该算法通常需要满足声源的来波方向位于网格点上，才能实现较高的方位估计精度。当入射声源的波达方向偏离预先划分好的网格，就导致方位估计精度的降低。因此，Yang等[13]提出了离网格的稀疏贝叶斯学习算法 (off-grid sparse Bayesian learning,OGSBL)，在离网格模型中，在信号真实到达角处采用一阶泰勒展开式近似表示，使得估计性能进一步提高。文献[14]中，提出了求根离格稀疏贝叶斯学习算法(root-off grid sparse Bayesian learning,Root-OGSBL)，降低了OGSBL方法的计算复杂性，在网格间距较小的情况下，提高了计算效率的同时保证了估计的精度。文献[15]中，使用平均场理论中的变分理论来估计参数的后验分布，并且提出了三层信号先验分布，这种分层先验进一步提升了信号的稀疏性，提升了方位估计的精度。文献[16]中，提出了稀疏信号新的先验分布框架—分层合成Lasso先验，相比于假设伽马先验分布，具有更高的稀疏性、更低的重建误差，提升了方位估计精度。

在实际的声呐测向系统中，不可避免的存在阵列误差。例如，水听器阵元间的耦合、水听器位置的偏差、水听器阵元的幅相通道不一致。导致现有常规算法的方位估计精度下降，甚至失效。在文献[17-18]中，提出了阵元间存在未知互耦下波达方向估计算法，解决了阵元间存在未知互耦情况下，稀疏信号的恢复问题。针对稀疏模型中入射信号方向与离散网格存在误差，且阵元间的未知互耦没有被同时考虑的问题，本文提出了一种水听器阵元间存在不确定互耦，基于稀疏贝叶斯学习的离网格DOA估计方法。

1 信号模型

考虑K个不同方向的远场窄带信号入射到如图1所示的均匀线列阵，该阵列包含M个全向阵元，阵元之间存在互耦，阵元间距为d，因此，该接收模型可以表示为

图1 M个水听器组成的均匀线列阵Fig.1 Uniform linear array composed of M-element hydrophones

n(t)

(1)

A=[a(θ1),a(θ2),…,a(θK)]

(2)

式中，a(θk)=[a1(θk),…,am(θk),…,aM(θk)]T为第k个入射信号的导向矢量，am(θk)=ej2π(m-1)d/λsin θk，λ为信号的波长。

式(1)为单快拍观测矢量模型，当水听器阵元接收多快拍数据时，多快拍观测矢量模型可以表示为

Y=CAS+N

(3)

式中，Y=[y(1),y(2),…,y(T)]∈CM×T为多快拍的接收数据矩阵，S=[s(1),s(2),…,s(T)]∈CK×T为多快拍的入射信号数据矩阵，N=[n(1),n(2),…,n(T)]∈CM×T为多快拍的阵列接收噪声矩阵，T为快拍数。

在本文中，我们研究水听器阵元间存在未知耦合情况下的DOA估计问题。由文献[19]可知，水听器阵元间的互耦可以用矩阵C表示

(4)

该矩阵为对称的托普利兹矩阵。根据文献[20]的引理三，我们将上述模型进行处理。首先，定义向量c=[1,c1,c2,…,cM-1]T，且满足C=toeplitz(c)，toeplitz(c)为使用向量c生成的循环对称矩阵 ，在t时刻的接收信号可以写为

yt=CAst+nt=Q(st⊗c)+nt

(5)

则对应水听器阵列的多快拍数据可表示为

Y=Q(S⊗c)+N

(6)

2 稀疏贝叶斯学习模型

2.1 稀疏信号模型

(7)

式(7)中的系统模型就可以转化为过完备的稀疏模型

(8)

式中，矩阵X是矩阵S的稀疏扩展

X=[x1,x2,…,xT]∈CN×T

(9)

图2 信号X结构稀疏化模型Fig.2 The structure of sparse matrix X

(10)

在实际的阵列测向系统中，接收信号并不是落在预设的网格点上，此时，采用离网格信号模型，使用离接收信号最近的网格矢量一阶泰勒展开式近似表示

(11)

图3 离格阵列空域模型Fig.3 Off-grid array spatial domain model

因此，由式(8)和(11)可以表示阵元间存在互耦下的离网格模型

Y≈Ω(δ)(X⊗c)+N

(12)

最终，水听器阵元间存在未知耦合的离网格稀疏模型建立如式(12)，接下来通过重建稀疏矩阵X来估计入射信号的波达方向，重建稀疏矩阵X中的非零行代表入射信号的实际波达方向。

2.2 模型参数分布假设

2.2.1 噪声及其精度分布

假设接收的噪声为独立同分布的复对称高斯白噪声，其分布为

(13)

p(αn)=Γ(αn|a,b)

(14)

稀疏模型的似然函数可表示为

(15)

2.2.2 信号及其精度分布

假设稀疏信号X的多快拍样本相互独立，各列服从均值为0、方差为，

(16)

(17)

其中，c和d为伽马分布的超参数。

2.2.3 互耦矢量及其精度分布

(18)

(19)

式中，e和f为伽马分布的超参数。

2.2.4 离网格矢量分布

p(δ)=U(-r/2,r/2)

(20)

通过组合式(14)～(20)，可得模型的联合概率密度

p(X,Y,δ,c,αx,αc,αn)=

p(Y|X,δ,c,αn)p(X|αx)p(c|αc)p(αn)

p(αx)p(αc)p(δ)

(21)

3 稀疏贝叶斯参数学习

利用上述的分布假设，我们提出了一种水听器阵元间存在互耦，基于稀疏贝叶斯学习的离网格波达方向估计算法。该算法中，我们使用EM算法学习更新离网格误差矢量和互耦矢量，直至收敛。

首先，利用上述接收的数据Y和参数先验分布，可以得出信号X的后验分布

p(X|Y,δ,c,αx,αc,αn)=

αn)p(X|αx)

(22)

式中，p(X|αx)的概率密度函数为

(23)

X后验分布也为高斯分布，满足

(24)

式中，μt为均值向量，Σx为协方差矩阵

(25)

Σx=[αnΥH(δ,c)Υ(δ,c)+diag(αx)]-1

(26)

其中，Υ(δ,c)=Ω(δ)(IN⊗c)，另外，定义所有快拍下均值矢量组成的矩阵

μ[μ1,…,μt,…,μT]

(27)

其次，我们对超参数进行学习，超参数的学习过程等价于最大化所有待估计参数的后验概率分布

(28)

最大化待估计参数的后验概率分布又等价于最大化似然函数

L(δ,c,αx,αc,αn)=

ε{lnp(Y|X,δ,c,αn)p(X|αx)p(c|αc)

p(αn)p(αx)p(αc)p(δ)}

(29)

式中，ε{·}表示条件后验期望EX|Y,δ,c,αx,αc,αn{·}。接下来，详细推导超参数迭代更新的表达式，下述推导中涉及到矩阵对向量变量的求导，参考附录B。

忽略式(29)中与互耦矢量无关的概率项，可得到关于互耦矢量的似然函数

L(c)=ε{lnp(Y|X,δ,c,αn)p(c|αc)}=

(30)

因为EM算法是不断收敛递增的，所以估计参数就等价于求似然函数的极值，有如下关系

(31)

因此，对式(31)求导，可得互耦矢量的迭代表达式(32)

(32)

令式(32)等于零，可得

(33)

式中，yt表示阵列接收数据Y的第t列，μt表示稀疏信号的均值矩阵μ的第t列。

同理，忽略与信号精度无关的概率分布，可得到如下似然函数表达式

L(αx)=ε{lnp(X|αx)p(αx)}

(34)

对似然函数求导，令其等于零，可得信号精度的迭代表达式

(35)

Σx(n,n)表示协方差矩阵的第n行n列的元素，μn,t表示稀疏信号的均值矩阵μ的第n行t列的元素。

同理，对于噪声精度，我们推导其似然函数和迭代表达式，如下

L(αn)=ε{lnp(Y|X,δ,c,αn)p(αn)}

(36)

(37)

对于互耦精度，给出似然函数和迭代表达式

L(αc)=ε{p(c|αc)p(αc)}

(38)

(39)

对于离网格间距，给出似然函数和迭代表达式

L(δ)=ε{lnp(Y|X,δ,c,αn)p(δ)}

(40)

(41)

其中，P∈RN×N，v∈RN×1下式(42)(43)给出

(42)

(43)

最后，基于上述分析算法总结如下：

第一步：输入水听器阵元接收数据Y；

第二步：初始化需要更新的参数X,δ,c,αx,αc,αn=mean(var(Y))，超参数a,b,c,d,e,f的设置[15]；

第三步：根据(25)(26)更新后验期望矢量和协方差矩阵；

第四步：根据(35)更新信号精度αx；

第五步：根据(37)更新噪声精度αn；

第六步：根据(33)和(39)更新互耦矢量及其精度；

第七步：根据(41)更新离网格间隔矢量；

4 仿真分析

4.1 空间谱分析

根据文献[21]，水听器间的耦合效应可以表示为

(44)

其中，互耦系数中的幅度ξ服从均匀分布，即ξ～U[-0.05,0.05]，互耦系数的相位φ也服从均匀分布，即φ～U[0,2π]，参数αc为耦合系数，表征相邻阵元之间的耦合效应。理论上，相邻阵元间的耦合效应是相同的，仿真参数的设置如表1所示。

表1 仿真参数Tab.1 Simulation parameters

阵元间的耦合系数αc为-15 dB，正则化参数为2时，四种方法归一化空间谱如图4所示；阵元间的耦合系数αc为-15 dB，正则化参数为3时，四种方法归一化空间谱如图5所示；阵元间的耦合系数αc为-5 dB，正则化参数为3时，四种方法归一化空间谱如图6所示。

从图4至图6可以看出MUSIC方法的空间谱主瓣宽度较宽，旁瓣较高，当增大耦合系数时，两个主瓣的峰值和真实的角度值相差较大；对比图4和图5，可以看出当L1-SVD方法正则化参数为2时，主瓣宽度较窄，旁瓣较低；改变正则化参数为3时，主瓣宽度变宽，且出现了许多旁瓣；对比图5和图6可以看出，当耦合系数增大时，OGSBL方法旁瓣数量增加，且旁瓣谱级升高；对比图4～图6，本文提出的方法主瓣宽度和旁瓣谱级基本不变。

图4 空间谱(正则化参数为2，αc=-15 dB)Fig.4 Spatial spectrum(Regularization parameters 2, αc=-15 dB)

图5 空间谱(正则化参数为3，αc=-15 dB)Fig.5 Spatial spectrum(Regularization parameters 3， αc=-15 dB)

图6 空间谱(正则化参数为3，αc=-5 dB)Fig.6 Spatial spectrum(Regularization parameters 3， αc=-5 dB)

四种方法的蒙特卡洛仿真结果如图7所示，其中L1-SVD方法的正则化参数设置为3。图7(a)是耦合系数为-15 dB时、图7(b)是耦合系数为-10 dB时、图7(c)是耦合系数为-5 dB时50次蒙特卡洛仿真的统计结果。从图7可以看出，水听器阵元间存在耦合时，本文方法的估计性能和稳健性是最优的。

(a) 耦合系数αc=-15 dB

L1-SVD方法在不同正则化参数时的蒙特卡洛仿真结果如图8所示。图8(a)为正则化参数为2时、图8(b)为正则化参数为3时，50次蒙特卡洛仿真统计结果。从图8可以看出，改变正则化参数对L1-SVD方法的估计性能影响较大，因而该方法的实用性较差。

(a) 正则化参数为2

综上所述，当水听器阵元间存在互耦时，相比于MUSIC、L1-SVD、OGSBL方法，本文提出方法的估计性能最优，稳健性最高。

阵元间耦合系数为-5 dB，正则化参数为3时，四种方法信号估计结果如表2所示，表中的估计结果为50次蒙特卡洛仿真试验的统计均值。从表中可以看出，本文方法的估计结果与真实角度的偏差最小，表明本文方法的估计性能优于MUSIC、L1-SVD、OGSBL方法。

表2 真实角度与本文方法估计DOA结果对比Tab.2 Comparison of the real angle and the DOA estimation result of the algorithm proposed in this paper

4.2 分辨力分析

不同方法的分辨能力不同，本小节从信噪比的高低、阵元间的耦合系数的大小研究本文提出方法的角度分辨能力。当满足下述的条件时，则认为能够正确分辨两个角度[22]。

(45)

仿真试验条件如表1所示，设定耦合系数αc=-7 dB，图9给出了分辨成功率随角度间隔的变化规律，其中(a)和(b)子图分别为信噪比为0和5 dB时的统计结果。从图中可以看出，信噪比为0时，三种方法成功估计出两个角度的间隔分别为9°(本文方法)，10°(OGSBL)，12°(MUSIC)。本文方法的角度分辨力优于OGSBL和MUSIC方法。从图9中可以看出，信噪比为0和5 dB时，可成功估计两个角度的间隔为9°和7°。随着信噪比的增加，本文方法的分辨力逐渐增强。

(a) SNR=0 dB 分辨成功率和角度间隔的关系

仿真试验条件如表1所示，设定信噪比为5 dB，图10给出了分辨成功率随角度间隔的变化规律，其中(a)和(b)子图分别为耦合系数为-12 dB和-3 dB时的统计结果。从图中可以看出耦合系数为-12 dB时，三种方法成功估计出两个角度的间隔分别为6°(本文方法)，7°(OGSBL)，8°(MUSIC)。本文方法的角度分辨力优于OGSBL和MUSIC方法。从图10中可以看出，耦合系数为-12 dB和-3 dB时，可成功估计两个角度的间隔为6°和11°。随着耦合系数的增加，本文方法的分辨力逐渐降低。

(a) 耦合系数αc=-12 dB，分辨成功率和角度间隔的关系

仿真试验条件如表1所示，设定入射角度间隔为8°和耦合系数为-7 dB时，图11给出了分辨成功率随信噪比的变化规律；设定入射角度间隔为11°和信噪比为5 dB时，图12给出了分辨成功率随耦合系数的变化规律。

图11 分辨成功率和信噪比的关系Fig.11 The relationship between resolution probability and SNR

从图11可以看出随着信噪比的增加，三种方法的分辨成功率逐渐增加。从图中可以看出，角度间隔为8°时，本文方法和OGSBL方法可成功估计出两个角度的信噪比分别为1 dB和4 dB，而MUSIC算法无法完全正确估计出两个真实的入射角度。本文方法的可成功估计两个信号的信噪比低于OGSBL方法，本文方法的分辨力优于OGSBL方法和MUSIC方法。

从图12可以看出随着耦合系数的增加，三种方法的分辨成功率逐渐降低。从图中可以看出，角度间隔为11°时，三种方法成功估计出两个角度的耦合系数分别为-3 dB(本文方法)，-7 dB(OGSBL)，-7 dB(MUSIC)。本文方法的可成功估计两个信号的耦合系数高于OGSBL方法和MUSIC方法，本文方法的分辨力优于OGSBL方法和MUSIC方法。

图12 分辨成功率和耦合系数的关系Fig.12 The relationship between resolution probability and coupling coefficient

综上所述，本文方法的角度分辨力优于OGSBL方法和MUSIC方法。当耦合系数一定时，分辨成功率随着信噪比的增加逐渐提升；当信噪比一定时，分辨成功率随着耦合系数的增加逐渐降低。

4.3 估计精度分析

本小节从信噪比、快拍数、网格间隔、阵元间的耦合系数研究不同方法的方位估计精度。把均方根误差(Root Mean Square Error，RMSE)作为衡量方位估计精度的标准，可表示为

(46)

仿真试验条件如表1所示，设定信噪比为-2 dB到10 dB变化，步长为2 dB，正则化参数为2，图13给出了均方根误差随信噪比的变化规律。从图中可以看出，当信噪比为10 dB时，本文方法的估计误差为0.4°低于0.45°(OGSBL)、0.6°(L1-SVD)、0.7°(MUSIC)。随着信噪比的增大均方根误差逐渐减小，且本文方法的估计精度高于OGSBL 方法、L1-SVD方法和MUSIC方法。

图13 不同信噪比，四种方法的方位估计精度Fig.13 The DOA estimation performance of the four algorithms with different SNR

仿真试验条件如表1所示，设定快拍数为50到300变化，步长为25，正则化参数为2，图14给出了均方根误差随快拍数的变化规律。从图中可以看出，当快拍数为300时，本文方法的估计误差为0.35°低于0.57°(OGSBL)、0.62°(L1-SVD)、0.9°(MUSIC)。随着快拍数的增大均方根误差逐渐减小，且本文方法的估计精度高于OGSBL 方法、L1-SVD方法和MUSIC方法。

图14 不同快拍数下，四种方法的方位估计精度Fig.14 The DOA estimation performance of the four algorithms with different snapshot numbers

仿真试验条件如表1所示，设定网格间隔为{1°,1.5°,2°,2.5°,3°,4°,5°,6°}变化，图15给均方根误差随网格间隔的变化规律。当网格间隔小于2°时，本文提出的方法与OGSBL方法相比，估计精度相差不大，当网格间隔增大到5°时，OGSBL估计误差将增大，高于本文方法。随着网格间隔增大，均方根误差增大，但本文方法的估计精度要优于OGSBL方法。

图15 不同信噪比，网格间隔与方位估计精度的关系Fig.15 Comparison of the relationship between grid spacing and DOA estimation performance with different SNR

仿真试验条件如表1所示，设定耦合系数为-15 dB到-3 dB，步长为2 dB，图16给出了均方根误差随耦合系数变化的规律。当信噪比一定时，增大耦合系数，估计误差逐渐增大，估计精度逐渐降低；当耦合系数为-3 dB时，增大信噪比为5 dB时，估计误差减小3°，继续增大信噪比为10 dB时，估计误差减小1°。当耦合系数一定时，随着信噪比的增大，本文方法估计误差逐渐减小，估计精度逐渐提升。

图16 不同信噪比条件下，耦合系数与方位估计精度的关系Fig.16 Compare the relationship between the adjacent coupling coefficient and the DOA estimation performance with different SNR

5 结 论

为了解决水听器阵元间存在未知耦合，现有方法波达方向估计精度较低、角度分辨力差的问题，本文提出了一种基于稀疏贝叶斯学习的DOA离网格估计方法。该方法通过建立新的贝叶斯学习模型，使用EM算法，推导离网格矢量和互耦矢量的分布，重新表示空间谱函数，最后通过谱峰搜索估计出波达方向。仿真试验结果表明：当水听器阵元间存在互耦时，本文方法的方位估计性能和稳健性优于MUSIC、OGSBL和L1-SVD方法；本文方法的角度分辨力优于MUSIC和OGSBL方法。鉴于上述试验使用的是两个不相关信号，在后续的研究中，将会对相关信号做进一步研究。

附录A

引理1：矩阵C和向量c定义如正文，对于任意向量a(θk)，我们有：

Ca=Qk(a(θk))c

式中，Qk(a(θk))∈CM×M是两个矩阵Qk1(a(θk))、Qk2(a(θk))的和：

因此，Qk(a(θk))=Qk1(a(θk))+Qk2(a(θk))，Qk1(a(θk))和Qk2(a(θk))分别为

根据上述引理1可得CAst=Q(st⊗c)，即

附录B

有复向量u∈CP×1，v∈CP×1，复矩阵A∈CM×P是复向量x∈CN×1的函数，可得