CN103954950B - 一种基于样本协方差矩阵稀疏性的波达方向估计方法 - Google Patents

Abstract

本发明公开了一种基于样本协方差矩阵稀疏性的波达方向估计方法，涉及阵列信号处理领域。其步骤为：步骤1，对雷达天线阵列接收信号的协方差矩阵进行降维，得到降维后的雷达天线阵列接收信号的协方差矩阵；步骤2，根据降维后的协方差矩阵建立基于稀疏约束的稀疏矢量的代价函数；步骤3，将基于稀疏约束的稀疏矢量的代价函数构造成适合于凸规划包求解的形式，根据凸规划包求解的形式求解稀疏矢量；步骤4，将稀疏矢量中的非零元素确定为目标的方位角，该目标的方位角即目标的波达方向。本发明主要解决现有技术中需要已知噪声功率，而且计算量很大的问题，本发明主要用于阵列信号处理的场景。

Description

一种基于样本协方差矩阵稀疏性的波达方向估计方法

技术领域

本发明属于雷达技术领域，涉及阵列信号处理领域，尤其涉及一种基于样本协方差矩阵稀疏性的波达方向估计方法。

背景技术

近几十年来，远场窄带信号的波达方向(Direction of Arrive，DOA)估计一直是阵列信号处理的热点问题，被广泛地应用于雷达、电磁场、无线通信、医学成像和地震勘探等诸多领域。DOA估计的主要目标是在噪声环境下，检测和估计多个信号的方位。针对DOA估计问题，人们提出了大量的DOA估计方法，包括：波束形成(Beam Forming,BF)法，基于子空间方法的多重信号分类(Multiple Signal Classifcation,MUSIC)法和基于最小方差无畸变(Minimum Variance Distortion Response,MVDR)的高分辨谱估计法(如CAPON法)，确定最大似然估计法和统计最大似然估计法等。然而这些方法大多只能处理非相关信号，对相关信号都需要进行去相关处理，因而限制了算法的应用。

针对传统DOA估计方法存在的不足，人们提出了一类新的DOA估计方法，该方法基于信号空间分布的稀疏性，用稀疏表示方法来解决DOA估计问题。例如Gorodnitsky等提出了一种加权迭代最小范数(FOCal-Underdetermined System Solver,FOCUSS)方法来求解稀疏表示问题，并在DOA估计上取得了较好的结果，但是该方法仅适用于单次快拍。传统的稀疏表示问题都是针对单帧接收到的回波信号进行的，而实际用于DOA估计的信号却有多帧。相较于单帧信号的DOA估计，多帧信号联合处理能够有效地提高算法的抗噪能力和可靠性，即使在低信噪比(Signal to Noise Ratio,SNR)的情况下也可以有较高的DOA估计精度。基于奇异值分解的多帧联合欠定系统正则化聚焦求解算法(Singular ValueDecomposition—Regularized FOCal-Underdetermined System Solver,SVD—RFOCUSS)能够实现DOA的超分辨估计，该算法可以准确地估计目标信号的方位，并且计算量小，但是它需要确定正则化参数，当正则化参数选择不准时，算法就会失效。l1-SRACV(l1-SparseRepresentation of Array Covariance Vectors)算法是一种基于l₁范数的样本协方差矩阵稀疏表示DOA估计方法，该算法在噪声环境下不需要直接确定正则化参数，但是需要已知噪声功率，而且计算量很大。

发明内容

针对现有DOA估计方法的不足，如只能处理非相关信号，对相关信号需要进行去相关处理，计算量大，需要确定正则化参数且正则化参数选取困难，需要已知噪声功率，本发明提出了一种新的基于样本协方差矩阵稀疏性的DOA估计方法，该方法对传统的基于协方差矩阵稀疏性的DOA估计模型进行改进，不需要已知或者估计噪声功率，稳健性和计算效率都得到了提升。

为达到上述目的，本发明采用以下技术方案预以实现。

一种基于样本协方差矩阵稀疏性的波达方向估计方法，其特征在于，包括以下步骤：

步骤1，雷达天线阵列接收目标的回波信号；对该回波信号的协方差矩阵R进行降维，得到降维后的回波信号的协方差矩阵

步骤2，根据降维后的协方差矩阵建立基于稀疏约束的稀疏矢量的代价函数；

步骤3，将基于稀疏约束的稀疏矢量的代价函数构造成适合于凸规划包求解的形式，根据凸规划包求解的形式求解稀疏矢量b₀；

步骤4，将稀疏矢量b₀中的非零元素确定为目标的方位角，该目标的方位角为目标的波达方向。

上述技术方案的特点和进一步改进在于：

(1)步骤1具体包括以下子步骤：

1a)设定雷达天线阵列为均匀线阵，阵元数为M，阵元间距为d，其中，d＝λ/2，λ为雷达工作波长，有P个随机分布的远场窄带平稳信号s_k(t)，分别以方向θ_k入射到M个阵元上，k＝1,2…,P，θ＝[θ₁,θ₂,…,θ_P]；设定雷达接收机噪声为加性高斯白噪声，则回波信号矢量表示为如下形式：

y(t)＝A(θ)S(t)+n(t),t＝1,…,L (1)

其中，t表示时刻，L为非负常数，y(t)表示t时刻M×1维的回波信号矢量，M表示雷达天线阵元数，n(t)为t时刻的雷达天线阵列接收噪声；S(t)＝[s₁(t),s₂(t),…,s_P(t)]^T为P×1维的平稳信号的矢量，且平稳信号s_k(t)为零均值的随机变量；A(θ)＝[a(θ₁),a(θ₂),…,a(θ_P)]表示M×P维的阵列流形矩阵并且其P个列矢量分别为： $a\left({\mathrm{\θ}}_k\right)={[1,e^{j2\mathrm{\π}\frac{d}{\mathrm{\λ}}\sin \left({\mathrm{\θ}}_k\right)},\·\·\·,e^{j2\mathrm{\π}\frac{d}{\mathrm{\λ}}\sin \left({\mathrm{\θ}}_k\right)(M-1)}]}^T,k=1,\·\·\·,P,$ θ_k表示目标的方位角，P表示随机分布的远场窄带平稳信号的个数，T表示转置操作；

1b)设定回波信号与噪声之间不相关，则回波信号的协方差矩阵R表示为如下形式：

R＝E[y(t)y^H(t)]＝A(θ)R_sA^H(θ)+σ²I (2)

其中，R表示回波信号的协方差矩阵，y(t)表示t时刻M×1维的回波信号矢量，A(θ)表示阵列流形矩阵，R_s＝E[S(t)S^H(t)]表示目标信号的协方差矩阵，σ²为噪声功率，I表示M×M维单位阵，M表示雷达天线阵元数，H表示共轭转置操作，E[·]表示求期望运算；

