CN112444773A - 基于空域融合的压缩感知二维doa估计方法 - Google Patents

Abstract

本发明公开了一种基于空域融合的压缩感知二维DOA估计方法，属于阵列信号处理技术领域。由于二维DOA估计算法目前存在着计算复杂度过高，估计精度较低的问题，本发明主要在压缩感知模型的基础上研究二维DOA估计方法。首先利用L阵的空间合成角进行降维，采用一种等角度结合等余弦的空域划分方式构建性能更优的阵列流型矩阵，以得到压缩感知框架下的DOA估计模型。再分别采用OMP算法进行角度重构并按照信源幅值大小进行配对，得到目标的二维方向信息。该发明具有估计精度高、解相干、抗噪性能好的特点，促进了压缩感知与DOA估计的进一步融合。

Description

基于空域融合的压缩感知二维DOA估计方法

技术领域

本发明是一种波达方向(DOA)估计算法，常应用于雷达探测、无线通信、地震勘探等任务，可自动完成二维波达方向角度的估计。该发明属于阵列信号处理技术领域。

背景技术作为阵列信号处理领域的关键问题，波达方向估计(DOA)在雷达、通讯、地震等众多领域得到了广泛的应用。由于一维方向估计无法描述空间特性，在许多方面不适用于工程实际，二维波达方向(2D-DOA)估计已经成为了研究热点。在二维信号子空间的估计算法中，最经典的估计算法是二维多重信号分类(two-dimensional multiple signalclassification，2D MUSIC)算法，它对接收数据的协方差矩阵进行特征值分解得到信号子空间和噪声子空间，再根据子空间的正交性构造出空间谱，然而该算法需要二维谱峰搜索因而带来了巨大的运算量，难以满足实际应用。为了解决上述缺陷，又有学者提出了一种无需谱峰搜索的二维旋转不变子空间(2-D ESPRIT)算法，但需要有移不变特性的数据以及特殊的阵列流型，并且需要特征值分解，算法的计算复杂度仍然较高。在此基础上有人提出了无需谱峰搜索以及特征值分解的二维传播算子(2D PM)算法，极大的提高了阵列信号的处理性能，且适用于不同阵列下的多维角度估计。但以上算法都无法突破子空间类算法的限制，在低信噪比、小快拍以及信源空间间距很小的时候通常不能准确的估计出波达方向，实用性受到一定的限制。二维波达方向估计的计算复杂度受到阵列几何结构的影响,与均匀矩形阵列或均匀圆形阵列等其它阵列结构相比，L型阵列结构下的二维DOA估计具有更高的估计精度和较低的计算复杂度，并且结构简单、易于实施，因此被广泛应用于工程应用中。

近年来Donoho等人提出的压缩感知(Comprehensive Sensing，CS)理论为现代信号处理带来了一种更高效，更精确的方法，它利用信号的稀疏性，在远小于奈奎斯特频率的采样率下，通过求解l₀最小范数优化问题得到高精度的重构信号。在实际空域中目标信号相比于整个空间角度是极少的，本身具有稀疏性，因此非常适合对空域进行角度划分，来构造压缩感知下的DOA估计模型进行求解。基于稀疏表示的DOA估计方法有效降低了信号的采样数、数据传输、存储及处理的成本，提高参数的估计性能，并且不受相干信源的影响。因此本发明将在压缩感知的基础上继续探索DOA的估计方法，进一步提高算法的实用性。

发明内容

传统二维DOA估计算法目前存在着计算复杂度过高，估计精度较低的问题，难以满足现代信号处理的实时性及准确性要求。为了解决上述缺陷，我们在压缩感知模型下进行二维DOA估计，并对空域的稀疏划分方式进行了研究，最终提出了一种基于空域融合划分的压缩感知二维DOA估计算法，大大减少了算法计算量，提高了DOA估计精度。

为达到上述目地，本发明包括以下步骤：

第一步：对二维阵列信号模型进行建模仿真，得到各天线阵元上的输出信号与目标信源的直接关系。

