侯进，陈鑫强,3
（1.西南交通大学信息科学与技术学院智能感知智慧运维实验室，四川 成都 611756；
2.西南交通大学综合交通大数据应用技术国家工程实验室，四川 成都 611756；
3.西南交通大学计算机与人工智能学院，四川 成都 611756）
随着通信、雷达等领域的发展，波达方向（DOA,direction of arrival）估计逐渐成为人们研究的热点课题。DOA 估计的基本问题是利用阵列对空间中某一区域内感兴趣的多个信号的空间位置进行估计，它面临的挑战主要集中在以下4 个方面：1) 在现实通信环境中，受高大建筑物和山丘的阻碍，电磁波在传播时将发生反射，从而产生多径效应，这将使部分DOA 估计算法失效；
2) 出于成本等因素的考虑，实际阵列的阵元数往往小于目标信源数，因此，如何在欠定情况下进行DOA 估计成为一个难题；
3) 在某些特定情况下，目标信源的来波方向可能相同，而目前许多DOA 估计算法无法对相同方向角的信源进行测向，因此，“角度兼并”问题成为DOA估计算法的另一挑战；
4) 多快拍数在一定程度上会降低DOA 估计算法的时效性，因此，如何在少快拍数下保证算法的高准确率和高分辨率，也成为DOA 估计算法待解决的一个问题。
目前，根据测向原理的不同，DOA 估计方法可分为3 类：基于特征结构子空间的方法、基于独立成分分析（ICA,independent component analysis）的方法和基于稀疏重构的方法。最早的超分辨率DOA估计方法是著名的多信号分类（MUSIC,multiple signal classification）算法和基于旋转不变技术的信号参数估计（ESPRIT,estimating signal parameter via rotational invariance technique）算法，它们都属于特征结构子空间方法，由于相干信源会使空间协方差矩阵降秩，因此传统的MUSIC 算法和ESPRIT 算法无法对相干信源进行测向。文献[1-2]和文献[3-4]分别提出了改进的MUSIC 算法和改进的ESPRIT算法用于对相干信号进行测向，但在欠定情况下该类算法性能将大幅度下降。受均匀线阵方向矢量结构的启发，Friedlander 等[5]提出了一种将均匀圆阵转换为虚拟均匀线阵的方法，并结合均匀线阵的解相干算法，如空间平滑算法[6-7]和矩阵重构算法[8-10]等，可以实现均匀圆阵对相干信号的测向，但在模式转换的过程中，噪声会进一步累积，因此该类算法在低信噪比（SNR,signal to noise ratio）下的性能较差。Ziskind 等[11]提出了一种基于确定性极大似然（DML,determined maximum likelihood）估计思想的DOA 估计方法，能够实现均匀圆阵对相干信号的测向，但在求解过程中需要进行高维搜索，因此计算量较大，且不适合欠定情况下的DOA 估计。Zamani等[12]提出了一种改进的ICA 方法，该方法可以对混合信源进行盲分离，但在欠定情况和信源中存在相干信号时将失效。基于稀疏重构的方法可以在低信噪比和信源中存在相干信号的情况下具有良好的性能，如Wang 等[13]和Trinh-Hoang 等[14]提出的正交匹配追踪（OMP,orthogonal matching pursuit）算法及其改进算法、Soubies 等[15]和Zheng 等[16]提出的基于L2 范数的方法和基于稀疏贝叶斯的方法等，但在信源数较多和入射信源示向度间隔较小的情况下，该类算法的测向性能会有所下降。
针对现有算法的不足，本文提出了一种基于几何序列分解和稀疏重构的方法用于DOA 估计。通过分析可以发现，阵列的等时间间隔采样使接收序列成为多个复数等比序列的叠加，各个等比序列的公比蕴含了各个信源的频率信息，而首项则蕴含了信源的幅值和相位信息。因此，本文首先利用几何序列分解算法对混合序列中各个序列的公比和首项进行估计，估计的各个公比可用于信号重建，而估计的首项则可以构建各个信源的实际方向向量，即相干组标记。然后利用稀疏重构方法，即可实现对每个相干组进行DOA 估计。
1.1 均匀圆阵模型
假设空间中均匀圆阵列的阵元数为M，半径为r，阵元1 的空间位置为 [r,0]T，信号的来波方向相对于x轴的逆时针旋转方向进行定义，如图1 所示。
图1 均匀圆阵示意
若从远场发射而来的复数独立信源数为k，则在t时刻，阵列的接收信号可表示为
其中，m=1,…,M，θi表示第i个独立信源的波达角，λ表示信号波长。
1.2 相干信号模型
对于阵列接收的数据，可以利用相关系数来衡量接收信号之间的相关程度，其定义为
式(4)表明，当2 个信号相干时，它们之间只相差一个常复数。现假设第i个独立信源中含有Pi个相干信号，则第i个独立信源的接收数据可表示为
其中，a(θip) ∈CM、Bip和φip分别表示在第i个独立信源中第p个相干信号的理论方向向量、幅值和相位。
令
此时，第i个独立信源的接收数据可表示为
1.3 相干组标记
1.4 网格划分
根据阵列相关理论可知，实际的阵列结构要求方向向量a(θ)与方向角一一对应，不能出现模糊现象，而入射信源的方向角相对于整个测向平面是稀疏的，因此，整个二维测向平面可以进行网格划分，以构成方向角参数空间。在实际应用中，通常以1°为间隔进行网格划分，改变方向角，使方向向量a(θ)在二维空间按方向角从小到大进行扫描，所形成的曲面称为阵列流形，所构成的矩阵称为过完备矩阵D∈CM×360，即
其中，Θ={1°,…,360°}表示方向角的参数空间。
现考虑如下情况：假设均匀圆阵列的阵元数为M，入射到阵列的独立信源数为k，第i个独立信源中含有Pi个相干信号，信号接收机的采样频率为fs，快拍数为L，并假设采样是从1 时刻开始，则阵列的接收数据可表示为
观察式(13)可以发现，均匀圆阵第q个阵元的接收数据（观测矩阵X的第q行数据）本质上是k个等比序列的叠加，其中，第i个等比序列的首项为，公比为，而首项和公比分别包含了第i个独立信源的波达角信息和频率信息。因此，若能估计出混合序列的k个公比，然后利用一个简单的矩阵除法，就能将k个独立信源的实际方向向量估计出来，从而用于DOA 估计。
Lee 等[17]提出了一种基于k维单纯形转换的几何序列分解（GSD-ST,geometric sequence decomposition withk-simplexes transform）方法，可用于对叠加几何序列的公比进行估计。k维单纯形是指包含了k+1 个节点的凸多面体，例如，一维单纯形就是包含了2 个节点的线段，二维单纯形就是三角形，三维单纯形就是四面体，以此类推。利用GSD-ST方法估计公比的思路主要分为4 个步骤：构造搜索空间；
构造k维单纯形序列；
构造联合单纯形；
构造多项式方程。而多项式方程的根即待求的各个序列公比的估计。
在构造k维单纯形之前，需要对原始的叠加序列按照一定方式进行采样，以构成k维单纯形的顶点集合，也称为搜索空间，其构造规则按以下方式进行。
给定自然数ic、k和φk，以φk为起始坐标、i Ick为索引变换按字典顺序对混合序列xq进行连续采样，所构成的集合称为k维搜索空间，其中，I k表示k维的全1 向量，xq表示第q个阵元接收的混合序列。例如，若ic=1，k=2，φk= [1,5]，则所构成的k维搜索空间为
由k维单纯形的定义可知，一个k维单纯形具有k+1 个顶点，现假设其中一个顶点为坐标原点，剩余k个顶点由搜索空间中的k个连续顶点构成，其定义为
需要指出的是，在构造本文要求的k维单纯形时，从搜索空间中选取的k个顶点必须是连续的，且在构造k维单纯形序列时，坐标j的取值也需要满足连续递增的条件。
接着，从构造的k维单纯形序列中挑选出任意2 个连续的单纯形，将其坐标按顺序排列，以构成一个具有k+2 个顶点的联合单纯形C∈Ck×(k+2)，定义为
从式(17)可以看出，联合单纯形C的k+1 个非原点坐标实际上由搜索空间的k+1 个连续顶点坐标构成；
然后，对这k+1 个连续顶点按组合的顺序进行采样，以构成一个含有k+1 个新的k维单纯形的序列。例如，若ic= 1,k= 2,φk= [1,5],j=1，则所构成的联合单纯形C和含有k+1 个新的k维单纯形的序列为
最后，需要构造一个k阶多项式，并且这个k阶多项式的k+1 个系数由序列中单纯形的体积通过运算得到。其中，k维单纯形的体积定义为
将k维单纯形序列的体积运算规则定义为
其中，表示公比的估计。
Lee 等[17]证明，给定混合序列xq和独立信源数k，当采样间隔ic=1时，对于任意的自然数φk和j，k阶多项式的系数集合是唯一的，并且由这个系数集合构成的k阶多项式的k个根即混合序列中k个公比的估计。由于系数集合在上述采样方式下与常数φk和j无关，因此，可以把系数集合记为v(k,xq)。需要注意的是，为了使构造的多项式是人们期望的结构，也就是要保证这个多项式的唯一性，在构造搜索空间时，一定要以1×Ik为索引变换对混合序列进行采样，即采样间隔一定要等于1，而采样起始坐标φk和联合多面体C的索引值j可以任取。同时，Lee 等[17]还证明，给定混合序列xq和独立信源数k，获取的最小样本集合由混合序列xq的2k个连续采样点构成，即在估计k个独立信源的实际方向向量时，所需要的最小快拍数L为2k，这也使该算法可以在极少快拍数下对DOA 进行有效估计。
通过以上分析可以发现，在对混合序列中各个序列的公比进行估计时，需要预先知道混合序列中的序列个数，也就是独立信源的个数k。Lee等[17]还提出了一种利用单纯形对独立信源数k进行估计的方法，但该方法需要进行多维搜索，计算量较大，因此会大幅度降低DOA 估计的时效性。Bazzi 等[18]提出了一种改进的最小描述长度（MDL,minimum description length）准则，可以实现在信噪比较低和快拍数较少的情况下对独立信源数k进行有效估计，且算法时效性较好，因此，可以使用该方法对k进行估计。
受现实扰动的影响，阵列在接收数据时会产生一定的噪声，因此，为了提升DOA 估计的性能，需要对接收数据进行去噪。Hansen[19]提出了一种基于k截断式奇异值分解（SVD,singular value decomposition）的去噪方法，该方法的主要思想是通过迭代的方法，不断提取原序列的k个主特征值和它们对应的特征向量，直到重构的序列收敛，从而达到滤除噪声的目的。该方法的去噪过程主要分为以下5 个步骤，这里假设去噪前和去噪后的序列分别为xq和，序列的索引值从1 开始，且快拍数为L。
其中，矩阵U和矩阵V分别表示左奇异矩阵和右奇异矩阵，矩阵S为对角矩阵。
步骤3利用矩阵Q的k个主特征值和对应的特征向量重构矩阵Q
其中，U[:,1:k]表示由矩阵U的所有行和前k列构成的矩阵。
步骤4通过平均值的方法，将矩阵Q转换为序列
其中，当l≤m时，z=l，否则，z=L-l+1。
步骤5重复以上4 个步骤，直到序列xq收敛到序列，收敛后的序列即滤除噪声后的序列。
以上就是估计k个实际方向向量的整体思路，现将算法流程总结如下。
算法1基于去噪GSD-ST 的实际方向向量估计算法
根据上述分析可知，通过网格划分构造了一个过完备矩阵D，并利用过完备矩阵中的稀疏基对第i个相干组的实际方向向量进行了稀疏表示，将第i个相干组中多径信号的DOA 估计问题转换为空间谱Vi的稀疏重构问题，若对每个相干组都进行稀疏重构，得到k张空间谱(V1,…,Vk)，即可估计出所有入射信源的波达角信息。
在求解第i个相干组的稀疏矩阵Vi时，由于稀疏度Pi事先未知，且实际方向向量的维度远小于矩阵Vi的维度，因此式(26)具有无数个解。根据压缩理论知识可知，如果过完备矩阵D满足约束等距性（RIP,restricted isometry property）条件，即可利用矩阵Vi的稀疏性，将式(26)转换为一个组合优化问题进行求解。现对该条件进行简单验证，由于
在求解式(28)时，可以通过正则化对这个非线性凸优化问题进行求解，其代价函数为
其中，λ为正则化系数。
在实际仿真时，为了高效地找到全局最优解，可以利用二阶锥规划（SOCP,second order cone programming）算法在内点法的框架下对式(29)中的优化表达式进行求解，基于此，可以通过引入辅助向量的方式，将式(29)转换为以下优化表达式进行求解
需要注意的是，正则化系数λ的选取对稀疏重构的效果有着举足轻重的作用，文献[14]给出了正则化系数λ的选取规则，由于本文通过k截断式SVD 算法对原观测矩阵进行了去噪处理，因此系数λ的值不应取得太大，在实际仿真时，系数λ在区间[0.1,3]内进行取值都能达到很好的效果。
以上就是本文所提算法的全部分析过程，现将算法流程总结如下。
算法2基于去噪GSD-ST 和稀疏重构的DOA估计算法
1) 利用算法1 得到独立信源数k和实际方向矩阵的估计
2) 以1°为间隔对空间进行网格划分，构造过完备矩阵D；
5) 将稀疏矩阵Vi的元素值画成空间谱，通过谱峰搜索返回第i个相干组的DOA 估计值；
6)i=i+1，如果i≤k，返回步骤4)，否则算法结束。
由于本文使用的信源数估计算法为改进的MDL 算法，因此，基于GSD-ST 的实际方向向量估计算法最多能够估计出M个独立信源的实际方向向量；
并且，一般情况下入射信源的幅值满足随机分布，对于每一个相干组，稀疏重构算法能够估计出M–1 个相干信号[20]，因此，本文所提算法最多能够估计出M(M–1)个入射信源，即能够实现在欠定情况下对相干信号的DOA 估计；
在第2 节的分析中，本文对相干组进行了拆分，即每个相干组的DOA 估计互不影响，因此，本文所提算法能够解决“角度兼并”问题；
此外，在估计实际方向向量时，所需要的最小快拍数L=2k，因此，本文所提算法可以在极少数快拍下实现DOA 估计。
为了验证本文所提算法的性能，本节将其与非线性快速独立成分分析（FastICA,fast independent component analysis）算法[20]、OMP 算法[13]、MUSIC算法[1]和DML 算法[11]进行对比，并通过5 个实验进行说明。
实验1相干组拆分性能分析
实验条件：均匀圆阵的阵元数M=3，半径r=0.58 m，信号的中心频率f=500 MHz，相干组个数k=3，正确示向度分别为[32°]、[210°]、[89°]，信噪比SNR=[-18,18] dB，快拍数L=200，蒙特卡罗实验200 次。相干组拆分效果通过改进的行元素优势指标[21]进行分析，定义为
图2 为非线性FastICA 算法与本文所提算法的盲分离效果对比。从图2 中可以看出，相比于非线性FastICA 算法，本文所提算法在低信噪比下的p指标系数都小于2，这表明本文所提算法可以在低信噪比下对相干组进行有效拆分，且拆分效果优于非线性FastICA 算法。
图2 2 种算法盲分离效果对比
为了进一步说明本文所提算法拆分相干组的效果，图3 和图4 给出了当信噪比为-18 dB 时相干组拆分前后的数据频谱曲线。
图3 拆分前的数据频谱曲线
图4 拆分后的数据频谱曲线
对比图3 和图4 可以看出，混合数据的频谱中有许多伪峰，而本文所提算法对混合数据进行盲分离后，其频谱曲线滤除了由噪声引起的伪峰，并且能够准确恢复出信号的频率信息，这也验证了算法1的正确性。
实验2 欠定情况下对相干信号的DOA 估计
实验条件：均匀圆阵的阵元数M=3，半径r=0.58 m，信号的中心频率f=500 MHz，相干组个数k=3，相干系数分别为=-0 .416+0.909i、= 0.921+0.389i 和= 0.853 - 0.522i，正确示向度分别为[32°,328°]、[30°,210°]、[89°,272°]，信噪比SNR=10 dB，快拍数L=200。
在上述实验条件下，改进的MDL 算法估计出的独立信源（相干组）数k=3，因此，MUSIC算法和DML 算法只能估计出3 个信号的示向度，为了验证本文所提算法在欠定情况下对相干信号的DOA 估计性能，将其与OMP 算法进行了对比。图5 和表1 分别为OMP 算法与本文所提算法在欠定情况下对相干信号的测向空间谱和对应的测向结果。
图5 2 种算法在欠定情况下对相干信号的测向空间谱
表1 2 种算法在欠定情况下对相干信号的测向结果
从图5 和表1 可以看出，在欠定情况下，OMP算法估计出的DOA 结果相较正确结果的误差最大达到了3°，并且由于30°和32°这2 个示向度较接近，OMP 算法不能有效分辨，从而出现了3°这个错误估计结果；
本文所提算法不仅能够有效拆分出3 个相干组，且DOA 估计结果全部正确，这也验证了本文所提算法所能准确估计的最大信源数为M(M–1)的结论。
实验3极少快拍数下的测向性能
实验条件均匀圆阵的阵元数M=9，半径r=0.58 m，信号的中心频率f=500 MHz，独立信源数k=3，正确示向度分别为[15°]、[167°]和[251°]，信噪比SNR=30 dB，快拍数L=6。
在极少快拍数的情况下，本节对比了OMP 算法、MUSIC 算法和DML 算法，由于DML 算法需要进行高维搜索，因此没有画出DML 算法的空间谱，只记录了对应的测向结果。图6 和表2 分别为对比算法与本文所提算法在极少快拍数下的测向空间谱和对应的测向结果。
图6 3 种算法在极少快拍数下的测向空间谱
表2 4 种算法在极少快拍数下的测向结果
从图6 和表2 可以看出，在极少快拍数的情况下，OMP 算法、MUSIC 算法和DML 算法在测向时会有2°～3°的测向误差，而本文所提算法能够精确地对入射信号进行DOA 估计。在实验过程中还发现，快拍数的减小会使噪声对测向的影响变得很大，因此，极少快拍数测向只适用于信噪比较高的情况，虽然现实测向情况较复杂，但本文所提算法为极少快拍数测向提供了一定的理论和仿真依据。
实验4入射信源角度差对测向成功率的影响
实验条件均匀圆阵的阵元数M=9，半径r=0.58 m，信号的中心频率f=500 MHz，独立信源数k=2，正确示向度分别为 [123°] 和[123°+Δθ(Δθ=0°,…,6°)]，信噪比SNR=10 dB，快拍数L=200，蒙特卡罗实验200 次。
图7 为4 种算法测向成功率随信号角度差变化的曲线。从图7 可以看出，OMP 算法、MUSIC 算法和DML 算法都无法分辨出角度差为0°的2 个入射信号，并且OMP 算法和MUSIC 算法在入射信源夹角超过4°时才能达到很好的分辨效果。本文所提算法由于对相干组进行了拆分，相干组之间的测向互不影响，因此可以对角度差较小的入射信号进行精确分辨。
图7 4 种算法测向成功率随信号角度差变化的曲线
实验5信噪比对测向成功率和测向精度的影响
实验条件均匀圆阵的阵元数M=9，半径r=0.58 m，信号的中心频率f=500 MHz，独立信源数k=3，正确示向度分别为[52°]、[189°]和[297°]，快拍数L=200，蒙特卡罗实验200 次。其中，测向精度采用均方误差根（RMSE,root mean square error）进行定义
其中，Q表示测向正确示向度的次数，即测向误差结果在3°范围内的个数；
表示在第q次测向时对第i个示向度的估计值。
图8 和图9 分别给出了4 种测向算法的成功率和RMSE 随信噪比变化的曲线。从图8 可以看出，当信噪比大于-16 dB 时，本文所提算法的测向成功率可达90%以上，且高于OMP 算法和MUSIC 算法；
当信噪比大于-14 dB 时，本文所提算法的测向成功率可达99%以上，这充分说明了本文所提算法在低信噪比下能够对多个信号进行有效的DOA 估计。从图9 可以看出，当信噪比大于-16 dB 时，本文所提算法的测向精度达到了0.5°以上，且对应的精度曲线收敛到了0°，而OMP算法、MUSIC 算法和DML 算法对应的测向精度曲线分别收敛到了0.8°、1.2°和0.7°。综合来看，相比于其他3 种测向算法，本文所提算法的抗噪能力和测向精度更优。
图8 4 种算法测向成功率随信噪比变化的曲线
图9 4 种算法测向精度随信噪比变化的曲线
本文以提升均匀圆阵对相干信号的测向性能为目标，提出了一种去噪GSD-ST 与稀疏重构相结合的DOA 估计算法，去噪GSD-ST 算法保证了对相干组拆分和估计实际方向矩阵的有效性和稳定性，稀疏重构算法则利用估计的实际方向矩阵对每个相干组进行DOA 估计。仿真实验表明，在欠定情况和入射信源存在相干信号时，本文所提算法的测向准确率和抗扰能力都更优。当入射信源角度差较小和信源数较多时，本文所提算法的测向优势更明显；
同时，本文所提算法为极少快拍数测向提供了理论和仿真依据。接下来，笔者将进一步探索如何在低信噪比下利用极少快拍数进行测向的问题，以进一步提高算法的测向性能。
——以鲁甸地震相关新浪微博为例湖南科技学院学报(2015年8期)2015-08-26