1c)给定超完备空间A(Θ)，回波信号的协方差矩阵R的第i列r_i用给定的超完备空间A(Θ)表示为：

$r_i=E\left[y\left(t\right)y_i^{*}\left(t\right)\right]=A\left(\mathrm{\Θ}\right)b_i+{\mathrm{\σ}}^2{e}_i,i=\mathrm{1,2}...,M---\left(3\right)$

其中，r_i表示回波信号的协方差矩阵R的第i列，E[·]表示求期望运算，y(t)表示t时刻M×1维的回波信号矢量，y_i(t)表示第i个雷达天线阵元t时刻接收到的信号矢量，*表示取共轭操作，A(Θ)＝{a(θ)|θ∈Θ}表示给定的超完备空间，表示目标信号可能的来波方向空域方位角集合，N为远大于M的正整数，b_i表示N×1维的系数矢量，σ²表示噪声功率，e_i表示误差项，e_i是第i个元素为1，其他元素均为零的M×1维的矢量列矢量，M表示雷达天线阵元数；

1d)用变换选择矩阵J_i左乘回波信号的协方差矩阵R的第i列r_i，得到降维后的回波信号的协方差矩阵的第i列的表达式为：

t_i＝J_ir_i＝J_iA(Θ)b_i+σ²J_ie_i＝J_iA(Θ)b_i,i＝1,2...,M (4)

其中，t_i表示降维后的回波信号的协方差矩阵的第i列，i＝1,2...,M，J_i表示变换选择矩阵且J_i为(M-1)×M维，J_i＝[e₁,…,e_i-1,e_i+1,…,e_M]^T，r_i表示回波信号的协方差矩阵R的第i列，A(Θ)表示给定的超完备空间，b_i表示系数矢量，σ²表示噪声功率，e_i表示误差项，e_i是第i个元素为1，其他元素均为零的M×1维的列矢量，M表示雷达天线阵元数，T表示转置操作；

1e)降维后的回波信号的协方差矩阵的表达式为：

其中，t_i表示降维后的回波信号的协方差矩阵的第i列，i＝1,2...,M，M表示雷达天线阵元数。

(2)步骤2具体包括以下子步骤：

2a)设定l₁范数为： ${\left|\right|b_0\left|\right|}_1={\mathrm{\Σ}}_{n=1}^N\left|b_0\left(n\right)\right|;$

其中，||·||₁表示取l₁范数操作，b₀(n)表示稀疏矢量b₀的第n个元素，n＝1,2...,N，N为远大于M的正整数，M表示雷达天线阵元数，|·|表示取绝对值操作；

2b)利用l₁范数作为稀疏矢量的代价函数，将降维后的回波信号协方差矩阵的第i列作为稀疏矢量的代价函数的稀疏约束条件，得到如下初始的基于稀疏约束的稀疏矢量的代价函数：

$\begin{array}{c}\min {\left|\right|b_0\left|\right|}_1\\ s.t.t_i=J_{i}A\left(\mathrm{\Θ}\right)b_i,i=\mathrm{1,2}...M\end{array}---\left(5\right)$

其中，min表示取最小值操作，s.t.表示取约束操作，b₀表示稀疏矢量，||·||₁表示取l₁范数操作，t_i表示降维后的回波信号协方差矩阵的第i列，i＝1,2...,M，J_i表示(M-1)×M维的变换选择矩阵，A(Θ)表示给定的超完备空间，b_i表示系数矢量，M表示雷达天线阵元数；

2c)用回波信号的采样数据样本对降维后的回波信号协方差矩阵的第i列t_i和回波信号的协方差矩阵R进行估计，获取基于稀疏约束的稀疏矢量的代价函数；

步骤2的子步骤2c)具体为：

2c1)利用回波信号的采样数据样本对降维后的回波信号协方差矩阵的第i列t_i进行估计，表示为下式：

${\hat{t}}_i=J_i(1/L{\mathrm{\Σ}}_{t=1}^{L}y\left(t\right){y_i}^H\left(t\right))=t_i+{\mathrm{\Δt}}_i$

其中，表示降维后的回波信号协方差矩阵的第i列t_i的样本估计值，J_i表示变换选择矩阵，且J_i为(M-1)×M维，y(t)表示t时刻M×1维的回波信号矢量，y_i(t)表示第i个雷达天线阵元t时刻接收到的信号矢量，t_i表示降维后的回波信号协方差矩阵的第i列，i＝1,2...,M，Δt_i为估计误差，L为样本的个数，M表示雷达天线阵元数，H表示共轭转置操作；

2c2)设定估计误差矢量 $\mathrm{\Δt}={[{\mathrm{\Δt}}_1^T,{\mathrm{\Δt}}_2^T,\·\·\·{\mathrm{\Δt}}_M^T]}^T=J\·\mathrm{vec}(R-\hat{R});$

其中，Δt_i为估计误差，i＝1,2...,M，J是由变换选择矩阵J_i构造的块对角矩阵，表示回波信号的协方差矩阵R的样本估计值，R表示回波信号的协方差矩阵，vec(·)表示列矢量化操作，它是将矩阵的各列按列排成一个列矢量，M表示雷达天线阵元数，T表示转置操作；

估计误差矢量Δt满足如下分布：

$\mathrm{\Δt}~\mathrm{AsN}(0,J(\frac{1}{L}{R}^T\⊗R)J^T)---\left(6\right)$

其中，Δt表示估计误差矢量，～表示服从某种分布，AsN(μ,C)表示均值矢量为μ，协方差矩阵为C的近似多维正态分布；L为样本的个数，R表示回波信号的协方差矩阵，T表示转置操作，表示克罗内克积，J表示由变换选择矩阵J_i构造的块对角矩阵，其形式如下：

2c3)根据估计误差矢量Δt的协方差矩阵的逆W^-1，以及估计误差矢量Δt满足的分布，得到W^-1/2Δt满足如下分布：

W^-1/2Δt～AsN(0,I_M(M-1))

其中，估计误差矢量Δt的协方差矩阵的逆W^-1为：J表示由变换选择矩阵J_i构造的块对角矩阵，L为样本的个数，R表示回波信号的协方差矩阵，T表示转置操作，表示克罗内克积，(·)^-1表示取逆操作，W^-1/2表示W^-1的平方根矩阵，～表示服从某种分布，AsN(0,I_M(M-1))表示均值矢量为0，协方差矩阵为I_M(M-1)的近似多维正态分布，I_M(M-1)表示M(M-1)维的单位阵，M表示雷达天线阵元数；

2c4)对W^-1/2Δt取2范数得到：

