超高斯噪声背景下的基于无网格压缩感知的DOA估计方法

摘要

本发明公开了一种超高斯噪声背景下的基于无网格压缩感知的DOA估计方法，首先确定天线阵信号模型，即通过一个阵元数为M，各阵元间距相等的线型天线阵列接收信号，并根据实际环境测定的噪声分布类型，选择lp的范数进行约束，然后确定求解无噪声信号x的函数表达式，无噪声信号x的恢复问题描述为原子范数最小化问题，采取半定规划理论方法就能求解；再通过ADMM算法求得原始解，并求得对偶解最后求解DOA，使之更加接近于原始信号的角度支撑集，从而完成DOA估计。本发明通过利用噪声的统计特性，对噪声项采用合适的范数进行有效的约束，最后利用原始解与对偶解的关系，求得原始信号的角度支撑集，实现高效精确的DOA估计。

说明书

技术领域

本发明属于无网格压缩感知(Grid-less Compressive Sensing,GCS)的波达方向(Direction-of-arrival,DOA)估计技术，特别是一种超高斯噪声背景下的基于无网格压缩感知的DOA估计方法。

背景技术

众所周知，波达方向估计在通信、语音、雷达、声纳等很多方面都有着不可替代的作用。压缩感知技术由于其独特的优势，在空间谱估计问题中得到了广泛的应用。对于DOA估计，压缩感知技术利用信号在空域中的稀疏性，只需要少量的观测数据，就能够实现对信号的重构，而且在信噪比(Signal to Noise Rate,SNR)较低，信源相关性较高的场合，算法仍然具有鲁棒性。然而，传统的超分辨算法，比如基于信号子空间的多重信号分类(Multiple Signal Classification,MUSIC)算法和旋转不变技术估计信号参数(Estimating Signal Parameter via Rotational Invariance Techniques,ESPRIT)算法等，并不具备这一优势。经过近年来的研究和发展，稀疏信号重构算法在小样本，甚至单次快拍的情况下也能实现超分辨的DOA估计。

在许多实际应用中，信号的参数空间是连续的。针对这一问题，传统压缩感知算法首先把连续参数空间进行离散化处理，将整个参数空间划分为有限个网格，并假设信号的参数能由某些网格点来表示。过于密集的网格虽然能提升估计的精度，但也会造成基字典中相邻原子之间的相关性太强，使得字典的约束等容性(RIP性)下降，从而降低压缩感知的重构性能。另外，信号的实际参数值很可能不会落在离散的网格上，假设的变换基则无法表示稀疏信号，研究者把这种现象称为基的失配(Basis Mismatch)。

针对该现象，研究者找到了解决欠定线性逆问题的最佳凸启发式方法——采用原子范数来刻画信号特征，即直接在连续参数空间中寻找最少的原子来表示给定的信号，而无需对稀疏域进行离散化处理，并提出了无网格压缩感知，或者称连续域压缩感知(Continuous Compressive Sensing,CCS)这一概念。参数空间不再被划分成有限个网格，而是直接被看作一个连续的变量。Tang等人在文献(1.Tang Gong-guo,Bhaskar B N,ShahP,et al.Compressed sensing off the grid[J].IEEE Transactions on InformationTheory,2013,59(11):7465-7490)中提出，当信号频率之间的间隔至少为时，只需O(K·logK·logM)个采样就能重构出信号的频率，其中ε为一个小常数，K为信号的个数，M为信号的长度。原子范数的求解依赖于半定规划(Semi-definite Programming,SDP)问题，而半定规划问题是一个非光滑的凸优化问题。当天线阵规模较大时，凸优化方法存在复杂度高、收敛速度慢等问题。为此，许多学者采用一种更好的近似求解算法——ADMM算法。

ADMM算法在压缩感知领域中取得了广泛的应用，比如图像处理，信道估计，无线传感器网络等。这些研究仍然依赖于格点划分。近些年来，相关研究已延伸到无网格压缩感知领域，比如雷达成像，加权原子范数最小化等。

