CN103983944A - 基于协方差矩阵稀疏表示的远场窄带doa估计方法 - Google Patents

Abstract

本发明提供一种高精度远场窄带DOA估计方法。在波达方向在空域上具有稀疏性的基础上，把协方差矩阵改写成稀疏表示的模型，在网格匹配的模型下，通过优化最小化方法求解出稀疏的空间功率谱，则该功率谱的支撑集在网格上对应的点即为估计得到的波达方向角。针对实际中真实的波达方向角不在网格上的情况，即网格失配的模型下，利用一阶泰勒展开来逼近真实波达方向的导向矢量，然后用最小二乘方法对估计得到的网格上的点进行修正，以达到更高的估计精度。本发明可以在粗糙的网格上到达高精度的DOA估计性能。

Description

基于协方差矩阵稀疏表示的远场窄带DOA估计方法

技术领域

本发明属于阵列信号处理领域，主要涉及远场窄带DOA估计。

背景技术

波达方向(DOA)估计一直是阵列信号处理中一个重要的研究领域，它在雷达、声纳、无线通信及电子对抗和侦查等领域中都有着广泛的应用。如何快速地，高精度地实现DOA估计一直是阵列信号处理不断研究和努力的方向。其中经典的算法有：多重信号分类(Multiple Signal Classification,MUSIC)算法、旋转不变子空间(Estimation of Signal Parameters via Rotational Invariance Technique,ESPRIT)算法等子空间类算法和最大似然估计类算法(Maximum Likelihood,ML)等。然而，基于子空间理论的DOA估计方法虽然实现了超分辨侧向，但是一旦阵列快拍数不足或者出现相干信号源时，这类方法不能有效地区分信号子空间和噪声子空间，其性能会急剧下降。而最大似然估计类算法由于要进行复杂的多维搜索而不具有实用性。

近年来，基于压缩感知的稀疏表示理论的兴起和发展为DOA估计问题提供了一种新思路。只要信号在某个变换域下具有稀疏性，就可以利用相应的算法以极高的概率精确地重构原始稀疏信号。在空间谱估计的阵列模型中，通常假设空域范围内只存在少数的点目标，因此，波达方向角在空域上具有稀疏性。标准的压缩感知理论中，使用l₀范数作为目标函数来保证重构信号的稀疏性，但是最小化l₀范数是一个组合优化问题，难以有效地求解，后来研究证明将l₀范数松弛为l₁范数依然能保证重构信号的稀疏性，因此如何用l₁范数来精确重构稀疏信号得到了广泛的研究。

基于稀疏重构进行DOA估计的研究工作一直非常活跃。Malioutov等基于稀疏重构思想提出了l₁-SVD算法，该算法最突出的贡献是在多快拍条件下通过奇异值分解(Singular Value Decomposition,SVD)来减小数据矩阵的规模以及降低噪声的影响，使得该算法的计算量不会随着快拍数的增加而增加，而且具有较高的估计精度。然而，这些基于稀疏表示的算法通常假设所有的真实波达方向角都位于预先设定的离散化网格上，即网格匹配的模型，这导致了这类算法在波达方向角不在网格上时估计性能急剧恶化。另一方面，虽然更密集的网格理论上可以减小重构误差，但是太过密集的离散化网格会使得过完备化字典原子间高度相关。

发明内容

本发明的目的在于提供一种基于协方差矩阵稀疏表示的远场窄带DOA估计方法。在网格匹配的情况下可以精确地估计波达方向角，在网格失配的情况下，能够对角度值进行修正，提高估计精度。

本发明的思路是：本发明基于空间协方差矩阵稀疏表示的模型，首先在网格匹配的假设下用最优最小化方法估计得到网格上的角度值，然后针对网格失配的情况，通过一阶泰勒展开来逼近真实的导向矢量，进而对估计得到的网格点进行修正。

本发明的目的通过如下步骤实现：