${\left|\right|W^{-1/2}\mathrm{\Δt}\left|\right|}_2^2~{\mathrm{As\χ}}^2\left(M(M-1)\right)$

其中，W^-1/2表示估计误差矢量Δt的协方差矩阵的逆W^-1的平方根矩阵，Δt表示估计误差矢量，||·||₂表示取2范数操作，～表示服从某种分布，Asχ²(M(M-1))表示自由度为M(M-1)的近似χ²分布，M表示雷达天线阵元数；

2c5)引入一个参数η，用来约束残差上限，使得以概率p成立；

其中，p大于等于0并且小于1；W^-1/2表示估计误差矢量Δt的协方差矩阵的逆W^-1的平方根矩阵，Δt表示估计误差矢量，||·||₂表示取2范数操作；

2c6)将概率p设定为p＝0.999，使得成立，则参数η根据χ²分布的概率累积函数以及p＝0.999求出，使得不等式以概率p＝0.999成立；

2c7)在不等式以概率p＝0.999成立的情况下，得到如下基于稀疏约束的稀疏矢量的代价函数，表示如下：

$\begin{array}{c}\min {\left|\right|b_0\left|\right|}_1\\ {\left|\right|W^{-1/2}(J\·\mathrm{vec}(\hat{R}-A\left(\mathrm{\Θ}\right)B))\left|\right|}_2^2\≤\mathrm{\η}\end{array}---\left(7\right)$

其中，min表示取最小值操作，b₀表示稀疏矢量，||·||₁表示取l₁范数操作，||·||₂表示取2范数操作，J表示由变换选择矩阵J_i构造的块对角矩阵，vec(·)表示列矢量化操作，A(Θ)表示给定的超完备空间，B是由系数矢量b_i构成的矩阵，B＝[b₁,b₂,…,b_M]，i＝1,2...,M，M表示雷达天线阵元数，η表示约束残差上限的参数，表示回波信号的协方差矩阵R的样本估计值，即x(t)表示t时刻的回波信号的采样数据样本，H表示共轭转置操作；用来逼近估计误差矢量Δt的协方差矩阵的逆W^-1的平方根矩阵W^-1/2，即 ${\hat{W}}^{-1/2}={\left(J(1/L{\hat{R}}^T\⊗\hat{R})J^T\right)}^{-1/2},$ L为样本的个数，表示克罗内克积。

(3)步骤3具体包括以下子步骤：

3a)将基于稀疏约束的稀疏矢量的代价函数式(7)构造成如下的适合于凸规划包求解的形式：

$\begin{array}{c}\min g\\ s.t.1^T\mathrm{\γ}\≤g,\\ {\left({\mathrm{\Σ}}_{m=1}^{M}B{(n,m)}^2\right)}^{1/2}=b_0\left(n\right)\≤{\mathrm{\γ}}_m,n=\mathrm{1,2}\·\·\·N\\ {\left|\right|{\hat{W}}^{-1/2}(J\·\mathrm{vec}(\hat{R}-A\left(\mathrm{\Θ}\right)B))\left|\right|}_2^2\≤\mathrm{\η}\end{array}---\left(8\right)$

其中，min表示取最小值操作，s.t.表示取约束操作，g表示松弛变量，1表示N×1维的全1矢量，中间变量γ＝[γ₁,γ₂,…γ_N]^T，B(n,m)表示矩阵B的第n行第m个元素，B是由系数矢量b_i构成的矩阵，B＝[b₁,b₂,…,b_M]，m＝1,2...,M，n＝1,2...,N，N为远大于M的正整数，b₀(n)表示稀疏矢量b₀的第n个元素，||·||₂表示取2范数操作，表示W^-1/2的逼近值，W^-1/2表示估计误差矢量Δt的协方差矩阵的逆W^-1的平方根矩阵，J表示由变换选择矩阵J_i构造的块对角矩阵，vec(·)表示列矢量化操作，表示回波信号的协方差矩阵R的样本估计值，A(Θ)表示给定的超完备空间，M表示雷达天线阵元数，η表示约束残差上限的参数，T表示转置操作；

3b)用凸规划包求解式(8)，得到稀疏矢量b₀。

与现有技术相比，本发明具有突出的实质性特点和显著的进步。本发明与现有方法相比，具有以下优点：

1、传统的角度高分辨DOA估计方法，例如MUSIC法、CAPON法等，只能处理非相关信号，对相关信号需要进行去相关处理，因而大大限制了这些算法的应用。而本发明方法是基于稀疏表示的DOA估计方法，对信号的相关性不敏感，如图3(a)和图3(b)，因而对具有任意相关性的信号源均能进行有效的DOA估计，不需要进行去相关处理，并且具有很高的分辨力和估计精度，应用范围比较广泛。

2、传统的稀疏表示DOA估计方法(如FOCUSS方法)大都是针对单帧接收到的回波信号进行的，但是实际用于DOA估计的信号一般都有多帧，因而导致算法的应用范围有限。而本发明方法适用于接收到的回波信号为多帧的情况，对多帧信号进行联合处理，有效地提高了本方法的抗噪能力和可靠性，而且在低信噪比情况下也有较高的DOA估计精度，如图4所示。

3、传统的基于协方差矩阵稀疏性的DOA估计方法，例如l1-SRACV算法，需要已知或者估计噪声功率，而且计算量很大。当噪声功率估计不准时，算法对目标的检测性能会下降，算法的性能很不稳定，如图4所示。而本发明方法不需要已知噪声功率，从而大大增加了该方法的可适应性和稳健性。此外，本发明方法通过采用降维措施，大大降低了估计目标源方位所需的计算量，以极低的性能损失获取了较高的计算效率。

本发明更进一步涉及阵列信号处理技术领域中的一种基于样本协方差矩阵稀疏性的波达方向估计方法。本发明可用于解决噪声功率未知、雷达接收到的信号具有相关性情况下的目标波达方向(DOA)估计问题。

附图说明

下面结合附图和具体实施方式对本发明做进一步说明。

图1是本发明方法流程图；

图2是仿真实验一中本发明方法和l1-SRACV算法的稀疏空域谱图；

图3是仿真实验一中得到的稀疏空域谱图：图3(a)是不相关信号的稀疏空域谱图；图3(b)是相关信号的稀疏空域谱图；

图4是仿真实验一中本发明方法、l1-SRACV算法和l1-SVD算法的检测概率随信噪比变化的曲线图；

图5是仿真实验二中本发明方法、l1-SRACV算法和l1-SVD算法在不同角度间隔下的检测概率图。

具体实施方式

参照图1，说明本发明一种基于样本协方差矩阵稀疏性的波达方向估计方法，其具体实施步骤如下：