目前，基于压缩感知稀疏重构的DOA估计算法大多是在高斯噪声背景下进行的研究。而自然界中的一些噪声往往具有非高斯性，无法用高斯分布来描述。而超高斯噪声是非高斯噪声当中的一类。低频大气噪声、语音信号、生物医学信号和其它脉冲噪声具有超高斯性。超高斯噪声的概率密度函数具有显著的尖峰脉冲和拖尾现象，噪声信号也会出现少数模值比较大的异常值。这种特性使这类非高斯过程的统计特性显著偏离高斯分布，如果仍然采用高斯分布模型来描述这些噪声，原来基于高斯噪声数学模型设计的DOA估计方法的性能会出现下降。

由于无网格压缩感知理论提出的时间比较晚，所以在现有的研究当中，只有在非高斯噪声背景下，未采用压缩感知理论进行DOA估计，比如文献(3.Wen-Jun Zeng,H.C.So,Lei Huang.lp-MUSIC:Robust Direction-of-Arrival Estimator for Impulsive NoiseEnvironments[J].IEEE Transactions on Signal Processing.Vol.61,No.17,Sep.2013,pp.4296-4308.)；和在非高斯噪声背景下，利用传统的压缩感知技术，通过划分格点进行DOA估计的文献，比如(4.孙飞，非高斯噪声环境下基于压缩感知的DOA估计[D]，大连海事大学，2015。)；以及利用无网格压缩感知理论，在高斯白噪声背景下进行的DOA估计，比如(4.Zai Yang,Lihua Xie.Enhancing Sparsity and Resolution via Reweighted AtomicNorm Minimization[J].IEEE Transactions on Signal Processing,Vol.64,No.4,Feb.2016,pp.995-1006.)。目前，尚未发现关于在超高斯噪声环境下，基于无网格压缩感知理论的DOA估计技术。因此，研究超高斯噪声环境下，基于无网格压缩感知理论的DOA估计方法具有非常重要的意义。

发明内容

本发明的目的在于提供一种超高斯噪声背景下的基于无网格压缩感知的DOA估计方法，实现高效精确的DOA估计。

实现本发明目的的技术解决方案为：一种超高斯噪声背景下的基于无网格压缩感知的DOA估计方法，步骤如下：

(1)确定天线阵信号模型，即通过一个阵元数为M，各阵元间距相等的线型天线阵列接收信号；天线阵接收到的信号模型表示为

其中y表示该天线阵接收到M个含噪声的信号；以M维观测向量表示无噪声信号x；是M维的导引矢量，则方向矩阵s代表K个源信号，w是M维的零均值，方差为σ²的超高斯噪声信号；

(2)根据实际环境测定的噪声分布类型，选择l_p的范数进行约束，1≤p＜2；

(3)确定求解无噪声信号x的函数表达式，无噪声信号x的恢复问题描述为原子范数最小化问题，采取半定规划理论方法就能求解，该原子范数最小化问题的函数表达式为：

其中s.t.为subject to的缩写，表示目标函数的约束条件；τ为正则项，取值为||·||_p表示向量或矩阵的l_p范数，为l_p范数的p次幂；另外u是M维的向量，T(u)表示u的Toeplitz矩阵，其维度是M×M；(·)^H表示向量或矩阵的共轭转置；

(4)ADMM算法求得原始解，利用凸优化工具包可以求解该原子范数最小化问题，或者为降低求解的复杂度，将步骤(3)当中的函数表达式写成拉格朗日函数的形式，并利用基于ADMM的一阶算法求解该问题，得到无噪声信号的估计，即该优化问题的原始解

(5)求得对偶解是x的对偶变量z的估计值；原子范数最小化问题的对偶问题的函数表达式为

其中为原子范数的对偶范数；由于原始和对偶问题之间存在强对偶性，分别是原始解和对偶解的最优解，二者之间的关系为

(6)求解DOA，在[-90°,90°]的范围内划分出格点，寻找的模为1，即的模为时的点所对应的角度作为DOA估计值；由于满足

看出的最大值等于用这种方法求解，使之更加接近于原始信号的角度支撑集，从而完成DOA估计。

本发明与现有技术相比，其显著优点：(1)通过利用噪声的统计特性，对噪声项采用合适的范数进行有效的约束，最后利用原始解与对偶解的关系，求得原始信号的角度支撑集。图3的仿真结果表明，由于成功利用了超噪声的统计特性，DOA估计的精度很高。(2)实现起来相对容易，具有很好的应用前景。

