应用双线性模型的离格波达角估计方法*

王友顺，高童迪，袁正道,2，丁永春

(1.河南开放大学 人工智能工程研究中心，郑州450002；2.郑州大学 信息工程学院，郑州450001；3.中国船舶集团有限公司第七一三研究所，郑州450001)

0 引 言

利用传感器阵列估计多个远场信号的波达角(Direction of Arrival,DOA)是阵列信号处理的基础问题，在雷达、声呐和地震预测等领域得到了广泛的应用[1-2]。针对DOA估计问题，国内外研究者提出了多种不同类型的方法，如子空间分解、支持向量机、矩阵特征空间分类等[3-4]。

随着压缩感知技术的出现[5]，研究者们将DOA估计归纳为稀疏恢复问题[6]，并提出了多种不同的估计算法。例如，基于正交匹配追踪算法[7]，将波达角范围划分为均匀网格，假定实际的波达角准确位于某个网格上。相比于子空间划分等其他方法，压缩感知类方法通常能够获得更好的性能，特别是在快拍数量较少的情况下。稀疏贝叶斯学习[8]是一类经典稀疏估计算法，通过给变量增加稀疏引导先验，能够以较少的观测实现更优的估计性能。但是当真实来波方向并不精确位于网格上时，上述方案会导致较严重的估计偏差。更精细的网格划分一方面会导致较高的运算复杂度，另一方面会使感知矩阵的列向量之间具有较高的相关度[9]，导致估计算法收敛性变差。

近年来国内外多个团队提出了基于离格信号的稀疏重构算法。文献[9]提出了一种基于稀疏最小二乘的离格算法，假定真实角度不必限于网格之上，能够计算出误差并对估计值进行修正。文献[10]则提出了一种求根稀疏贝叶斯学习方法，在线性阵列场景下能够取得很好的估计性能。由于复杂度较低，经典的近似消息传递算法(Approximate Message Passing,AMP)在稀疏贝叶斯学习框架下得到了广泛应用[11]，但是普通的AMP算法在遇到较“坏”的观测矩阵时(如病态条件数、列相关、低秩等)会导致迭代过程中出现发散[12]。特别是在DOA场景下，阵列流形具有较强的相关性，普通的AMP算法不稳定。为了解决上述问题，多种不同的改进算法相继被提出：文献[13]提出了一种向量化AMP算法(Vector-AMP)，通过两次标量化-期望传播方法，大大提高了鲁棒性；在此基础上，文献[14]归纳了双线性VAMP算法(Bilinear Adaptive-VAMP,BAd-VAMP)，将VAMP算法扩展到了双线性场景，在字典矩阵学习、矩阵补全等多个场景下得到了应用。

为解决离格波达角估计中复杂度和鲁棒性的矛盾，本文提出将DOA估计模型转换为双线性问题，利用文献[14]中所提BAd-VAMP算法解决该问题。本文的创新点如下：第一，提出将离格DOA估计模型转换为双线性问题，并利用BAd-VAMP算法对双线性模型进行推导和运算；第二，充分利用离格DOA模型中待估计向量的相同稀疏特征减小迭代中矩阵的维度，降低了迭代算法的复杂度；第三，提出一种基于残差的冷/热重启的迭代机制，避免了双线性算法容易陷入鞍点的问题。在仿真部分，本文建立DOA估计模型对所提算法与文献中已有方法进行对比，结果表明所提算法在机载雷达等快拍数较少的场景下具有明显的性能优势。

本文符号说明：Y∈N×L表示维度为N×L的复矩阵，IN表示N维单位阵，(·)T和(·)H分别代表转置和共轭转置，CN(x;μ,ν)表示均值为μ方差为ν的复高斯分布，x∈U(a,b)表示变量x服从区间(a,b)的均匀分布，对向量x取平均表示为〈x〉，矩阵或向量的2范数表示为‖‖2，tr(A)表示矩阵的求迹运算，diag(a)表示利用向量a构建对角阵。

1 离格DOA模型

(1)

式中：yt∈M×1为第t个快拍下M个阵元的接收向量；K×1表示第t个快拍下的K个来波信号；w表示均值为0、方差为1/λ的高斯白噪声，其中λ表示噪声精度。由于噪声方差不随时间变化，本文中省略w的上标t。阵列流形表示为