步骤1，雷达天线阵列接收目标的回波信号；对该回波信号的协方差矩阵R进行降维，得到降维后的回波信号的协方差矩阵

1a)设定雷达天线阵列为均匀线阵，阵元数为M，阵元间距为d，其中，d＝λ/2，λ为雷达工作波长，有P个随机分布的远场窄带平稳信号s_k(t)，分别以方向θ_k入射到M个阵元上，k＝1,2…,P，θ＝[θ₁,θ₂,…,θ_P]；设定雷达接收机噪声为加性高斯白噪声，则回波信号矢量表示为如下形式：

y(t)＝A(θ)S(t)+n(t),t＝1,…,L (1)

其中，t表示时刻，L为非负常数，y(t)表示t时刻M×1维的回波信号矢量，M表示雷达天线阵元数，n(t)为t时刻的雷达天线阵列接收噪声；S(t)＝[s₁(t),s₂(t),…,s_P(t)]^T为P×1维的平稳信号的矢量，且平稳信号s_k(t)为零均值的随机变量；A(θ)＝[a(θ₁),a(θ₂),…,a(θ_P)]表示M×P维的阵列流形矩阵并且其P个列矢量分别为：

$a\left({\mathrm{\θ}}_k\right)={[1,e^{j2\mathrm{\π}\frac{d}{\mathrm{\λ}}\sin \left({\mathrm{\θ}}_k\right)},\·\·\·,e^{j2\mathrm{\π}\frac{d}{\mathrm{\λ}}\sin \left({\mathrm{\θ}}_k\right)(M-1)}]}^T,k=1,\·\·\·,P,$ θ_k表示目标的方位角，P表示随机分布的远场窄带平稳信号的个数，T表示转置操作。

1b)设定回波信号与噪声之间不相关，则回波信号的协方差矩阵R表示为如下形式：

R＝E[y(t)y^H(t)]＝A(θ)R_sA^H(θ)+σ²I (2)

其中，R表示回波信号的协方差矩阵，y(t)表示t时刻M×1维的回波信号矢量，A(θ)表示阵列流形矩阵，R_s＝E[S(t)S^H(t)]表示目标信号的协方差矩阵，σ²为噪声功率，I表示M×M维单位阵，M表示雷达天线阵元数，H表示共轭转置操作，E[·]表示求期望运算；

1c)给定超完备空间A(Θ)，回波信号的协方差矩阵R的第i列r_i用给定的超完备空间A(Θ)表示为：

$r_i=E\left[y\left(t\right)y_i^{*}\left(t\right)\right]=A\left(\mathrm{\Θ}\right)b_i+{\mathrm{\σ}}^2{e}_i,i=\mathrm{1,2}...,M---\left(3\right)$

其中，r_i表示回波信号的协方差矩阵R的第i列，E[·]表示求期望运算，y(t)表示t时刻M×1维的回波信号矢量，y_i(t)表示第i个雷达天线阵元t时刻接收到的信号矢量，*表示取共轭操作，A(Θ)＝{a(θ)|θ∈Θ}表示给定的超完备空间，表示目标信号可能的来波方向空域方位角集合，N为远大于M的正整数，b_i表示N×1维的系数矢量，σ²表示噪声功率，e_i表示误差项，e_i是第i个元素为1，其他元素均为零的M×1维的矢量列矢量，M表示雷达天线阵元数。

在本发明中，由于A(Θ)是超完备空间，因而满足方程式(3)的b_i有无数个，可以利用信号相对于超完备空间A(Θ)的稀疏性来约束式(3)，使式(3)有唯一解。考虑理想情况下，即Θ在整个角度空间采样足够稠密，那么A(Θ)的某些列将会很接近其中，P表示随机分布的远场窄带平稳信号的个数。所以，理想情况下的系数矢量b_i具有稀疏结构，即b_i中除对应于信号方向的元素非零外，其余元素均为零。根据b_i的稀疏性，应用稀疏表示理论，可以求解欠定方程组(3)的唯一解，从而估计出目标的方位角。

由式(3)可知：噪声功率σ²仅影响r_i的第i个元素，通过步骤1d)对回波信号的协方差矩阵R的各列进行降维处理，即实现对R的降维，以达到避免估计噪声功率的目的。

1d)用变换选择矩阵J_i左乘回波信号的协方差矩阵R的第i列r_i，得到降维后的回波信号的协方差矩阵的第i列的表达式为：

t_i＝J_ir_i＝J_iA(Θ)b_i+σ²J_ie_i＝J_iA(Θ)b_i,i＝1,2...,M (4)

其中，t_i表示降维后的回波信号的协方差矩阵的第i列，i＝1,2...,M，J_i表示变换选择矩阵且J_i为(M-1)×M维，J_i＝[e₁,…,e_i-1,e_i+1,…,e_M]^T，r_i表示回波信号的协方差矩阵R的第i列，A(Θ)表示给定的超完备空间，b_i表示系数矢量，σ²表示噪声功率，e_i表示误差项，e_i是第i个元素为1，其他元素均为零的M×1维的列矢量，M表示雷达天线阵元数，T表示转置操作。

1e)降维后的回波信号的协方差矩阵的表达式为：

其中，t_i表示降维后的回波信号的协方差矩阵的第i列，i＝1,2...,M，M表示雷达天线阵元数。

雷达天线阵列接收到的回波信号的协方差矩阵R为M×M维矩阵，降维之后为(M-1)×M维矩阵，虽然只降了一维，但仍然可以大大减少计算量。此外，降维措施也避免了对噪声功率的估计，大大增加了本发明方法的可适应性和稳健性。

步骤2，根据降维后的协方差矩阵建立基于稀疏约束的稀疏矢量的代价函数；

显然，降维处理对系数矢量b_i的稀疏性没有影响。令B＝[b₁,b₂,…,b_M]，可以看出，B的每一列都满足相同的稀疏结构，即非零元素在每一列的位置相同。对B的每行取2-范数，得到一个新的N×1维的列矢量b₀，其第n个元素b₀(n)为：其中，B(n,m)表示矩阵B的第n行第m个元素，m＝1,2...,M，n＝1,2...,N，N为远大于M的正整数，M表示雷达天线阵元数。理想的稀疏矩阵B的联合稀疏结构可以由b₀的稀疏性来描述，即只需找到一个足够稀疏的矢量b₀，就会使B很好地满足约束，同时B的每一列都具有与b₀相同的稀疏结构。那么，基于式(4)的DOA估计就等价于找到一个稀疏矢量b₀。

稀疏矢量b₀的稀疏程度的计算，直接影响到对于式(4)的求解。衡量矢量稀疏度的最直接的方法是计算该矢量中非零元素的个数，即计算矢量的l₀范数，但是这种方法在数学上难以求解。

