CN104539340A - 一种基于稀疏表示和协方差拟合的稳健波达角估计方法 - Google Patents

Abstract

本发明提出了一种基于稀疏表示和协方差拟合的稳健波达角估计方法。主要包括根据天线阵列接收信号建立稀疏空间谱表示模型，并将模型误差进行参数化表示，根据协方差拟合准则建立最优化问题。因为所得问题为非凸优化，因此通过等价转换、增加参数、分步求解对问题进行转化和求解。首先不考虑模型误差，原问题可简化为凸优化问题。利用现有方法快速求解该凸优化问题，得到初始解。再迭代求解原问题，估计模型误差参数并更新初始估计。本发明能够以低复杂度获得精确的DOA估计。

Description

一种基于稀疏表示和协方差拟合的稳健波达角估计方法

技术领域

本发明属于信号处理技术领域，尤其涉及达波角方向估计，可用于无源定位和目标检测等。

背景技术

波达方向(DOA)估计是利用处于空间不同位置的天线阵列接收多个不同方向的信号源发出的信号，运用现代信号处理方法快速准确的获得信号源的方向，在雷达、声纳、无线通信等领域具有重要应用价值。

传统的DOA估计技术包括基于子空间的方法，例如MUSIC方法、最大似然估计。近年来，信号稀疏恢复和稀疏表示的发展为DOA估计技术提供了新思路。信号稀疏表示的主要思路是利用精细网格覆盖待估参数的取值空间，再找出距离参数真实值最近的格点，从而实现参数估计。目前，基于稀疏表示的DOA估计方法中较常用的是L1_SVD方法。L1_SVD方法对阵列接收信号进行稀疏空间谱表示，得到含有L1范数约束的二阶锥规划(SOCP)问题，并利用奇异值分解(SVD)来降低所求问题的维数。但当方向估计精度要求较高时，这种方法仍然面临求解高维SOCP、L1范数罚因子选择等问题。为了避开这些问题，出现了基于稀疏自相关拟合的迭代估计方法(Sparse Iterative Co-variance-based Estimation,SPICE)。此类方法避免了求解SOCP问题，复杂度低。但是该方法要求信号源位于稀疏空间谱表示时所采用的网格上，其估计精确度受到了网格密度的限制。

发明内容

技术问题：本发明针对已有技术的不足，将信号源真实位置和网格格点间的误差参数化、线性化，并结合协方差拟合准则提出一种迭代的方向估计算法。此方法是现有技术的改进和扩展，能够以低复杂度获得DOA在连续取值空间上的估计值，提高了估计精度。

技术方案：为实现上述目标，本发明的主要步骤如下：

一种基于稀疏表示和协方差拟合准则的稳健波达角估计方法，包括以下步骤：

1)采用N个全向天线形成均匀线性阵列，对空间中K个窄带信号源进行方向估计。信号源可能方向集合为Ω，其上覆盖候选方向网格其中θ_k表示通用方向参数。根据阵列天线的输入信号建立接收信号模型：

$y\left(t\right)=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}(a_k+e_k\left({\widetilde{\mathrm{\θ}}}_k\right))s_k\left(t\right)+n\left(t\right),t=\mathrm{1,2,3},....T$

其中a_k是对应于θ_k的导向矢量，表示建模误差，为相对于θ_k的实际位置，s_k(t)是未知的源信号，n(t)是噪声，T为总的信号数目。当稀疏空间谱表示采用的网格的密度足够大时e_k(θ_k)可以线性近似为：

$e_k\left({\widetilde{\mathrm{\θ}}}_k\right)=e_k{\∈}_k$

其中为网格方向偏差，其均匀的分布在[-l_g,l_g]中，2l_g为网格的步长。

2)根据天线阵列的输出信号计算信号协方差矩阵R:

$R=E\left[y\left(t\right)y{\left(t\right)}^H\right]=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{p}_k{(a_k+e_k{\∈}_k)(a_k+e_k{\∈}_k)}^H+\mathrm{\σI}$

其中，E(·)表示求数学期望,(·)^H表示共轭转置运算，其中，E＝[e₁,e₂,…e_k]，A＝[a₁,a₂,…a_k]，P＝diag([p^T,σ]^T)，p＝[p₁,p₂,...p₂]^T。

3)基于参数化的协方差矩阵构造参数估计的最优化问题：

3a)利用步骤2)中协方差基于p,σ,∈的参数化表示得到方向估计的加权协方差拟合准则：

