CN104020438B - 基于稀疏表示的波达方向角估计方法 - Google Patents

Abstract

本发明公开了一种基于稀疏表示的波达方向角估计方法，主要解决现有同类技术运算量大和角度分辨率低，造成目标侦察和无源定位估计速度慢误差大的问题。其实现步骤是：1)获取阵列输出信号，并计算信号协方差矩阵R；2)利用协方差矩阵R下三角的元素构建稀疏重构向量；3)将空域进行网格划分处理，构造超完备基；4)根据稀疏重构向量和超完备基的稀疏表示关系，将波达方向角估计问题转化为求解约束优化方程的问题；5)采用凸优化方法求解约束优化方程得到最优估计向量；6)利用最优估计向量与空域角度的一一对应关系绘制幅度谱图，得到波达方向角度值。本发明降低了目标侦察和无源定位的运算速度和估计误差，可用于目标侦察和无源定位。

Description

基于稀疏表示的波达方向角估计方法

技术领域

本发明属于信号处理技术领域，特别涉及一种波达方向角估计方法，可用于目标侦察与无源定位。

背景技术

波达方向角DOA估计是利用处于空间不同位置的天线阵列接收多个不同方向的信号源发出的信号，运用现代信号处理方法快速准确的获得信号源的方向，在雷达、声纳、无线通信等领域具有重要应用价值。针对该问题构建的模型中，出现较早且应用较为广泛的是基于子空间的模型，之前的大部分DOA估计方法都是利用该模型生成的。近年来，由Donoho等提出的压缩感知理论为DOA估计问题提供了一种新思路，从而产生出一类基于稀疏表示模型的DOA估计方法。

目前，基于稀疏信号表示的DOA估计方法中最经典的是L1_SVD方法。L1_SVD方法利用阵列接收数据奇异值分解(SVD)得到的信号子空间构造稀疏表示模型，然后通过二阶锥规划对L1范数约束模型进行求解。但这种方法运算量较大，且在低信噪比情况下，角度分辨率不理想。在实际应用中，由于目标侦察与无源定位均需要在角度估计的基础上进行，使用以上方法将造成目标侦察和无源定位反应速度慢和估计误差较大的问题。

发明内容

本发明的目的在于针对上述已有技术的不足，提出一种基于稀疏表示的波达方向角估计方法，以在降低运算量的情况下，提高目标侦察和无源定位的发现概率和角度分辨率，避免因角度估计误差引起目标侦察失误。

为实现上述目的，本发明的实现步骤包括如下：

一种基于稀疏表示的波达方向角度估计方法，包括以下步骤：

1)采用M个天线接收机形成均匀线性阵列，并假设有K个信号入射到该均匀线性阵列，各天线接收机间距均为d，M≥2，K≥1，0＜d≤λ/2，λ为入射窄带信号波长；

2)根据阵列天线的输出信号X(t)，计算信号协方差矩阵R：

R＝E{X(t)X^H(t)}，

其中，E(·)表示求数学期望，(·)^H表示共轭转置运算；

3)利用协方差矩阵构造稀疏重构向量

根据协方差矩阵R中对角线方向上同一条斜线上的各元素的值近似相等，且协方差矩阵R是共轭对角矩阵的特点，对协方差矩阵下半角的元素按对角线求平均，得到一个M-1维稀疏重构向量 其中，(·)^T表示转置运算，表示稀疏重构向量中的元素，其表达式为：

$\hat{z}\left(m\right)=\frac{1}{M-m}\underset{p=m+1}{\overset{M}{\mathrm{\Σ}}}R(p,p-m),m=\mathrm{1,2},...,M-1,$

其中，R(p,p-m)表示协方差矩阵R的第p行，p-m列的元素；

4)构造超完备基并定义一个空域稀疏向量u：

4a)根据信号源的空域稀疏特性，采用空间网格划分方法，将观测空域[-90°,90°]等间隔划分成Q个角,定义为波达方向角范围θ＝[θ₁,θ₂,…,θ_q,…,θ_Q]，Q＞＞M；

4b)构造一个信号稀疏化后对应的(M-1)×Q维的超完备基

$\tilde{A}\left(\mathrm{\θ}\right)=[\mathrm{\α}\left({\mathrm{\θ}}_1\right),...,\mathrm{\α}\left({\mathrm{\θ}}_q\right),...,\mathrm{\α}\left({\mathrm{\θ}}_Q\right)],$

