CN105303009A - 基于压缩感知与正则mfocuss的超分辨谱估计方法 - Google Patents

Abstract

本发明属于信号处理技术领域，具体提供了一种基于压缩感知与正则MFOCUSS的超分辨谱估计方法，包括如下步骤：步骤1.构造接收的阵列信号，步骤2.构造稀疏字典和一个稀疏向量，步骤3.初始化迭代参数和正则化参数，步骤4.由正则化的MFOCUSS算法求解，当满足误差终止条件时，结束循环，否则，更新权值参数，步骤5.根据角度谱的计算方法得出最终的DOA估计。本发明在天线阵元数较少的情况下，相比传统的MUSIC算法和MFOCUSS算法取得了较高分辨率，同时正则化参数的引入使得算法在解的稀疏度和误差之间有一个合理的权衡，而自适应参数的引入进一步提高了解的精度。

Description

基于压缩感知与正则MFOCUSS的超分辨谱估计方法

技术领域

本发明属于信号处理技术领域，具体涉及一种基于压缩感知与正则MFOCUSS的超分辨谱估计方法，可应用于天线阵元数目较少时的高分辨率DOA估计中。

背景技术

压缩感知是近几年在信号处理领域中发展起来的一种关于信号传输与存储的新理论，它突破了解决传统奈奎斯特采样中关于采样速率的限制，可在低采样率下实现信息的精确感知。

波达方向(DOA)估计是高分辨率阵列信号处理的重要内容之一，广泛应用于雷达、导航、声呐等诸多领域。目前有效的DOA估计算法主要有Capon提出的最小方差谱估计(MVDR)和Schimidt提出的多重信号分类法(MUSIC)为代表的子空间算法。但是，MVDR方法虽然突破了阵列孔径的内在限制以及不需要事先预知信号源个数，但是由于受到瑞利限的制约，无法分辨一个波束宽度内的两个信号，而子空间算法虽然具有高分辨率的特性，但是由于相干信号的信号子空间和噪声子空间的相互渗透，使得子空间算法不能识别出相干信号源。实际上，在DOA估计中，感兴趣的目标(即信源)仅占据少量的角度分辨单元，因而信源在空域中是稀疏的，从而可以将压缩感知应用于DOA估计中。

多测量矢量欠定系统聚焦求解(MFOCUSS)算法实现了DOA的高分辨率估计，但是该方法在快拍数较多的情况下运算量急剧增加，而且在低信噪比，即噪声很大的情况下无法适用。因此，在天线单元较少的情况下，如何使DOA估计在取得了较高分辨率的同时具有更好的抗噪性能以及较高的精度是重要的研究方向。

发明内容

本发明的目的是克服上述现有技术的缺点，提出了一种基于压缩感知与正则MFOCUSS的超分辨谱估计方法，以在解的稀疏性和误差之间有一个合理的权衡,从而得到高分辨率的，对噪声比较鲁棒的DOA估计。

为实现上述目的，本发明提供了一种基于压缩感知与正则MFOCUSS的超分辨谱估计方法，包括以下步骤：

步骤1.构造接收的阵列信号，具体过程如下：

1a)设定天线阵元数目N,快拍数M；

1b)第j(j＝1,2,…p)个信源的信号其中u_j(t)和为第j个接收信号的幅度和相位，ω₀为接收信号的频率，则第i个阵元接收到的信号为其中g_ij为第i个阵元对第j个信号的增益，n_i(t)为第i个阵元在t时刻的噪声，τ_ij表示第j个信号到达第i个阵元时相对于参考阵元的延迟；

1c)构造接收信号组成的矩阵Y＝[y₁(k)y₂(k)…y_N(k)]^T,作为N个天线阵元的输出，其中y_i(k)＝[y_i(1×t/M)y_i(2×t/M)…y_i(M×t/M)]，k＝1,2,…M；

步骤2.构造稀疏字典和一个稀疏向量，具体过程如下：

2a)构造N×N_s的角度扫描矩阵A，即稀疏字典 $A=\left[\begin{array}{cccc}a\left({\mathrm{\θ}}_1\right)& a\left({\mathrm{\θ}}_2\right)& \·\·\·& a\left({\mathrm{\θ}}_{N_s}\right)\end{array}\right],$ 其中 $a\left({\mathrm{\θ}}_i\right)={\left[\begin{array}{cccc}1& e^{-{\mathrm{j\ω}}_i}& \·\·\·& e^{-j(N-1){\mathrm{\ω}}_i}\end{array}\right]}^T,$ ω_i＝2πsin(θ_i)/λ,i＝1,2,…N_s，λ为接收信号对应的波长，N_s为要观测的全部角度θ_i的数目；