本发明将l₁范数作为评价矢量稀疏度的代价函数，这是因为：一、l₁范数已经被证明在稀疏矢量b₀足够稀疏的条件下高概率等价于l₀范数；二、l₁范数是凸函数，能保证代价函数收敛到全局最优解。

2a)设定l₁范数为： ${\left|\right|b_0\left|\right|}_1={\mathrm{\Σ}}_{n=1}^N\left|b_0\left(n\right)\right|;$

其中，||·||₁表示取l₁范数操作，b₀(n)表示稀疏矢量b₀的第n个元素，n＝1,2...,N，N为远大于M的正整数，M表示雷达天线阵元数，|·|表示取绝对值操作；

2b)利用l₁范数作为稀疏矢量的代价函数，将降维后的回波信号协方差矩阵的第i列作为稀疏矢量的代价函数的稀疏约束条件，得到如下初始的基于稀疏约束的稀疏矢量的代价函数：

$\begin{array}{c}\min {\left|\right|b_0\left|\right|}_1\\ s.t.t_i=J_{i}A\left(\mathrm{\Θ}\right)b_i,i=\mathrm{1,2}...M\end{array}---\left(5\right)$

其中，min表示取最小值操作，s.t.表示取约束操作，b₀表示稀疏矢量，||·||₁表示取l₁范数操作，t_i表示降维后的回波信号协方差矩阵的第i列，i＝1,2...,M，J_i表示(M-1)×M维的变换选择矩阵，A(Θ)表示给定的超完备空间，b_i表示系数矢量，M表示雷达天线阵元数。

2c)用回波信号的采样数据样本对降维后的回波信号协方差矩阵的第i列t_i和回波信号的协方差矩阵R进行估计，获取基于稀疏约束的稀疏矢量的代价函数；

子步骤2c)具体为：

2c1)利用回波信号的采样数据样本对降维后的回波信号协方差矩阵的第i列t_i进行估计，表示为下式：

${\hat{t}}_i=J_i(1/L{\mathrm{\Σ}}_{t=1}^{L}y\left(t\right){y_i}^H\left(t\right))=t_i+{\mathrm{\Δt}}_i$

其中，表示降维后的回波信号协方差矩阵的第i列t_i的样本估计值，J_i表示变换选择矩阵，且J_i为(M-1)×M维，y(t)表示t时刻M×1维的回波信号矢量，y_i(t)表示第i个雷达天线阵元t时刻接收到的信号矢量，t_i表示降维后的回波信号协方差矩阵的第i列，i＝1,2...,M，Δt_i为估计误差，L为样本的个数，M表示雷达天线阵元数，H表示共轭转置操作。

2c2)设定估计误差矢量 $\mathrm{\Δt}={[{\mathrm{\Δt}}_1^T,{\mathrm{\Δt}}_2^T,\·\·\·{\mathrm{\Δt}}_M^T]}^T=J\·\mathrm{vec}(R-\hat{R});$

其中，Δt_i为估计误差，i＝1,2...,M，J是由变换选择矩阵J_i构造的块对角矩阵，表示回波信号的协方差矩阵R的样本估计值，R表示回波信号的协方差矩阵，vec(·)表示列矢量化操作，它是将矩阵的各列按列排成一个列矢量，M表示雷达天线阵元数，T表示转置操作；

估计误差矢量Δt满足如下分布：

$\mathrm{\Δt}~\mathrm{AsN}(0,J(\frac{1}{L}{R}^T\⊗R)J^T)---\left(6\right)$

其中，Δt表示估计误差矢量，～表示服从某种分布，AsN(μ,C)表示均值矢量为μ，协方差矩阵为C的近似多维正态分布；L为样本的个数，R表示回波信号的协方差矩阵，T表示转置操作，表示克罗内克积，J表示由变换选择矩阵J_i构造的块对角矩阵，其形式如下：

在本发明中，当在式(5)中用降维后的回波信号协方差矩阵的第i列t_i的样本估计值代替t_i时，由于估计误差的存在，式(5)的等式约束将不再成立。

将对目标的DOA估计看作是求解最稀疏的稀疏矢量b₀，同时保证最小，把初始的基于稀疏约束的稀疏矢量的代价函数式(5)改写为如下形式：

$\min {\left|\right|b_0\left|\right|}_1+\mathrm{\λ}{\left|\right|\mathrm{\Δt}\left|\right|}_2^2$

其中，min表示取最小值操作，b₀表示稀疏矢量，||·||₁表示取l₁范数操作，Δt表示估计误差矢量，||·||₂表示取2范数操作，λ为正则化参数，用来平衡稀疏矢量b₀的稀疏程度和约束残差，对于正确估计目标信号的波达方向具有重要的作用。但是，要确定一个合适的λ非常困难。因此，为了避免直接确定正则化参数λ，本发明给出一种足够多样本情况下的参数确定标准来控制约束残差上限。执行以下步骤2c3)：

2c3)根据估计误差矢量Δt的协方差矩阵的逆W^-1，以及估计误差矢量Δt满足的分布，得到W^-1/2Δt满足如下分布：

W^-1/2Δt～AsN(0,I_M(M-1))

其中，估计误差矢量Δt的协方差矩阵的逆W^-1为：J表示由变换选择矩阵J_i构造的块对角矩阵，L为样本的个数，R表示回波信号的协方差矩阵，T表示转置操作，表示克罗内克积，(·)^-1表示取逆操作，W^-1/2表示W^-1的平方根矩阵，～表示服从某种分布，AsN(0,I_M(M-1))表示均值矢量为0，协方差矩阵为I_M(M-1)的近似多维正态分布，I_M(M-1)表示M(M-1)维的单位阵，M表示雷达天线阵元数；

2c4)对W^-1/2Δt取2范数得到：

${\left|\right|W^{-1/2}\mathrm{\Δt}\left|\right|}_2^2~{\mathrm{As\χ}}^2\left(M(M-1)\right)$

其中，W^-1/2表示估计误差矢量Δt的协方差矩阵的逆W^-1的平方根矩阵，Δt表示估计误差矢量，||·||₂表示取2范数操作，～表示服从某种分布，Asχ²(M(M-1))表示自由度为M(M-1)的近似χ²分布，M表示雷达天线阵元数。

2c5)引入一个参数η，用来约束残差上限，使得以概率p成立。

其中，p大于等于0并且小于1；W^-1/2表示估计误差矢量Δt的协方差矩阵的逆W^-1的平方根矩阵，Δt表示估计误差矢量，||·||₂表示取2范数操作。