其中，α(θ_q)＝(exp(-j2πdsinθ_q/λ),…,exp[-j2πd(M-1)sinθ_q/λ])^T，q＝1,2,…,Q；

4c)定义一个Q×1维的空域稀疏向量u，u＝[u₁,u₂…,u_q,…,u_Q]^T，u中元素均为未知变量；

5)通过稀疏重构获得空域稀疏向量u的最优估计

5a)根据稀疏表示的思想，将波达方向角估计问题转化为求解约束优化方程式：

$\underset{\hat{u}}{\min }{\left|\right|\hat{u}\left|\right|}_{1}s.t.{\left|\right|\hat{z}-\tilde{A}\left(\mathrm{\θ}\right)\hat{u}\left|\right|}_2\≤\mathrm{\ϵ},$

其中，||·||₁和||·||₂分别表示求向量1-范数和向量2-范数，s.t.表示约束关系，ε表示误差的允许值，ε的取值为现实环境中的噪声标准差；

5b)利用凸优化方法求解上述优化方程，得到稀疏向量的最优解

6)以波达方向角范围θ＝[θ₁,θ₂,…,θ_q,…,θ_Q]的值为x轴坐标，以稀疏向量的最优解的幅度值为y轴坐标，绘制幅度谱图，从该幅度谱图中按照从高到低的顺序寻找幅值较大的前K个谱峰，这些谱峰的峰值点所对应的x轴坐标即为所求的波达方向角度值。

本发明与现有技术相比具有以下优点：

1)本发明采用稀疏表示的思想将波达方向角度估计问题转化成稀疏重构问题，是新理论技术与传统问题的结合，利用信号源在空域稀疏的特点进行建模，避免了传统算法的角度搜索或角度匹配过程，同时突破了阵列分辨率的瑞丽限，提高了测角精度。

2)本发明根据协方差矩阵的结构特点，采用其下三角元素构建稀疏表示模型，通过该过程将多测量矢量MMV问题转化为单测量矢量SMV模型，在稀疏重构过程中大大降低了运算量；该模型的构建方法起到抑制噪声信息的作用，因而提高了低信噪比下目标侦察和无源定位的发现概率和精度。

3)本发明无需事先估计入射信号的个数，同时可用于处理相干信号源和非相干信号源，在现实环境中具有更实际的应用价值。

附图说明

图1是本发明的实现流程图；

图2是本发明与现有两种波达方向角估计方法的运算速度对比图；

图3是本发明与现有两种波达方向角估计方法的角度分辨率对比图。

具体实施方式

参照图1，本发明的实现步骤如下：

步骤1：利用天线接收机形成均匀线阵。

每隔间距d放置1个天线接收机，共放置M个，形成一个均匀线性天线阵列，假设有K个信号源入射到该天线阵列，其中，M≥2，K≥1，0＜d≤λ/2，λ为入射窄带信号波长。

步骤2：根据阵列天线的输出信号X(t)，计算其协方差矩阵R:

R＝E{X(t)X^H(t)}，

其中，E(·)表示求数学期望，(·)^H表示矩阵共轭转置运算。

步骤3：利用协方差矩阵R构造稀疏重构向量

将协方差矩阵R第p行，第q列的元素R(p,q)表示为：

$R(p,q)=\left\{\begin{array}{ccc}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{\mathrm{\σ}}_k^2+{\mathrm{\σ}}_n^2& \mathrm{if}& p=q\\ \underset{k=1}{\overset{K}{\mathrm{\Σ}}}{\mathrm{\σ}}_k^2\exp [j(p-q)2\mathrm{\π}d\sin {\mathrm{\θ}}_k/\mathrm{\λ}]& \mathrm{if}& p\≠q\end{array}\right.,$

其中，p,q∈[1,2,…,M]，分别表示信号能量和噪声能量，θ_k表示入射信号的方向角，根据协方差矩阵R中对角线方向上同一条斜线上的各元素的值近似相等，且协方差矩阵R是共轭对角矩阵的特点，对协方差矩阵下半角的元素按对角线求平均，得到一个M-1维稀疏重构向量

$\hat{z}={[\hat{z}\left(1\right),\hat{z}\left(2\right),...,\hat{z}\left(m\right),...,\hat{z}(M-1)]}^T,$

其中，(·)^T表示转置运算，表示稀疏重构向量中的元素，其表达式为：

$\hat{z}\left(m\right)=\frac{1}{M-m}\underset{p=m+1}{\overset{M}{\mathrm{\Σ}}}R(p,p-m),m=\mathrm{1,2},...,M-1,$

其中，R(p,p-m)表示协方差矩阵R的第p行，p-m列的元素。

步骤4：构造超完备基

4a)根据信号源所具有的空域稀疏特性，对观测空域进行空间网格划分处理，即将观测空域[-90°,90°]等间隔划分成Q个角，θ＝[θ₁,θ₂,…,θ_q,…,θ_Q]，θ表示波达方向角范围，θ_q为第q个角度区间，q＝1,2,…,Q，Q＞＞M；网格划分间隔的取值根据期望达到的角度估计精度进行设定，网格划分间隔越小，则最终得到的角度估计值精度越高。