$f(p,\mathrm{\σ},\∈)=\mathrm{tr}\left({\hat{R}R}^{-1}\right)+\mathrm{tr}\left({\hat{R}}^{-1}R\right)+{\left|\right|\∈\left|\right|}^2/{\mathrm{\ω}}_e$

其中σ表示噪声功率，∈表示网格误差，ω_e为∈的协方差，∈＝[∈₁,∈₂,…∈_k]^T。

令并假设和R的逆都存在，并结合参数自身约束，协方差拟合准则可等价转换为如下优化问题：

$\underset{p_k,\mathrm{\σ},{\∈}_k}{\min }\mathrm{tr}\left({\hat{R}R}^{-1}\right)+{\left|\right|\∈\left|\right|}^2/{\mathrm{\ω}}_e$

$s.t.\mathrm{tr}\left({\hat{R}}^{-1}R\right)=N$

p_K≥0,|∈_k|≤l_g k＝1,2,…K

σ≥0

其中l_g为网格步长的1/2，N为天线数，tr(·)表示求矩阵的迹。

3b)由于上一步中建立的最优化问题为非凸优化问题，不易求解。因此对问题进行转换，引入变量C，构造如下最优化问题：

$s.t.B^{H}C={\hat{R}}^{1/2}$

其中B＝[A,E,I]^H， $D=\left[\begin{array}{ccc}I& \mathrm{\γ}& 0\\ 0& 0& I\end{array}\right],$ γ＝diag(∈)，表示矩阵的违逆。

3c)求解上述二次凸优化问题，并带回目标函数可得最优值为：此值恰好为3a)中的最优化问题的目标函数。所以可以通过求解如下最优化问题来得到波达角估计的最优值：

$s.t.B^{H}C={\hat{R}}^{1/2}$

$\mathrm{tr}\left({\hat{R}}^{-1}R\right)=N$

p_K≥0,|∈_k|≤l_g,k＝1,2,…K

σ≥0

4)上述3c)中的最优化问题为非凸优化问题，求解复杂度极高。因此我们将问题的求解分为两步：首先假设模型误差为0，从而将该问题转化为易于求解的凸优化问题，并利用闭式解迭代的方法求得最优解。然后以该最优解为初始点，求解模型误差参数并迭代更新初始估计。

4a)在不考虑建模误差的情况下利用闭式解迭代的方法求解p,σ,C的初始值。将初始值代入原优化问题得到如下关于∈估计的SOCP问题：

$\underset{\∈,t,f}{\min }f$

$\begin{array}{cc}s.t.f=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{t}_k+l& g\≥\left|\right|[{\∈}^T/\sqrt{{\mathrm{\ω}}_e},f]\left|\right|\end{array}$

$t_k\≥\frac{\left|\right|c_k\left|\right|\left|\right|{\hat{a}}_k+{\∈}_k{\tilde{e}}_k\left|\right|}{\sqrt{N}}$

-l_g≤∈_k≤l_g k＝1,…K

其中 $l=\sqrt{\mathrm{tr}\left({\hat{R}}^{-1}\right)c_{\mathrm{\σ}}/N},{\tilde{a}}_k={\hat{R}}^{-1/2}{a}_k,{\tilde{e}}_k={\tilde{R}}^{-1/2}{e}_k.$

4b)根据上一步求解得到的∈更新。p,σ将C，∈代入原优化问题得到包含未知数p,σ的优化问题，并求解得：

$p_k=\frac{\sqrt{N}(1+{\widehat{\∈}}_k{\∈}_k)\left|\right|c_k\left|\right|}{(1+{\∈}_k^2)\left|\right|{\hat{a}}_k+{\∈}_k{\tilde{e}}_k\left|\right|\mathrm{\ρ}},k=\mathrm{1,2},...K$

$\mathrm{\σ}=\sqrt{N{c}_{\mathrm{\σ}}/\mathrm{tr}\left({\hat{R}}^{-1}\right)/\mathrm{\ρ}}$

$\mathrm{\ρ}=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\frac{(1+{\widehat{\∈}}_k{\∈}_k)\left|\right|c_k\left|\right|\left|\right|{\hat{a}}_k+{\∈}_k{\tilde{e}}_k\left|\right|}{\sqrt{N}(1+{{\∈}_k}^2)}+l$

根据上式更新p_k，σ的值。

4c)根据上述步骤中求得的∈,p,σ可进一步更新C：

