基于可分离替代函数算法的DOA估计方法

郑文康, 魏志晴, 白艳萍, 黄嘉俊, 禹秀梅, 谭秀辉, 王 鹏

(中北大学 数学学院, 山西 太原 030051)

0 引言

波达方向(Direction of Arrival,DOA)估计的应用非常广泛,其在航海、军事、航空和日常生活中发挥着重要的作用[1-3].最早的DOA估计算法是用空域阵元采集的信息代替传统时域处理中的时域信息,本质是在空域中使用传统的时域傅立叶谱估计方法.常规波束形成(Conventional Beamforming,CBF)[4]就是这类算法,然而这类算法对处在一个波束宽度内的目标无法做出准确分辨,因此存在很大局限性.

传统的DOA估计算法主要是子空间分解类算法和子空间拟合类算法.子空间分解算法主要分为信号子空间类算法和噪声子空间类算法.多重信号分类(Multiple Signal Classification,MUSIC)[5,6]算法是经典的噪声子空间类算法,该算法利用信号子空间和噪声子空间的正交性来估计信号的入射方向.旋转不变子空间(Estimation of Signal Parameters using Rotational Invariance Techniques,ESPRIT)[7]算法是经典的信号子空间类算法,该算法利用的是信号子空间的旋转不变性求出入射角方向信息.子空间分解类算法在小快拍或者入射信号为相干源的情况下性能会受到严重影响.子空间拟合类算法中具有代表性的算法是最大似然(Maximum Likelihood,ML)[8]算法,该算法利用了阵列流形矩阵和阵列输出数据子空间之间的拟合关系,思想较为简单,然而该算法复杂度较高,运算量大.

近年来,许多学者将阵列接收数据转化到空域上进行稀疏表示,利用阵列接收信号的空域稀疏性将DOA估计问题转化为阵列信号在过完备基上的稀疏系数求解问题,这些研究源自于Donoho[9]提出的压缩感知理论.过完备基上的稀疏系数求解问题是一个非凸问题,本质是求解欠定方程,因此很多学者使用贪婪类算法(Greedy Algorithm,GA)进行信号重构,如匹配追踪(Matching Pursuit,MP)算法、正交匹配追踪(Orthogonal-Matching-Pursuit,OMP)算法、迂回式匹配追踪(Detouring Matching Pursuit,DMP)算法等[10].贪婪类算法需要预先确定入射信号个数,在实际使用中具有局限性.在一定条件下,l0-范数非凸问题可以转化为l1-范数的凸优化问题,通过凸优化算法来求解,目前已经有许多凸优化算法应用在信号重构中,如基追踪降噪(Basis Pursuit Denoising,BPDN)[11]和套索回归(Least Absolute Shrinkage and Selection Operator,LASSO)[12],然而这些算法在DOA估计中性能较差,如果划分空间网格较细会导致计算耗时过长,如果网格划分较粗则无法保证计算精度,在低信噪比环境中重构精度低,并且随着入射信号快拍数增多算法计算效率会大大降低.

迭代收缩算法也是最近几年出现的一类用于求解l1-优化问题的算法,这类算法对于优化问题的求解更加高效,在迭代过程中会对解进行标量收缩操作,过去几年这类算法得到了严格的理论证明,其收敛速度也有充分的研究[13].基于Daubechies等[14]的工作,Michael Elad[15]提出了可分离替代函数(Separable Surrogate Functions,SSF)算法,该算法是迭代收缩算法的一种.

为了解决现有算法在低信噪比条件下估计精度差、执行效率低等问题,本文首次将迭代收缩算法思想应用到DOA估计中,提出了一种基于可分离替代函数算法的DOA估计方法,并在单快拍DOA估计模型的基础上提出了基于信号多快拍数据的DOA估计模型,该方法可利用入射信号的多快拍数据进行DOA估计,无需预先设定入射信号个数,运算速度快,在低信噪比的场景下具有优良的性能.仿真实验部分验证了算法的性能.

1 基本理论

1.1 矢量水听器阵列信号模型