第二步：利用L阵的空间几何关系构造空间合成角，将二维角度降维成两个一维空间角进行估计。

第三步：根据目标信号在整个空域中的稀疏性，将信号存在的所有可能的空间角度采用等角度划分结合等余弦划分的方式得到大量离散化的角度方向。在确定了空域的划分方式后，包含信号特征信息的阵列导向矩阵A(θ)也会随之扩展，从而得到基于压缩感知的DOA估计模型。

第四步：利用经典的OMP算法分别求解上述稀疏重构问题得到两组空间合成角的DOA估计结果。

第五步：将两组角度对应的稀疏系数h_x和h_y按幅值大小排序后使其成组配对，最后按照配对结果反解得到各信源的方位角θ和俯仰角

第六步：实验验证本文基于空域融合的压缩感知二维DOA估计方法的有效性，将均方误差作为指标验证算法性能优劣。

与现有技术相比，本发明具有以下有益效果：

(1)利用L阵模型的空间合成角将二维信号降维为两个一维信号，大大减小了阵列流型矩阵的维度，进而有效降低了二维DOA估计的计算复杂度以及空间的复杂度。

(2)利用压缩感知理论构建稀疏重构模型，使得本文算法在少快拍低信噪比的情况下也能够高精度的还原出信号的二维角度信息，并且不涉及矩阵特征值分解，具有天然解相干的优良特性。

(3)通过对等角划分与等弦划分的正交性进行分析，找到各划分方式占优的角度空间，并将两者进行结合构成新的空域划分方式，相比于等角划分和等弦划分具有更高的估计性能。

附图说明

图1为L型阵列接收模型

图2入射信号的空域稀疏化示意图

图3两种划分方式下阵列流型矩阵随网格数变化的正交性曲线

图4方案一下两种划分方式随信噪比变化时的均方根误差曲线

图5方案二下两种划分方式随信噪比变化时的均方根误差曲线

图6 2D-MUSIC算法的DOA估计散点图

图7 2D-PM算法的DOA估计散点图

图8本文算法的DOA估计散点图

图9各算法随信噪比变化时的均方根误差曲线

具体实施方式

以下结合具体实施例，并参照附图，对本发明进一步详细说明。

步骤一：本文采用L型均匀阵列，该模型中含有互相垂直的两个子阵X、Y，且分别位于x轴和y轴上，原点由两子阵共用，每个子阵中均含有M个阵元。假设存在k个远场窄带信号(k＝1，2，...，K)入射到该阵列上，每个阵元间的距离为半个波长，入射信号与x轴和y轴的夹角分别为α，β，入射方向角θ，

分别为信号的方位角和俯仰角，阵列接收信号的示意图如图1所示。在时刻t两子阵接收信号的单快拍矢量模型为：

其中s(t)＝[s₁(t)，s₂(t)…，s_k(t)]^T表示K个信源t时刻信号矢量，x(t)＝[x₁(t)，x₂(t)，…，x_M(t)]^T和y(t)＝[y₁(t)，y₂(t)，…，y_M(t)]^T分别为子阵X、Y的阵列接收向量，

以及

均为方差等于

均值为零的加性高斯白噪声，且该噪声信号与原始接收信号相互独立。M×K阶阵列流型矩阵可表示为：

矩阵中第m行k列的元素表示第m个天线阵元接收到第k个远场信号的增益和相位延迟信息，其中：

公式(5)和公式(6)分别表示第k个信号在子阵X和子阵Y的导向矢量，

λ表示信号的波长，d代表两子阵阵元之间的距离。

步骤二：如果一维DOA估计中线阵的空域角度离散化个数为n，那么相同精度的平面阵同时估计方位角和俯仰角时，其空域角度离散化个数为n²。这会大大提高算法的时间复杂度和空间复杂度，因此本文选择使用L阵降维来解决这个问题。L阵是由两个具有空间关系的线阵组成。