2b)构造一个N_s×1的稀疏向量 $x\left(t\right)={\left[\begin{array}{cccc}{x}_{{\mathrm{\θ}}_1}\left(t\right)& x_{{\mathrm{\θ}}_2}\left(t\right)& \·\·\·& x_{{\mathrm{\θ}}_{\mathrm{Ns}}}\left(t\right)\end{array}\right]}^T,$ 其中有p个非零的系数x_θ(t)＝s_p(t)在对应于信号源角度的位置，剩下的N_s-p个系数为零，进而构造M个快拍组成的N_s×M的矩阵X＝[x(1)x(2)…x(M)]；

步骤3.初始化迭代参数和正则化参数，具体过程如下：

3a)设定初始解X⁽⁰⁾＝0；

3b)初始化迭代次数k＝0，自适应参数p_k＝1,误差阈值ε＝1e-8，总的迭代限制次数Iter＝800，正则化参数γ＝σ²，σ²为噪声方差；

步骤4.由正则化的MFOCUSS算法求解，当满足误差终止条件时，结束循环，否则，更新权值参数，具体过程如下：

4a)利用式X^(k+1)＝W^(k+1)Q^(k+1)计算第k+1次迭代求解出的目标矩阵，其中

$W^{(k+1)}=\mathrm{diag}\left({\left(c_i^{\left(k\right)}\right)}^{1-p/2}\right),p\∈\left[\mathrm{0,1}\right],$

$c_i^{\left(k\right)}=\left|\right|X_i^{\left(k\right)}\left|\right|={\left({\mathrm{\Σ}}_{i=1}^L{\left|X_i^{\left(k\right)}\left(l\right)\right|}^2\right)}^{1/2},$

$Q^{(k+1)}={\left(A^{(k+1)}\right)}^{+}{(A^{(k+1)}{\left(A^{(k+1)}\right)}^{+}+\mathrm{\γI})}^{-1}Y,\mathrm{\γ}\≥0,$

I为N阶单位矩阵，A^(k+1)＝AW^(k+1)

4b)计算相对误差rRMSE_k＝||X^(k)-X^(k-1)||_F/||X^(k)||_F，||·||_F表示先对矩阵各元素求平方和，再开2次方，若rRMSE_k＜ε或k≥Iter,结束循环，否则k＝k+1，p_k＝min(rRMSE_k,1)，重复4a)和4b)；

步骤5.根据角度谱的计算方法得出最终的DOA估计，具体过程如下：

5a)根据RMSE＝||X^(k)-X||_F计算最终的均方误差；

5b)根据角度谱的计算公式求出角度谱，其中 $\mathrm{\θ}={\mathrm{\θ}}_1,{\mathrm{\θ}}_2,...{\mathrm{\θ}}_{N_s},$ 为x(t)的估计值。

本发明的有益效果：本发明在传统的MFOCUSS算法基础上，提出了一种基于压缩感知与正则MFOCUSS的超分辨谱估计方法，在天线阵元数较少的情况下，相比传统的MUSIC算法和MFOCUSS算法取得了较高分辨率，同时正则化参数的引入使得算法在解的稀疏度和误差之间有一个合理的权衡，而自适应参数的引入进一步提高了解的精度。

以下将结合附图对本发明做进一步详细说明。

附图说明

图1是本发明的总流程示意图。

图2(a)是三个信源，信噪比为80dB，天线阵元数为32时，标准角度谱图。

图2(b)是三个信源，信噪比为80dB，天线阵元数为32时，传统的MUSIC算法的恢复效果图。

图2(c)是三个信源，信噪比为80dB，天线阵元数为32时，MFOCUSS算法的恢复效果图。

图2(d)是三个信源，信噪比为80dB，天线阵元数为32时，本发明算法的恢复效果图。

图3(a)是三个信源，信噪比为80dB，天线阵元数为16时，标准角度谱图。

图3(b)是三个信源，信噪比为80dB，天线阵元数为16时，传统的MUSIC算法的恢复效果图。

图3(c)是三个信源，信噪比为80dB，天线阵元数为16时，MFOCUSS算法的恢复效果图。

图3(d)是三个信源，信噪比为80dB，天线阵元数为16时，本发明算法的恢复效果图。