S1、由阵列接收的K个信号源的数据 $x\left(t\right)=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}a\left({\mathrm{\θ}}_k\right)s_k\left(t\right)+n\left(t\right)=A\left(\mathrm{\θ}\right)s\left(t\right)+n\left(t\right),$ 得到空间协方差矩阵R＝E[x(t)x^H(t)]＝A(θ)R_sA^H(θ)+σ²I_M，其中，x(t)＝[x₁(t),x₂(t),...,x_M(t)]^T表示各个阵元接收信号构成的矩阵，M为阵元数目，K为远场窄带信号源个数，θ_k为第k个信号源入射到阵列的角度，为第k个信号源的导向矢量，为第k个信号源入射到第m个阵元与所述第k个信号源入射到参考阵元的相位差，λ为入射信号的波长，d为相邻两个阵元的间距，A(θ)＝[a(θ₁),a(θ₂),...,a(θ_K)]为阵列流行矩阵，s(t)＝[s₁(t),s₂(t),...,s_K(t)]^T为入射信号，附加噪声n(t)为与各个信号源不相关的加性零均值高斯白噪声，协方差矩阵R中，R_s＝diag(r₁,r₂,...,r_K)为入射信号的功率，σ²为噪声功率，I_M为M阶的单位矩阵，E[·]表示期望，(·)^H表示矩阵的转置，(·)^T表示矩阵的共轭转置，k＝1,2,...,K，m＝1,2,...,M；

S2、对S1所述空间协方差矩阵R进行向量化操作，写成稀疏表示的模型r，具体如下：

S21、将S1所述空间协方差矩阵R依次按列排列，写成向量的形式，即 $\begin{array}{c}r=\mathrm{vec}\left(R\right)\\ =\mathrm{vec}(A\left(\mathrm{\θ}\right)R_s{A}^H\left(\mathrm{\θ}\right)+{\mathrm{\σ}}^2{I}_M)\\ =G\left(\mathrm{\θ}\right)r_s+{\mathrm{\σ}}^2{I}_v\end{array},$ 其中，g(θ_k)＝vec(a(θ_k)a^H(θk))，vec(·)表示向量化操作，r_s为矩阵R_s对角线元素构成的向量，I_v为单位矩阵I_M按列排列得到的向量，表示M²行K列的复矩阵；

S22、把S1所述角度θ_k在[-90°,90°)的空间范围上过完备化为一个离散的网格即将向量写成稀疏表示的形式其中，过完备矩阵 是一个K稀疏向量，N＞＞K，表示N行1列的向量；

S3、对噪声功率σ²进行估计，得到噪声功率的估计值其中，λ_i为S1所述协空间方差矩阵R中M-K+1个最小的特征值；

S4、利用S3所述噪声功率的估计值通过优化最小化(majorization-minimization，MM)的方法求解出网络上真实波达方向最近的点，具体为：

S41、用优化最小化方法求解优化式得出稀疏空间功率谱 ${\tilde{r}}_s^{l+1}=\arg \min \underset{i}{\overset{N}{\mathrm{\Σ}}}\frac{\left|{\tilde{r}}_i\right|}{\left|{\tilde{r}}_i^l\right|+\mathrm{\ξ}}+\mathrm{\β}{\left|\right|r-{\widehat{\mathrm{\σ}}}^2{I}_v-G\left(\widetilde{\mathrm{\θ}}\right){\tilde{r}}_s\left|\right|}_2,$ 其中，符号≥表示对向量的每个元素进行操作，0为全零列向量，β为正则化参数，代表第l次迭代得到的空间谱，0＜ξ且ξ为常数；

S42、估计得到的波达方向角，具体如下：

S421、当真实的波达方向角θ_k在预先设定的离散化网格上时，则当满足迭代停止条件或者达到最大迭代次数时，得到稀疏的空间谱则最大的K个元素在网格上对应的点就是估计得到的波达方向角，其中，所述网格上对应的位置记作支撑集Λ；

S422、当真实的波达方向角不在预先设定的离散化网格上时，即需要对S421得到的结果进行修正，得到修正后的波达方向角 ${\widehat{\mathrm{\θ}}}_k=\arcsin [\sin {\widetilde{\mathrm{\θ}}}_{i_k}+{\mathrm{\δ}}_{i_k}].$

进一步地，S41所述β＝0.5。

进一步地，S41所述ξ＝0.1。

进一步地，S422所述对S421得到的结果进行修正，具体方法如下：

