自动化学报  2018, Vol. 44 Issue (3): 443-452   PDF    
一种针对单快拍DOA估计的子空间搜索近似消息传递算法
曾令豪1, 刘静1, 韩崇昭1     
1. 西安交通大学智能网络与网络安全教育部重点实验室 西安 710049
摘要: 压缩感知(Compressed sensing,CS)技术应用于单快拍波达方向(Direction of arrival,DOA)估计中可以实现相关信号的超分辨估计,但会遇到感知矩阵高相干性以及对噪声敏感的问题.本文提出一种基于近似消息传递的子空间搜索算法以解决上述问题.该算法首先通过近似消息传递算法得到一个粗解,随后利用该粗解划分子空间,最后在子空间中寻找精确解.仿真结果验证了所提算法的有效性.文章最后通过理论分析了该算法性能并讨论了算法在信号数未知时的扩展应用.
关键词: DOA估计     压缩感知     近似消息传递算法     子空间搜索    
A Subspace Searching Approximation Message Passing Algorithm for Single Snapshot DOA Estimation
ZENG Ling-Hao1, LIU Jing1, HAN Chong-Zhao1     
1. Ministry of Education Key Laboratory for Intelligent Networks and Network Security, Xi'an Jiaotong University, Xi'an 710049
Manuscript received : July 21, 2016, accepted: May 6, 2017.
Foundation Item: Supported by National Basic Research Program of China (973 Program) (2013CB329405), National Natural Science Foundation of China (61221063, 61573271, 61573276, 61370037)
Author brief: ZENG Ling-Hao Ph. D. candidate at the Institute of Integrated Automation, School of Electronic and Information Engineering, Xi'an Jiaotong University. His research interest covers compressed sensing and target tracking;
HAN Chong-Zhao Professor at the School of Electronic and Information Engineering, Xi'an Jiaotong University. His research interest covers multi-source information fusion, stochastic control and adaptive control, and nonlinear spectral analysis
Corresponding author. LIU Jing Associate professor at the School of Electronics and Information Engineering, Xi'an Jiaotong University. Her research interest covers compressed sensing and information fusion. Corresponding author of this paper
Recommended by Associate Editor XIN Jing-Min
Abstract: When compressed sensing (CS) is applied to single snapshot direction of arrival (DOA) estimation, super-resolution reconstruction of the correlated signals becomes possible. However, problems such as high coherent sensing matrix and noise sensitivity come into being at the same time. In order to solve these problems, an approximate message passing based subspace search algorithm is proposed. Firstly, a rough solution is obtained by approximate message passing algorithm. Secondly, subspaces are divided according to the rough solution. Thirdly, an exact solution is found by searching the subspaces. Simulation results show the effectiveness of the proposed approach. Finally, the performance of the proposed algorithm is analzed, and the application in which the number of signals is unknown is discussed as well the performance.
Key words: Direction of arrival (DOA) estimation     compressed sensing (CS)     approximate message passing (AMP)     subspace searching    

波达方向(Direction of arrival, DOA)估计是阵列信号处理的关键问题, 在雷达、红外、声呐和地震等领域有着广泛的应用.以多重信号分类(Multiple signal classification, MUSIC)算法[1]和信号参数旋转不变技术(Estimating signal parameters via rotational invariance techniques, ESPRIT)[2]为代表的空间谱估计方法突破了瑞利限, 实现了目标DOA的超分辨估计.但这些方法是根据阵列接收信号的统计特性来估计目标的到达角, 因此需要大量独立同分布的量测数据.为了减少所需的测量数据, 甚至实现单快拍DOA估计, 近年来出现的压缩感知(Compressive sensing, CS)技术[3]带来了新的解决方案.当原始信号满足一定稀疏性条件时, 压缩感知可以实现以远少于经典采样理论所要求的采样数来精确重构信号.压缩感知技术带来的另一好处是可以实现相干信号的超分辨估计.目前, 国内外已有不少基于压缩感知的DOA估计的研究[4-8].

基于压缩感知的DOA估计算法, 大致可以分为以下两类:

一类是将CS与传统的空间谱估计相结合.文献[4]首先将稀疏性引入DOA估计, 提出了$l_1$范数奇异值分解算法, 该算法使用奇异值分解来降低重构的计算复杂度.文献[5]提出了一种CS-MUSIC算法, 利用多快拍数据与压缩感知中多测量矢量模型结合进行DOA估计, 该算法是先利用压缩感知技术估计噪声子空间, 再使用MUSIC算法实现DOA的估计.由于该算法利用了压缩感知技术减少了空间谱估计所需的量测数量, 因此其计算复杂度相对较高, 是一种在直接压缩感知重构与空间谱估计之间权衡的算法.

另一类是将DOA估计看作压缩感知的重构问题, 直接使用压缩感知重构算法来得到DOA估计.文献[6]是在单快拍情况下直接使用压缩感知重构算法得到DOA估计, 并对几种不同的压缩感知重构算法的结果进行了比较.这类算法充分利用了压缩感知技术能够实现稀疏重构的特点, 实现了单快拍DOA估计.但这类算法也存在固有的缺陷: 1) DOA估计所使用的感知矩阵具有较高的相干性; 2)现有大多数的压缩感知算法对噪声敏感.针对感知矩阵高相干性的问题, 文献[7-8]提出了一种广义相似感知矩阵匹配追踪算法, 该算法利用构造相似感知矩阵来降低原感知矩阵的高相干性.然而由于相似感知矩阵的构造是预先确定的, 其估计结果在有噪声的情况下会出现一定的偏差.考虑到噪声在DOA估计中是不可忽视的, 算法的抗噪声能力也是十分重要的的性能指标, 因此需要对该算法进行改进.

一种被称为近似消息传递(Approximate message passing, AMP)的算法由于抗噪声能力强, 且计算复杂度低而受到了关注. AMP算法是由Donoho等首先提出来的[9], 它是由和积信度传播算法推导而来, 并被证明与基追踪降噪(Basis pursuit denoising, BPDN)算法是等价的[10].文献[11]分析了一大类被称之为广义AMP的算法, 使用一个简单的状态演化方程描述广义AMP算法在大的独立同分布高斯变换矩阵下的渐进行为.而文献[12]提出了一种针对原信号非零元素先验分布未知的期望最大化高斯混合AMP算法. AMP算法在图像处理[13-14]、医学图像处理[15]等领域已经得到了应用, 并被证明是简单有效的.然而将AMP算法用于DOA估计, 也会受到感知矩阵高相干性的影响, 无法得到分辨率较高的结果.