下面结合附图对本发明作进一步详细描述。

附图说明

图1是均匀线型阵列接收信号模型。

图2是本发明的流程图。

图3是本发明超高斯噪声背景下，采用不同范数对噪声项约束的DOA估计效果图。

具体实施方式

结合图1，对本发明的天线阵列信号模型进行说明：将无网格压缩感知技术应用到DOA估计中。考虑这样一个问题：有K个角度未知的目标随机地分布在θ∈[-90°,90°]区域中。为了确定它们的角度，本发明通过一个阵元数为M的线型天线阵列(Uniform LinearArray,ULA)接收信号，进而实现DOA估计。

本发明是基于均匀直线阵列，其中每个阵元都是全向天线，假设阵元数目为M，阵元之间的距离为d。假设K个窄带远场点源信号分别从方向θ_k，k＝1,…,K，入射到天线阵，即θ_k为信号入射角。则在某一时刻，天线阵接收到的信号表示为

其中y＝[y1,…,y_M]^T表示在该时刻M个天线阵元接收到的数据组成的M×1维列向量，x＝[x₁,…,x_M]^T表示无噪声的信号，(·)^T表示向量或矩阵的转置。为M×K维方向矩阵，K为信源数，是对应于相位差ω_k的导引矢量，ω_k＝2πdsin(θ_k)/λ表示两个相邻阵元间的相位差。λ为信号波长，S＝[s₁,…,s_k]^T为K个信号s₁,…,s_k组成的K×1维信号向量，w＝[w₁,…,w_M]^T为各个天线阵元上噪声组成的M×1维向量。假定信号与阵元噪声统计独立，各阵元噪声之间同样相互独立。

结合图2，本发明超高斯噪声背景下的基于无网格压缩感知的DOA估计方法，以无网格压缩感知理论为基础，采用l_p范数对接收信号中的超高斯噪声进行约束，通过求解原子范数最小化问题及其对偶问题，从而实现DOA的估计。总之，在无网格压缩感知领域中，DOA估计问题可分解为两个子问题：求解原子范数最小化问题，得到原始信号的估计值求解原子范数最小化问题的对偶问题，得到DOA估计值。这两个子问题一共可分解为以下六个步骤。

具体步骤如下：

步骤一，确定天线阵信号模型。采用一个均匀直线型阵列(Uniform lineararray,ULA)，其中每个阵元都是全向天线，假设阵元数目为M，阵元之间的距离为d。假设K个窄带远场点源信号分别从方向θ_k∈[-90°,90°]，k＝1,…,K，入射到天线阵，即θ_k为信号入射角。则在某一时刻，天线阵接收到的信号表示为

其中y＝[y1,…,y_M]^T表示在该时刻M个天线阵元接收到的数据组成的M×1维列向量，x＝[x₁,…,x_M]^T表示无噪声的信号，(·)^T表示向量或矩阵的转置。是对应于相位差ω_k的导引矢量，ω_k＝2πdsin(θ_k)/λ表示两个相邻阵元间的相位差，A＝[a(θ₁),a(θ₂),…,a(θ_K)]是M×K维的方向矩阵。λ为信号波长，S＝[s₁,…,s_k]^T为K个信号s₁,…,s_k组成的K×1维信号向量，w＝[w₁,…,w_M]^T为各个天线阵元上噪声组成的M×1维向量。假定信号与阵元噪声统计独立，各阵元噪声之间同样相互独立；

为σ²的超高斯噪声信号。

步骤一中的超高斯噪声信号服从广义高斯分布，该广义高斯分布的概率密度函数为：

其中w是函数变量，m表示其均值，Γ(·)为gamma函数，α反映广义高斯概率密度函数峰的宽度；若限制x的方差为1，α表示为β控制广义高斯概率密度函数的形状，它与概率密度函数的尖锐程度相关，β值越小则密度函数越尖锐；

β的值确定了源信号的类型：在x均值为0，方差为1时，若β＝2，信号x为标准高斯分布的信号；若0＜β＜2，x为超高斯信号，其中β＝1时，信号x呈Laplacian分布的信号；若β＞2，x为亚高斯信号；若β＝+∞，信号x依概率收敛于均匀分布U(-α,α)。