步骤1、利用一阶泰勒展开对真实的导向矢量进行逼近， $a\left(\sin {\mathrm{\θ}}_k\right)\≈a\left(\sin {\widetilde{\mathrm{\θ}}}_i\right)+a^{\′}\left(\sin {\widetilde{\mathrm{\θ}}}_i\right)(\sin {\widetilde{\mathrm{\θ}}}_i-{\sin \mathrm{\θ}}_k),$ 其中，为网格上离θ_k最近的点，为导向矢量的一阶导数；

步骤2、记代入空间协方差矩阵向量化的模型r中，即其中，Δ＝diag(δ)，δ＝[δ₁,δ₂,...,δ_N]^T，符号⊙表示向量对应元素分别相乘；

步骤3、根据S421所得支撑集Λ，则通过最小二乘法得到则δ的支撑集上的元素其中，表示对r_s的每一个元素分别取倒数构成的向量；

步骤4、假设步骤3所述δ_Λ在网格上的索引值为i₁,...,i_K，则修正后的波达方向角的估计值为 ${\widehat{\mathrm{\θ}}}_k=\arcsin [\sin {\widetilde{\mathrm{\θ}}}_{i_k}+{\mathrm{\δ}}_{i_k}],k=\mathrm{1,2},...,K.$

本发明的有益效果是：

可在粗糙的网格上进行波达方向的精确估计，避免了密集的网格带来的高计算量，提高了估计精度。利用本发明的估计方法，修正后的估计性能有明显的提高，即使在粗糙的网格，且快拍数为100时都能达到0.1°以内的估计精度。本发明的估计方法与其他估计方法相比，估计性能有明显的提高。

附图说明

图1是本发明方法的流程图。

图2是远场窄带信号接收阵列模型图。

图3是网格匹配时，本发明方法与其他方法DOA估计的均方根误差随信噪比变化曲线图。

图4是网格匹配时，本发明方法与其他方法DOA估计的均方根误差随快拍数变化曲线图。

图5是网格失配时，本发明方法与其他方法DOA估计的均方根误差随信噪比变化曲线图。

图6是网格失配时，本发明方法与其他方法DOA估计的均方根误差随快拍数变化曲线图。

具体实施方式

下面结合实施例和附图，详细说明本发明的技术方案。

图1是本发明基于协方差矩阵稀疏表示的高精度DOA估计方法的一种具体实施方式流程图。如图1所示，本发明基于协方差矩阵稀疏表示的高精度DOA估计方法包括以下步骤：

S1、由阵列接收的K个信号源的数据 $x\left(t\right)=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}a\left({\mathrm{\θ}}_k\right)s_k\left(t\right)+n\left(t\right)=A\left(\mathrm{\θ}\right)s\left(t\right)+n\left(t\right),$ 得到空间协方差矩阵R＝E[x(t)x^H(t)]＝A(θ)R_sA^H(θ)+σ²I_M，其中，x(t)＝[x₁(t),x₂(t),...,x_M(t)]^T表示各个阵元接收信号构成的矩阵，M为阵元数目，K为远场窄带信号源个数，θ_k为第k个信号源入射到阵列的角度，为第k个信号源的导向矢量，为第k个信号源入射到第m个阵元与所述第k个信号源入射到参考阵元的相位差，λ为入射信号的波长，d为相邻两个阵元的间距，A(θ)＝[a(θ₁),a(θ₂),...,a(θ_K)]为阵列流行矩阵，s(t)＝[s₁(t),s₂(t),...,s_K(t)]^T为入射信号，附加噪声n(t)为与各个信号源不相关的加性零均值高斯白噪声，协方差矩阵R中，R_s＝diag(r₁,r₂,...,r_K)为入射信号的功率，σ²为噪声功率，I_M为M阶的单位矩阵，E[·]表示期望，(·)^H表示矩阵的转置，(·)^T表示矩阵的共轭转置，k＝1,2,...,K，m＝1,2,...,M；

S2、对S1所述空间协方差矩阵R进行向量化操作，写成稀疏表示的模型r，具体如下：