2c6)将概率p设定为p＝0.999，使得成立，则参数η根据χ²分布的概率累积函数以及p＝0.999求出，使得不等式以概率p＝0.999成立。

2c7)在不等式以概率p＝0.999成立的情况下，得到如下基于稀疏约束的稀疏矢量的代价函数，表示如下：

$\begin{array}{c}\min {\left|\right|b_0\left|\right|}_1\\ {\left|\right|W^{-1/2}(J\·\mathrm{vec}(\hat{R}-A\left(\mathrm{\Θ}\right)B))\left|\right|}_2^2\≤\mathrm{\η}\end{array}---\left(7\right)$

其中，min表示取最小值操作，b₀表示稀疏矢量，||·||₁表示取l₁范数操作，||·||₂表示取2范数操作，J表示由变换选择矩阵J_i构造的块对角矩阵，vec(·)表示列矢量化操作，A(Θ)表示给定的超完备空间，B是由系数矢量b_i构成的矩阵，B＝[b₁,b₂,…,b_M]，i＝1,2...,M，M表示雷达天线阵元数，η表示约束残差上限的参数，表示回波信号的协方差矩阵R的样本估计值，即x(t)表示t时刻的回波信号的采样数据样本，H表示共轭转置操作；用来逼近估计误差矢量Δt的协方差矩阵的逆W^-1的平方根矩阵W^-1/2，即 ${\hat{W}}^{-1/2}={\left(J(1/L{\hat{R}}^T\⊗\hat{R})J^T\right)}^{-1/2},$ L为样本的个数，表示克罗内克积。

在本发明中，当样本数足够大时，估计误差能得到很好的抑制。

步骤3，将基于稀疏约束的稀疏矢量的代价函数构造成适合于凸规划包求解的形式，根据凸规划包求解的形式求解稀疏矢量b₀；

基于稀疏约束的稀疏矢量的代价函数式(7)中，由于||b₀||₁为凸函数，而二次约束也为凸约束，因此式(7)可以由凸规划包进行求解。

3a)将基于稀疏约束的稀疏矢量的代价函数式(7)构造成如下的适合于凸规划包求解的形式：

$\begin{array}{c}\min g\\ s.t.1^T\mathrm{\γ}\≤g,\\ {\left({\mathrm{\Σ}}_{m=1}^{M}B{(n,m)}^2\right)}^{1/2}=b_0\left(n\right)\≤{\mathrm{\γ}}_m,n=\mathrm{1,2}\·\·\·N\\ {\left|\right|{\hat{W}}^{-1/2}(J\·\mathrm{vec}(\hat{R}-A\left(\mathrm{\Θ}\right)B))\left|\right|}_2^2\≤\mathrm{\η}\end{array}---\left(8\right)$

其中，min表示取最小值操作，s.t.表示取约束操作，g表示松弛变量，1表示N×1维的全1矢量，中间变量γ＝[γ₁,γ₂,…γ_N]^T，B(n,m)表示矩阵B的第n行第m个元素，B是由系数矢量b_i构成的矩阵，B＝[b₁,b₂,…,b_M]，m＝1,2...,M，n＝1,2...,N，N为远大于M的正整数，b₀(n)表示稀疏矢量b₀的第n个元素，||·||₂表示取2范数操作，表示W^-1/2的逼近值，W^-1/2表示估计误差矢量Δt的协方差矩阵的逆W^-1的平方根矩阵，J表示由变换选择矩阵J_i构造的块对角矩阵，vec(·)表示列矢量化操作，表示回波信号的协方差矩阵R的样本估计值，A(Θ)表示给定的超完备空间，M表示雷达天线阵元数，η表示约束残差上限的参数，T表示转置操作。

3b)用凸规划包求解式(8)，得到稀疏矢量b₀。

在本发明，凸规划包为本领域技术人员所公知的软件处理程序，例如SeDuMi和CVX等。

步骤4，将稀疏矢量b₀中的非零元素确定为目标的方位角，该目标的方位角为目标的波达方向。

下面结合仿真实验对本发明的效果做进一步说明。

(1)实验条件：

为了进一步说明本发明的基于样本协方差矩阵稀疏性的DOA估计方法较传统DOA估计方法(如l1-SVD算法和l1-SRACV算法)的优越性，做如下两个仿真实验。

系统模型：采用阵元数M＝9的均匀线阵，阵元间距为半波长，平稳快拍观测数L＝300。将方位角按-90°～90°等间隔地分为181份，角度间隔为1°来构造超完备阵列流形矩阵A(Θ)。本发明采用常用的针对基于稀疏表示的DOA估计方法的性能评测指标，即正确检测出信号的概率来评价不同方法的性能，仿真实验中的每一幅图中的每一个点均由200次独立实验得到。当估计得到的信号源方位角度与信号源实际角度相差的绝对值和小于2°时，即定义为一次正确检测。

(2)实验结果分析

实验一：假设有2个远场信号源，也就是远场窄带平稳信号，其信号到达角为：[2°，11°]。设通道噪声为零均值的复高斯白噪声，信噪比设为SNR＝0dB。

图2给出了本发明方法和l1-SRACV算法的稀疏空域谱。如图2所示，横坐标表示角度，纵坐标表示归一化幅度。本发明方法所获得的空域谱与l1-SRACV算法的空域谱性能相似，但是，本发明方法无需已知或者估计噪声功率。

如图3(a)和图3(b)所示，横坐标均表示角度，纵坐标均表示归一化幅度。图3(a)和图3(b)分别给出了不相关信号和相关信号的稀疏空域谱，从图3可以看出，基于稀疏表示的DOA估计方法在两种情况下得到的空域谱接近一致，从而从实验上证明了该类DOA估计方法对信号的相关性不敏感。

图4的横坐标表示信噪比，纵坐标表示检测概率。如图4所示。在低信噪比情况下，本发明方法在检测概率上要高于l1-SRACV算法，但是在信噪比较高时，本发明方法与l1-SRACV算法性能相似，而要优于l1-SVD(l1-Singular Value Decomposition)算法。噪声功率精确已知的情况下，理论上，l1-SRACV算法在检测概率上应优于本发明方法，但是在实验中结果却相反，这是由于在l1-SRACV算法中噪声功率采用最小特征值，噪声功率估计不准所致。当噪声功率估计不准确时，会造成l1-SRACV算法性能的不稳定，而本发明方法不需要估计噪声功率，也能达到与理想性能相当的性能。

