基于二阶统计特性的方向向量估计算法的DOA估计

侯 进 盛尧宝 张 波

(西南交通大学信息科学与技术学院智能感知智慧运维实验室 成都 611756)

1 引言

作为阵列信号处理领域的核心研究内容之一，波达方向(Direction Of Arrival, DOA)估计已被广泛应用到了各个领域[1-3]。在DOA估计的实际应用场景中，由于天线阵列不可避免地存在内部干扰、阵元互耦等因素的影响，天线阵列对信号的幅值和相位的实际响应往往会偏离理论响应，即天线的实际阵列流形与理论阵列流形往往会存在差别[4]。多信号分类算法(MUltiple SIgnal Classification, MUSIC)[5]、旋转不变子空间算法(Estimating Signal Parameter via Rotational Invariance Techiques,ESPRIT)[6]、确定性极大似然估计算法(Deterministic Maximum likelihood, DML)[7]、正交匹配追踪算法(Orthogonal Matching Pursuit, OMP)[8]等经典空间谱算法都需要用到天线的阵列流形。在阵列流形存在误差时，上述空间谱算法的谱峰的峰宽比较大，DOA估计结果容易出现较大的误差，尤其在两个信号之间的夹角较小时，两个信号的谱峰很容易相互影响造成DOA估计出错。因此，在阵列流形存在误差的应用背景下，如何提升DOA估计算法的性能是值得研究的问题。

针对这一问题，国内外很多学者采取的方法是互耦矫正，这类方法的研究热点是在线自校正方法[9,10]，这类方法将互耦系数、DOA等参数转化为多参数的联合估计问题，存在计算量庞大、不易工程实施等缺点。此外，另一种有效的方法是先采用盲源分离算法估计出各个单信号的实际方向向量，然后再以实际方向向量为基础使用空间谱算法估计各个单信号的DOA[11,12]。传统盲源分离算法主要以分离后的信号的互信息最小、非高斯性最强等作为分离准则，常见的经典算法有快速不动点算法(Fast Fixed-point Algorithm, FastICA)[13]和特征矩阵联合相似对角化算法(Joint Approximative Diagonalization of Eigenmatrix, JADE)[14]。

在DOA估计的实际应用场景中，由于通道接收机体积庞大和价格昂贵等原因，测向设备一般为少通道测向设备，即测向设备的通道接收机数小于天线阵列阵元数的测向设备[15]，例如单通道九阵元圆阵测向设备、三通道九阵元圆阵测向设备。当测向设备为少通道测向设备时，通道接收机会以时间片轮转的方式在不同的阵元处接收信号数据，由于传统盲源分离算法需要使用信号的高阶统计信息，因此传统盲源分离算法仅能获得同一个时间片中的某几个阵元的实际方向矩阵，而不能获得阵列中所有阵元实际方向矩阵，因此基于传统盲源分离算法的DOA估计算法不能应用于少通道测向设备。因此，探索一种适用于少通道测向设备的基于方向向量估计算法的DOA估计算法有重要意义。

针对以上问题，本文提出了一种基于2阶统计特性的方向向量估计算法的DOA估计算法。首先，根据DML估计算法谱函数的特征，构造关于协方差矩阵的酉约束下的优化问题；然后，采用酉约束下的粒子群算法和酉约束下的梯度下降算法对问题进行优化，求得所需酉矩阵；随后，使用所求酉矩阵求得各个单信号的实际方向向量；最后，将各个单信号的实际方向向量输入到空间谱算法中实现DOA估计。由于将多信号的DOA估计转化为多个单信号的DOA估计，因此在天线阵列流形存在误差时，该算法比传统的DOA方法具有更好的DOA估计性能。由于所提算法仅需使用协方差矩阵，因此所提算法可应用于少通道测向设备。由仿真实验结果可知，在阵列流形存在误差以及测向设备为少通道测向设备时，与传统DOA方法相比，所提算法的DOA估计的准确度、抗扰度以及分辨率更高。

2 系统模型

2.1 天线阵列模型

假设空间中均匀圆阵天线阵列(Uniform Circular Array, UCA)的阵元数为M，通道接收机数为m。假设存在多个信号在空间传播，它们均为窄带远场信号，且与天线阵列处于同一平面。均匀圆阵示意图如图1所示。

图1 均匀圆阵示意图

当天线阵列对信号幅值和相位的响应不偏离理论值时，若从远场发射而来的独立信源数为N，在t时刻，阵列接收数据可表示为