S21、将S1所述空间协方差矩阵R依次按列排列，写成向量的形式，即 $\begin{array}{c}r=\mathrm{vec}\left(R\right)\\ =\mathrm{vec}(A\left(\mathrm{\θ}\right)R_s{A}^H\left(\mathrm{\θ}\right)+{\mathrm{\σ}}^2{I}_M)\\ =G\left(\mathrm{\θ}\right)r_s+{\mathrm{\σ}}^2{I}_v\end{array},$ 其中，g(θ_k)＝vec(a(θ_k)a^H(θ_k))，vec(·)表示向量化操作，r_s为矩阵R_s对角线元素构成的向量，I_v为单位矩阵I_M按列排列得到的向量；

S22、把S1所述角度θ_k在[-90°,90°)的空间范围上过完备化为一个离散的网格即将向量写成稀疏表示的形式其中，过完备矩阵 是一个K稀疏向量，N＞＞K；

S3、对噪声功率σ²进行估计，得到噪声功率的估计值其中，λ_i为S1所述协空间方差矩阵R中M-K+1个最小的特征值；

S4、利用S3所述噪声功率的估计值通过优化最小化(majorization-minimization，MM)的方法求解出网络上真实波达方向最近的点，具体为：

S41、用优化最小化方法求解优化式得出稀疏空间功率谱 ${\tilde{r}}_s^{l+1}=\arg \min \underset{i}{\overset{N}{\mathrm{\Σ}}}\frac{\left|{\tilde{r}}_i\right|}{\left|{\tilde{r}}_i^l\right|+\mathrm{\ξ}}+\mathrm{\β}{\left|\right|r-{\widehat{\mathrm{\σ}}}^2{I}_v-G\left(\widetilde{\mathrm{\θ}}\right){\tilde{r}}_s\left|\right|}_2,$ 其中，符号≥表示对向量的每个元素进行操作，0为全零列向量，β为正则化参数，β＝0.5，代表第l次迭代得到的空间谱，ξ为大于零的正常数，ξ＝0.1；

S42、估计得到的波达方向角，具体如下：

S421、当真实的波达方向角θ_k在预先设定的离散化网格上时，则当满足迭代停止条件或者达到最大迭代次数时，得到稀疏的空间谱则最大的K个元素在网格上对应的点就是估计得到的波达方向角，其中，所述网格上对应的位置记作支撑集Λ；

S422、当真实的波达方向角不在预先设定的离散化网格上时，即需要对S421得到的结果进行修正，得到修正后的波达方向角 ${\widehat{\mathrm{\θ}}}_k=\arcsin [\sin {\widetilde{\mathrm{\θ}}}_{i_k}+{\mathrm{\δ}}_{i_k}],$ 具体为：

步骤1、利用一阶泰勒展开对真实的导向矢量进行逼近， $a\left(\sin {\mathrm{\θ}}_k\right)\≈a\left(\sin {\widetilde{\mathrm{\θ}}}_i\right)+a^{\′}\left(\sin {\widetilde{\mathrm{\θ}}}_i\right)(\sin {\widetilde{\mathrm{\θ}}}_i-{\sin \mathrm{\θ}}_k),$ 其中，为网格上离θ_k最近的点，为导向矢量的一阶导数；

步骤2、记代入空间协方差矩阵向量化的模型r中，即其中，Δ＝diag(δ)，δ＝[δ₁,δ₂,...,δ_N]^T，则δ与联合稀疏，即非零元素的位置相同，符号⊙表示向量对应元素分别相乘；

步骤3、根据S421所得支撑集Λ，则和也即已知，可以通过最小二乘法得到则δ的支撑集上的元素其中，表示对r_s的每一个元素分别取倒数构成的向量；

步骤4、假设步骤3所述δ_Λ在网格上的索引值为i₁,...,i_K，则修正后的波达方向角的估计值为 ${\widehat{\mathrm{\θ}}}_k=\arcsin [\sin {\widetilde{\mathrm{\θ}}}_{i_k}+{\mathrm{\δ}}_{i_k}],k=\mathrm{1,2},...,---K.$

实施例1和实施例2的仿真条件为网格匹配的模型，因此估计值即为最终的波达方向估计值，无需进行修正。实施例3和实施例4的仿真条件为网格失配的模型。实施例中用均方根误差(RMSE)来评估各算法的性能，其定义为： $\mathrm{RMSE}=\frac{1}{\mathrm{Mon}}\underset{m=1}{\overset{\mathrm{Mon}}{\mathrm{\Σ}}}\sqrt{\frac{1}{K}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{({\widehat{\mathrm{\θ}}}_{k,m}-{\mathrm{\θ}}_k)}^2},$ 其中，Mon为蒙特卡洛实验次数，和θ_k分别代表第m次蒙特卡洛实验估计得到的第k个角度和第k个真实角度。