图4(a)是三个信源，信噪比为30dB，天线阵元数为32时，标准角度谱图。

图4(b)是三个信源，信噪比为30dB，天线阵元数为32时，传统的MUSIC算法的恢复效果图。

图4(c)是三个信源，信噪比为30dB，天线阵元数为32时，MFOCUSS算法的恢复效果图。

图4(d)是三个信源，信噪比为30dB，天线阵元数为32时，本发明算法的恢复效果图。

图5(a)是三个信源，信噪比为30dB，天线阵元数为16时，标准角度谱图。

图5(b)是三个信源，信噪比为30dB，天线阵元数为16时，传统的MUSIC算法的恢复效果图。

图5(c)是三个信源，信噪比为30dB，天线阵元数为16时，MFOCUSS算法的恢复效果图。

图5(d)是三个信源，信噪比为30dB，天线阵元数为16时，本发明算法的恢复效果图。

图6是不同信噪比时，天线阵元数为16时，MFOCUSS算法和本发明算法的均方误差图。

图7是不同信噪比时，天线阵元数为8时，MFOCUSS算法和本发明算法的均方误差图。

具体实施方式

实施例1：

本发明提出了一种基于压缩感知与正则MFOCUSS的超分辨谱估计方法，如图1所示，本发明的实现过程包含如下步骤：

步骤1.构造接收的阵列信号，具体过程如下：

1a)设定天线阵元数目N,快拍数M；

1b)第j(j＝1,2,…p)个信源(假定入射信源是远场窄带信号)的信号其中u_j(t)和为第j个接收信号的幅度和相位，ω₀为接收信号的频率，则第i个阵元接收到的信号为其中g_ij为第i个阵元对第j个信号的增益，n_i(t)为第i个阵元在t时刻的噪声，τ_ij表示第j个信号到达第i个阵元时相对于参考阵元的延迟；

1c)构造接收信号组成的矩阵Y＝[y₁(k)y₂(k)…y_N(k)]^T,作为N个天线阵元的输出，其中y_i(k)＝[y_i(1×t/M)y_i(2×t/M)…y_i(M×t/M)]，k＝1,2,…M；

步骤2.构造稀疏字典和一个稀疏向量，具体过程如下：

2a)构造N×N_s的角度扫描矩阵A，即稀疏字典 $A=\left[\begin{array}{cccc}a\left({\mathrm{\θ}}_1\right)& a\left({\mathrm{\θ}}_2\right)& \·\·\·& a\left({\mathrm{\θ}}_{N_s}\right)\end{array}\right],$ 其中 $a\left({\mathrm{\θ}}_i\right)={\left[\begin{array}{cccc}1& e^{-{\mathrm{j\ω}}_i}& \·\·\·& e^{-j(N-1){\mathrm{\ω}}_i}\end{array}\right]}^T,$ ω_i＝2πsin(θ_i)/λ,i＝1,2,…N_s，λ为接收信号对应的波长，N_s为要观测的全部角度θ_i的数目；

2b)构造一个N_s×1的稀疏向量 $x\left(t\right)={\left[\begin{array}{cccc}{x}_{{\mathrm{\θ}}_1}\left(t\right)& x_{{\mathrm{\θ}}_2}\left(t\right)& \·\·\·& x_{{\mathrm{\θ}}_{\mathrm{Ns}}}\left(t\right)\end{array}\right]}^T,$ 其中有p个非零的系数x_θ(t)＝s_p(t)在对应于信号源角度的位置，剩下的N_s-p个系数为零，进而构造M个快拍组成的N_s×M的矩阵X＝[x(1)x(2)…x(M)]；

步骤3.初始化迭代参数和正则化参数，具体过程如下：

3a)设定初始解X⁽⁰⁾＝0；

3b)初始化迭代次数k＝0，自适应参数p_k＝1,误差阈值ε＝1e-8，总的迭代限制次数Iter＝800，正则化参数γ＝σ²，σ²为噪声方差；

步骤4.由正则化的MFOCUSS算法求解，当满足误差终止条件时，结束循环，否则，更新权值参数，具体过程如下：

4a)利用式X^(k+1)＝W^(k+1)Q^(k+1)计算第k+1次迭代求解出的目标矩阵，其中

$W^{(k+1)}=\mathrm{diag}\left({\left(c_i^{\left(k\right)}\right)}^{1-p/2}\right),p\∈\left[\mathrm{0,1}\right],$