在仿真实验中，我们根据上述的广义高斯概率密度函数，利用文献(6.朱晓玲，姜浩，任意概率分布的伪随机数研究和实现[J]，计算机技术与发展，Vol.17，No.12，Dec.2007，pp.116-118,168.)提出的舍选法，产生超高斯噪声信号。其中，通过设置β＝0.5产生超高斯噪声。

步骤二，根据实际环境测定的噪声分布类型，选择l_p的范数进行约束，1≤p＜2。

针对超高斯噪声环境下基于无网格压缩感知的DOA估计问题，本发明提出对于优化目标函数的噪声项，采用l_p范数约束的方法。传统的DOA估计方法假定背景噪声为高斯白噪声，所以采用l₂范数约束噪声项；而在超高斯噪声背景下，同样的方法存在性能下降的问题；对此，本发明提出的一种应用于超高斯噪声背景下的无网格压缩感知DOA估计方法，通过利用噪声的统计特性，对噪声项采用合适的l_p范数进行有效的约束；超高斯噪声适合采用l_p(0＜p＜2)范数进行约束。

由于一些噪声信号具有很强的超高斯性，需要l_p(0＜p＜1)的范数进行约束，原问题不再是一个凸优化问题，而是一个非凸的l_p范数最小化问题，难以获得全局极小值。因此，在本发明中，将不会考虑l_p(0＜p＜1)范数这一情形，对超高斯噪声只采用l_p(1≤p＜2)的范数进行约束。

步骤三，确定求解无噪声信号x的函数表达式，无噪声信号x的恢复问题描述为原子范数最小化问题，采取半定规划理论方法就能求解。根据文献(2.B.N.Bhaskar,G.Tang,B.Recht.Atomic norm denoising with applications to line spectral estimation[J].IEEE Transactions on Signal Processing,2013,61(23):5987-5999)，该原子范数最小化问题可表示为

其中tr(·)表示矩阵的迹。

因此，该原子范数最小化的函数表达式为：

其中s.t.为subject to的缩写，表示目标函数的约束条件；τ为正则项，取值为||·||_p表示向量或矩阵的l_p范数，为l_p范数的p次幂；另外u是M维的向量，T(u)表示u的Toeplitz矩阵，其维度是M×M；(·)^H表示向量或矩阵的共轭转置。T(u)表示u的Toeplitz矩阵：

u_j表示u的第j个值。该函数表达式可以通过SeDuMi工具包求得参数x，u，t的解(8.J.F.Sturm,“Using>

在无网格压缩感知技术中，参数空间不再被划分成有限个网格，而是直接被看作一个连续的变量。令为一个原子集合，在DOA估计问题中，该集合可以写成其凸包相对于原点是一个中心对称的紧集，且包含原点作为内点。a为集合中的原子，由此即中的元素都是的极值点。此时由凸包的尺度函数定义的范数为原子范数，用表示，则有

其中inf为函数的下确界(或最大下界)。对于原子范数来说，原子集合中的原子就是用于构造任何一个信号的基本单元，那么凸包的低维平面同样对应的是由很少的相应原子所构成的信号，因此原子范数实际上是给集合增加了稀疏约束。这种约束方式将集合看作一个描述连续变化参数的无限字典，同时增加稀疏约束时没有引入离散化表示，避免了基的失配问题的产生。x则是相应原子的线性组合。原子范数具有半定规划(Semi-definite programming,SDP)性质，采取线性半定规划方法就能求解原子范数。根据Caratheodory引理，任何半正定Toeplitz矩阵都能进行Vandermonde分解，从而将原子范数最小化问题转化为半定规划问题。因此，需要运用半定规划理论来求解。我们利用现有的一些凸优化工具，如SeDuMi等，即可求解上述问题，得到接收到的无噪声信号的估计

在文献(2.B.N.Bhaskar,G.Tang,B.Recht.Atomic norm denoising withapplications to line spectral estimation[J].IEEE Transactions on SignalProcessing,2013,61(23):5987-5999)当中，原始解和对偶解之间的关系，以及正则项τ的上界和下界，是基于l₂范数约束的优化模型的推导得来的结果。但这一结果在采用l_p(1≤p＜2)范数约束的优化模型中不再适用。

在采用l₂范数约束的条件下，当时，即在波达方向上，的值始终等于常数τ。然而，在采用l_p(1≤p＜2)范数约束下，在波达方向上虽然会出现峰值，但这个值不再是一个常数。在实际情况下，信号源的数目往往是未知的。当信噪比较低时，DOA估计的空间谱上会出现杂峰，杂峰有时甚至会高于在波达方向上的峰值。此时，信号源和杂峰将难以区分。我们对正则项的上界和下界，以及原始解和对偶解之间的关系进行推导，使之适用于超高斯噪声背景和l_p(1≤p＜2)范数约束的优化模型。

步骤三中正则项τ的上界和下界的求解过程如下：

1)正则项τ的上界