本文提出了一种基于AMP算法的子空间搜索近似消息传递(Subspace searching AMP, SSAMP)算法来解决感知矩阵高相干性与噪声同时存在而产生的问题.该算法分为两步:第一步, 由AMP算法得到一个较粗的解.由分析可知, 在噪声较小时真实信号必然位于该粗解之中.第二步, 在该粗解对应的子空间中进一步寻找精确解, 从而实现超分辨估计.本文所提算法对粗解所对应的子空间的划分, 是以真实非零元素为中心的自适应的划分, 避免了预先划分不当造成的损失. SSAMP算法继承了AMP算法计算复杂度低、抗噪声能力强的特点, 同时又解决了AMP算法在感知矩相干性较高时所产生的问题.

1 基于压缩感知的单快拍DOA估计 1.1 单快拍DOA估计的稀疏化描述

本文考虑单快拍条件下的DOA估计, 即在已知信号来波频率的情况下, 每一时刻仅利用单个量测数据进行DOA估计.假设有$K$个远场窄带信号, 其载波波长为$\lambda$.接收端是由$M$个阵元组成的均匀线阵沿$x$轴排列, 它们的间距为$d = \lambda /{2}$, 每个阵元的接收噪声是相互独立的零均值高斯噪声.定义入射角为波信号与$y$的夹角.在某一时刻, 接收端得到的信号为

$ {\pmb{y}}(t) = \sum\limits_{k = 1}^K {{\pmb{a}}({\theta _k}){\beta_k}(t)} + {\pmb{e}}(t) $ (1)

其中, ${\pmb{e}}(t)$表示接收端的噪声, ${\beta_k}(t)$表示第$k$个目标的复振幅, ${\pmb{a}}({\theta _k})$表示相应来波角度$\theta _k$的导向矢量, 可表示为[16]

$ \begin{array}{l} {{\pmb a}({\theta _k})} = {\left[{1, {{\rm{e}}^{{\rm{-j}}\frac{{2\pi d}}{\lambda }\sin {\theta _k}}}, \cdots, {{\rm{e}}^{{\rm{-j}}\left( {M-1} \right)\frac{{2\pi d}}{\lambda }\sin {\theta _k}}}} \right]^{\rm{T}}}\\ \end{array} $ (2)

将整个监视区域离散化, 并按角度等分为$N$份, 用符号$\theta_i, i=1, \cdots, N$来表示.可以将式(1)写成如下形式:

$ {\pmb y}(t) = {\Phi}{\pmb \beta }(t) + {\pmb e}(t) $ (3)

其中, ${\pmb{\beta}} = {\left[{{\beta_1}, \cdots, {\beta_N}} \right]^{\rm{T}}}$表示相应角度上的信号复振幅, $\Phi=[{\pmb a}(\theta_1), \cdots, {\pmb a}(\theta_N)]$$M \times N$维阵列流型矩阵, 也称感知矩阵.由于后续压缩感知算法需要使用列归一化的感知矩阵, 因此将式(3)改写为

$ {\pmb y}(t) = {A}{\pmb x }(t) + {\pmb e}(t) $ (4)

其中按列归一化的感知矩阵$A$的每一列为

$ \begin{array}{l} {A_i} = \dfrac{1}{{\sqrt M }}{{\pmb a}({\theta _i})}\\ \end{array} $ (5)

在本文中称${\pmb x}$为信号的幅值, 有${\pmb{x}} = \sqrt M {\pmb{\beta}}$, 并且假设${\pmb e}$$M$维独立同分布的高斯白噪声, 其每一维的均值为$0$, 方差为$\sigma^2$.本文中对信噪比的定义为[17]

$ {\rm{SNR}} = 20{\rm lg}\left( \frac{\left\| {\pmb{\beta}} \right\| }{\sqrt{M} \sigma} \right) = 20{\rm lg}\left( \frac{\left\| {\pmb{x}} \right\| }{M \sigma } \right) $ (6)

根据信号的稀疏模型可知, $K$个来波信号相对于整个监视区域可以认为是稀疏的, 即$K \ll N$.这满足了压缩感知的使用条件, 因此可以使用压缩感知技术重构${\pmb{x}}$, 完成单快拍DOA的估计.将上述单快拍DOA估计描述为如下压缩感知重构问题:

$ \min{\left\| {\pmb x} \right\|_0}, {\rm{s.t.}}{\pmb y} = A{\pmb x} + {\pmb e} $ (7)

其中, ${\left\| {\pmb x} \right\|_0}$表示向量${\pmb x}$的非零元素数目.

1.2 AMP算法

在众多压缩感知重构算法中, 一种被称作近似消息传递的算法因为计算复杂度低以及具有良好的去噪声能力而引起了人们的关注. AMP算法是一种软阈值迭代算法, 其迭代的初值选取为${\pmb x}^0 = 0, {\pmb z}^0={\pmb y}$, 迭代公式如下所示[9]:

$ {{\pmb x}^{t + 1}} = \eta \left( {{{\pmb x}^t} + {A^*}{{\pmb z}^t}, {{\tau} ^t}} \right) $ (8)
$ {{\pmb z}^t} = {\pmb y} - A{{\pmb x}^t} + \frac{{{{\pmb z}^{t - 1}}}}{\delta }\left\langle {\eta^\prime\left( {{{\pmb x}^{t - 1}} + {A^*}{{\pmb z}^{t - 1}}, {{ \tau} ^{t - 1}}} \right)} \right\rangle $ (9)

其中, $A^*$表示$A$的共轭转置, ${\pmb x}^t$表示第$t$次迭代时原信号的估计值, ${\pmb z}^t$表示第$t$次迭代的残差. $\eta \left(\cdot \right)$${\pmb{x}}$的每一维而言是一个标量非线性的阈值函数, 它起到了在每次迭代中将结果导向更加稀疏的方向的作用.若去掉阈值函数该算法将收敛于式(8)的最小$l_2$范数解[18]. ${{\tau} ^t}$是第$t$次迭代时的阈值, $\eta ^\prime\left(\cdot \right)$表示$\eta \left( \cdot \right)$的导数, $\left\langle {\pmb x} \right\rangle $表示向量${\pmb x}$的平均值, $\delta = {M}/{N}$.阈值函数的具体形式为[9]