$c_k=p_k{\hat{R}}^{1/2}{R}^{-1}(a_k+{\∈}_k{e}_k),k=1,...K$

$[c_{2K+1},...,c_{2K+N}]=\mathrm{\σ}{\hat{R}}^{1/2}{R}^{-1}$

其中C＝[c₁,c₂,...c_2K+N]^H

5)由求得的C值重复上述4)中的三个步骤直到最大迭代次数或者满足迭代收敛条件最终求得最优解。

有益效果：本发明采用的技术方案通过将稀疏空间谱表示误差参数化、线性化和协方差拟合准则相结合，得到稳健的基于协方差的连续稀疏迭代估计方法。由于本发明中的方案考虑到了建模误差，能够获得DOA在连续空间上的估计值，以提高估计精度。

附图说明

图1为本发明的流程图。

图2为是本发明与现有波达方向估计方法的性能对比图。

具体实施方式

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

步骤1：采用N个全向天线接收机形成均匀线性阵列，对空间中K个窄带信号源进行方向估计。信号源可能方向集合为Ω，其上覆盖候选方向网其中θ_k表示通用方向参数。

步骤2根据天线阵列的输出信号计算信号协方差矩阵R的理论表达式子和实际采样值

$R=E\left[y\left(t\right)y{\left(t\right)}^H\right]=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{p}_k{(a_k+e_k{\∈}_k)(a_k+e_k{\∈}_k)}^H+\mathrm{\σI}$

$\hat{R}=\underset{t=1}{\overset{T}{\mathrm{\Σ}}}y\left(t\right)y^H\left(t\right)/T$

步骤3：利用协方差矩阵构造参数估计的优化问题：

利用步骤2)中协方差关于p,σ,∈的参数化表示得到方向估计的加权协方差拟合准则并结合参数自身约束建立优化问题。由于该优化问题为非凸优化，因此引入参数C构造如下最优化问题：

$s.t.B^{H}C={\hat{R}}^{1/2},\mathrm{tr}\left({\hat{R}}^{-1}R\right)=N,$

p_K≥0,|∈_k|≤l_g,k＝1,2,…K

σ≥0

步骤4：上述步骤3中的最优化问题为非凸优化问题，求解复杂度极高。因此我们将问题的求解分为两步：首先假设模型误差为0，从而将该问题转化为易于求解的凸优化问题，并利用闭式解迭代的方法求得最优解；然后以该最优解为初始点，求解模型误差参数并迭代更新初始估计值。

4a)在不考虑建模误差的情况下利用闭式解迭代的方法求解p,σ,C的初始值。随机产生p,σ的初始值，令

$c_{\mathrm{\σ}}={\left|\right|{\mathrm{\σ}}^2{\hat{R}}^{1/2}{R}^{-1}\left|\right|}^2,l=\sqrt{\mathrm{tr}\left({\hat{R}}^{-1}\right)c_{\mathrm{\σ}}/N}.$

给定∈_k＝0，迭代求解p_k,σ

$c_k=p_k{\hat{R}}^{1/2}{R}^{-1}{a}_k,\mathrm{\ρ}=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\frac{\left|\right|c_k\left|\right|\left|\right|{\tilde{a}}_k\left|\right|}{N}+l$

$p_k=\sqrt{N}\left|\right|c_k\left|\right|/\left(\right|\left|{\tilde{a}}_k\right|\left|\mathrm{\ρ}\right),\mathrm{\σ}=\sqrt{c_{\mathrm{\σ}}/\mathrm{tr}\left({\hat{R}}^{-1}\right)}/\mathrm{\ρ}$

4b)将4a)中求解的p,σ,C的初始值代入原优化问题，得到如下关于∈估计的SOCP问题：

$\underset{\∈,t,f}{\min }f$

$\begin{array}{cc}s.t.f=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{t}_k+l& g\≥\left|\right|[{\∈}^T/\sqrt{{\mathrm{\ω}}_e},f]\left|\right|\end{array}$

$t_k\≥\frac{\left|\right|c_k\left|\right|\left|\right|{\hat{a}}_k+{\∈}_k{\tilde{e}}_k\left|\right|}{\sqrt{N}}$

-l_g≤∈_k≤l_g k＝1,…K

4c)根据上一步求得的∈更新p,σ。将C，∈代入原优化问题得到包含未知数p,σ的优化问题，并根据柯西-施瓦茨不等式只有当下式成立时取得等号：