由图1的空间几何位置关系可知：

通过式(7)的空间合成角将二维空间角估计降维成两个一维空间角估计。

将式(7)分别带入式(5)，式(6)得到空间合成角的导向量：

步骤三：式(7)定义的空间合成角α_k和β_k已包含了二维信号的方位角和俯仰角等信息，并且在整个空域中具有稀疏性。因此可以将信号存在的所有可能的空间角度θ进行空域划分得到N个离散化的角度方向，用h_i表示第i个方向的信号，则空域信号可表示为h＝[h₁，h₂，…，h_N]^T，如果空间中有K个目标信源，显然h中这K个对应方向的信号有非零元素，其余方向都为零。目前常见的空间划分方式有等弦划分和等角划分两类，不同的划分方法对应得到的阵列流形矩阵是不一样的，它会直接影响稀疏重构的性能，进而影响DOA估计性能。因此，为了选取一种好的空域划分方法，需要分析两种划分方法下阵列流型矩阵A的RIP(有限等距性)条件。

但是，直接验证给定矩阵的RIP条件是一个NP-hard问题，于是学者Candès给出了RIP条件的等价描述：若从感知矩阵中任意抽取2K列矢量之间几乎是正交的，那么可以保证稀疏信号的精确重构。于是，不同划分方法下阵列流形矩阵中各列之间的正交性大小决定着能否精确恢复信号，进而得到准确的信源DOA估计值。根据向量内积的定义可知，阵列流形矩阵A中任意向量之间的正交性可以用下式表示：

式中θ_p、θ_q表示空域划分后的不同角度值，其中p，q＝1，2，…，N。显然当两者角度最相近时其正交性最大。

将观测空间[-90°，90°]分别按等角和等弦均匀地分成N个角度

假设μ(A₁)表示等角划分下流型矩阵的正交性，则

假设μ(A₂)表示等弦划分下流型矩阵的正交性，则

对两者的正交性进行仿真比较，采用16阵元的均匀线阵，分别按等弦划分与等角划分将空间均匀划分180份。图为两种划分方式的正交性随网格变化的结果图。可以看出，当目标角度在两端时等正弦方法下阵列流形矩阵具有更好的近似正交性，更易满足RIP条件；而当角度在中间范围内时，等角度划分方法下阵列流形矩阵的近似正交性更好，具有更显著的RIP性质。那么必然存在一个临界角度满足等角划分与等弦划分非正交性相同，即：

μ(A₁)＝μ(A₂) (13)

将式(11)式(12)代入式(13)化简得：

cosθ_P-cosθ_P-1＝2/N (14)

由上式可推出

对上式左边进行和差化积，可得

则上式关于θ_P的两个解可分别表示为θ_P1、θ_P2

即在(-90°，θ_P1)和(θ_P2，90°)等弦划分的方式更优，在(θ_P1，θ_P2)等角度划分的方式更优。因此我们可以根据角度范围选择空域的划分方式，再将其对应的阵列流型矩阵进行融合，得到最优的阵列流型矩阵。假设空域网格总数为n，则等角度划分网格数为(θ_P2-θ_P1)*n/180，剩余网格数为等正弦划分网格数。合并后的矩阵可表示为：

于是稀疏信号表示的DOA估计数据模型为：

其中X、Y分别为子阵X、Y的接收信号向量，N_x、N_y分别表示两子阵的噪声信号向量。扩张后的稀疏字典矩阵为：

式(20)式(21)是典型的压缩感知模型，扩展后的流形矩阵

扮演着测量矩阵φ的角色。因为目标信号h_x、h_y相对于整个空域是稀疏的，不需要稀疏基矩阵来稀疏表示，所以

既是观测矩阵也是感知矩阵。到此基于压缩感知的DOA估计模型就构造完了。

步骤三：通过上述估计模型的构造，求解式(24)(25)的稀疏重构问题即可得到两组空间合成角α_k和β_k。即已知观测信号X、Y和感知矩阵