$ \eta \left( {{\pmb x}, { \tau} } \right) = {\mathop{\rm sgn}} \left( {\pmb x} \right){\left( {\left| {\pmb x} \right| - { \tau} } \right)_ + } $ (10)

其中, ${\left( {\pmb z} \right)_ + }$表示${\pmb z}$大于$0$的部分.文献[10]中给出了参数${ \tau}$的各种选择方法.若参数${ \tau}$由迭代式(11)表示:

$ {{ \tau} ^{t + 1}} = \rho + \frac{{{ \tau} ^t}}{\delta }\left\langle {\eta '\left( {{{\pmb x}^t} + A{{\pmb z}^t}, {{ \tau} ^t}} \right)} \right\rangle $ (11)

则AMP算法的迭代结果将收敛于如下被称为BPDN问题的解[10].

$ \min \left( {\frac{1}{2}\left\| {{\pmb y} - A{\pmb x}} \right\|_2^2 + \rho {{\left\| {\pmb x} \right\|}_1}} \right) $ (12)

其中, $\rho$是平衡参数.

2 感知矩阵的高相干问题

直接将AMP算法应用于单快拍DOA估计, 无法得到一个令人满意的结果.如图 1所示, AMP算法将原信号中的点状非零元素重构成在一个区域, 就如同在重构过程中将原信号的能量向其邻域内“泄漏”.这样, 我们无法直接得到一个精确的DOA估计结果.本节将分析造成该现象的原因.

图 1 AMP算法重构结果 Figure 1 Reconstruction result of AMP algorithm

定义矩阵的列相干系数:

$ {g_{ij}} :=\frac{{\left\| {A_i^*{A_j}} \right\|}}{{\left\| {{A_i}} \right\|\left\| {{A_j}} \right\|}} $ (13)

其中, $\left\| {\cdot} \right\|$表示向量的模, “$*$”表示共轭转置.当$i=j$时有$g_{ii}=1$, 此时被称为自相干; 当$i\neq j$时有${g_{ij}} = g_{ji} < 1$表示矩阵$A$的列互相干系数.令$\mu \left( A \right) = {\max }_{i \ne j} \left( {{g_{ij}}} \right)$表示最大列互相干系数.

文献[19]给出了无噪声时保证阈值算法性能的充分条件:

定理1[19].对于式(4)定义的线性系统方程${\pmb y} = A{\pmb x}$, 如果存在一个解${\pmb x}$, 其最大非零值为$\left\| x_{\max} \right\|$, 最小非零值为$\left\| x_{\min} \right\|$满足:

$ {\left\| {\pmb x} \right\|_0} < 1 + \frac{\left\| x_{\min} \right\|}{\left\| x_{\max} \right\|} \frac{1}{{\mu \left( A \right)}} $ (14)

则可以保证阈值算法能找到这个解.

该定理使用了最大互相干系数描述了感知矩阵与算法性能之间的关系.因此, 需要分析单快拍DOA估计所知用的感矩阵的相干性.

与压缩感知常用的随机矩阵不同, 单快拍DOA感知矩阵是一个确定性的矩阵.在DOA估计中, 将相同的监视区域划分得越细($N$越大)估计结果就越精确.但与此同时, 就会带来的感知矩阵的高相干性问题.若用$G=A^*A$表示矩阵$A$的格拉姆矩阵, 则定义感知矩阵$A$的列相干系数矩阵$C = {\left\| {A^*{A}} \right\|}$, 即$C$的元素$C_{i, j}$就表示矩阵的第$i$列与第$j$列的相干系数.

图 2图 3分别展示了DOA感知矩阵与随机感知矩阵的相干系数矩阵.图中可以看出随机感知矩阵列的最大互相干系数一般较小, 并且最大互相干列的位置是随机的.而DOA感知矩阵的最大互相干系数与自相干系数十分接近, 且相距越近的列其互相干系数越大, 相距较远的列相干程度较低.

图 2 DOA感知矩阵的列相干系数矩阵 Figure 2 Coherent coefficient matrix of DOA sensing matrix
图 3 随机感知矩阵的列相干系数矩 Figure 3 Coherent coefficient matrix of random sensing matrix

研究AMP算法的迭代公式(式(8)和(9)), 利用感知矩阵$A$的格拉姆矩阵$G$, 可以将式(8)改写为

$ {{\pmb x}^{t + 1}} = \eta \left( {{\pmb x} + \left( {I - G} \right){{\widetilde {\pmb x}}^t} + {\delta ^{t}}{A^*}{{\pmb z}^{t - 1}}+A^*{\pmb e}, {\tau ^t}} \right) $ (15)

其中, $I$表示$N$维的单位矩阵, ${\pmb x}$表示原信号, ${\widetilde {\pmb x}^t} = {{\pmb x}^t} - {\pmb x}$, ${\mathit{\boldsymbol{\tilde x}}^t} = {\mathit{\boldsymbol{x}}^t} - \mathit{\boldsymbol{x}},\;{\delta ^t} = \left\langle {\eta '\left( {{\mathit{\boldsymbol{x}}^{t - 1}} + {A^*}{\mathit{\boldsymbol{z}}^{t - 1}},\;{\tau ^{t - 1}}} \right)} \right\rangle /\delta $.

文献[9-10]推导AMP算法迭代公式是基于假设$M \to \infty $.当其他条件不变时, 若$M \to \infty $$G \to I $, 此时由式(15)可以看出${\pmb x}^{t} \to {\pmb{x}}$.文献[9]又指出当$M$较大($10^{2} \sim 10^{3}$)时, AMP算法依然有效.然而在本文的应用场景中$M$表示阵元数量, 这是客观限定的, 一般不能满足上述假设, 因此AMP算法不能得到较好的结果.

3 SSAMP算法 3.1 算法思路

研究图 4所示的DOA感知矩阵列的相干系数, 发现该相干系数可以被分为两个部分:一部分是相干系数较大且较为集中的列, 被称为高相干列, 即图 4中“主瓣”相应的部分, 且相距越近相干系数就越大; 另一部分是相干系数较小且更为分散的列, 被称为低相干列.通过对比图 1所示的AMP算法重构结果与图 4所示的列相干系数, 可以发现AMP算法的重构结果对应于原信号支撑集的高相干列.

图 4 DOA感知矩阵某一列的相干系数 Figure 4 The coherence coefficient of one column in sensing matrix

虽然AMP算法在感知矩阵相干程度较高时不能得到一个精确的重构结果, 但是其重构结果的支撑集包含了原信号的支撑集及其高相干列.从这一点出发, 可以考虑利用AMP算法的重构结果进行更进一步的精确求解.

因此我们提出了一种基于AMP算法的子空间搜索近似消息传递(SSAMP)算法. SSAMP算法的求解过程总体上分为两步:第一步由AMP算法求得粗解; 第二步在粗解的支撑集中找出精确解.在不考虑噪声以及低相干列的影响时, 每一步迭代结果中的局部最大值就对应原信号的支撑集.而噪声会使迭代中局部最大值偏离原信号的支撑集.但是可以发现当信噪比大于某一水平时, 就能保证粗解的支撑集包含原信号的支撑集.将AMP算法得到的粗解记为${\widehat {\pmb x}_{\rm{rough}}} $, 由其支撑集${Supp}$划分出若干子空间.进而在这些子空间中, 利用残差最小原则寻找原信号的支撑集, 并求得精确解.由于AMP算法本身的计算复杂度较低, 同时这些子空间相对较小, 因此该算法的计算效率也较好.值得说明的是这种子空间的划分不是预先设定的, 而是由量测驱动的自适应划分.

3.2 阈值的选择

AMP算法是一种迭代阈值算法, 阈值初值设置过大会导致幅值较小的信号被去除, 而阈值初值设置较小会使收敛速度降低.阈值初值的选取原则应当是在保证算法速度的前提下尽可能的小.我们通过仿真的方式来分析不同的阈值初值的选择对AMP算法收敛速度的影响.假设信号的幅值$x=1$, 阈值初值选择从0.001到0.3, 图 5显示了1 000次蒙特卡洛仿真中AMP算法收敛的平均迭代次数.从图 5中可以看出阈值初值选择在大于$0.05$时, 迭代次数几乎不变.

图 5 阈值选择与迭代次数关系 Figure 5 The iteration number versus threshold value

因此, 若已知信号的幅值的水平为$x$, 则建议阈值初值设置为略大于$0.05x$.若信号幅值水平未知, 已知信号个数$K$, 则可以估计信号幅值水平$x={{\left\| y \right\|}}/{K}$.本文仿真实验中阈值初值的选择为$\tau_0=0.1x$.

3.3 算法描述

为了解决高相干性对AMP算法造成的影响并保留AMP算法在计算复杂度和抗噪声能力上的优势, 本文提出了一种基于子空间搜索近似消息传递(SSAMP)算法. SSAMP算法需要已知感知矩阵$A$、量测${\pmb y}$、初始阈值${\tau_0}$和信号数$K$, 最终得到重构结果${\widehat {\pmb x}}$.共分为如下$5$个步骤:

步骤1. 求粗解

根据AMP算法式(8)~(-11)迭代计算直到得到一个收敛的结果或者得到一个是残差最小的结果, 并将其记为粗解${\widehat {\pmb x}_{\rm{rough}}}$.

步骤2. 子空间划分

按顺序遍历${\widehat {\pmb x}}_{\rm{rough}}$中的元素, 子集数$r$初始设置为$0$.如果第一个元素${\widehat {\pmb x}}_{{\rm {rough}}, 1} \ne 0$, 子集数$r=1$, 并将其指标$1$放入子集${Supp}_1$.

从第二个元素开始, 若某一元素${\widehat {\pmb x}}_{{\rm {rough}}, i} \ne 0$${\widehat {\pmb x}}_{{\rm {rough}}, i-1} = 0$, 子集数$r=r+1$, 并将指标$i$作为第一个元素放入新子集${Supp}_r$.若某${\widehat {\pmb x}}_{{\rm {rough}}, i} \ne 0$${\widehat {\pmb x}}_{{\rm {rough}}, i-1} \ne 0$, 子集数不变, 将指标$i$并入子集${Supp}_r$.

这样可以提取出$N_R$个子集(子空间), 分别记为$ {{Supp}}_{r}, r = 1, \cdots, {N_R}$.将各子集的中位元素作为该子集对应的尝试解支撑元素位置的初值, 并将它们合并为初始尝试解的支撑集$S$.

步骤3. 求精确解1

按顺序选择每一个子集, 保持尝试解支撑集${S}$中其他子集对应的元素不变, 用该子集中的每一个元素替换$S$中该子集相应的元素$S_r$.用$A_S$表示$A$中与$S$相应的列组成的矩阵, 根据

$ {\widehat {\pmb x}}_{t} ={A}^{*}_{S} {\pmb y} $ (16)

求出尝试解, 并找出使残差

$ {{res}}_{t} = {\left\| {{\pmb y} - {A}_{S}{{\widehat {\pmb x}}_{t}}} \right\|} $ (17)

最小的支撑集$S_r$作为该子集的一元最优解.当$N_R = K$时, 说明一个子集只包含一个非零元素即一元解$S_r^{(1)}$.此时所有子集的一元最优解的集合就组成了精确解的支撑集$S$.

步骤4. 求精确解2

$N_R < K$时, 说明在某些子集中存在多个非零元素.因此在每个子集中寻找最优解就需要考虑最多存在$N_{sol}=k-N_R+1$个非零元素的情况.假设第$r$个子集中非零元素的个数$p$, 遍历该子集中所有$p$元组合构成的备选集合$T$.用$p$个元素替换$S_r$, 按式(16)和式(17)计算得到$p$元最优解$S_r^{(p)}$和相应的残差$res^{(p)}$.从$p=2$开始, 如果$res^{(p)}>res^{(p-1)}$则说明该子集中只有$p-1$个非零元素, 并用$S_r^{(p-1)}$替换$S_r$得到新的$S$.更新$N_{sol}=K-{\rm size}(S)+1$.直到size$(S)=K$.

步骤5. 求最终解

根据步骤3和步骤4得到的支撑集$S$和式(16)求出最终解${\widehat {\pmb x}}$.

SSAMP算法的伪代码描述见表 1.

表 1 SSAMP算法伪代码 Table 1 SSAMP algorithm pseudocode
4 单快拍DOA估计应用的仿真

本节将利用仿真实验来分析讨论SSAMP算法在单快拍DOA估计应用中的表现.由于DOA估计对算法的抗噪声性能和超分辨性能有着特殊的要求, 因此将通过两组仿真来验证SSAMP算法在单快拍DOA估计应用中的性能, 并与经典的DOA估计MUSIC算法[1]、ESPRIT算法[2], 压缩感知类的SP[20]、OMP算法[6]、SSDOA-L1算法[6]、SSDOA-RFOCUSS算法[6]以及GSSMP算法[8]进行了比较.值得注意的是, MUSIC算法以及ESPRIT算法是无法实现单快拍DOA估计的, 因此在仿真2中这两种算法均使用了$10$个快拍的数据来得到估计的结果.

4.1 仿真1

仿真场景1的设置如下:接收阵列有$M=30$个阵元沿$x$轴呈直线排列, 其阵元间距为半波长.监控角度的范围为$\theta \in \left[{{{15}^\circ}, {{60}^\circ}} \right]$, 将上述监控区域按等角度分为$N=200$个分辨单元.设信号数目为$K=1$, 其幅值设为$x_i=1$, 其方位角$\theta_i$在上述监控区域内随机均匀分布.将信噪比范围设置为从$-20\, \rm{dB}$$70\, \rm{dB}$.一共进行了$500$次蒙特卡洛仿真实验.仿真硬件环境为: Matlab R2011b, Windows 7 64 bit, Intel Core i5-4570 CPU 3.20 GHz, RAM 4.00 GB.

设置该场景的目的是比较不同压缩感知算法在单快拍DOA估计中的性能, 即在不同的信噪比条件下对比SP、OMP、SSDOA-L1、SSDOA-RFOCUSS、GSSMP与SSAMP算法的单信号DOA估计误差和运行时间.

DOA估计误差将由$500$次蒙特卡洛仿真实验的均方根误差来表示.

$ {\rm RMSE }= \sqrt {\sum\limits_{mc = 1}^{500} \dfrac{{\left| {\theta - \widehat \theta_{mc} } \right|}^2 } {500} } $ (18)

其中, $\theta, \widehat \theta$分别表示真实来波方向与估计来波方向.

图 6展示了几种压缩感知算法来波方向估计均方根误差的比较结果.特别指出的是, 为了结果在对数坐标显示的方便, 图 6$10^{-9}$即对应RMSE为$0$.

图 6 角度估计的RMSE与信噪比关系 Figure 6 RMSE in angle estimation versus SNR

当信噪比$ >20\, {\rm{dB}}$时, SSAMP、OMP、SSDOA-L1与SSDOA-RFOCUSS算法单信号DOA估计的RMSE为$0$, 意味着这些算法此时能完全抵消噪声的影响实现精确的DOA估计. SSDOA-RFOCUSS算法在小于$ 20\, {\rm{dB}}$时, SSAMP算法在小于$ 5\, {\rm{dB}}$时, OMP算法在小于$ 3\, {\rm{dB}}$时以及SSDOA-L1算法在小于$ -16\, {\rm{dB}}$时, 不能保证精确的DOA估计结果.其余几种算法则一直不能得到精确的DOA估计结果, 但当信噪比$ >5\, {\rm{db}}$时RMSE一直在$10^{-2}$的水平左右.大多数算法在信噪比从$ 5\, {\rm{dB}}$下降到$ -15\, {\rm{dB}}$时, RMSE (Root mean suare error)从$10^{-2}$增加到$10^{-1}$.当信噪比小于$ -19\, {\rm{dB}}$左右时, 全部算法的RMSE都跳变到了$10^1$的水平.

图 7展示了SSAMP与其他几种算法在不同信噪比条件下运行时间的比较结果, 结果是500次蒙特卡罗仿真结果的平均.特别指出的是SSDOA-L1算法使用的时CVX工具包求解.

图 7 运行时间与信噪比关系 Figure 7 Execute time versus SNR

可以看出SSAMP、OMP、SP三种算法的运行时间较小, 而其余三种算法的运行时间要大得多. SSAMP算法的运行时间在大部分情况下优于其他几种算法的, 仅略高于OMP算法.当信噪比小于$ -15\, {\rm{dB}}$时, SSAMP算法运行时间会略为增加.

4.2 仿真2

仿真2的目的是为了对比几种经典算法与压缩感知类算法在不同信噪比条件下可以分辨的最小角度.

仿真2的设置与仿真1的主要区别在于仿真2中设置的信号数$K=2$, 其幅值均设为$1$.信号1的方位角$\theta_{i1}$为随机变量, 在监控区域内均匀分布, 信号2的方位角$\theta_{i2}$为在监控区域内与$\theta_{i1}$相差一个固定的角度间隔$\Delta \theta$.

在每一次重复的仿真实验中, $\Delta \theta$$22.5^\circ$逐渐减小到$0.225^\circ$, 即间隔的分辨单元从100依次减小到1.当估计结果中两个信源的估计误差分别都小于$2^\circ$, 则认为分辨成功; 否则, 认为该分辨失败.当某一算法不能成功分辨一个角度时, 则认为上一个角度间隔为其最小可分辨角度.记录下在不同信噪比条件下, 各个算法能分辨的最小角度.仿真结果是500次蒙特卡罗仿真中最小可分辨角度的平均值.特别说明当在最大角度间隔时, 某一算法仍无法有效分辨两个来波信号, 则将其最小分可辨角度记为$22.5^\circ$.

图 8展示了信噪比与该算法能分辨的最小间隔角度的关系.可以看出SSAMP算法与GSSMP算法在信噪比大于-2 dB时, 以最小可分辨角度定义的分辨性能要好于其他算法. ESPRIT算法认为不能分辨相干信号. MUSIC算法与OMP算法最小可分辨角度大约等于$3^\circ$. MUSIC算法分辨性能几乎不随信噪比改变. OMP算法信噪比大于$ -5\, {\rm dB}$时, 分辨性能几乎不随信噪比改变. SSAMP、GSSMP、OMP以及SSDOA-RFOCUSS算法在$ -5\, {\rm dB}$$0\, {\rm dB}$时, 分辨性能急剧下降.

图 8 最小分辨角度与信噪比的关系 Figure 8 Minimum resolution angle versus SN
4.3 仿真结果分析

从单信号DOA估计的均方根误差角度看SSDOA-L1算法是最优秀的, 然而其使用二阶锥规划方法求解使其运行时间较高.同时其最小可分辨角度仿真结果也较差.这可能是由于其分辨两个信号时得到的结果具有随机性, 即使相距角度较大也不能保证$100\, \%$的概率.而仿真2的结果说明的是一个算法能够稳定得到的最小可分辨角度, 在这个意义上SSDOA-L1算法性能较差.

而在最小可分辨角度性能中较好的GSSMP算法, 但单信号DOA估计的RMSE较大. GSSMP算法也是一种针对感知矩阵高相关问题的算法.与本文所提SSAMP算法类似, GSSMP算法也是先求粗解再进行精细解的寻找的两步式算法.但GSSMP算法处理高相干感知矩阵问题与本文所提算法的不同. GSSMP算法先将感知矩阵的高相干列进行合并形成一个互相干度较低的相似感知矩阵.而后通过相似感知矩阵求出粗解, 进而通过粗解对应的高相干列寻找出精确解.然而由于其对高相干列的划分是预先确定的, 不能够保证非零元素位于高相干列组成的子空间的中心.若非零元素恰好位于子空间的边缘, 由于噪声的影响会出现粗解错误的指向与非零元素相邻的子空间, 因而无法找出正确的解.因此GSSMP算法无法得到如本文所提算法在高信噪比条件下RMSE$=0$的情况.

SSDOA-RFOCUSS算法仿真中表现的性能与文献[6]中的结论有较大差距.在本文给出的仿真条件下, SSDOA-RFOCUSS算法迭代公式中的$AA^*$这一项是不满秩的, 其逆不存在, 只能使用伪逆代替.这可能是SSDOA-RFOCUSS在本文仿真中性能较差的原因.

OMP算法在运行时间、单信号DOA估计RMSE和最小可分辨角度三项性能指标的比较中综合性能较好.本文所提SSAMP算法与其比较.在运行时间和单信号DOA估计RMSE的比较中, OMP算法略微好于SSAMP算法.而在最小可分辨角度的比较中, 当信噪比大于$ -5 \, {\rm dB}$时SSAMP要明显优于OMP算法. SP算法在各种性能指标上均要差于SSAMP算法.因此SSAMP算法在综合性能上具有一定的优势.

5 SSAMP算法性能分析

由仿真结果图 6图 8可以看出SSAMP算法有3个特点: 1)在SNR<19 dB时, 其单信号估计的RMSE和运行时间恶化; 2)在SNR>5 dB时, 其单信号估计的RMSE为$0$; 3)在SNR<2 dB时, 其分辨能力大幅下降.下面将通过分析上述现象产生的原因进而分析SSAMP算法的性能.