$c_i^{\left(k\right)}=\left|\right|X_i^{\left(k\right)}\left|\right|={\left({\mathrm{\Σ}}_{i=1}^L{\left|X_i^{\left(k\right)}\left(l\right)\right|}^2\right)}^{1/2},$

$Q^{(k+1)}={\left(A^{(k+1)}\right)}^{+}{(A^{(k+1)}{\left(A^{(k+1)}\right)}^{+}+\mathrm{\γI})}^{-1}Y,\mathrm{\γ}\≥0,$

I为N阶单位矩阵，A^(k+1)＝AW^(k+1)

4b)计算相对误差rRMSE_k＝||X^(k)-X^(k-1)||_F/||X^(k)||_F，||·||_F表示先对矩阵各元素求平方和，再开2次方，若rRMSE_k＜ε或k≥Iter,结束循环，否则k＝k+1，p_k＝min(rRMSE_k,1)，重复4a)和4b)；

步骤5.根据角度谱的计算方法得出最终的DOA估计，具体过程如下：

5a)根据RMSE＝||X^(k)-X||_F计算最终的均方误差；

5b)根据角度谱的计算公式求出角度谱，其中 $\mathrm{\θ}={\mathrm{\θ}}_1,{\mathrm{\θ}}_2,...{\mathrm{\θ}}_{N_s},$ 为x(t)的估计值。

实施例2：

基于压缩感知与正则MFOCUSS的超分辨谱估计方法同实施例1，本发明的效果可以通过以下实验进一步说明。

实验条件：本实验所用数据为仿真的阵列信号，采用软件MATLAB7.0作为仿真工具，计算机配置为IntelCore2/2.33G/2G。

实验采用传统的MUSIC算法，MFOCUSS算法和基于压缩感知与正则MFOCUSS的算法(本发明)进行对比实验；实验中，信号源的个数均为3个(两个相干信源)。

首先在信噪比较高(SNR＝80dB)，噪声方差σ²＝9.99e-009，在天线阵元数目为32的情况下，传统的MUSIC算法，MFOCUSS算法和本发明的对比实验，结果如图2(a)，图2(b)，图2(c)和图2(d)所示。在天线阵元数目为16的情况下比较传统的MUSIC算法，MFOCUSS算法和本发明三种算法，结果如图3(a)，图3(b)，图3(c)和图3(d)所示。结果表明，在较高的信噪比下，天线阵元数足够多时，三种方法都能准确估计出信号源的角度谱，而在天线阵元数减为原天线阵元数的一半时，MUSIC算法失效，而MFOCUSS算法和本发明仍然可以准确估计出信号的角度谱。

然后在信噪比较低(SNR＝30dB)，噪声方差σ²＝0.001，在天线阵元数目为32的情况下比较传统的MUSIC算法，MFOCUSS算法和本发明三种算法，结果如图4(a)，图4(b)，图4(c)和图4(d)所示。在天线阵元数目为16的情况下比较传统的MUSIC算法，MFOCUSS算法和本发明三种算法，结果如图5(a)，图5(b)，图5(c)和图5(d)所示。结果表明，在信噪比较低的情况下，天线阵元数足够多时，三种方法都能估计出相干信号源的角度谱，但由于噪声的影响，MUSIC算法和MFOCUSS算法不够准确，而本发明能够准确地估计出相干信号源的角度谱；天线阵元数较少时，MUSIC算法和MFOCUSS算法已经不能估计出相干信号源的角度谱，而本发明仍然能够准确地估计出相干信号源的角度谱，具有更强的抗噪性能。

另外，如图6所示，在不同的信噪比(5、10、20、30、40、50、60、70、80dB)，天线阵元数目为16时，比较MFOCUSS算法和本发明的均方误差(50次平均值)；如图7所示，在不同的信噪比(5、10、20、30、40、50、60、70、80dB)，天线阵元数目为8时，比较MFOCUSS算法和本发明的均方误差(50次平均值)。结果表明，在信噪比较低的情况下，本发明均方误差较小，抗噪性能强，能够取得更好的效果。

综上所述，在天线单元较少的情况下，本发明比传统MUSIC算法具有更高的角度分辨率，而且在信噪比较低的情况下，本发明比MFOCUSS的DOA估计具有更高的鲁棒性，解的精度也有所提高。

以上例举仅仅是对本发明的举例说明，并不构成对本发明的保护范围的限制，凡是与本发明相同或相似的设计均属于本发明的保护范围之内。

Claims (1)

1.基于压缩感知与正则MFOCUSS的超分辨谱估计方法，其特征在于：包括以下步骤：