A(θ)=[a(θ1),a(θ2),…,a(θN)]∈M×N。

要提高估计精度，最直接的方法是将信号空间划分为更小的网格，即加大N。但更大的N会导致复杂度提升，且阵列流形A(θ)的列相关性变大。相关性高的阵列流形A(θ)会导致估计性能变差，对于消息传递等迭代类算法甚至会引起发散，导致估计失败。

(2)

式中：

e(θn)=a′(θn)为误差矢量，其中a′(θn)表示对导向矢量a(θn)依θn求导。式(2)所示的接收模型可以重写为

yt=(A(θ)+E(θ)diag(β))st+w。

(3)

式中：误差矩阵和误差向量分别定义为E(θ)=[e(θ1),e(θ2),…,e(θN)]和β=[β1,β2,…βN]T。当DOA估计模型中有L个快拍时，上式可以写为多重观测形式：

Y=(A(θ)+E(θ)diag(β))S+W。

(4)

式中：矩阵Y=[y1,y2,…,yL]和S=[s1,s2,…,sL]分别表示观测矩阵和待估计矩阵，且S具有相同列稀疏特征，即每列非零元素位置相同。由于上式中向量θ为固定值，则可简记为

Y=(A+Ediag(β))S+W。

(5)

值得注意的是，上述模型中未知参数有S和β，在信号处理领域，通常将利用一组观测同时估计两组独立参数的模型称为双线性模型，在下一节中将式(5)重新建模为双线性形式。

2 双线性系统建模

通用双线性模型可以表示为[14-15]

Y=Φ{β}C+W。

(6)

式中：Y=[y1,y2,…,yL]和C=[c1,c2,…,cN]分别表示观测和待估计矩阵。目前针对上述双线性系统的方法有俄亥俄州立大学Schniter教授提出的PBiGAMP[15]和BAd-VAMP[14]两种算法。根据文献中的仿真结果可知BAd-VAMP从复杂度和性能上全面超越PBiGAMP算法，所以本文选择BAd-VAMP算法作为算法框架。

通过对式(5)和式(6)的类比，可以将式(5)所示的DOA估计转换为双线性模型。定义式(5)中的流形和误差矩阵为

Φ{β}=A+∑nβnEn。

式中：En=[0,…,e(θn),…,0]∈M×N。则DOA估计模型可以写为双线性形式：

Y=(A+∑nβnEn)S+W=Φ{β}S+W。

(7)

根据式(7)所示的数据之间的约束关系，可以对系统中所有已知和未知变量的联合分布函数利用贝叶斯公式进行因式分解[12-13]：

p(Y,S,β,λ)=p(Y|S,β,λ)p(S|γ)p(γ)p(β)p(λ)=

(8)

3 基于BAd-VAMP的DOA估计算法

BAd-VAMP算法可以分为两部分，分别针对参数S和β估计，本节也按照以上划分进行推导。

3.1 参数S估计部分

定义

(9)

定义

(10)

由于函数节点fδ为等效节点，则消息

其中均值方差分别为

(11)

然后再利用MF规则估计超先验参数γ为

(12)

其中期望和方差分别计算为

(13)

从而完成S估计部分的迭代计算。

3.2 参数β估计部分

E[lnp(Y|S2,β,λ)]=

Const-λ(βHHβ-βHα-αHβ)。

(14)

式中：H∈N×N，α∈N×1，Const表示常数。式(14)中矩阵H和向量α中每个元素计算为

(15)

(16)

其中：i,j∈[1:N]，

(17)

(18)

式中：β为前一次迭代中求得的估计值，据此可得到Φ{β}。对式(14)所示期望依β求梯度，并令其等于零，可得到使期望取最大值的β值，即

βmax=H-1α。

3.3 算法改进

(1)迭代阻尼

(2)复杂度降低

(3)局部极值点判决

由于双线性问题的“非凸”特性，使得迭代算法可能会陷入局部极值点(鞍点)，给整体的估计精度带来较大损失。鞍点在双线性问题中容易识别，即判定收敛后的残差r=‖Y-Φ{β}S‖2是否大于阈值，如果是则可判定为鞍点。此外为提升算法收敛速度，本文提出了一种冷/热重启的解决方案，实现方式是最初迭代时仅迭代较少次数，选择初始化NOuter=10，每次迭代进行收敛性判断，当识别为收敛即非鞍点时，采用热重启方式，将本次所估计的矩阵S作为入口参数进入重启迭代；反之识别为发散或鞍点，则采用冷重启方式，即将矩阵S设置为随机。