5.1 SSAMP算法有效的条件

本节通过分析SNR<-19 dB时DOA估计的RMSE恶化与运行时间增加的原因来找出SSAMP算法有效的信噪比条件.

由AMP算法迭代公式(式(8)~式(11))可以看出所有包含噪声的项都位于阈值函数之中, 可以考虑分析在${\pmb x}$的每一维噪声对阈值函数的影响来解释算法失效的原因.假设噪声为高斯白噪声, 并记为$\varepsilon $.由于$\varepsilon $是一个随机变量, 因此我们分析它对阈值函数的期望造成的影响.设阈值$\tau > 0$, 那么根据式(10)则阈值函数的期望可以表示为

$ {\rm{E}}(\eta \left( {x + \varepsilon, \tau } \right)) = \int_{ - \infty }^\infty {\eta \left( {x + \varepsilon, \tau } \right)f\left( \varepsilon \right){\rm{d}}\varepsilon } $ (19)

其中, ${f\left( \varepsilon \right)}$表示噪声的分布函数, 其累积分布函数用${F\left( {\varepsilon} \right)}$表示.将式(19)积分展开可得

$ \begin{array}{l} {\rm{E}}(\eta \left( {x + \varepsilon, \tau } \right)) = 2x + \frac{\sigma }{{\sqrt {2\pi } }}\left( {{{\rm{e}}^{\frac{{ - {{(\tau - x)}^2}}}{{2{\sigma ^2}}}}} - {{\rm{e}}^{\frac{{ - {{(\tau + x)}^2}}}{{2{\sigma ^2}}}}}} \right)+ \\ \qquad \left( {\tau - x} \right)F\left( {\tau - x} \right) - \left( {\tau + x} \right)F\left( {\tau + x} \right) \end{array} $ (20)