其中，s(t)=[s1(t),s2(t),...,sN(t)]T为未知源信号，x(t)=[x1(t),x2(t),...,xM(t)]T为阵列接收数据，M为天线阵列阵元数，n(t) 为功率为σ2的加性高斯白噪声，A=[a(θ1),a(θ2),...,a(θN)]为天线阵列的理论方向矩阵，第n个信源的理论方向向量a(θn)定义为

其中，k=1,2,...,M，k为阵元号，θn表示第n个独立信源的波达方向角，λ是信号的波长，r为圆阵的半径。

信号的DOA相对于整个估计空间来说是符合空域稀疏的，因此可对DOA估计空间进行网格划分。在实际应用中，可以以1°为间隔进行网格划分。改变方向角θ，让方向向量a(θ)在2维空间按方向角从小到大进行扫描，所有网格间隔方向的方向向量a(θ)构成的集合D为理论阵列流形。

其中，Θ={1◦,2◦,...,360◦}为方向角的参数空间。

当天线阵列对信号幅值和相位的响应偏离理论值时，阵列接收数据x(t)可表示为

其中，W=[w1,w2,...,wN]为实际方向矩阵，wn为第n个信源的实际方向向量。阵列接收数据的协方差矩阵为

在实际应用中，常用样本协方差矩阵Rˆ来代替协方差矩阵。当M= m时，样本协方差矩阵Rˆ为

其中，L为阵列接收数据的快拍数。

当M > m时，即测向设备为少通道测向设备时，阵列响应输出是m个通道接收机以时间片轮转的方式，按照某一特定的顺序在不同的阵元处接收信号数据，此时整个阵列的样本协方差矩阵不能直接由式(6)得到。在空间信号是平稳过程的情况下，可以用矩阵拼凑的方法获取整个阵列的样本协方差矩阵。矩阵拼凑的方法是：在某一轮次的通道轮转中，若有两个通道接收机在阵元i和阵元j处同时接收到了信号数据，则可计算获得样本协方差rij，然后将rij填在M行M列的矩阵的第i行、第j列，经过多次时间片轮转，可将1个M行M列的矩阵填满，填满后的矩阵即为整个阵列的样本协方差矩阵Rˆ。

2.2 DML伪谱

引理1在上述数据模型中，假设噪声为平稳高斯零均值随机白噪声，且是循环对称的，其实部和虚部为同一分布，并且有1个反对称的互协方差。在上述统计假设下，信号波达方向角的最大似然估计是下列最大化问题的解[7]

因此，当只有1个信号入射到阵列时，信号波达方向角最大似然估计的解等价于式(8)的解。 证毕

定理3当只有1个信号入射到阵列，其实际方向向量为w时，定义DML伪谱的谱函数为

证明根据定理2，单个信号的DML谱峰搜索函数为PDML(θ)=(aH(θ)a(θ))-1aH(θ)Rˆa(θ)，根据式(5)，可得

其中， E{s(t)s∗(t)}的大小仅与源信号本身幅值大小有关，而与θ和w无关，σ2为噪声功率，也与θ和w无关。所以在实际应用中，可忽略 E{s(t)s∗(t)}和σ2，得到式(10)。 证毕

DML伪谱的谱函数是由DML算法所推导得到的，其反映的是信号实际方向向量与理论方向向量的相关关系。当只有1个信号入射到阵列时，该信号的实际方向向量与其波达方向所对应的理论方向向量的相关性最大，而与其他方向所对应的理论方向向量相关性较小，因此DML伪谱只在信号的波达方向处有1个大的峰值。

3 基于2阶统计特性的方向向量估计算法的DOA估计

3.1 理论推导

假设有N个不同方向的独立未知源信号入射到天线阵列。阵列接收数据的样本协方差矩阵Rˆ的特征值分解公式为

其中，ξi表示Rˆ 的 第i大的特征值，vi为其对应的特征向量。选取前N个大特征值所对应的特征向量构成信号子空间。结合式(5)可知，实际方向向量张成的子空间与信号子空间是同一个子空间，所以可以用N个大特征值对应的特征向量来初始化实际方向向量。在实际应用中，可使用改进的最小长度算法(MMDL, Modified Minimum Description Length)[16]，实现对独立信源数N的有效估计。

令实际方向矩阵W的初始化值为W0，W0=[v1,v2,...,vN]。W0的列向量张成的子空间与W的列向量张成的子空间是同一个子空间，故W0到W的过渡矩阵是1个酉矩阵U，即W=W0U。此时，求实际方向矩阵W的问题转化为求合适的酉矩阵U的问题。