步骤1.构造接收的阵列信号，具体过程如下：

1a)设定天线阵元数目N,快拍数M；

1b)第j(j＝1,2,…p)个信源的信号其中u_j(t)和为第j个接收信号的幅度和相位，ω₀为接收信号的频率，则第i个阵元接收到的信号为其中g_ij为第i个阵元对第j个信号的增益，n_i(t)为第i个阵元在t时刻的噪声，τ_ij表示第j个信号到达第i个阵元时相对于参考阵元的延迟；

1c)构造接收信号组成的矩阵Y＝[y₁(k)y₂(k)…y_N(k)]^T,作为N个天线阵元的输出，其中y_i(k)＝[y_i(1×t/M)y_i(2×t/M)…y_i(M×t/M)]，k＝1,2,…M；

步骤2.构造稀疏字典和一个稀疏向量，具体过程如下：

2a)构造N×N_s的角度扫描矩阵A，即稀疏字典 $A=\left[\begin{array}{cccc}a\left({\mathrm{\θ}}_1\right)& a\left({\mathrm{\θ}}_2\right)& \·\·\·& a\left({\mathrm{\θ}}_{N_s}\right)\end{array}\right],$ 其中 $a\left({\mathrm{\θ}}_i\right)={\left[\begin{array}{cccc}1& e^{-{\mathrm{j\ω}}_i}& \·\·\·& e^{-j(N-1){\mathrm{\ω}}_i}\end{array}\right]}^T,$ ω_i＝2πsin(θ_i)λ,i＝1,2,…N_s，λ为接收信号对应的波长，N_s为要观测的全部角度θ_i的数目；

2b)构造一个N_s×1的稀疏向量 $x\left(t\right)={\left[\begin{array}{cccc}{x}_{{\mathrm{\θ}}_1}\left(t\right)& x_{{\mathrm{\θ}}_2}\left(t\right)& \·\·\·& x_{{\mathrm{\θ}}_{\mathrm{Ns}}}\left(t\right)\end{array}\right]}^T,$ 其中有p个非零的系数x_θ(t)＝s_p(t)在对应于信号源角度的位置，剩下的N_s-p个系数为零，进而构造M个快拍组成的N_s×M的矩阵X＝[x(1)x(2)…x(M)]；

步骤3.初始化迭代参数和正则化参数，具体过程如下：

3a)设定初始解X⁽⁰⁾＝0；

3b)初始化迭代次数k＝0，自适应参数p_k＝1,误差阈值ε＝1e-8，总的迭代限制次数Iter＝800，正则化参数γ＝σ²，σ²为噪声方差；

步骤4.由正则化的MFOCUSS算法求解，当满足误差终止条件时，结束循环，否则，更新权值参数，具体过程如下：

4a)利用式X^(k+1)＝W^(k+1)Q^(k+1)计算第k+1次迭代求解出的目标矩阵，其中

$W^{(k+1)}=\mathrm{diag}\left({\left(c_i^{\left(k\right)}\right)}^{1-p/2}\right),p\∈\left[\mathrm{0,1}\right],$

$c_i^{\left(k\right)}=\left|\right|X_i^{\left(k\right)}\left|\right|={\left({\mathrm{\Σ}}_{i=1}^L{\left|X_i^{\left(k\right)}\left(l\right)\right|}^2\right)}^{1/2},$

$Q^{(k+1)}={\left(A^{(k+1)}\right)}^{+}{(A^{(k+1)}{\left(A^{(k+1)}\right)}^{+}+\mathrm{\γI})}^{-1}Y,\mathrm{\γ}\≥0,$

I为N阶单位矩阵，A^(k+1)＝AW^(k+1)

4b)计算相对误差rRMSE_k＝||X^(k)-X^(k-1)||_F/||X^(k)||_F，||·||_F表示先对矩阵各元素求平方和，再开2次方，若rRMSE_k＜ε或k≥Iter,结束循环，否则k＝k+1，p_k＝min(rRMSE_k,1)，重复4a)和4b)；

步骤5.根据角度谱的计算方法得出最终的DOA估计，具体过程如下：

5a)根据RMSE＝||X^(k)-X||_F计算最终的均方误差；

5b)根据角度谱的计算公式求出角度谱，其中 $\mathrm{\θ}={\mathrm{\θ}}_1,{\mathrm{\θ}}_2,...{\mathrm{\θ}}_{N_s}$ ，为x(t)的估计值。