为了进一步量化分析噪声对阈值函数的影响, 我们定义一个等效阈值的概念.

等效阈值. $\bar \tau = \sup \left( {{\rm{E}}\left( {\eta \left( {x, \tau } \right)} \right) < \epsilon } \right)$.其中, $\epsilon$表示一个接近零的较小的正数($10^{-4}$).

图 9展示了在不同信噪比条件下的等效阈值变化的情况.阈值函数的作用就是将小于阈值的量直接设置为$0$, 从而实现噪声的消除和迭代结果的稀疏化, 而噪声的出现使得阈值函数的这种作用变弱.当信噪比在$-20\, \rm{dB}$$-15\, \rm{dB}$时, 阈值函数作用下降明显.因此, 此时SSAMP算法的性能会明显下降.当信噪比小于$-19\, \rm{dB}$时, 利用AMP算法得到的子空间未必包含真实的非零元素解, 此时SSAMP算法得到的DOA估计结果误差将会大幅增加.同时根据文献[9]的结论, 此时AMP算法的结果趋近于最小$l_2$范数解.粗解的稀疏性大幅降低, 由此得到的子空间数量会大幅增加, 此时SSAMP算法的运行时间也会增加.

图 9 不同噪声水平下的等效阈值 Figure 9 Equivalent threshold under different noise levels
5.2 SSAMP算法精确重构的条件