重构出目标稀疏信号h_x、h_y，可表示为：

其中||h_x||₀、||h_y||₀是h_x、h_y的0范数，表示h_x、h_y中非零元素的数量，ε为观测信号中所含噪声的能量大小，式(24)(25)在本质上是具有非凸稀疏约束的优化问题。本发明采用最经典的OMP算法对其进行求解得到两组DOA估计结果。OMP的算法流程如下：

输入：M×N维的测量矩阵

M维的测量值X；稀疏度K

输出；N×1维稀疏信号

K×1维DOA估计结果

(1)初始化。残差r₀＝X，支撑集

支撑矩阵

迭代次数t＝1；

(2)扩充候选集。计算测量矩阵

中各原子与当前残差r_t-1的相关性：

并从测量矩阵中选取和r_t-1相关性最大的一个原子；

(3)更新支撑集。当前支撑集Λ_i由上一次迭代得到的支撑集Λ_t-1和上一步中选出的一个原子下标的并集组成；

(4)更新支撑矩阵。当前支撑矩阵

通过将本次选取的一个原子添加到上一次迭代的支撑矩阵

中得到；

(5)更新信号逼近。求

的最小二乘解：

(6)更新残差值。

(7)判断迭代终止条件，如达到终止条件t≥K则停止迭代并输出信号的K稀疏估计

否则更新迭代次数t＝t+1并转入步骤(2)。

稀疏度和前面信源数的个数应该是相等的，都用K。

步骤四：通过上述重构算法即可求出两个子线阵对应的一组夹角α和β。若来波信号数量为1，那么将夹角α以及夹角β按照式(19)直接换算即可得到信源的二维空间角度，但实际应用中，空间往往存在多个信号源，两个子线阵分别估计出的角度结果如何配对，确定哪两个夹角属于同一个信源就是一个必须解决的问题。

因为h_x和h_y中的非零元素大小对应着信号的幅值，所以可以采用信号幅度排序的方法来实现角度的配对，即对式(15)式(16)估计出的稀疏系数h_x和h_y中的非零元素按幅值大小排序后使其一一对应，这种角度配对方案的前提是各信号间的幅度存在一定的差异。最后将配对好的一组角度按照式(26)反解得到各信源的方位角θ和俯仰角

步骤五：仿真验证本发明的有效性

实验一：验证空域融合算法对DOA估计的有效性，仿真比较不同信噪比下两种角度划分方式下的DOA估计性能。实验设置阵元数为16，离散网格数为360，信噪比从-10dB以2dB为间隔逐渐提升到10dB，信源数为16，按照上述理论分析得到的临界值将测试角度分为两种方案进行设置(1)信源网格编号为160，200，230的角度(2)信源网格编号为10，20，30的角度。本次实验采用均方根误差作为衡量DOA估计算法性能好坏的指标，DOA估计的均方根误差定义为：

其中K为信源数，N设为100是蒙特卡罗实验次数，

为第k个信源第n次试验的角度估计值，θ_k为第k个信源真实的波达角。图4图5分别为方案一和方案二下两种划分方式的均方根误差随信噪比变化的结果。

可以看到，随着信噪比的提升，均方根误差都在逐渐减小，在方案一下等角度划分的均方根误差低于等正弦划分，可以得知在中间角度区域等角度划分性能更好，在方案二下等正弦划分的均方根误差低于等角度划分，可以得知在两端区域等弦划分得到的DOA估计更准确。实验结果表明不同的角度空间所适合的角度划分方式是不一样的，也说明了本文空域融合方法的可行性。

实验二：比较2-D MUSIC算法、2-D PM算法以及本文算法的2维DOA估计性能。假设有3个非相关远场信号，分别从(60°，120°)、(100°，100°)、(140°，80°)的方向入射到阵元总数为16的L阵中。设定仿真试验的快拍数等于100，信噪比为-10dB，本文算法的角度网格为0.5°。分别对经典2-D MUSIC算法、2-D PM算法以及本文算法进行仿真试验，并重复绘制100次散点。得到的仿真结果如图6-图8所示。