变量z原子范数的对偶范数记为：

多项式Z_n的最大模值定义为

假定变量w₁,…,w_N满足广义高斯分布，其广义高斯分布的概率密度函数为：

参数的取值如下：β＝0.5，m＝0，

变量z满足z＝|w|^p-2·w，即

由此z的概率密度函数为

由此可得

其中，的详细推导过程如下：

其中Γ(v,z)是不完全伽玛函数，定义如下：

设定参数δ＝2log(N)，得

其中噪声变量w的标准差为σ，则变量z的标准差为σ^p-1，由式(3.1)、(3.7)得

设定N＝4πnlog(n)，则的上界写作

当n≥3时，根据权利要求1可知，正则项τ的上界是上界的倍；

其中，当p＝1时，即正则项τ的上界为：

2)正则项τ的下界

令F为变量u₁,…,u_N的分布函数，同时令M_n＝max(u₁,…,u_N)，则零均值超高斯变量的分布函数为

当n≥5时，正则项τ的下界为

其中，当p＝1时，即正则项τ的上界为

步骤四，ADMM算法求得原始解，利用凸优化工具(如SeDuMi工具包)来求解该原子范数最小化问题，或者为降低求解的复杂度，将步骤三当中的函数表达式写成拉格朗日函数的形式，并利用基于ADMM的一阶算法求解，得到无噪声信号的估计，即该优化问题的原始解(估计值)

该算法引入拉格朗日乘子Λ，同时令根据(9.孙志强，白圣建，郑永斌，刘伟(译)，最优化导论(第四版)[M]，电子工业出版社，pp.380-392)，步骤三当中的优化目标函数表达式，可以转化为拉格朗日函数的形式：

其中||·||_F表示向量或矩阵的Frobenius范数，<·>表示向量或者矩阵的内积，则(·)^H表示向量或矩阵的共轭转置；ρ表示拉格朗日函数的惩罚项，能够控制函数的收敛速率。通过对矩阵或向量求导的方法，对参数(t,u,x,Z,Λ)反复迭代更新，直至算法收敛。

l表示第l次迭代。

步骤五，求得对偶解是x的对偶变量z的估计值。对于每种范数来说都有相应的对偶范数，与原始范数相比，对偶范数通常具有一些有用的结构和性质，因此被广泛用于许多具体问题的分析和应用。原子范数的对偶范数等价于原子集合的支撑集函数，即

其中表示内积z^Hx的实部，sup为函数的上确界(或最小上界)，为原子范数的对偶范数。由于原始和对偶问题之间存在强对偶性，这意味着分别是原始和对偶问题的最优解。

原子范数最小化问题的对偶问题的函数表达式为

其中为原子范数的对偶范数；由于原始和对偶问题之间存在强对偶性，分别是原始解和对偶解的最优解，二者之间的关系为

步骤五中对偶函数表达式的求解过程如下：

原子范数最小化问题的函数表达式重新写作：

并将其转化为拉格朗日函数的形式：

用g(z)表示对偶函数：

令g(z)最大化，则第三项等于零，如式(5.5)所示。

步骤五中的原始解和对偶解之间关系的求解过程如下：

函数在处取得最小值，若令常数α∈(0,1)，对于变量x而言，即

是凸函数，根据凸函数的性质所以

为表示方便，令则(1)式的右式化为

令f(α)在α→0条件下满足Taylor展开的条件，故

令α→0，将上式代入式(1)、(2)，得