精确重构是指算法重构的RMSE$=0$.通过寻找仿真1中SNR>5 dB时的DOA估计RMSE为$0$的原因来得到确重构的条件.

在搜索过程中我们使用残差${{res}} = {\pmb y}-A {\widehat {\pmb x}}$来判别最优解.根据仿真1的设定, 原信号中只有一个非零元素元素.假设非零元素真实位置为$i$, 那么有$A{\pmb x} = {x_i}{{A}_i}$, 其中, $x_i$表示了非零元素真实的幅值.而一元尝试解由式(21)求得, 其残差由式(22)求得.

$ \widehat x_t = {A}_t^*{\pmb y} $ (21)
$ {{res}}_{t} = \left\| {{\pmb y} - \widehat {x}_t{{A}_t} } \right\| $ (22)

${{\pmb \zeta }} = {x_i}({{A}_i} - {g_{ti}}{{A}_t})$表示了向量$ {\pmb x}_i={x_i}{{A}_i}$与向量${{A}_t}$正交的分量.同理设${{\pmb e}_{e \bot t}} = {\pmb e} - {{A}_t}({A}_t^*{\pmb e})$表示了噪声向量${\pmb e}$与向量${{A}_t}$正交的分量.结合式(21)与式(22)可得:

$ res_t = \left\{ \begin{array}{l} \left\| {{{\pmb \zeta}} + {{\pmb e}_{e \bot t}}} \right\|, \\ \left\| {{{\pmb e}_{e \bot i}}} \right\|, \end{array} \right.\begin{array}{*{20}{c}} {t \ne i}\\ {t = i} \end{array} $ (23)

${{res}_{t \ne i}}$根据三角不等式有:

$ res_{t \ne i} \ge \left\| {{{\pmb \zeta}}} \right\| - \left\| {{{\pmb e}_{e \bot t}}} \right\| $ (24)

若要RMSE$=0$必须有$res_{t \ne i} > {{res}_{t=i}}$

$ \left\| {{{\pmb \zeta}_{i \bot t}}} \right\| > \left\| {{{\pmb e}_{e \bot t}}} \right\| + \left\| {{{\pmb e}_{e \bot i}}} \right\| $ (25)

假设出现极端情况即${{\pmb e}}=-{\pmb \zeta}$时, 有$\left\| {{{\pmb e}_{e \bot t}}} \right\| = \left\| {\pmb e} \right\|$, $\left\| {{{\pmb e}_{e \bot i}}} \right\| = g_{ti} \left\| {\pmb e} \right\|$.此时式(25)变为

$ \frac{{\left| {{x_i}} \right|}}{{\left\| {\pmb e} \right\|}} > \frac{1+g_{ti}}{{\left\| {{{A}_i} - {g_{ti}}{{A}_t}} \right\|}} \approx \frac{1+g_{ti}}{{\sqrt {1 - g_{ti}^2} }} $ (26)

由此得到SSAMP算法RMSE$=0$的充分条件为

$ {\rm SNR} > 20 \lg \left(\frac{1+g_{ti}}{\sqrt{{M} ({1-g_{ti}^2})}}\right) $ (27)

代入具体的仿真参数$\max \{ {g_{ti}}\} = 0.9964$, 得到的信噪比应大于$12.67\, \rm{dB}$.

然而由仿真得到的结果为SNR>5 dB倒推得到的${\left\| {\pmb e} \right\|}$要比理论值大$2.24$倍.必须要考虑到${\pmb e}$是一个每一维独立同分布且均值为$0$、方差为$\sigma^2$$M$维高斯随机向量.因此${\rm E} \left({\pmb \zeta}^*{{\pmb e}}\right)=0$, 说明${\pmb \zeta}$${\pmb e}$在大多数情况下接近正交.而在上述假设出现的极端情况下, ${\pmb \zeta}$${\pmb e}$是平行的.同时又因为${\pmb \zeta}$${\pmb e}$的维数$M=30$, 所以在有限次的仿真实验中很难取得这种极端的情况.因此可以认为实验的结果与理论分析是一致的.

5.3 SSAMP算法信噪比与分辨能力的关系

分辨能力的下降主要是因为当噪声水平增加时, 二元真实解的残差大于一元解的残差.设真实的非零元素位置为${i_1}, {i_2}$, 有$A{\pmb x} = {x_{i_1}}{{A}_{i_1}}+{x_{i_2}}{{A}_{i_2}}$.一元解及其残差由式(21)和式(22)求得.根据仿真2的设置, 有${x_{i_1}}={x_{i_2}}=x$.用一个合成向量来表示这个二元真实解, 即${\pmb v} = x({{A}_{i_1}}+{{A}_{i_2}})$.则真实解的残差可以表示为$res_{v}=\left\| {{{\pmb e}_{e \bot v}}} \right\|$.这样处理后, 该问题变得与上一节的情况相似, 可以由式(26)直接得到一个相似的结果:

$ \frac{{\left\| {\pmb v} \right\|}}{{\left\| {\pmb e} \right\|}} \ge \frac{2}{{\left\| {{\pmb v} - {{A}_t}{g_{tv}}} \right\|}} $ (28)

其中, ${g_{tv}} = {\left\| {A}_t^*({{A}_{i_1}} + {{A}_{i_2}})\right\|}/{\left\| {{{A}_{i_1}} + {{A}_{i_2}}} \right\|}$.考虑到出现该现象时, $i_1, i_2$相距较近, 有${\left\| {\pmb v} \right\|} = \frac{1+g_{tv}}{\sqrt{2}}{\left\| {\pmb x} \right\|}$.因此式(28)可以变为

$ {\rm SNR} \ge 20{\rm lg} \left( \frac{1+g_{tv}}{\sqrt {2M(1-g_{tv}^2)}} \right) $ (29)

将最小可分辨角度为$1.125^\circ$时相应的参数$g_{tv}=0.9675$代入式(29), 可得对应的信噪比的下界是$0.04\, {\rm dB}$.而图 8中与该最小分辨角度对应的信噪比为$-1\, {\rm{dB}}$.因此可以认为上述关于信噪比与最小可分辨角度的分析是有效的.

5.4 信号数未知时SSAMP算法的扩展

本文所提的SSAMP算法是在信号数$K$已知的基础上设计的.然而在实际的DOA估计应用中, 信号数已知的前提有时难以满足.因此本节将讨论SSAMP算法能否扩展应用于信号数未知的情况.

本文所提算法使用信号数的地方有两处:确定初始阈值和求精确解.在确定初始阈值时需要使用来波信号的幅值水平.当来波信号的幅值水平未知时, 使用$x={\left\| y \right\|}/{K}$进行估计以便确定初始阈值.若信号数未知时, 没有办法按照第3.2节中的建议选择初始阈值.为了保证幅值较小的信号能被正确重构, 需要选择足够小的初始阈值.根据图 5所示, 这样做会大幅降低AMP算法收敛的速度, 增加程序运行时间.在求取精确解的过程中若信号数$K$已知, 每个子集中最多可能存在的信号数$N_{sol}=K-N_R+1$就已知, 这会使算法计算量较低.当信号数$K$未知时, 如果存在多个信号位于一个子空间中的情况, 为了在所有子集中遍历全部可能的解, 会使算法计算量大大增加.

综上所述, 本文所提SSAMP算法在信号数未知时也可以扩展使用.但此时算法的运算量会大幅增加, 这意味着算法可以通过增加计算量来弥补信号数未知的信息损失.

6 结论

本文通过分析AMP算法得到的粗解与感知矩阵高相干列的关系, 提出了一种SSAMP算法.该算法利用分步求解的思想, 首先利用AMP算法求出粗解, 然后在粗解所对应的子空间中进一步寻找精确解.该算法保留了AMP算法计算复杂度低以及相对于其他压缩感知类算法具有抗噪声能力强的优势, 同时又解决了AMP算法在感知矩阵相干系数较高时无法精确求解的问题.仿真实验通过对比SSAMP算法与一些经典的压缩感知算法在单信号RMSE、运行时间和最小可分辨角度上的性能, 验证了SSAMP算法在综合性能上要优于其他算法.最后本文在理论上分析了SSAMP算法有效范围、精确重构条件以及信噪比与最小分辨角度的关系, 并讨论了信号数未知时算法的扩展.在下一步的工作中, 可以考虑将目标状态估计与本文所提算法相结合, 使用目标状态估计的结果作为先验信息来进一步降低单快拍DOA估计算法的计算复杂度.

参考文献
1
Schmidt R. Multiple emitter location and signal parameter estimation. IEEE Transactions on Antennas and Propagation, 1986, 34(3): 276-280. DOI:10.1109/TAP.1986.1143830
2
Roy R, Paulraj A, Kailath T. ESPRIT-a subspace rotation approach to estimation of parameters of cisoids in noise. IEEE Transactions on Acoustics, Speech, and Signal Processing, 1986, 34(5): 1340-1342. DOI:10.1109/TASSP.1986.1164935
3
Donoho D L. Compressed sensing. IEEE Transactions on Information Theory, 2006, 52(4): 1289-1306. DOI:10.1109/TIT.2006.871582
4
Malioutov D, Cetin M, Willsky A S. A sparse signal reconstruction perspective for source localization with sensor arrays. IEEE Transactions on Signal Processing, 2005, 53(8): 3010-3022. DOI:10.1109/TSP.2005.850882
5
Kim J M, Lee O K, and Ye J C. Compressive MUSIC:revisiting the link between compressive sensing and array signal processing. IEEE Transactions on Information Theory, 2012, 58(1): 278-301. DOI:10.1109/TIT.2011.2171529
6
Wang Xiu-Hong, Mao Xing-Peng, Zhang Nai-Tong. Single-snap DOA estimation based on compressed sensing in pulse compression radar system. Systems Engineering and Electronics, 2014, 36(9): 1737-1743.
( 王秀红, 毛兴鹏, 张乃通. 基于CS的脉冲压缩雷达单快拍DOA估计. 系统工程与电子技术, 2014, 36(9): 1737-1743. DOI:10.3969/j.issn.1001-506X.2014.09.11)
7
Liu J, Mallick M, Han C Z, Yao X H, Lian F. Similar sensing matrix pursuit:an efficient reconstruction algorithm to cope with deterministic sensing matrix. Signal Processing, 2014, 95: 101-110. DOI:10.1016/j.sigpro.2013.08.009
8
Liu J, Mallick M, Lian F, Han C Z, Sheng M X, Yao X H. General similar sensing matrix pursuit:An efficient and rigorous reconstruction algorithm to cope with deterministic sensing matrix with high coherence. Signal Processing, 2015, 114: 150-163. DOI:10.1016/j.sigpro.2015.03.002
9
Donoho D L, Maleki A, Montanari A. Message passing algorithms for compressed sensing: Ⅰ. Motivation and construction. In: Proceedings of the 2010 Information Theory Workshop on Information Theory. Cairo, Egypt: IEEE, 2010. 1-5
10
Donoho D L, Maleki A, Montanari A. Message passing algorithms for compressed sensing: Ⅱ. Analysis and Validation. In: Proceedings of the 2010 IEEE Information Theory Workshop on Information Theory. Cairo, Egypt: IEEE, 2010. 6-10 http://ieeexplore.ieee.org/xpls/icp.jsp?arnumber=5503228
11
Rangan S. Generalized approximate message passing for estimation with random linear mixing. In: Proceedings of the 2011 IEEE International Symposium on Information Theory Proceedings (ISIT). St Petersburg, Russia: IEEE, 2011. 2168-2172 http://arxiv.org/abs/1010.5141
12
Vila J P, Schniter P. Expectation-maximization Gaussian-mixture approximate message passing. IEEE Transactions on Signal Processing, 2013, 61(19): 4658-4672. DOI:10.1109/TSP.2013.2272287
13
Tan J, Ma Y T, Baron D. Compressive imaging via approximate message passing with image denoising. IEEE Transactions on Signal Processing, 2015, 63(8): 2085-2092. DOI:10.1109/TSP.2015.2408558
14
Ren Yue-Mei, Zhang Yan-Ning, Li Ying. Advances and perspective on compressed sensing and application on image processing. Acta Automatica Sinica, 2014, 40(8): 1563-1575.
( 任越美, 张艳宁, 李映. 压缩感知及其图像处理应用研究进展与展望. 自动化学报, 2014, 40(8): 1563-1575.)
15
Ziniel J, Schniter P. Dynamic compressive sensing of time-varying signals via approximate message passing. IEEE Transactions on Signal Processing, 2013, 61(21): 5270-5284. DOI:10.1109/TSP.2013.2273196
16
Orlando D, Venturino L, Lops M, Ricci G. Track-before-detect strategies for STAP radar. IEEE Transactions on Signal Processing, 2010, 58(2): 933-938. DOI:10.1109/TSP.2009.2032991
17
Richards M A. Fundamentals of Radar Signal Processing. New York: McGraw-Hill, 2005, 88-92.
18
Maleki A, Montanari A. Analysis of approximate message passing algorithm. In: Proceedings of Information Sciences and Systems. Princeton, USA: IEEE, 2010. 1-7 http://ieeexplore.ieee.org/xpls/icp.jsp?arnumber=5464887
19
Elad M. Sparse and Redundant Representations. New York: Springer, 2010, 65-68.
20
Wei D, Milenkovic O. Subspace pursuit for compressive sensing signal reconstruction. IEEE Transactions on Information Theory, 2009, 55(5): 2230-2249. DOI:10.1109/TIT.2009.2016006