比较MUSIC算法和PM算法以及本文算法在噪声环境下的估计性能。本次实验采用均方根误差作为衡量DOA估计算法性能好坏的指标，二维DOA估计的均方根误差定义为：

表示第m次蒙特卡罗仿真实验得到的第K个信源的方位角和俯仰角。(α_k，β_k)为第k个信源真实的角度，K为信源数，N设为100是蒙特卡罗实验次数。仿真中加入高斯白噪声干扰，信噪比从-10dB逐步提升到10dB，其它仿真参数同上个实验。实验对比了2-DMUSIC算法、2-D PM算法以及本文算法的均方根误差随信噪比变化的曲线，仿真结果如图9所示：

由图9可知，随着SNR增大，3种算法的均方根误差逐渐减小，估计精度变高。在SNR处于-10dB到2dB区间时本文算法的均方根误差低于另外两种算法，2db以上时三种算法的均方根误差几乎相同。从整体来看，本文算法的均方根误差都处于较低的位置，估计精度更高，这是由于本文算法利用压缩感知理论得到了更多信号信息，提高了阵列信号的利用率，大大提升了少快拍数低信噪比条件下的估计性能。因此可得，改进后的算法在同样的仿真条件下估计精度更高，且在信噪比较低的情况下相比其它算法具有更明显的优势。

以上实施例仅为本发明的示例性实施例，不用于限制本发明，本发明的保护范围由权利要求书限定。本领域技术人员可以在本发明的实质和保护范围内，对本发明做出各种修改或等同替换，这种修改或等同替换也应视为落在本发明的保护范围内。

Claims (2)

1.一种基于空域融合的压缩感知二维DOA估计方法，其特征在于，具体如下：

第一步：对二维阵列信号模型进行建模仿真，得到各天线阵元上的输出信号与目标信源的直接关系；

第二步：利用L阵的空间几何关系构造空间合成角，将二维角度降维成两个一维空间角进行估计；

第三步：根据目标信号在整个空域中的稀疏性，将信号存在的所有可能的空间角度采用等角度划分结合等余弦划分的方式得到离散化的角度方向；在确定了空域的划分方式后，包含信号特征信息的阵列导向矩阵A(θ)也会随之扩展，从而得到基于压缩感知的DOA估计模型；

第四步：利用OMP算法分别求解上述稀疏重构问题得到两组空间合成角的DOA估计结果；

第五步：将两组角度对应的稀疏系数h_x和h_y按幅值大小排序后使其成组配对，最后按照配对结果反解得到各信源的方位角θ和俯仰角

2.根据权利要求1所述的方法，其特征在于：

步骤一：采用L型均匀阵列，含有互相垂直的两个子阵X、Y，且分别位于x轴和y轴上，原点由两子阵共用，每个子阵中均含有M个阵元；假设存在k个远场窄带信号入射到该阵列上，k＝1，2，...，K表示K个信源；每个阵元间的距离为半个波长，入射信号与x轴和y轴的夹角分别为α，β，入射方向角θ，

分别为信号的方位角和俯仰角，在时刻t两子阵接收信号的单快拍矢量模型为：

其中s(t)＝[s₁(t)，s₂(t)…，s_k(t)]T表示K个信源t时刻信号矢量，x(t)＝[x₁(t)，x₂(t)，...，x_M(t)]^T和y(t)＝[y₁(t)，y₂(t)，...，y_M(t)]^T分别为子阵X、Y的阵列接收向量，

以及

均为方差等于

均值为零的加性高斯白噪声，且该噪声信号与原始接收信号相互独立；M×K阶阵列流型矩阵表示为：

矩阵中第m行k列的元素表示第m个天线阵元接收到第k个远场信号的增益和相位延迟信息，其中：

公式(5)和公式(6)分别表示第k个信号在子阵X和子阵Y的导向矢量，