实施例1

网格匹配的情况下，本发明估计值的均方根误差随信噪比变化仿真：

实施例1采用的接收阵列如图2所示的由8个阵元组成的半波长均匀线性阵列。参考阵元为编号1的阵元天线。四个相同功率的信号源按入射方向[-35°,-10°,15°,40°]入射到阵列。为了使入射方向角落在网格上，取离散化网格为{-90°,-89°,...,89°}，间隔1°。采样快拍数为200。参考信噪比SNR从-10dB到10dB变化，间隔为4dB，每个信噪比进行1000次蒙特卡洛实验。

实施例1中DOA估计方法包括以下步骤：

根据不同信噪比下的阵列接收信号x(t)得到协方差矩阵R；

对R进行特征值分解得到噪声功率的估计值然后对R进行向量化得到r，其中M＝8，K＝4；

通过优化最小化算法求解得到稀疏空间功率谱并找到其4个峰值对应的索引值，得到其支撑集Λ并得到网格上的点即估计得到的角度：

按照本发明的方法估计得到的波达方向角的均方根误差随信噪比变化的曲线如图3所示。图3可以看到，利用本发明的估计方法，估计性能会随着信噪比的增加而显著提高，当信噪比大于0dB时，估计精度可以达到0.1°以内。本发明的估计方法与其他估计方法相比，估计性能有明显的提高，说明了本发明的估计方法是有效的。

实施例2

网格匹配的情况下，本发明估计值的均方根误差随快拍数变化仿真：

实施例1采用的接收阵列为如附图2所示的由8个阵元组成的半波长均匀线性阵列。参考阵元为编号1的阵元天线。四个相同功率的信号源按入射方向[-35°,-10°,15°,40°]入射到阵列。为了使入射方向角落在网格上，取离散化网格为{-90°,-89°,...,89°}，间隔1°。参考信噪比SNR固定为10dB。快拍数从100到600，间隔100，每个快拍数进行1000次蒙特卡洛实验。

实施例2中DOA估计方法包括以下步骤：

根据不同快拍数下的阵列接收信号x(t)得到协方差矩阵R；

对R进行特征值分解得到噪声功率的估计值然后对R进行向量化得到r，其中M＝8，K＝4；

通过优化最小化算法求解得到稀疏空间功率谱并找到其4个峰值对应的索引值，得到其支撑集Λ并得到网格上的点即估计得到的角度：

按照本发明的方法估计得到的波达方向角的均方根误差随快拍数变化的曲线如图4所示。图4可以看到，利用本发明的估计方法，估计性能会随着快拍数的增加而提高。当快拍数为100时，估计精度已经达到0.2°，当快拍数大于300时，估计均方根误差基本保持不变。本发明的估计方法与其他估计方法相比，估计性能有明显的提高，说明了本发明的估计方法是有效的。

实施例3

网格失配的情况下，本发明估计值的均方根误差随信噪比变化仿真：

实施例3采用的接收阵列为如附图2所示的由8个阵元组成的半波长均匀线性阵列。参考阵元为编号1的阵元天线。两个相同功率的信号源按入射方向[-14.5°,36.3°]入射到阵列。为了使入射方向角不落在网格上，取离散化网格为{-90°,-88°,...,88°}，间隔2°。采样快拍数为200。参考信噪比SNR从-4dB到20dB变化，间隔为4dB，每个信噪比进行1000次蒙特卡洛实验。

实施例3中DOA估计方法包括以下步骤：

根据不同信噪比下的阵列接收信号x(t)得到协方差矩阵R；

对R进行特征值分解得到噪声功率的估计值然后对R进行向量化得到r，其中M＝8，K＝2；

通过优化最小化算法求解得到稀疏空间功率谱并找到其4个峰值对应的索引值，得到其支撑集Λ并得到网格上的点：

通过最小二乘方法求解上一步骤波达方向角的修正值