实验二：考察三种DOA估计方法(本发明方法、l1-SVD算法和l1-SRACV算法)在不同角度间隔下的检测概率，来评估本发明方法的角度分辨能力。如图5所示，横坐标表示角度差，纵坐标表示检测概率。图5比较了本发明方法、l1-SVD算法和l1-SRACV算法对于两个信号在不同角度间隔下的DOA检测概率，其中，第一个信号的角度为θ₁＝2°，第二个信号的角度为θ₂＝θ₁+Δθ，Δθ为第一个信号与第二个信号的角度差，从2°～18°以1°为步长，SNR＝0dB。图5表明，在角度分辨性能方面，本发明方法和l1-SRACV算法要优于传统的l1-SVD方法，同时，本发明方法略优于l1-SRACV算法。

Claims (4)

1.一种基于样本协方差矩阵稀疏性的波达方向估计方法，其特征在于，包括以下步骤：

步骤1，雷达天线阵列接收目标的回波信号；对该回波信号的协方差矩阵R进行降维，得到降维后的回波信号的协方差矩阵

步骤1包括以下子步骤：

1a)设定雷达天线阵列为均匀线阵，阵元数为M，阵元间距为d，其中，d＝λ/2，λ为雷达工作波长，有P个随机分布的远场窄带平稳信号s_k(t)，分别以方向θ_k入射到M个阵元上，k＝1，2…，P，θ＝[θ₁，θ₂，…，θ_P]；设定雷达接收机噪声为加性高斯白噪声，则回波信号矢量表示为如下形式：

y(t)＝A(θ)S(t)+n(t)，t＝1，…，L (1)

其中，t表示时刻，L为非负常数，y(t)表示t时刻M×1维的回波信号矢量，M表示雷达天线阵元数，n(t)为t时刻的雷达天线阵列接收噪声；S(t)＝[s₁(t)，s₂(t)，…，s_P(t)]^T为P×1维的平稳信号的矢量，且平稳信号s_k(t)为零均值的随机变量；A(θ)＝[a(θ₁)，a(θ₂)，…，a(θ_P)]表示M×P维的阵列流形矩阵并且其P个列矢量分别为：θ_k表示目标的方位角，P表示随机分布的远场窄带平稳信号的个数，T表示转置操作；

1b)设定回波信号与噪声之间不相关，则回波信号的协方差矩阵R表示为如下形式：

R＝E[y(t)y^H(t)]＝A(θ)R_sA^H(θ)+σ²I (2)

其中，R表示回波信号的协方差矩阵，y(t)表示t时刻M×1维的回波信号矢量，A(θ)表示阵列流形矩阵，R_s＝E[S(t)S^H(t)]表示目标信号的协方差矩阵，σ²为噪声功率，I表示M×M维单位阵，M表示雷达天线阵元数，H表示共轭转置操作，E[·]表示求期望运算；

1c)给定超完备空间A(Θ)，回波信号的协方差矩阵R的第i列r_i用给定的超完备空间A(Θ)表示为：

$r_i=E\[y\left(t\right)y_i^{*}\left(t\right)\]=A\left(\mathrm{\Θ}\right)b_i+{\mathrm{\σ}}^2{e}_i,i=1,2...,M---\left(3\right)$

其中，r_i表示回波信号的协方差矩阵R的第i列，E[·]表示求期望运算，y(t)表示t时刻M×1维的回波信号矢量，y_i(t)表示第i个雷达天线阵元t时刻接收到的信号矢量，*表示取共轭操作，A(Θ)＝{a(θ)|θ∈Θ}表示给定的超完备空间，表示目标信号可能的来波方向空域方位角集合，N为远大于M的正整数，b_i表示N×1维的系数矢量，σ²表示噪声功率，e_i表示误差项，e_i是第i个元素为1，其他元素均为零的M×1维的矢量列矢量，M表示雷达天线阵元数；

1d)用变换选择矩阵J_i左乘回波信号的协方差矩阵R的第i列r_i，得到降维后的回波信号的协方差矩阵的第i列的表达式为：

t_i＝J_ir_i＝J_iA(Θ)b_i+σ²J_ie_i＝J_iA(Θ)b_i，i＝1，2...，M (4)

其中，t_i表示降维后的回波信号的协方差矩阵的第i列，i＝1，2...，M，J_i表示变换选择矩阵且J_i为(M-1)×M维，J_i＝[e₁，…，e_i-1，e_i+1，…，e_M]^T，r_i表示回波信号的协方差矩阵R的第i列，A(Θ)表示给定的超完备空间，b_i表示系数矢量，σ²表示噪声功率，e_i表示误差项，e_i是第i个元素为1，其他元素均为零的M×1维的列矢量，M表示雷达天线阵元数，T表示转置操作；

1e)降维后的回波信号的协方差矩阵的表达式为：

其中，t_i表示降维后的回波信号的协方差矩阵的第i列，i＝1，2...，M，M表示雷达天线阵元数；

步骤2，根据降维后的协方差矩阵建立基于稀疏约束的稀疏矢量的代价函数；

步骤3，将基于稀疏约束的稀疏矢量的代价函数构造成适合于凸规划包求解的形式，根据凸规划包求解的形式求解稀疏矢量b₀；

步骤4，将稀疏矢量b₀中的非零元素确定为目标的方位角，该目标的方位角为目标的波达方向。

2.根据权利要求1所述的一种基于样本协方差矩阵稀疏性的波达方向估计方法，其特征在于，步骤2包括以下子步骤：

2a)设定l₁范数为：

其中，||·||₁表示取l₁范数操作，b₀(n)表示稀疏矢量b₀的第n个元素，n＝1，2...，N，N为远大于M的正整数，M表示雷达天线阵元数，|·|表示取绝对值操作；

2b)利用l₁范数作为稀疏矢量的代价函数，将降维后的回波信号协方差矩阵的第i列作为稀疏矢量的代价函数的稀疏约束条件，得到如下初始的基于稀疏约束的稀疏矢量的代价函数：

$\begin{array}{c}\min \left|\right|b_0|{|}_1\\ \begin{array}{cc}s.t.& t_i=J_{i}A\left(\mathrm{\Θ}\right)b_i,i=1,\mathrm{2...}M\end{array}\end{array}---\left(5\right)$

其中，min表示取最小值操作，s.t.表示取约束操作，b₀表示稀疏矢量，||·||₁表示取l₁范数操作，t_i表示降维后的回波信号协方差矩阵的第i列，i＝1，2...，M，J_i表示(M-1)×M维的变换选择矩阵，A(Θ)表示给定的超完备空间，b_i表示系数矢量，M表示雷达天线阵元数；

2c)用回波信号的采样数据样本对降维后的回波信号协方差矩阵的第i列t_i和回波信号的协方差矩阵R进行估计，获取基于稀疏约束的稀疏矢量的代价函数。

3.根据权利要求2所述的一种基于样本协方差矩阵稀疏性的波达方向估计方法，其特征在于，步骤2的子步骤2c)具体为：