4 数值仿真和复杂度分析

本文选择均匀线性阵列，阵元个数M=30，则虚拟网格个数N=30，即网格宽度为6°。设置快拍数为L=2～12，设置阈值γmax=103。为衡量本文所提算法的有效性，本节选择文献中已有的网格化稀疏贝叶斯学习(Off-Grid Sparse Bayesian Learning,OGSBL)、求根SBL(Root SBL,RTSBL)和高斯-赛德尔根(Gauss-Seidel Root,GSROOT)作为对比，此外还选择克拉美罗界(Cramer-Rao Bound,CRB)作为参考。

4.1 性能分析

图2给出了一次蒙特卡洛仿真中估计值和真实值的对比，仿真中选择真实角度为[-28.6°,-18.6°,3.5°,15.6°,31.7°]。本仿真中设置信噪比为30 dB。从图2中可以看出，本文所提算法能够精确捕捉到5个来波信号的离格角度和准确功率，在无来波信号的角度估计结果趋近于0。

图2 真实角度/功率与估计值对比

图3对比了BAd-VAMP与已有算法归一化均方误差(Normalized Mean Square Error,NMSE)随快拍次数的变化曲线，数值仿真结果为500次蒙特卡洛仿真的平均值。每次仿真中波达方向以[-30°,-18°,6°,18°,30°]为基础，叠加U(-3°,3°)的随机角度偏移产生。图3中设置信噪比SNR=20 dB，来波信号个数K=5。从仿真结果中可以看出，本文所提算法在快拍数较多时性能与现有的RTSBL和GSROOT较为接近。但是随快拍数的减少，RTSBL和GSROOT方法性能明显恶化，使得本文所提算法表现出较大的性能增益。从图中还能看出，由于不考虑网格偏差，OGSBL算法性能较差，从而证明了离格估计算法的必要性。

图3 估计性能与快拍个数关系图(SNR=20 dB)

图4对比了在不同快拍个数下估计性能随信噪比的变化曲线，从图中可以得出与图3相同的结论，即在快拍个数较少时本文所提算法具有明显的性能优势。要达到相同的估计性能，本文所提算法需要的快拍个数较少，更适合于机载雷达等快时变场景。

图4 估计性能随信噪比变化曲线

4.2 复杂度分析

本文所提算法可以划分为S估计和β估计两部分。S估计在式(9)需要矩阵乘法ΦT{β}Φ{β}，具有复杂度O(MN2)，其余运算均为标量形式，复杂度为O(MNL)。β估计部分复杂度集中于式(15)～(19)，其中式(15)中利用Ei矩阵的特征和矩阵求迹的特点，可知只需要计算i=j情况，所以式(15)具有复杂度O(MN)，式(17)需要计算L个矩阵求逆，具有极高的复杂度。为此,文献[14]提出了一种奇异值分解(Singular Value Decomposition,SVD)的替代方法。通过SVD分解得

Udiag(η)V=SVD(ΦT{β}Φ{β})，

则式(17)可等效计算为

从而将式(17)复杂度降低为O(MN2)。式(19)需要矩阵求逆操作，具有复杂度O(N3)，所以本文所提算法整体复杂度可写为O(MN2+MNL+N3)。根据3.2节所提方法，通过设置阈值可以随迭代降低式(19)中的矩阵H维度，从而使得本文所提算法复杂度进一步降低。参考文献[10]中的复杂度分析，表1给出了各算法的复杂度对比，可以看出本文方法与文献中现有方法保持同阶复杂度。

表1 复杂度对比

5 结 论

本文将离格波达角估计问题建模为典型双线性问题，通过最新的BAd-VAMP算法对其进行估计。利用离格DOA估计模型的共同稀疏性等特点，本文还对BAd-VAMP算法进行了改进，降低了算法复杂度，并且引入迭代阻尼提升了鲁棒性。仿真结果表明，相比于现有的RTSBL等方法，本文所提算法在机载雷达等快拍数较少时表现出较大的性能优势，具有较高的应用价值。