得到修正后的角度 ${\widehat{\mathrm{\θ}}}_k=\arcsin [\sin {\widetilde{\mathrm{\θ}}}_{i_k}+{\mathrm{\δ}}_{i_k}],k=\mathrm{1,2}.$

按照本发明的方法估计得到的波达方向角的均方根误差随信噪比变化的曲线如图5所示。图5可以看到，利用本发明的估计方法，修正后的估计性能有明显的提高，即使在粗糙的网格，且信噪比为0dB时都能达到0.15°以内的估计精度。本发明的估计方法与其他估计方法相比，估计性能有明显的提高，说明了本发明的估计方法是有效的。

实施例4

网格失配的情况下，本发明估计值的均方根误差随快拍数变化仿真：

实施例4采用的接收阵列为如附图2所示的由8个阵元组成的半波长均匀线性阵列。参考阵元为编号1的阵元天线。两个相同功率的信号源按入射方向[-14.5°°,36.3°]入射到阵列。为了使入射方向角不落在网格上，取离散化网格为{-90°,-88°,...,88°}，间隔2°。参考信噪比SNR固定为10dB。快拍数从100到400，间隔50，每个快拍数进行1000次蒙特卡洛实验。

实施例4中DOA估计方法包括以下步骤：

根据不同快拍数下的阵列接收信号x(t)得到协方差矩阵R；

对R进行特征值分解得到噪声功率的估计值然后对R进行向量化得到r，其中M＝8，K＝2；

通过优化最小化算法求解得到稀疏空间功率谱并找到其4个峰值对应的索引值，得到其支撑集Λ并得到网格上的点：

通过最小二乘方法求解上一步骤波达方向角的修正值

得到修正后的角度 ${\widehat{\mathrm{\θ}}}_k=\arcsin [\sin {\widetilde{\mathrm{\θ}}}_{i_k}+{\mathrm{\δ}}_{i_k}],k=\mathrm{1,2}.$

按照本发明的方法估计得到的波达方向角的均方根误差随快拍数变化的曲线如图6所示。图6可以看到，利用本发明的估计方法，修正后的估计性能有明显的提高，即使在粗糙的网格，且快拍数为100时都能达到0.1°以内的估计精度。本发明的估计方法与其他估计方法相比，估计性能有明显的提高，说明了本发明的估计方法是有效的。

Claims (4)

1.基于协方差矩阵稀疏表示的远场窄带DOA估计方法，其特征在于，包括如下步骤：

S1、由阵列接收的K个信号源的数据 $x\left(t\right)=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}a\left({\mathrm{\θ}}_k\right)s_k\left(t\right)+n\left(t\right)=A\left(\mathrm{\θ}\right)s\left(t\right)+n\left(t\right),$ 得到空间协方差矩阵R＝E[x(t)x^H(t)]＝A(θ)R_sA^H(θ)+σ²I_M，其中，x(t)＝[x₁(t),x₂(t),...,x_M(t)]^T表示各个阵元接收信号构成的矩阵，M为阵元数目，K为远场窄带信号源个数，θ_k为第k个信号源入射到阵列的角度，为第k个信号源的导向矢量，为第k个信号源入射到第m个阵元与所述第k个信号源入射到参考阵元的相位差，λ为入射信号的波长，d为相邻两个阵元的间距，A(θ)＝[a(θ₁),a(θ₂),...,a(θ_K)]为阵列流行矩阵，s(t)＝[s₁(t),s₂(t),...,s_K(t)]^T为入射信号，附加噪声n(t)为与各个信号源不相关的加性零均值高斯白噪声，协方差矩阵R中，R_s＝diag(r₁,r₂,...,r_K)为入射信号的功率，σ²为噪声功率，I_M为M阶的单位矩阵，E[·]表示期望，(·)^H表示矩阵的转置，(·)^T表示矩阵的共轭转置，k＝1,2,...,K，m＝1,2,...,M；

S2、对S1所述空间协方差矩阵R进行向量化操作，写成稀疏表示的模型r，具体如下：