$\mathrm{\ρ}=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\frac{(1+{\widehat{\∈}}_k{\∈}_k)\left|\right|c_k\left|\right|\left|\right|{\hat{a}}_k+{\∈}_k{\tilde{e}}_k\left|\right|}{\sqrt{N}(1+{{\∈}_k}^2)}+l$

$\mathrm{\σ}=\sqrt{{\mathrm{Nc}}_{\mathrm{\σ}}/\mathrm{tr}\left({\hat{R}}^{-1}\right)}/\mathrm{\ρ},p_k=\frac{\sqrt{N}(1+{\widehat{\∈}}_k{\∈}_k)\left|\right|c_k\left|\right|}{(1+{\∈}_k^2)\left|\right|{\hat{a}}_k+{\∈}_k{\tilde{e}}_k\left|\right|\mathrm{\ρ}},k=\mathrm{1,2},...K,$

根据上式更新p^k,σ的值。

4d)根据上述步骤中求得的∈,p,σ可进一步更新C：

$c_k=p_k{\hat{R}}^{1/2}{R}^{-1}(a_k+{\∈}_k{e}_k),k=1,...K,[c_{2K+1},...,c_{2K+N}]=\mathrm{\σ}{\hat{R}}^{1/2}{R}^{-1}$

其中C＝[c₁,c₂,...c_2K+N]^H。

5)由求得的C值重复上述4)中的三个步骤直到最大迭代次数

或者满足迭代收敛条件最终求得最优解。

仿真结果：下面结合仿真分析本发明的性能。

仿真条件与方法：

天线阵列为均匀线性阵列且天线数N为20，各天线接收机的间距d均等于半波长。源信号个数M＝13且相互独立系数稳定，其方向均匀分布在-60～60度间。源信号的模型分别为其中分布在[-π,π]。噪声为均值为零方差为σ＝10^(-SNR/20)的高斯白噪声，SNR定义为s₁(t)功率(设a₁＝1)和σ的比值。采样数T等于200。

2.仿真内容与结果

仿真将本专利所提算法与不考虑稀疏表示误差时方向估计的均方根误差，以及方向估计的克拉美罗界相比，Monte Carlo仿真次数为1000。仿真结果如图2所示。图中，60、180、720分别表示用3、1、0.25度为网格密度覆盖波达方向的取值空间。由图可见，不考虑格点误差时，方向估计精度受到密度的限制，而本专利算法性能逼近克拉美罗界，并且所需的网格密度大大降低。显然，当网络密度降低时覆盖波达方向取值空间的格点数大大降低，因此也大大降低了算法的求解复杂度。

Claims (3)

1.一种基于稀疏表示和协方差拟合的稳健波达角估计方法，其特征在于包括以下步骤：

步骤1：通过稀疏空间谱表示对天线阵列接收信号建立模型，并将模型误差进行参数化表示，根据协方差拟合准则得到参数估计的优化问题；

步骤2：将步骤1中得到的非凸优化问题进行转化；首先不考虑模型误差，得到能够快速求解的凸优化问题，从而求得初始解；然后再利用迭代方法更新初始估计和模型误差。

2.根据权利要求1所述的一种基于稀疏表示和协方差拟合的稳健波达角估计方法，其特征在于步骤1中所述的针对天线阵列接收信号建立基于协方差拟合准则的稀疏优化问题，具体为：采用N个全向天线接收机形成均匀线性阵列，对空间中K个窄带信号源进行方向估计；基于阵列输出信号的稀疏空间谱表示，同时将稀疏空间谱表示的误差进行参数化表示，根据天线阵列的输出信号计算信号协方差矩阵R：

根据协方差拟合准则并且引入未知数C建立最优化问题：

p_K≥0,|∈_k|≤l_g,k＝1,2,...K

σ≥0 。

3.根据权利要求1所述的一种基于稀疏表示和协方差拟合的稳健波达角估计方法，其特征在于步骤2中所述的求解步骤1中建立的优化问题，具体为：步骤1中建立的优化问题为非凸优化问题，因此将问题的求解分为两步：首先假设模型误差为零，在不考虑模型误差的情况下利用闭式解迭代的方法求解p,σ,C的初始值；利用求解的p,σ,C代入原优化问题，得到关于∈的SOCP问题：

其中根据求解的∈更新p,σ：

根据上述步骤中求解的∈,p,σ可进一步更新C：

其中C＝[c₁,c₂,...c_2K+N]^H；重复上述步骤直到满足收敛条件求得最优解。