λ表示信号的波长，d代表两子阵阵元之间的距离；

步骤二：如果一维DOA估计中线阵的空域角度离散化个数为n，那么相同精度的平面阵同时估计方位角和俯仰角时，其空域角度离散化个数为n²；

由空间几何位置关系知：

通过式(7)的空间合成角将二维空间角估计降维成两个一维空间角估计；

将式(7)分别带入式(5)，式(6)得到空间合成角的导向量：

步骤三：式(7)定义的空间合成角α_k和β_k已包含了二维信号的方位角和俯仰角等信息，并且在整个空域中具有稀疏性；因此将信号存在的所有可能的空间角度θ进行空域划分得到N个离散化的角度方向，用h_i表示第i个方向的信号，则空域信号可表示为h＝[h₁，h₂，...，h_N]^T，如果空间中有K个目标信源，显然h中这K个对应方向的信号有非零元素，其余方向都为零；

根据向量内积的定义可知，阵列流形矩阵A中任意向量之间的正交性用下式表示：

式中θ_p、θ_q表示空域划分后的不同角度值，其中p，q＝1，2，...，N；显然当两者角度最相近时其正交性最大；

将观测空间[-90°，90°]分别按等角和等弦均匀地分成N个角度

假设μ(A₁)表示等角划分下流型矩阵的正交性，则

假设μ(A₂)表示等弦划分下流型矩阵的正交性，则

对两者的正交性进行仿真比较，采用16阵元的均匀线阵，分别按等弦划分与等角划分将空间均匀划分180份；

存在一个临界角度满足等角划分与等弦划分非正交性相同，即：

μ(A₁)＝μ(A₂) (13)

将式(11)式(12)代入式(13)化简得：

cosθ_P-cosθ_P-1＝2/N (14)

由上式推出

对上式左边进行和差化积，可得

则上式关于θ_P的两个解可分别表示为θ_P1、θ_P2

即在(-90°，θ_P1)和(θ_P2，90°)等弦划分的方式更优，在(θ_P1，θ_P2)等角度划分的方式更优；根据角度范围选择空域的划分方式，再将其对应的阵列流型矩阵进行融合，得到最优的阵列流型矩阵；假设空域网格总数为n，则等角度划分网格数为(θ_P2-θ_P1)*n/180，剩余网格数为等正弦划分网格数；合并后的矩阵可表示为：

于是稀疏信号表示的DOA估计数据模型为：

其中X、Y分别为子阵X、Y的接收信号向量，N_x、N_y分别表示两子阵的噪声信号向量；扩张后的稀疏字典矩阵为：

式(20)式(21)是典型的压缩感知模型，扩展后的流形矩阵

扮演着测量矩阵φ的角色；因为目标信号h_x、h_y相对于整个空域是稀疏的，不需要稀疏基矩阵来稀疏表示，所以

既是观测矩阵也是感知矩阵；到此基于压缩感知的DOA估计模型就构造完了；

步骤三：通过上述估计模型的构造，求解式(24)(25)的稀疏重构问题即得到两组空间合成角α_k和β_k；即已知观测信号X、Y和感知矩阵

重构出目标稀疏信号h_x、h_y，表示为：

其中||h_x||₀、||h_y||₀是h_x、h_y的0范数，表示h_x、h_y中非零元素的数量，ε为观测信号中所含噪声的能量大小；

采用OMP算法对其进行求解得到两组DOA估计结果；

步骤四：通过上述重构算法求出两个子线阵对应的一组夹角α和β；若来波信号数量为1，那么将夹角α以及夹角β按照式(19)直接换算得到信源的二维空间角度；

若来波信号数量为多个，因为h_x和h_y中的非零元素大小对应着信号的幅值，所以采用信号幅度排序的方法来实现角度的配对，即对式(15)式(16)估计出的稀疏系数h_x和h_y中的非零元素按幅值大小排序后使其一一对应，最后将配对好的一组角度按照式(26)反解得到各信源的方位角θ和俯仰角

Images