4b)构建一个信号稀疏化后对应的(M-1)×Q维的超完备基

$\tilde{A}\left(\mathrm{\θ}\right)=[\mathrm{\α}\left({\mathrm{\θ}}_1\right),\mathrm{\α}\left({\mathrm{\θ}}_2\right),...,\mathrm{\α}\left({\mathrm{\θ}}_q\right),...,\mathrm{\α}\left({\mathrm{\θ}}_Q\right)],$

其中，α(θ_q)表示超完备基的第q列，α(θ_q)为一个M-1维的列向量：α(θ_q)＝(exp(-j2πdsinθ_q/λ),…,exp[-j2πd(M-1)sinθ_q/λ])^T，q＝1,2,…,Q。

步骤5：定义一个Q×1维的空域稀疏向量u。

假设每个角度区间θ_q都对应一个入射信号，这样就构造一个Q×1维的空域稀疏向量u＝[u₁,u₂…,u_q,…,u_Q]^T，u中元素均为未知变量。

步骤6：稀疏重构获得空域稀疏向量u的最优估计

6a)根据稀疏表示的思想，稀疏重构向量为：

$\hat{z}=\tilde{A}\left(\mathrm{\θ}\right)\hat{u}+\mathrm{\δ}$

其中，δ表示噪声修正量；

进一步将波达方向角估计问题可转化为求解以下的约束优化方程式：

$\underset{\hat{u}}{\min }{\left|\right|\hat{u}\left|\right|}_{1}s.t.{\left|\right|\hat{z}-\tilde{A}\left(\mathrm{\θ}\right)\hat{u}\left|\right|}_2\≤\mathrm{\ϵ}$

其中，||·||₁和||·||₂分别表示求向量的1-范数和2-范数，s.t.表示约束关系，ε表示误差的允许值，ε的取值为现实环境中的噪声标准差。

6b)采用凸优化方法求解上述约束优化方程式，凸优化是一种比较特殊的优化，是指目标函数和约束函数均为凸函数的优化问题，凸优化问题有一套非常完备的解决算法，在此采用现有针对凸优化问题的软件包CVX(Grant M,Boyd S.CVX:Matlabsoftware for disciplined convex programming[J].2008[Online]Available:http://stanfordedu/～boyd/cvx)来求解，通过该方法能够快速地得到空域稀疏向量的最优估计

步骤7：根据得到的最优估计绘制幅度谱图，得到波达方向角度值。

以波达方向角范围θ＝[θ₁,θ₂,…,θ_q,…,θ_Q]的值为x轴坐标，以最优估计的幅度值为y轴坐标，绘制幅度谱图，从该幅度谱图中按照从高到低的顺序寻找幅值较大的前K个谱峰，这些谱峰的峰值点所对应的x轴坐标即为所求的波达方向角度值。

本发明的效果可通过以下仿真说明：

1.仿真条件与方法：

天线阵列为均匀线性阵列，各天线接收机的间距d均等于半波长。信噪比SNR为0dB，采样快拍数L等于500，观测空域角度范围为[-90°，90°]，空间网格划分间隔为1°。

2.仿真内容与结果

仿真1，假设有两个窄带信号分别以角度15°和25°入射到均匀线性阵列，均匀线性阵列的天线接收机个数由3个增加到12个，采用本发明和现有的L1_SVD、L1_SRACV算法在不同实验条件下分别进行100次独立的波达方向角估计实验，并统计运算时间，三种方法的运算时间随天线接收机个数变化的对比图如图2所示。图2中横坐标表示天线接收机个数，纵坐标表示运算时间。

从图2可以看出，本发明与现有L1_SVD、L1_SRACV算法相比，大幅度降低了角度估计的运算量，即使天线接收机个数增加，运算时间仍基本保持稳定不变。

仿真2，采用8个天线接收机的均匀线性阵列，假设有两个窄带信号分别以角度15°和15°+Δθ入射到均匀线性阵列，其中入射信号角度间隔Δθ取值由1°增加到15°。采用本发明和现有的L1_SVD、L1_SRACV算法分别进行100次独立的波达方向角估计实验，分别计算角度估计正确的发现概率，三种方法的发现概率对比图如图3所示，横坐标表示入射信号的角度间隔，纵坐标表示发现概率。

从图3中可以看出，本发明比现有的L1_SVD、L1_SRACV算法均具有更高的角度分辨率。

综上，本发明在降低了角度估计运算量的同时，还提高了算法的角度分辨率，保证了目标侦察和无源定位的快速反应，避免了因信息滞后引起的参数估计误差。

Claims (1)

1.一种基于稀疏表示的波达方向角估计方法，包括以下步骤：

1)采用M个天线接收机形成均匀线性阵列，并假设有K个信号入射到该均匀线性阵列，各天线接收机间距均为d，M≥2，K≥1，0<d≤λ/2，λ为入射窄带信号波长；