令变量其中·表示向量的点积，即对向量中每一个元素取p-2次幂的绝对值后得到的新向量，与的点积，等于故式(3)可化为

根据原子范数和对偶范数的定义，

其中inf为函数的下确界或最大下界，sup为函数的上确界或最小上界，根据式(5)和函数上下界的关系，式(4)可整理为

通过(6)得出，当且仅当时，取得最小值；

最后，在p＝2的条件下，τ满足的关系依然成立。

步骤六，求解DOA，在[-90°,90°]的范围内划分出格点，寻找的模为1，即的模为时的点所对应的角度作为DOA估计值；由于满足

的最大值等于用这种方法求解，使之更加接近于原始信号的角度支撑集，从而完成DOA估计。

在步骤六中DOA求解具体过程如下：

在DOA估计问题中，θ_i∈[-90°,-90°+Δθ,…,90°]，Δθ表示角度的分辨率；向量a的导引矢量形式为由于满足

也就是说，在此问题中的最大值等于可以在[-90°,90°]的范围内进行搜索，即求得与每一个向量a内积的实部，寻找到的模为时的点所对应的角度作为DOA估计值。角度的分辨率Δθ的值完全可以根据需要来设定，在仿真实验中，可以将Δθ精确到0.0001°甚至更高，如此高的分辨率使得实际参数和估计值之间的误差很小，足以实现超分辨DOA估计。通过这种方法求解，使之更加接近于原始信号的角度支撑集，完成DOA估计。

在步骤六中，我们利用推导的原始解和对偶解之间的关系，求得原始信号的角度支撑集。本发明成功利用了超高斯噪声的统计特性，相较于传统的采用l₂范数约束的方案，减少了DOA估计的误差，增大了波达方向角估计成功的概率，在相同的条件下获得了更高的DOA估计精度。概括地说，本发明首先通过求解优化目标函数，得到原始解和对偶解然后求解优化目标函数的对偶问题，得到DOA估计值。

值得一提的是，根据传统的压缩感知理论进行DOA估计，需要预先在参数空间内划分网格，并假设信号参数恰好在网格点上，同时还要选择合适的稀疏变换基。实际上，信号参数通常不会恰好在预先划分的网格点上，这会造成基的失配现象；如果网格划分过于精细，使得信号参数完全落在网格点上，又会造成基字典相邻的原子之间的相关性太强。以上两个因素都会导压缩感知重构性能降低的问题。

然而，应用无网格压缩感知理论进行DOA估计，在建模过程中没有涉及到网格划分和选择稀疏变换基的问题，不存在由以上的两个因素，而导致压缩感知重构性能降低的问题，体现出了无网格压缩感知理论的优越性。

结合图3对本发明进行性能分析。为了检验本发明的性能，我们进行了以下仿真实验，首先对仿真环境进行说明：选择一个线型天线阵列，阵元数M＝64，信源数K＝4。在仿真过程中，取ρ＝2。所有仿真实验重复100次，并且对结果取平均值。

下面对本发明性能进行检验：在超高斯噪声背景下，运用不同的l_p范数对噪声项进行约束。如前所述，超高斯噪声适合采用l_p(1≤p＜2)范数进行约束。然后，在相同信噪比条件下，通过对比采用不同范数约束下DOA的估计值与真实值，说明采用合适的l_p范数，能够实现波达方向的高精度估计。在此，验证在超高斯噪声背景下，采用l_p(1≤p＜2)范数对噪声项约束的DOA估计方法的优越性；

本方案分别采用l_1.2范数，l_1.5范数。输入信号的信噪比为5dB，β＝1.5，随机假设四个波达方向角分别为-59.35°,12.70°,23.52°,68.70°。

图3是第一次仿真实验中，DOA估计的空间谱，体现了l_p(1≤p＜2)范数约束的有效性。采用l_1.2，l_1.5两种范数约束的DOA估计方法，结果分别为-59.57°,12.55°,23.89°,68.59°，-59.23°,12.46°,23.88°,68.54°，平均误差分别为0.21°,0.22°。空间谱都能在相应的波达方向角上达到峰值。因此，在超高斯噪声背景下，采用l_1.2，l_1.5两种范数都能对噪声项进行有效约束，成功估计出DOA。