从DML的性质可知，当只有1个信号入射到阵列时，DML伪谱只在信号的波达方向处有1个较大的峰值。所以当有多个不同方向的信号入射到阵列时，求酉矩阵U，使得W的每个列向量代入式(10)后所得到的每个DML伪谱都只有1个较大的峰值时，所求酉矩阵U满足要求。

当J(U) 取得最小值时，W的每个列向量wn代入式(10)后所得到的每个DML伪谱都只有1个较大的峰值，即此时酉矩阵U满足要求。

3.2 酉约束下的粒子群算法

粒子群算法是种群中的粒子以某一速度从当前位置进行迭代寻优的一种全局优化算法[17]。需要注意的是：在粒子迭代寻优的过程中，当粒子的位置不为酉矩阵时，需要将其修正为酉矩阵。式(13)中的待求酉矩阵属于复数Stiefel流形矩阵，将矩阵映射到距离其欧几里得距离最小的复数Stiefel流形矩阵的方法是 π变换[18]。对粒子的位置X做变换 π的公式为

其中，V1和V2为对X做奇异值分解后所得到的酉矩阵，即V1ΣV2=svd(X)。

3.3 酉约束下的梯度下降算法

待求酉矩阵U可以看作损失函数J(U)的最佳权系数，可采用梯度下降算法来求解最佳酉矩阵[18-20]。梯度下降算法的基本迭代步骤是权系数U(t+1)等于权系数U(t) 加上负梯度Z(t)的比例项，即

算法1 酉约束下的粒子群算法

其中，γ为收敛步长，用来控制算法的收敛速度与稳定性。在迭代过程中，当权系数不为复数Stiefel流形矩阵时，可用 π变换将其修正[18]，故上式变换为

在待求酉矩阵属于复数Stiefel流形矩阵的约束条件下，损失函数式(13)的负梯度Z的计算公式为

计算出负梯度Z后，需要选择1个正的步长γ。根据Armijo步长选择准则，步长γ应该满足下列不等式

其中， 表示负梯度Z与其本身的内积，计算公式为

其中，I为单位阵。求合适γ的一种简单的方法是对γ连续加倍直到不等式(23)不成立，然后对γ连续减半直到不等式(22)成立。

3.4 基于2阶统计特性的方向向量估计算法的DOA估计算法

求出酉矩阵U后，令W=W0U，即可得到实际方向矩阵W。实际方向矩阵W的每一列wn为单个信号的实际方向向量，将其输入到式(10)中可得到的DML伪谱，通过寻找峰值即可得到波达方向的估计值。

算法2 酉约束下的梯度下降算法

4 仿真分析

仿真1酉约束下的粒子群算法以及酉约束下的梯度下降算法收敛速度测试。设置均匀圆阵的阵元数M = 9，独立信源数N = 3，快拍数L = 1 024，粒子群算法粒子个数设置为50，两个算法的最大迭代次数都设置为200 次。同一台计算机、同一软件平台下运行两个算法，蒙特卡罗实验200 次。

在上述实验中，粒子群算法平均收敛次数较少，平均收敛时间却较长，梯度下降算法平均收敛次数较多，平均收敛时间却较短，这是因为粒子群算法的每一次优化中需要优化全部粒子的位置和速度。粒子群算法和梯度下降算法都可以得到正确的实验结果，粒子群算法原理较简单，但验证了算法的正确性，而梯度下降算法则提升了算法的运行速度，两种算法收敛速度对比如表1。

仿真2在阵列流形存在误差以及测向设备为少通道测向设备时，不同算法的空间谱测试。设置均匀圆阵的阵元数M = 9，通道接收机数m = 3，半径0.58 m，信号载波频率500 MHz，独立信源数N = 3，正确的DOA分别为[120°]、[140°]、[160°]，快拍数L = 1 024，向天线的理论阵列流形加20 dB的高斯白噪声模拟构造天线的实际阵列流形。

在上述实验条件下，对比测试MUSIC[5]算法、OMP[8]算法以及本文算法的空间谱。3种算法的空间谱图如图2、图3、图4所示。

图2 MUSIC算法空间谱

图3 OMP算法空间谱

图4 本文算法空间谱

3种算法的DOA估计结果如表2所示。

表2 3种算法DOA估计结果