2c1)利用回波信号的采样数据样本对降维后的回波信号协方差矩阵的第i列t_i进行估计，表示为下式：

${\hat{t}}_i=J_i(1/{\mathrm{L\Σ}}_{t=1}^{L}y(t){y_i}^H\left(t\right))=t_i+{\mathrm{\Δt}}_i$

其中，表示降维后的回波信号协方差矩阵的第i列t_i的样本估计值，J_i表示变换选择矩阵，且J_i为(M-1)×M维，y(t)表示t时刻M×1维的回波信号矢量，y_i(t)表示第i个雷达天线阵元t时刻接收到的信号矢量，t_i表示降维后的回波信号协方差矩阵的第i列，i＝1，2...，M，Δt_i为估计误差，L为样本的个数，M表示雷达天线阵元数，H表示共轭转置操作；

2c2)设定估计误差矢量

其中，Δt_i为估计误差，i＝1，2...，M，J是由变换选择矩阵J_i构造的块对角矩阵，表示回波信号的协方差矩阵R的样本估计值，R表示回波信号的协方差矩阵，vec(·)表示列矢量化操作，它是将矩阵的各列按列排成一个列矢量，M表示雷达天线阵元数，T表示转置操作；

估计误差矢量Δt满足如下分布：

$\mathrm{\Δ}t~AsN(0,J(\frac{1}{L}{R}^T\⊗R\left)J^T\right)---\left(6\right)$

其中，Δt表示估计误差矢量，～表示服从某种分布，AsN(μ，C)表示均值矢量为μ，协方差矩阵为C的近似多维正态分布；L为样本的个数，R表示回波信号的协方差矩阵，T表示转置操作，表示克罗内克积，J表示由变换选择矩阵J_i构造的块对角矩阵，其形式如下：

2c3)根据估计误差矢量Δt的协方差矩阵的逆W^-1，以及估计误差矢量Δt满足的分布，得到W^-1/2Δt满足如下分布：

W^-1/2Δt～AsN(0，I_M(M-1))

其中，估计误差矢量Δt的协方差矩阵的逆W^-1为：J表示由变换选择矩阵J_i构造的块对角矩阵，L为样本的个数，R表示回波信号的协方差矩阵，T表示转置操作，表示克罗内克积，(·)^-1表示取逆操作，W^-1/2表示W^-1的平方根矩阵，～表示服从某种分布，AsN(0，I_M(M-1))表示均值矢量为0，协方差矩阵为I_M(M-1)的近似多维正态分布，I_M(M-1)表示M×(M-1)维的单位阵，M表示雷达天线阵元数；

2c4)对W^-1/2Δt取2范数得到：

$\left|\right|W^{-1/2}\mathrm{\Δ}t|{|}_2^2~{\mathrm{As\χ}}^2\left(M\right(M-1\left)\right)$

其中，W^-1/2表示估计误差矢量Δt的协方差矩阵的逆W^-1的平方根矩阵，Δt表示估计误差矢量，||·||₂表示取2范数操作，～表示服从某种分布，Asχ²(M(M-1))表示自由度为M(M-1)的近似χ²分布，M表示雷达天线阵元数；

2c5)引入一个参数η，用来约束残差上限，使得以概率p成立；

其中，p大于等于0并且小于1；W^-1/2表示估计误差矢量Δt的协方差矩阵的逆W^-1的平方根矩阵，Δt表示估计误差矢量，||·||₂表示取2范数操作；

2c6)将概率p设定为p＝0.999，使得成立，则参数η根据χ²分布的概率累积函数以及p＝0.999求出，使得不等式以概率p＝0.999成立；

2c7)在不等式以概率p＝0.999成立的情况下，得到如下基于稀疏约束的稀疏矢量的代价函数，表示如下：

$\begin{array}{c}\min \left|\right|b_0|{|}_1\\ \left|\right|W^{-1/2}(J\·vec(\hat{R}-A\left(\mathrm{\Θ}\right)B\left)\right)|{|}_2^2\≤\mathrm{\η}\end{array}---\left(7\right)$

其中，min表示取最小值操作，b₀表示稀疏矢量，||·||₁表示取l₁范数操作，||·||₂表示取2范数操作，J表示由变换选择矩阵J_i构造的块对角矩阵，vec(·)表示列矢量化操作，A(Θ)表示给定的超完备空间，B是由系数矢量b_i构成的矩阵，B＝[b₁，b₂，…，b_M]，i＝1，2...，M，M表示雷达天线阵元数，η表示约束残差上限的参数，表示回波信号的协方差矩阵R的样本估计值，即x(t)表示t时刻的回波信号的采样数据样本，H表示共轭转置操作；用来逼近估计误差矢量Δt的协方差矩阵的逆W^-1的平方根矩阵W^-1/2，即L为样本的个数，表示克罗内克积。

4.根据权利要求3所述的一种基于样本协方差矩阵稀疏性的波达方向估计方法，其特征在于，步骤3包括以下子步骤：

3a)将基于稀疏约束的稀疏矢量的代价函数式(7)构造成如下的适合于凸规划包求解的形式：

$\begin{array}{cc}& \min g\\ s.t.& 1^T\mathrm{\γ}\≤g,\\ & {\left({\Σ}_{m=1}^{M}B{(n,m)}^2\right)}^{1/2}=b_0\left(n\right)\≤{\mathrm{\γ}}_m,n=1,\mathrm{2...}N\\ & \left|\right|{\hat{W}}^{-1/2}(J\·vec(\hat{R}-A\left(\mathrm{\Θ}\right)B\left)\right)|{|}_2^2\≤\mathrm{\η}\end{array}---\left(8\right)$

其中，min表示取最小值操作，s.t.表示取约束操作，g表示松弛变量，1表示N×1维的全1矢量，中间变量γ＝[γ₁，γ₂，…γ_N]^T，B(n，m)表示矩阵B的第n行第m个元素，B是由系数矢量b_i构成的矩阵，B＝[b₁，b₂，…，b_M]，m＝1，2...，M，n＝1，2...，N，N为远大于M的正整数，b₀(n)表示稀疏矢量b₀的第n个元素，||·||₂表示取2范数操作，表示W^-1/2的逼近值，W^-1/2表示估计误差矢量Δt的协方差矩阵的逆W^-1的平方根矩阵，J表示由变换选择矩阵J_i构造的块对角矩阵，vec(·)表示列矢量化操作，表示回波信号的协方差矩阵R的样本估计值，A(Θ)表示给定的超完备空间，M表示雷达天线阵元数，η表示约束残差上限的参数，T表示转置操作；

3b)用凸规划包求解式(8)，得到稀疏矢量b₀。