2)根据阵列天线的输出信号X(t)，计算信号协方差矩阵R：

R＝E{X(t)X^H(t)}，

其中，E(·)表示求数学期望，(·)^H表示共轭转置运算；

3)利用协方差矩阵构造稀疏重构向量

根据协方差矩阵R中对角线方向上同一条斜线上的各元素的值近似相等，且协方差矩阵R是共轭对角矩阵的特点，对协方差矩阵下半角的元素按对角线求平均，得到一个M-1维稀疏重构向量 其中，(·)^T表示转置运算，表示稀疏重构向量中的元素，其表达式为：

$\hat{z}\left(m\right)=\frac{1}{M-m}\underset{p=m+1}{\overset{M}{\&Sigma;}}R(p,p-m),m=1,2,...,M-1,$

其中，R(p,p-m)表示协方差矩阵R的第p行，p-m列的元素；

4)构造超完备基并定义一个空域稀疏向量u：

4a)根据信号源的空域稀疏特性，采用空间网格划分方法，将观测空域[-90°，90°]等间隔划分成Q个角,定义为波达方向角范围θ＝[θ₁,θ₂,…,θ_q,…,θ_Q]，Q>>M；

4b)构造一个信号稀疏化后对应的(M-1)×Q维的超完备基

$\tilde{A}\left(\mathrm{\&theta;}\right)=\&lsqb;\mathrm{\&alpha;}\left({\mathrm{\&theta;}}_1\right),...,\mathrm{\&alpha;}\left({\mathrm{\&theta;}}_q\right),...,\mathrm{\&alpha;}\left({\mathrm{\&theta;}}_Q\right)\&rsqb;,$

其中，α(θ_q)＝(exp(-j2πdsinθ_q/λ),…,exp[-j2πd(M-1)sinθ_q/λ])^T，q＝1,2,…,Q；

4c)定义一个Q×1维的空域稀疏向量u，u＝[u₁,u₂…,u_q,…,u_Q]^T，u中元素均为未知变量；

5)通过稀疏重构获得空域稀疏向量u的最优解

5a)根据稀疏表示的思想，将波达方向角估计问题转化为求解约束优化方程式：

$\begin{array}{ccc}\underset{\hat{u}}{min}{||\hat{u}||}_1& s.t.& {||\hat{z}-\tilde{A}\left(\mathrm{\&theta;}\right)\hat{u}||}_2\&le;\mathrm{\&epsiv;}\end{array},$

其中，||·||₁和||·||₂分别表示求向量1-范数和向量2-范数，s.t.表示约束关系，ε表示误差的允许值，ε的取值为现实环境中的噪声标准差；

5b)利用凸优化方法求解上述优化方程，得到稀疏向量的最优解

6)以波达方向角范围θ＝[θ₁,θ₂,…,θ_q,…,θ_Q]的值为x轴坐标，以稀疏向量的最优解的幅度值为y轴坐标，绘制幅度谱图，从该幅度谱图中按照从高到低的顺序寻找幅值较大的前K个谱峰，这些谱峰的峰值点所对应的x轴坐标即为所求的波达方向角度值。