如图2、图3、图4以及表2所示，在上述实验条件下，本文算法基本上能完成3个信号的分辨，而MUSIC算法和OMP算法的分辨结果出错较大。MUSIC算法和OMP算法等传统空间谱算法，需要使用天线的阵列流形。在阵列流形存在误差时，传统空间谱算法的谱峰的峰宽会比较大，尤其在两个信号之间的夹角较小时，两个信号的谱峰很容易相互影响造成DOA估计出错。本文算法采用了将多信号的DOA估计转化为多个单信号的DOA估计的方法，故本文算法避免了信号间谱峰相互影响的问题，故在当前实验条件下，与MUSIC算法和OMP算法相比，本文算法的分辨性能更优。

此外，当前的实验条件为三通道九阵元的均匀圆阵，如图4以及表2所示，本文算法能正确完成DOA估计任务。当测向设备为少通道设备时，通道接收机会以时间片轮转的方式在不同的阵元处接收信号数据，由于传统盲源分离算法[13,14]需要使用信号的高阶统计信息，因此传统盲源分离算法仅能获得同一个时间片中的某几个阵元的实际方向矩阵，而不能获得阵列中所有阵元实际方向矩阵，故传统盲源分离算法不能完成此实验条件下的DOA估计任务。相比基于传统盲源分离算法的DOA估计算法，本文算法能应用于少通道测向设备，这是本文算法的一大优势。

仿真3 测试信噪比对DOA估计准确度的影响，如图5所示。均匀圆阵的阵元数M = 9，通道接收机数m = 3，半径0.58 m，信号载波频率500 MHz，独立信源数N = 2，正确的DOA分别为[113°]、[130°]，快拍数L = 1 024，蒙特卡罗实验200 次。在上述实验条件下，测试MUSIC算法、OMP算法、DML[7]算法以及本文算法的DOA估计均方根误差(RootMeanSquareError,RMSE)，其定义为

图5 不同信噪比下，DOA估计准确度测试

其中，K为蒙特卡罗次数，N为信号源个数，θˆn(k)为第k次蒙特卡罗实验中算法对第n个信号源DOA估计结果，θn为第n个信号源的真实方位角。

在对多个信号进行DOA估计时，本文算法采用了将多信号的DOA估计转化为多个单信号的DOA估计的方法，这种将多信号分离成单信号再进行DOA估计的方法会提升DOA估计算法的抗扰度。所以在上述实验条件下，在信噪比较小时，本文算法的DOA估计RMSE低于MUSIC算法、OMP算法和DML算法的DOA估计RMSE。而在信噪比大于–5 dB时，本文算法的DOA估计RMSE比MUSIC算法的DOA估计RMSE高1°至2°。这是由于上述实验未对天线阵列流形添加误差，在天线的实际阵列流形与理论阵列流形不存在差别且噪声干扰较小时，MUSIC算法可以直接使用理论阵列流形以实现较高精度的分辨，而本文算法中估计实际方向向量的优势将得不到体现。

仿真4在阵列流形存在误差时，测试两个信号间夹角大小对DOA估计准确度的影响。设置均匀圆阵的阵元数M = 9，通道接收机数m = 3，半径0.58 m，信号载波频率500 MHz，独立信源数N = 2，正确的D O A 分别为[1 2 0°]、[120◦+∆θ(∆θ=3◦,6◦,...,30◦)]，快拍数L = 1 024，向天线的理论阵列流形加20 dB的高斯白噪声模拟构造天线的实际阵列流形，蒙特卡罗实验200 次。在上述实验条件下，测试MUSIC算法、OMP算法、DML算法以及本文算法的RMSE。

如图6所示，随着两个信号间夹角的变小，本文算法的DOA估计RMSE低于MUSIC算法、OMP算法和DML算法的DOA估计RMSE。在阵列流形存在误差时且两个信号间夹角较小时，本文算法中将多信号的DOA估计转化为多个单信号的DOA估计的方法优势得到体现，本文算法DOA估计准确度高于MUSIC算法、OMP算法和DML算法，也就是说本文算法的DOA估计分辨率更高。

图6 不同夹角下，DOA估计准确度测试

5 结束语

在DOA估计的实际应用场景中，阵列流形往往会存在误差以及测向设备常常为少通道测向设备，对此，本文提出了一种基于2阶统计特性的方向向量估计算法的DOA估计算法。由于将多信号的DOA估计转化为多个单信号的DOA估计，因此在天线阵列流形存在误差时，该算法比传统的DOA方法具有更好的DOA估计性能。此外，由于仅需使用阵列接收数据的协方差矩阵，因此所提算法克服了基于传统盲源分离算法的DOA估计算法不能应用于少通道测向设备的不足。由仿真实验可知，在阵列流形存在误差以及测向设备为少通道测向设备时，与传统DOA方法相比，所提算法的DOA估计准确度、抗扰度以及分辨率更高。