传统的声压水听器阵元只能接收到入射信源的声压信息,而矢量水听器阵元可以空间共点同步拾取声场信息中的声压与质点的振速信息,因此可以更准确地探测水中目标的角度、位置信息.阵列信号处理指的是将多个阵元按照一定的顺序摆放在空间中的不同位置组成一个阵列,然后使用阵列对声场中的信号进行接收和处理.均匀线性阵列指的是等间距线性排列的多个阵元构成的阵列.只考虑水平方位角,矢量水听器均匀线性阵列接收信号模型如图1所示.

图1 矢量水听器阵列信号接收模型

在图1中:θ是入射信源的方位角,ux,uy是沿x轴和y轴的矢量分量,每个矢量水听器由一个声压传感器和两个振速传感器组成,且在空间中位于同一点,t时刻声压和振速有如下关系[16]:

(1)

(2)

因此单个矢量水听器的输出为:

z(t)=[p(t),vx(t),vy(t)]T

(3)

下面考虑矢量水听器阵列接收信号模型,假设K个远场窄带信号s1(t),s2(t),…,sK(t)入射到由M个间距为d的矢量水听器阵元构成的均匀线性阵列中,假设在t时刻第k个信号的方位角为θk,k=1,2,…,K,将阵列中的第一个阵元作为参考阵元,下面计算t时刻第m个阵元的输出.

对于处于远场的信号,同一个信号到达不同阵元时会产生波程差,这个波程差又会导致各个接收阵元间产生相位差,利用各阵元间的相位差就可以得到信号的方位信息.如图2所示,考虑两个阵元间的波程差,假设c为水下声速,θ为远场信号入射角度,ψ为阵元间的相位延迟,则两阵元接收信号波程差为:

图2 波程差示意图

(4)

从而可得相位差为:

(5)

式(5)中:f0是中心频率,λ是信号波长.因为窄带信号中心频率远大于其带宽,所以相位差可表示为:

(6)

于是可得t时刻第m个阵元的输出为:

(7)

uk=[1,cosθk,sinθk]T,k=1,2,…,K

(8)

(9)

将阵列输出写为矩阵形式即为:

X(t)=A(θ)S(t)+N(t)

(10)

1.2 压缩感知理论

压缩感知理论的核心思想是:如果信号是可压缩的或者在某个变换域是稀疏的,那么就能够通过某种重构算法从部分样本中重构信号.该理论主要分为两个部分:感知测量和重构恢复.

假设空间中的一个有限长的离散时间信号x=[x1,x2,…,xN]T,可以看作RN空间中的N×1维的列向量,x可以用一组线性组合的标准正交基来表示,即:

(11)

为方便表示和描述,记:Ψ=[ψ1,ψ2,…,ψN],a=[a1,a2,…,aN]T.Ψ是N×N的正交基矩阵,a是信号x在标准正交基ψ1,ψ2,…,ψN下的线性表示系数向量.若其系数向量a是稀疏的,即‖a‖0=K≪N,则称信号x在Ψ域上是K-稀疏的.

选择合适的M×N(M≪N)维测量矩阵[17]Ф,对稀疏信号x进行线性投影测量,得到测量值y,即:

y=Фx=ФΨα

(12)

定义矩阵Θ=ΦΨ,Θ是一个M×N维的矩阵,用该矩阵来表示推广之后的观测矩阵,称为感知矩阵(Sensing Matrix).测量值y可以表示为:

y=Θα.

(13)

式(13)中:y表示的是信号x在感知矩阵Θ下的线性投影.

Candes等[18]证明,如果感知矩阵Θ满足有限等距性质(Restrited Isomely Property,RIP),那么就可以利用重构算法通过求解以下l0-范数的优化问题得到稀疏信号:

min‖a‖0s.t.y=Θ·a

(14)

1.3 基于压缩感知的多快拍DOA估计模型

假设空间中有K个远场窄带信号入射到M个阵元组成的均匀线阵上,第k个信号入射角度为θk,一段时间内阵列接收到的L个快拍数据矢量可以表示为:

(15)

同时,将原有的方向矩阵A(θ)扩充为包含划分好的所有方向角的3M×N阶超完备冗余字典G(θ),即:

G(θ)=[a(θ1),a(θ2),…,a(θN)]

(16)

综合式(15)和式(16)可以得到压缩感知框架下的多快拍DOA估计模型,即:

(17)

因此在得到水听器接收到的信号后,通过求解式(17)所述模型就可以估计入射信号的角度,而式(17)是一个欠定方程,无法直接求解,需要转化为最小化l0-范数的优化问题进行求解,即可以转化成如下式子求解:

(18)

通常使用贪婪类算法求解最小化l0-范数的优化问题,如正交匹配追踪算法(OMP)、匹配追踪算法(MP)等.然而贪婪类算法在解决高维问题时效率并不高,往往需要过多的迭代过程和计算步骤,于是为了保证信号的稀疏度和信号重构的精度,一般将式(15)转化为带有不等式约束的优化问题:

(19)

式(19)中:γ2代表一个可能的噪声标准差.于是高度不连续l0-范数问题就变成了l1-范数优化问题,然后可以使用优化算法来求解.求解式(19)可以转化求解如下无约束优化问题:

(20)

式(20)中:λ是一个与噪声有关的正则化参数.

最小化l1-范数优化问题是一个广义的线性规划问题,需要求解二阶锥规划[19],以牺牲计算复杂度来实现准确重构,计算复杂度较高,难以适用于大规模信号快速DOA估计.利用Matlab中的CVX工具箱求解最小化l1-范数优化问题是目前比较常用的方法,但是CVX工具箱的求解速度较慢,不适合实际应用.针对这些问题,本文使用了一种迭代收缩算法来求解压缩感知DOA估计模型中的信号重构问题.

2 基于可分离替代函数的DOA估计方法

2.1 可分离替代函数算法

可分离替代函数(Separable Surrogate Functionals,SSF)算法是迭代收缩算法的一种.考虑函数:

(21)

接下来往式(21)中增加下面一项:

(22)

式(22)中的函数d(·)本质上是x和x0之间的距离度量函数,为了确保它是严格凸的,应该使函数d(·)的黑塞(Hessian)矩阵为正定矩阵,即cI-ATA>0,当c>‖ATA‖2=λmax(ATA)满足该条件.

此时可以得到新的目标函数为:

(23)

式(23)就是算法中要使用的替代函数,将式子中的可变项展开,并重新组织,可以获得如下形式的表示:

(24)

式(24)中的常数项Const1包含了只依赖于b和x0的项.定义:

(25)

(26)

(27)

对于1维优化问题的处理更加容易,于是将G(x)转化成由m个独立并且类型相同的1维优化问题来解决,可以简化总的优化任务,即:

(28)

式(28)中:x[k],v0[k]分别指的是向量x,v0的第k(k=1,2,…,n)个元素,v0是常数项.

(29)

(30)

下面将SSF算法与基于压缩感知的DOA估计模型相结合.为了保证收敛,选择ρ(x)=‖x‖1,即ρ(x[i])=|x[i]|,于是可以得到:

(31)

根据上面的推导,此时一定能够找到公式(31)所表示的问题的最优解.将公式(31)中的向量x替换为矩阵,于是可以将该问题转化为2.3节中得到的基于压缩感知的多快拍DOA估计模型中要解决的凸优化问题,因此可以使用SSF算法的思想来求解公式(20)中的凸优化问题.

2.2 基于可分离替代函数的稀疏重构算法

输入:3M×L维观测信号Y(t),稀疏信号维度为N,稀疏度为K,重构信号的快拍数L,3M×N维传感矩阵G(θ).

(1)初始化:初始设置k=0;

初始残差为r0=Y-Gxk=Y;

(2)主要迭代步骤:k=k+1,

执行以下步骤

步骤1：向后投影,计算e=GTrk-1;

步骤2：:收缩操作,c为参数,利用阈值λ计算eS=Shrink(xk-1+e/c);

步骤3：更新解,计算xk=eS

步骤4：更新残差,rk=Y-Gxk;

2.3 获得DOA估计值

2.4 算法时间复杂度分析