S21、将S1所述空间协方差矩阵R依次按列排列，写成向量的形式，即 $\begin{array}{c}r=\mathrm{vec}\left(R\right)\\ =\mathrm{vec}(A\left(\mathrm{\θ}\right)R_s{A}^H\left(\mathrm{\θ}\right)+{\mathrm{\σ}}^2{I}_M)\\ =G\left(\mathrm{\θ}\right)r_s+{\mathrm{\σ}}^2{I}_v\end{array},$ 其中，g(θ_k)＝vec(a(θ_k)a^H(θ_k))，vec(·)表示向量化操作，r_s为矩阵R_s对角线元素构成的向量，I_v为单位矩阵I_M按列排列得到的向量，表示维数为M²×K的复矩阵；

S22、把S1所述角度θ_k在[-90°,90°)的空间范围上过完备化为一个离散的网格即将向量写成稀疏表示的形式其中，过完备矩阵 是一个K稀疏向量，N＞＞K；

S3、对噪声功率σ²进行估计，得到噪声功率的估计值其中，λ_i为S1所述协空间方差矩阵R中M-K个最小的特征值；

S4、利用S3所述噪声功率的估计值通过优化最小化(majorization-minimization，MM)的方法求解出网络上真实波达方向最近的点，具体为：

S41、用优化最小化方法求解优化式得出稀疏空间功率谱 ${\tilde{r}}_s^{l+1}=\arg \min \underset{i}{\overset{N}{\mathrm{\Σ}}}\frac{\left|{\tilde{r}}_i\right|}{\left|{\tilde{r}}_i^l\right|+\mathrm{\ξ}}+\mathrm{\β}{\left|\right|r-{\widehat{\mathrm{\σ}}}^2{I}_v-G\left(\widetilde{\mathrm{\θ}}\right){\tilde{r}}_s\left|\right|}_2,$ 其中，符号≥表示对向量的每个元素进行操作，0为全零列向量，β为正则化参数，代表第l次迭代得到的空间谱，0＜ξ且ξ为常数；

S42、估计得到的波达方向角，具体如下：

S421、当真实的波达方向角θ_k在预先设定的离散化网格上时，则当满足迭代停止条件或者达到最大迭代次数时，得到稀疏的空间谱则最大的K个元素在网格上对应的点就是估计得到的波达方向角，其中，所述网格上对应的位置记作支撑集Λ；

S422、当真实的波达方向角不在预先设定的离散化网格上时，即需要对S421得到的结果进行修正，得到修正后的波达方向角 ${\widehat{\mathrm{\θ}}}_k=\arcsin [\sin {\widetilde{\mathrm{\θ}}}_{i_k}+{\mathrm{\δ}}_{i_k}].$

2.根据权利要求1所述的基于协方差矩阵稀疏表示的远场窄带DOA估计方法，其特征在于：S41所述β＝0.5。

3.根据权利要求1或2所述的基于协方差矩阵稀疏表示的远场窄带DOA估计方法，其特征在于：ξ＝0.1。

4.根据权利要求1所述的基于协方差矩阵稀疏表示的远场窄带DOA估计方法，其特征在于：S422所述对S421得到的结果进行修正，具体方法如下：

步骤1、利用一阶泰勒展开对真实的导向矢量进行逼近， $a\left(\sin {\mathrm{\θ}}_k\right)\≈a\left(\sin {\widetilde{\mathrm{\θ}}}_i\right)+a^{\′}\left(\sin {\widetilde{\mathrm{\θ}}}_i\right)(\sin {\widetilde{\mathrm{\θ}}}_i-{\sin \mathrm{\θ}}_k),$ 其中，为网格上离θ_k最近的点，为导向矢量的一阶导数；

步骤2、记代入空间协方差矩阵向量化的模型r中，即其中，Δ＝diag(δ)，δ＝[δ₁,δ₂,...,δ_N]^T，符号⊙表示向量对应元素分别相乘；

步骤3、根据S421所得支撑集Λ，则通过最小二乘法得到则δ的支撑集上的元素其中，表示对r_s的每一个元素分别取倒数构成的向量；

步骤4、假设步骤3所述δ_Λ在网格上的索引值为i₁,...,i_K，则修正后的波达方向角的估计值为 ${\widehat{\mathrm{\θ}}}_k=\arcsin [\sin {\widetilde{\mathrm{\θ}}}_{i_k}+{\mathrm{\δ}}_{i_k}],k=\mathrm{1,2},...,K.$