SSF算法时间复杂度的主要影响因素为迭代次数I,每次迭代需要执行O(2×MN)次乘法运算和O(2N)次减法运算,收缩步骤的执行次数为常数级,因此复杂度为O(1).综合每次迭代的复杂度和迭代次数,可以得到SSF算法的复杂度为O(I(2M+2)N).贪婪算法的复杂度通常为O(KMN),其中K迭代次数,SSF算法与贪婪类算法的复杂度处于同一量级,迭代过程中的收缩操作会减少SSF算法的迭代次数,因此可以加快收敛速度.

3 仿真与实验

在本节中,为了全面地验证本文所提算法的性能,选取OMP、MUSIC、CBF、LASSO算法进行对比.仿真实验在Matlab平台进行,实验中均采用阵元间距为半波长的均匀线阵,信号入射角属于[0°,180°],信号中所含噪声为高斯白噪声.使用DOA估计的均方根误差(Root Mean Square Error,RMSE)来评价算法的性能,DOA估计的均方根误差定义为:

(32)

3.1 SSF算法与其它算法DOA估计精度对比

本实验中,矢量均匀线阵的阵元数为20个,阵元间距0.5 m,设置信噪比为-15 dB,入射信号个数为3,入射角度分别为30.5°、60.1°、120.5°,信号频率分别为1 000 Hz、1 200 Hz、1 400 Hz,信号快拍数为300个.

在图3中,绘制出了每个算法得到的系数解向量的图像,峰值越高则其对应的角度为入射角度的可能性越大.可以看出只有需要预先设置入射信号个数的OMP算法没有伪峰,然而实际应用中不可能预先知道入射信号个数,所以通常OMP算法需要先假设一个信号个数,最后再通过算法逼近真实的个数;尽管其余算法都有一定伪峰存在,但其主要的峰值都远大于伪峰峰值,这些主峰对应的横坐标是入射信号的角度,因此很容易就能确定入射信号个数和角度.

图3 各算法估计结果

表1展示了各算法的DOA估计值,可以看出本文所提出的SSF算法求得的主峰最精确,准确估计出了两个信号的入射角度,其余算法只有OMP算法准确估计出了一个信号的入射角度,实验结果可以证明SSF算法具有更高的DOA估计精度.

表1 各算法估计角度

3.2 不同快拍数下各算法DOA估计误差比较

本实验中,设置矢量均匀线阵的阵元数为20个,阵元间距为0.5 m,信噪比为0 dB,入射信源快拍数从5增加到305个,步长为30个快拍,在每个快拍数下进行100次蒙特卡洛实验.图3展示了1个、3个、5个信号源条件下SSF算法、OMP算法、MUSIC算法、CBF算法和LASSO算法DOA估计的均方根误差随快拍数变化的曲线,1个信源条件下入射角为10°,3个信源条件下入射角度为10°、30°和80°,5个信源条件下入射角度为10°、30°、60°、80°和150°.实验结果见图4所示.

图4 不同信源数下均方根误差随快拍数变化曲线

由图4(a)可以看出,单信源条件下五种算法的DOA估计性能都很好,LASSO算法和SSF算法性能明显强于其余算法,本文提出的SSF算法最优;从图4(b)可以看出,各算法的性能相较单信源条件都出现了下降,CBF算法的性能下降最多;从图4(c)可以看出,在5个信源的情况下五种算法DOA估计值的均方根误差都有增大,各算法稳定性也有所下降,这是因为随着信源增多会出现某一个信源入射角度估计出现较大误差,从而导致均方误根差的增大.综合来看信噪比为0 dB,单信源和多信源条件下,本文提出的SSF算法的DOA估计性能都是最佳的,性能随信源数增多影响最小,并且随着信号源快拍数增加DOA估计性能越来越强.

3.3 不同信噪比下各算法DOA估计误差比较

本实验中,将矢量均匀线阵设置为20个,阵元间距为0.5 m.设置快拍数为10个,信噪比以2 dB为步长从-10 dB增长到20 dB.图5展示了1个、3个、5个信号源条件下SSF算法、OMP算法、MUSIC算法、CBF算法和LASSO算法DOA估计的均方根误差随信噪比变化的曲线,1个信号源条件下入射角度为10°,3个信源条件下入射角度为10°、30°和80°,5个信号源条件下入射角度为10°、30°、60°、80°和150°.

图5 不同信源数均方根误差随信噪比变化曲线

从图5(a)可以看出,在单信源条件下MUSIC、CBF这两种传统算法和OMP算法的均方根误差接近,LASSO算法在信噪比大于0 dB时均方根误差小于MUSIC、OMP、CBF算法,在信噪比为-10 dB～20 dB范围内本文提出的SSF算法性能最佳.从图5(a)、(b)、(c)可以看出,随着信号源增加五种算法均方根误差都有所增长,CBF算法在信源增多时均方根误差明显增大,MUSIC算法和OMP算法均方根误差接近,LASSO算法在信噪比大于0 dB时表现好于前三个算法,但是LASSO算法的性能很不稳定,受信噪比影响最大.在信噪比为-10～20 dB范围内SSF算法的均方根误差都小于其他算法,并且SSF算法的DOA估计性能稳定.随着信噪比增大五种算法的DOA估计值的均方根误差都在逐步减少,信源数增多会造成算法性能下降.综合来看在10个快拍条件下,SSF算法在单信源和多信源条件下DOA估计性能都强于其他算法,具有更好的鲁棒性.

3.4 各算法在不同信噪比下DOA估计成功率

首先规定DOA估计成功的定义:DOA估计角度与真实角度相差在±0.5°以内视为估计成功.假设设置入射信号源个数为3,角度分别为10°、60°和120°,快拍数为10个,矢量水听器阵元个数为20个.实验中信噪比从-10 dB到20 dB以2 dB为步长变化,为了保证实验的准确性每个信噪比条件下做100次蒙特卡罗实验,实验结果如图6所示.

图6 各算法随信噪比变化的成功率比较

从图6可以看出,随着信噪比增大各算法的DOA估计成功率都均有所提升,在相同信噪比下本文提出的SSF算法的成功率都高于其它四个算法,在信噪比为10 dB时就达到了100%的成功率,与其余算法相比表现更优.

3.5 算法执行效率对比

上述实验中所对比的算法中OMP和LASSO算法与本文所提SSF算法为同类算法,因此主要对比SSF、OMP、LASSO算法的执行效率.本实验从DOA估计误差和运行时间两个方面综合考虑,实验条件为:信噪比0 dB,信源数为5,快拍数分别为10、30、50、100、200、300,每个快拍数下重复实验10次,计算各算法的DOA估计均方误差和算法平均执行时间.本实验在Windows10操作系统下进行,处理器信息为:Intel(R) Core(TM) i5-7360U CPU @ 2.30GHz,RAM为8.00 GB,仿真软件为Matlab R2019b.实验结果如表2所示.

表2 算法均方根误差及执行时间

由表2可以看出,在各个快拍数下SSF算法的均方根误差均最低;SSF算法和OMP算法的平均执行时间接近,SSF算法执行速度稍快,而LASSO算法由于使用了Matlab中的CVX工具箱,执行效率较低,耗费时间更长,随着信号快拍数增加耗费时间大幅增加.综合来看,本文提出的SSF算法的执行效率高于同类算法.

4 结论

针对低信噪比、多信源情况下DOA估计不准的问题,本文在稀疏重构和压缩感知理论的基础上,提出了一种基于可分离替代函数算法的DOA估计算法,该算法在稀疏重构问题中加入了与上一步稀疏解之间的距离度量函数,通过求解修改后的问题得到稀疏重构问题的解,并且在迭代过程中加入了对解向量的收缩操作,加快了算法的收敛速度,算法中还提出了多快拍的DOA估计模型,进一步提高了估计精度.该算法的优势和特点为:

(1)算法中使用了多快拍的DOA估计模型,因此在低信噪比、少量信源、多信源等条件下该算法都具有很高的DOA估计精度,并且该算法不涉及对协方差矩阵进行空间分解,拥有更强的鲁棒性.

(2)该算法在DOA估计时无需预先设定信源个数,有更高的实际应用价值.

(3)得益于迭代过程中的收缩操作,该算法较同类算法有更快的DOA估计速度.