CN103901396B - 相干信号源亚分辨率超分辨到达角估计方法 - Google Patents

Abstract

本发明属于相干信号源到达角估计技术领域，特别涉及相干信号源亚分辨率超分辨到达角估计方法。该相干信号源亚分辨率超分辨到达角估计方法包括以下步骤：利用雷达的接收阵列接收相干信号源的回波信号Y，所述雷达的接收阵列为具有M个阵元的均匀线阵；M>K，K为回波信号Y对应的相干信号源的个数；建立关于回波信号Y的信号模型，所述关于回波信号Y的信号模型包括两个待求解的矩阵：方向矩阵D和信号稀疏矩阵S；由关于回波信号Y的信号模型，建立关于方向矩阵D和信号稀疏矩阵S的代价函数；求解步骤S4中的代价函数，得出方向矩阵D和信号稀疏矩阵S；根据方向矩阵D中主对角线上的非零元素、以及每个相干信号源的到达角的初始估计值，得出各个相干信号源的到达角。

Description

相干信号源亚分辨率超分辨到达角估计方法

技术领域

本发明属于相干信号源到达角估计技术领域，特别涉及相干信号源亚分辨率超分辨到达角估计方法，能够在低信噪比低快拍相干源的到达角估计中提高小于分辨率的到达角的估计精度。

背景技术

空间信号的到达角估计与时间信号的频率估计十分相似。在理论上，他们均可表述为基本的非线性参数问题。因此，时域非线性信号处理技术在空域处理上的推广曾成为热点问题，许多时域非线性谱估计方法推广成为空域谱估计方法，于是便产生了所谓的高分辨谱估计技术，即超分辨到达角估计技术，它能突破并进一步改善一个波束宽度内的空间不同来向的分辨能力。然而在实际工程运用中，例如在米波雷达低仰角估计时，就遇到在低快拍，低信噪比和相干信源的情形下对仰角估计的问题。

在实际中提出了很多关于求解到达角的方法，如多重信号分类（MUSIC），旋转不变子空间（ESPRIT）等方法，但对于相干信号源的到达角估计，相干信号源信号子空间与噪声子空间相互渗透，故不能对相干信号源进行有效分辨或测向。当前对相干信号源的处理方法，如空间平滑技术，其对相干信号源到达角估计有理想的效果，但它通过牺牲有效阵元数来换取的，对阵列孔径有一定的损失，且在低信噪比时此方法的性能较差；Toeplitz技术，具有估计的到达角偏差较大的缺点。因此，相干信号源的估计问题始终是一个十分棘手的问题，它也始终是空间谱估计中研究的热点问题。

对于低快拍的问题，通常的到达角估计是基于理想情况下进行的，但在工程应用场合中各种各样的误差不可避免，因此，当实际模型与假设模型不符合时，通常会出现有限数据长度快拍的问题，许多基于理想模型基础上的到达角方法的性能将严重下降。有限数据长度快拍对模型的影响在于，当快拍数小时，信号与噪声没有足够的时间解相关，噪声协方差矩阵也没有收敛，只有当快拍数大或信噪比高时，有限数据长度的影响才可以忽略。如何消除这些误差的影响，使估计到达角方法在应用中具有较强的稳健性是到达角方法实用化的一个关键环节。

亚分辨率问题，这个问题描述为当使用搜索算法来求解到达角时，如果一味地对细化搜索间隔，就会造成运量大的问题；而当把搜索间隔变粗时，又会造成到达角估计的不精确，亚分辨率问题就是当实际的到达角在两次相邻搜索样点之间时，再不增加运算量的前提下，求出真实的到达角，就是亚分辨率问题。

一类基于稀疏信号重构的方法被提出，这类方法快速的被应用到到达角估计的问题上，在低信噪比和快拍数较少的环境下提供了一种有效的解决方法。该算法在低快拍下更具有优势，该算法随着快拍数的增多可以提高估计效果，此外建立模型时并未要求信号源是相互独立的，因此对空域信号的相关性不敏感，可以直接用于相干信号的到达角估计，且具有更高的角度分辨力。但当实际的到达角在两次相邻搜索样点之间时，此时需要解决到达角的亚分辨率的问题。这样一个问题，有学者提出了一种基于稀疏贝叶斯表示（SBI）的亚分辨率到达角估计的方法，该方法能够有效、准确的估计当到达角在两次相邻搜索样点之间的到达角，但这个方法需较高的信噪比，且其运算量也较大。

发明内容

本发明的目的在于针对稀疏贝叶斯表示（SBI）的亚分辨率到达角估计技术的不足，提出相干信号源亚分辨率超分辨到达角估计方法，在低信噪比、低快拍情况下，就能够很好的解决亚分辨率的相干源到达角估计问题。

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

相干信号源亚分辨率超分辨到达角估计方法包括以下步骤：

S1：利用雷达的接收阵列接收相干信号源的回波信号Y，所述雷达的接收阵列为具有M个阵元的均匀线阵；M>K，K为回波信号Y对应的相干信号源的个数；

S2：建立关于回波信号Y的信号模型，所述关于回波信号Y的信号模型包括两个待求解的矩阵：方向矩阵D和信号稀疏矩阵S；

S3：由关于回波信号Y的信号模型，建立关于方向矩阵D和信号稀疏矩阵S的代价函数；

S4：求解步骤S4中的代价函数，得出方向矩阵D和信号稀疏矩阵S；

S5：根据方向矩阵D中主对角线上的非零元素，得出各个相干信号源的到达角。

本发明的特点和进一步改进在于：

所述步骤S2具体包括以下步骤：

设落在雷达搜索样点上的采样到达角共有N个，N>M；将落在雷达搜索样点上的第n个采样到达角表示为则将落在雷达搜索样点上的所有N个采样到达角表示为集合 n取1至N；

当雷达的接收阵列的阵元个数M为偶数时，其中，T表示矩阵或向量的转置；当雷达的接收阵列的阵元个数M为奇数时，然后得出常量矩阵

定义其中，⊙表示向量的点积；当雷达的接收阵列的阵元个数M为偶数时，

$p={[-\frac{M-1}{2},-\frac{M-3}{2},...,-\frac{1}{2},\frac{1}{2},...\frac{M-3}{2},\frac{M-1}{2}]}^T;$

当雷达的接收阵列的阵元个数M为奇数时，

$p={[-\frac{M-1}{2},-\frac{M-3}{2},...,-\mathrm{1,0,1},...\frac{M-3}{2},\frac{M-1}{2}]}^T;$

然后得出系数矩阵

然后，建立如下关于回波信号Y的信号模型：

其中，E为设定的噪声矩阵；方向矩阵D和信号稀疏矩阵S为两个待求解的矩阵，方向矩阵D为N阶对角矩阵。

所述步骤S3具体包括以下步骤：

将回波信号Y中K个最大的奇异值对应的左特征向量按列排列，形成左信号矩阵U_s；设立K阶对角矩阵Λ_s，回波信号Y的K个最大的奇异值按照从大到小的顺序依次排列在对角矩阵Λ_s的主对角线上；令

${\mathrm{\σ}}_e^2=\frac{1}{M-K}\underset{g=K+1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\λ}}_g$

其中，λ_g表示回波信号Y的第g个奇异值，g取K+1至M；

令X=U_sW^1/2，其中，I_K为K阶单位矩阵，W^1/2表示矩阵W的1/2次方，上标-1表示求矩阵的逆；

然后，建立如下关于方向矩阵D和信号稀疏矩阵S的代价函数：

其中，||·||_F表示求F范数，λ为设定的归一化参数，λ>0，||S||_∞,0为：

${\left|\right|S\left|\right|}_{\∞,0}={\left|\right|\underset{q}{\mathrm{\Σ}}\left|S_{\mathrm{pq}}\right|\left|\right|}_0$

其中，S_pq表示信号稀疏矩阵S中第p行第q列的元素；|·|表示取绝对值，||·||₀表示求0范数。

所述步骤S4具体包括以下步骤：

S41：将方向矩阵D的初始取值设置为D(0)；设置迭代参数i，i=1,2,3...；当i=1时，执行步骤S42；

S42：令选取Φ_i中不全为零的列，利用Φ_i中选取的各列的编号组成集合Ω_i；然后通过求解得出Ω_i；

所述步骤S42中求解得出Ω_i的过程包括以下步骤：

S421：将集合Ω_i的初始取值设为Ω_i,0，令R_i,0=X，设置迭代参数j，j=1,2,3...；当j=1时，执行步骤S422；

S422：根据下式求解β_i：

${\mathrm{\β}}_i=\underset{p\left(i\right)}{\arg \min }\left[\underset{q\left(i\right)}{\mathrm{\Σ}}\right|{({\mathrm{\Φ}}_{i,{\mathrm{\Ω}}_{i,j-1}}^H{R}_{i,j-1}-{\left|\right|{\mathrm{\Φ}}_{i,{\mathrm{\Ω}}_{i,j-1}}^H{R}_{i,j-1}\left|\right|}_{\∞})}_{p\left(i\right)q\left(i\right)}\left|\right]$

其中，p(i)表示第p(i)行，q(i)表示第q(i)列，(·)_p(i)q(i)表示矩阵的第p(i)行第q(i)列的元素，|·|表示取绝对值；||·||_∞表示求无穷范数，矩阵的求取过程为：在Φ_i中，筛选出与Ω_i,j-1B中元素对应的列向量，然后将Φ_i中筛选出的列向量按照列顺序组成矩阵H表示矩阵的共轭转置；

S423：通过以下迭代关系式得出Ω_i,j和R_i,j：

Ω_i,j=Ω_i,j-1∪β_i， $R_{i,j}=X-{\mathrm{\Φ}}_{{\mathrm{\Ω}}_{i,j}}{\left({\mathrm{\Φ}}_{{\mathrm{\Ω}}_{i,j}}^H{\mathrm{\Φ}}_{{\mathrm{\Ω}}_{i,j}}\right)}^{-1}{\mathrm{\Φ}}_{{\mathrm{\Ω}}_{i,j}}^{H}X$

其中，∪表示集合并运算，矩阵的求取过程为：在Φ_i中，筛选出与Ω_i,j中元素对应的列向量，然后将Φ_i中筛选出的列向量按照列顺序组成矩阵 ${\mathrm{\Φ}}_{i,{\mathrm{\Ω}}_{i,j}};$

S424：令j的值自增1，判断j是否等于K+1，K为回波信号Y对应的相干信号源的个数；如果j=K+1，则求解出Ω_i：Ω_i=Ω_i,K；然后，执行步骤S43；否则，则返回至步骤S422；

S43：将矩阵定义如下：在矩阵S(i)中筛选出不全为零的列，将矩阵S(i)中筛选出的列按照列顺序组合成矩阵 的计算公式为：

${S\left(i\right)}_{{\mathrm{\Ω}}_i}={\mathrm{\Φ}}_{i,{\mathrm{\Ω}}_i}{\left({\mathrm{\Φ}}_{i,{\mathrm{\Ω}}_i}^H{\mathrm{\Φ}}_{i,{\mathrm{\Ω}}_i}\right)}^{-1}{\mathrm{\Φ}}_{i,{\mathrm{\Ω}}_i}^{H}X$

然后，将矩阵的各个元素还原至S(i)中，在矩阵S(i)的其余位置补零，从而求解出矩阵S(i)；

S44：根据以下公式得出方向矩阵D的第i次迭代取值D(i)：

s.t.||S||_∞，0=K

其中，||·||_F表示求F范数；

S45：对角矩阵D(i-1)中主对角线上的元素按顺序排列形成的列向量记为d(i-1)，对角矩阵D(i)中主对角线上的元素按顺序排列形成的列向量记为d(i)；判断||d(i)-d(i-1)||₂是否小于ε，其中，||·||₂表示求2-范数，ε为设定的迭代终止值且ε>0；如果||d(i)-d(i-1)||₂小于ε，则得出方向矩阵D和信号稀疏矩阵S，D=D(i)，S=S(i)；否则，令i的值自增1，返回至步骤S42。

在步骤S44中，根据以下公式得出方向矩阵D的第i次迭代取值D(i)：

D(i)=diag(d(i))

其中，||·||₂表示求2-范数，vec(·)表示矩阵的列向量化处理，T表示向量或矩阵的转置，表示K-R积，在得出的D(i)中，d(i)中的元素按照顺序排列在D(i)的主对角线上。

在步骤S44中，根据以下公式得出方向矩阵D的第i次迭代取值D(i)：

D(i)=diag(d(i))

其中，||·||₁表示求1-范数，vec(·)表示矩阵的列向量化处理，T表示向量或矩阵的转置，表示K-R积， $\stackrel{\‾}{W}=I_K\⊗W,$ 表示Kronecker积，I_K为K阶单位矩阵，在得出的D(i)中，d(i)中的元素按照顺序排列在D(i)的主对角线上。

所述步骤S5具体包括以下步骤：方向矩阵D中主对角线上的非零元素共有K个，将方向矩阵D中主对角线上的非零元素筛选出来，方向矩阵D中主对角线上的第k个非零元素表示为D_k，k取1至K；D_k是方向矩阵D中主对角线上的第k'个元素，1≤k'≤N；则根据以下公式求解得出第k个相干信号源的到达角θ_k：

本发明的有益效果为：1）本发明使用子空间拟合和加权1范数逼近的方法，可以有效的估计不在搜索样点上的相干信源到达角估计问题。2）本发明估计到达角的方法可以克服现有方法（例如基于稀疏贝叶斯表示的亚分辨率到达角估计的方法）运算量大、信噪比要求较高的缺点。

附图说明

图1为本发明的相干信号源亚分辨率超分辨到达角估计方法的流程图；

图2为雷达阵地的举例示意图；

图3为仿真实验1中采用本发明得出的两个相干信号源的到达角的估计示意图；

图4为仿真实验2中几种方法得出的两个相干信号源的到达角的估计的均方根误差的对比示意图；

图5为仿真实验3中几种方法得出的阵元个数与运算时间的关系曲线的对比示意图。

具体实施方式

下面结合附图对本发明作进一步说明：

参照图1，为本发明的相干信号源亚分辨率超分辨到达角估计方法的流程图。该相干信号源亚分辨率超分辨到达角估计方法包括以下步骤：

S1：利用雷达的接收阵列接收相干信号源的回波信号Y，所述雷达的接收阵列为具有M个阵元的均匀线阵；M>K，K为回波信号Y对应的相干信号源的个数。

S2：建立关于回波信号Y的信号模型，所述关于回波信号Y的信号模型包括两个待求解的矩阵：方向矩阵D和信号稀疏矩阵S。具体说明如下：

步骤S2具体包括以下步骤：

设落在雷达搜索样点上的采样到达角共有N个，N>M，例如雷达搜索样点上的采样到达角为：从-90度至90度之间的整数角度，此时，N=181。

此时，无论每个相干信号源的到达角的真实值位于雷达搜索样点上，还是没有落在雷达搜索样点上，均可以按照以下步骤得出每个相干信号源的到达角。

将落在雷达搜索样点上的第n个采样到达角表示为（例如为-16度），则将落在雷达搜索样点上的所有N个采样到达角表示为集合 n取1至N。

在本发明实施例中，将第k个相干信号源的到达角的真实值表示为θ_k，将第k个相干信号源的到达角的导向矢量表示为a(θ_k)。将a(θ_k)按照如下公式进行一阶泰勒展开：

$a\left({\mathrm{\θ}}_k\right)\≈a\left({\widetilde{\mathrm{\θ}}}_k\right)+b\left({\widetilde{\mathrm{\θ}}}_k\right)({\mathrm{\θ}}_k-{\widetilde{\mathrm{\θ}}}_k)$

其中，表示第k个相干信号源的到达角的估计值的导向矢量，表示对应的一阶泰勒系数。

当雷达的接收阵列的阵元个数M为偶数时，其中，T表示矩阵或向量的转置；当雷达的接收阵列的阵元个数M为奇数时，无论雷达的接收阵列的阵元个数M是偶数还是奇数，列向量表示的导向矢量，列向量的元素个数均为M。然后得出常量矩阵

则为M×N维矩阵。

定义其中，⊙表示向量的点积（Hadamard积）；当雷达的接收阵列的阵元个数M为偶数时，

$p={[-\frac{M-1}{2},-\frac{M-3}{2},...,-\frac{1}{2},\frac{1}{2},...\frac{M-3}{2},\frac{M-1}{2}]}^T;$

当雷达的接收阵列的阵元个数M为奇数时，

$p={[-\frac{M-1}{2},-\frac{M-3}{2},...,-\mathrm{1,0,1},...\frac{M-3}{2},\frac{M-1}{2}]}^T;$

无论雷达的接收阵列的阵元个数M是偶数还是奇数，列向量p的元素个数均为M。然后得出系数矩阵

则为M×N维矩阵。

然后，根据a(θ_k)的一阶泰勒展开式，可以建立如下关于回波信号Y的信号模型：

其中，E为设定的噪声矩阵，E为M×L维矩阵；方向矩阵D为N阶对角矩阵，信号稀疏矩阵S为N×L维矩阵；方向矩阵D和信号稀疏矩阵S为两个待求解的矩阵。L为雷达接收阵列中每个阵元的采样数，则回波信号Y为M×L维矩阵。

S3：由关于回波信号Y的信号模型，建立关于方向矩阵D和信号稀疏矩阵S的代价函数。具体说明如下：

步骤S3具体包括以下步骤：

将回波信号Y中K个最大的奇异值对应的左特征向量按列排列，形成左信号矩阵U_s；设立K阶对角矩阵Λ_s，回波信号Y的K个最大的奇异值按照从大到小的顺序依次排列在对角矩阵Λ_s的主对角线上；令

${\mathrm{\σ}}_e^2=\frac{1}{M-K}\underset{g=K+1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\λ}}_g$

其中，λ_g表示回波信号Y的第g个奇异值，g取K+1至M。

令X=U_sW^1/2，其中，I_K为K阶单位矩阵，W表示渐近优化加权量，W^1/2表示矩阵W的1/2次方，上标-1表示求矩阵的逆；

然后，建立如下关于方向矩阵D和信号稀疏矩阵S的代价函数：

其中，的含义为：以D和S为变量求解括号内的矩阵的最小值。||·||_F表示求F范数，λ为设定的归一化参数，λ>0，||S||_∞,0为：

${\left|\right|S\left|\right|}_{\∞,0}={\left|\right|\underset{q}{\mathrm{\Σ}}\left|S_{\mathrm{pq}}\right|\left|\right|}_0$

其中，S_pq表示信号稀疏矩阵S中第p行第q列的元素；|·|表示取绝对值，||·||₀表示求0范数。

S4：求解步骤S4中的代价函数，得出方向矩阵D和信号稀疏矩阵S。具体说明如下：

对于步骤S3的代价函数求解问题，从表达式可以看出，是属于二元优化问题，优化变量为S和D，如果直接求解，其求解过程相当复杂，且运算量大。对此本发明采用迭代计算S和D，其策略为先固定D的前一次迭代值，再求解D的本次迭代值和S的本次迭代值。

步骤S4具体包括以下步骤：

S41：将方向矩阵D的初始取值设置为D(0)；设置迭代参数i，i=1,2,3...；当i=1时，执行步骤S42。

S42：令选取Φ_i中不全为零的列（该列元素不全是零或该列元素全不为零），利用Φ_i中选取的各列的编号组成集合Ω_i；然后通过求解得出Ω_i。Ω_i表示Φ_i的支撑集。

所述步骤S42中求解得出Ω_i的过程包括以下步骤：

S421：将集合Ω_i的初始取值设为Ω_i,0，令R_i,0=X，设置迭代参数j，j=1,2,3...；当j=1时，执行步骤S422。

S422：根据下式求解β_i：

${\mathrm{\β}}_i=\underset{p\left(i\right)}{\arg \min }\left[\underset{q\left(i\right)}{\mathrm{\Σ}}\right|{({\mathrm{\Φ}}_{i,{\mathrm{\Ω}}_{i,j-1}}^H{R}_{i,j-1}-{\left|\right|{\mathrm{\Φ}}_{i,{\mathrm{\Ω}}_{i,j-1}}^H{R}_{i,j-1}\left|\right|}_{\∞})}_{p\left(i\right)q\left(i\right)}\left|\right]$

其中，p(i)表示第p(i)行，q(i)表示第q(i)列，(·)_p(i)q(i)表示矩阵的第p(i)行第q(i)列的元素，|·|表示取绝对值；||·||_∞表示求无穷范数，矩阵的求取过程为：在Φ_i中，筛选出与Ω_i,j-1中元素对应的列向量，然后将Φ_i中筛选出的列向量按照列顺序组成矩阵H表示矩阵的共轭转置。

S423：通过以下迭代关系式得出Ω_i,j和R_i,j：

Ω_i,j=Ω_i,j-1∪β_i， $R_{i,j}=X-{\mathrm{\Φ}}_{{\mathrm{\Ω}}_{i,j}}{\left({\mathrm{\Φ}}_{{\mathrm{\Ω}}_{i,j}}^H{\mathrm{\Φ}}_{{\mathrm{\Ω}}_{i,j}}\right)}^{-1}{\mathrm{\Φ}}_{{\mathrm{\Ω}}_{i,j}}^{H}X$

其中，∪表示集合并运算，矩阵的求取过程为：在Φ_i中，筛选出与Ω_i,j中元素对应的列向量，然后将Φ_i中筛选出的列向量按照列顺序组成矩阵 ${\mathrm{\Φ}}_{i,{\mathrm{\Ω}}_{i,j}}.$

S424：令j的值自增1，判断j是否等于K+1，K为回波信号Y对应的相干信号源的个数；如果j=K+1，则求解出Ω_i：Ω_i=Ω_i,K；然后，执行步骤S43；否则，则返回至步骤S422。

S43：将矩阵定义如下：在矩阵S(i)中筛选出不全为零的列（该列元素不全是零或该列元素全不为零），将矩阵S(i)中筛选出的列按照列顺序组合成矩阵 的计算公式为：

${S\left(i\right)}_{{\mathrm{\Ω}}_i}={\mathrm{\Φ}}_{i,{\mathrm{\Ω}}_i}{\left({\mathrm{\Φ}}_{i,{\mathrm{\Ω}}_i}^H{\mathrm{\Φ}}_{i,{\mathrm{\Ω}}_i}\right)}^{-1}{\mathrm{\Φ}}_{i,{\mathrm{\Ω}}_i}^{H}X$

然后，将矩阵的各个元素还原至S(i)中，在矩阵S(i)的其余位置补零，从而求解出矩阵S(i)；

S44：得出方向矩阵D的第i次迭代取值D(i)可以有如下几种方法：

D(i)的第一种求解方法：

根据以下公式得出方向矩阵D的第i次迭代取值D(i)：

s.t.||S||_∞，0=K

其中，||·||_F表示求F范数。

D(i)的第二种求解方法：

在D(i)的第一种求解方法中，D(i)的闭式解为：

$D\left(i\right)=B^{+}\left(\widetilde{\mathrm{\θ}}\right)(X-A\left(\widetilde{\mathrm{\θ}}\right))S^{+}\left(i\right)$

其中，[·]⁺表示Moore-Penrose逆运算符。但上述闭式解在进行Moore-Penrose逆运算时，如果对应的条件数较大，则会造成较大的相对误差。因此，为了克服该缺点，在D(i)的第二种求解方法中，根据以下公式得出方向矩阵D的第i次迭代取值D(i)：

D(i)=diag(d(i))

其中，||·||₂表示求2-范数，vec(·)表示矩阵的列向量化处理，T表示向量或矩阵的转置，表示K-R积，在得出的D(i)中，d(i)中的元素按照顺序排列在D(i)的主对角线上。

D(i)的第三种求解方法：

当雷达的接收阵列中每个阵元的采样数较少（L<30），且设定的噪声矩阵E不是白噪声矩阵时，使用D(i)的第二种求解方法得出的结果就容易出现不稳定的情形，因此，根据以下公式得出方向矩阵D的第i次迭代取值D(i)：

D(i)=diag(d(i))

其中，||·||₁表示求1-范数，vec(·)表示矩阵的列向量化处理，T表示向量或矩阵的转置，表示K-R积， $\stackrel{\&OverBar;}{W}=I_K\&CircleTimes;W,$ 表示Kronecker内积，I_K为K阶单位矩阵，在得出的D(i)中，d(i)中的元素按照顺序排列在D(i)的主对角线上。用I_M表示M阶单位矩阵，当W≠I_M时，D(i)的第三种求解方法称为子空间拟合和加权1范数法（简称L1-WSSF）；当W=I_M时，D(i)的第三种求解方法称为子空间拟合和1范数（简称L1-SSF）。

S45：对角矩阵D(i-1)中主对角线上的元素按顺序排列形成的列向量记为d(i-1)，对角矩阵D(i)中主对角线上的元素按顺序排列形成的列向量记为d(i)；判断||d(i)-d(i-1)||₂是否小于ε，其中，||·||₂表示求2-范数，ε为设定的迭代终止值且ε>0；如果||d(i)-d(i-1)||₂小于ε，则得出方向矩阵D和信号稀疏矩阵S，D=D(i)，S=S(i)；否则，令i的值自增1，返回至步骤S42。

S5：根据方向矩阵D中主对角线上的非零元素、以及每个相干信号源的到达角的初始估计值，得出各个相干信号源的到达角。具体说明如下：

步骤S5具体包括以下步骤：方向矩阵D中主对角线上的非零元素共有K个，将方向矩阵D中主对角线上的非零元素筛选出来，方向矩阵D中主对角线上的第k个非零元素表示为D_k，k取1至K；D_k是方向矩阵D中主对角线上的第k'个元素，1≤k'≤N；则根据以下公式求解得出第k个相干信号源的到达角θ_k：

本发明的效果可以通过以下仿真实验进一步说明。

仿真条件：

参照图2，为雷达阵地的举例示意图，在图2中，目标1至目标4为举例说明的4个目标，1至M表示雷达的接收阵列中的M个阵元。在仿真环境中，采用8根均匀线性阵列架设阵列天线。设有两个相干信号源目标回波信号，其相干系数为0.96，到达均匀线性阵列的到达角（指两个相干信号源的到达角）分别为θ₁=-10.3°和θ₁=-15.6°，雷达的接收阵列中每个阵元的采样数为11次，对环境中[-90°,90°]的角度按照间隔1°进行均匀划分，得出雷达搜索样点上的181个采样到达角，则θ₁和θ₂均不落在雷达搜索样点上。

仿真实验1：采用本发明对不落在雷达搜索样点上的两个相干信号源的到达角进行估计。参照图3，为仿真实验1中采用本发明得出的两个相干信号源的到达角的估计示意图。

仿真实验2：分别采用几种方法对两个相干信号源的到达角的估计进行蒙特卡洛实验。几种方法包括：本发明、对比的方法空间平滑MUSIC方法（简称SS-MUSIC）、空间平滑求根MUSIC（简称SS root-MUSIC）、稀疏贝叶斯表示法（简称SBI）、子空间拟合和加权L1范数法（简称L1-WSSF）、以及子空间拟合和L1范数（简称L1-SSF）。在仿真实验2中，蒙特卡洛实验的次数为200次，则根据下式得出两个相干信号源的到达角的估计的均方根误差RMSE：

$\mathrm{RMSE}=\sqrt{E{({\widehat{\mathrm{\&theta;}}}_1-{\mathrm{\&theta;}}_1)}^2+E{({\widehat{\mathrm{\&theta;}}}_2-{\mathrm{\&theta;}}_2)}^2}$

其中，表示θ₁的估计值，表示θ₂的估计值，E(·)表示求期望。参照图4，为仿真实验2中几种方法得出的两个相干信号源的到达角的估计的均方根误差的对比示意图。图4中，横轴表示信噪比，纵轴表示10log₁₀(RMSE)。

仿真实验3：分别采用几种方法对两个相干信号源的到达角进行估计。这几种方法包括：本发明、空间平滑求根MUSIC（简称SS root-MUSIC）、稀疏贝叶斯表示法（简称SBI）、以及子空间拟合和加权L1范数法（简称L1-WSSF）。当雷达的接收阵列的阵元数发生变化时，得出对应的运算时间。参照图5，为仿真实验3中几种方法得出的阵元个数与运算时间的关系曲线的对比示意图。图5中，横轴表示雷达的接收阵列中的阵元个数，纵轴表示几种方法得出两个相干信号源的到达角的运算时间。

仿真分析：

对于仿真实验1，从图3中可以看出，本发明可以很好的估计相干源在低采样数、低信噪比下的到达角。

对于仿真实验2，从图4中可以看出，在信噪比低于1dB时，相比其他方法，本发明在对于不在雷达搜索样点上的的到达角估计具有更低的均方根误差。子空间拟合和加权L1范数法（L1-WSSF）比子空间拟合和L1范数（L1-SSF）稍微好一些。

对于仿真实验3，从图5可以看出，本发明比稀疏贝叶斯表示法的计算时间少很多，也就是说比稀疏贝叶斯表示法具有更高的运算效率。

显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。

Claims (6)

1.相干信号源亚分辨率超分辨到达角估计方法，其特征在于，包括以下步骤：

S1：利用雷达的接收阵列接收相干信号源的回波信号Y，所述雷达的接收阵列为具有M个阵元的均匀线阵；M＞K，K为回波信号Y对应的相干信号源的个数；

S2：建立关于回波信号Y的信号模型，所述关于回波信号Y的信号模型包括两个待求解的矩阵：方向矩阵D和信号稀疏矩阵S；

步骤S2具体包括以下步骤：

设落在雷达搜索样点上的采样到达角共有N个，N＞M；将落在雷达搜索样点上的第n个采样到达角表示为则将落在雷达搜索样点上的所有N个采样到达角表示为集合 n取1至N；

当雷达的接收阵列的阵元个数M为偶数时，

其中，T表示矩阵或向量的转置；当雷达的接收阵列的阵元个数M为奇数时，

然后得出常量矩阵

定义其中，⊙表示向量的点积；当雷达的接收阵列的阵元个数M为偶数时，

$p={\&lsqb;-\frac{M-1}{2},-\frac{M-3}{2},...,-\frac{1}{2},\frac{1}{2},...\frac{M-3}{2},\frac{M-1}{2}\&rsqb;}^T;$

当雷达的接收阵列的阵元个数M为奇数时，

$p={\&lsqb;-\frac{M-1}{2},-\frac{M-3}{2},...,-1,0,1,...\frac{M-3}{2},\frac{M-1}{2}\&rsqb;}^T;$

然后得出系数矩阵

然后，建立如下关于回波信号Y的信号模型：

其中，E为设定的噪声矩阵；方向矩阵D和信号稀疏矩阵S为两个待求解的矩阵，方向矩阵D为N阶对角矩阵；

S3：由关于回波信号Y的信号模型，建立关于方向矩阵D和信号稀疏矩阵S的代价函数；

S4：求解步骤S4中的代价函数，得出方向矩阵D和信号稀疏矩阵S；

S5：根据方向矩阵D中主对角线上的非零元素，得出各个相干信号源的到达角。

2.如权利要求1所述的相干信号源亚分辨率超分辨到达角估计方法，其特征在于，所述步骤S3具体包括以下步骤：

将回波信号Y中K个最大的奇异值对应的左特征向量按列排列，形成左信号矩阵U_s；设立K阶对角矩阵Λ_s，回波信号Y的K个最大的奇异值按照从大到小的顺序依次排列在对角矩阵Λ_s的主对角线上；令

${\mathrm{\&sigma;}}_e^2=\frac{1}{M-K}\underset{g=K+1}{\overset{M}{\&Sigma;}}{\mathrm{\&lambda;}}_g$

其中，λ_g表示回波信号Y的第g个奇异值，g取K+1至M；

令X＝U_sW^1/2，其中，I_K为K阶单位矩阵，W^1/2表示矩阵W的1/2次方，上标-1表示求矩阵的逆；

然后，建立如下关于方向矩阵D和信号稀疏矩阵S的代价函数：

其中，||·||_F表示求F范数，λ为设定的归一化参数，λ＞0，||S||_∞，0为：

${\left|\right|S\left|\right|}_{\&infin;,0}={\left|\right|\underset{q}{\mathrm{\&Sigma;}}\left|S_{\mathrm{pq}}\right|\left|\right|}_0$

其中，S_pq表示信号稀疏矩阵S中第p行第q列的元素；|·|表示取绝对值，||·||₀表示求0范数。

3.如权利要求2所述的相干信号源亚分辨率超分辨到达角估计方法，其特征在于，所述步骤S4具体包括以下步骤：

S41：将方向矩阵D的初始取值设置为D(0)；设置迭代参数i，i＝1，2，3...；当i＝1时，执行步骤S42；

S42：令选取Φ_i中不全为零的列，利用Φ_i中选取的各列的编号组成集合Ω_i；然后通过求解得出Ω_i；

所述步骤S42中求解得出Ω_i的过程包括以下步骤：

S421：将集合Ω_i的初始取值设为Ω_i，0，令R_i，0＝X，设置迭代参数j，j＝1，2，3...；当j＝1时，执行步骤S422；

S422：根据下式求解β_i：

${\mathrm{\&beta;}}_i=\underset{p\left(i\right)}{\arg \min }\&lsqb;\underset{q\left(i\right)}{\&Sigma;}\left|{({\mathrm{\&Phi;}}_{i,{\mathrm{\&Omega;}}_{i,j-1}}^H{R}_{i,j-1}-|\left|{\mathrm{\&Phi;}}_{i,{\mathrm{\&Omega;}}_{i,j-1}}^H{R}_{i,j-1}\right|{|}_{\mathrm{\&infin;}})}_{p\left(i\right)q\left(i\right)}\right|\&rsqb;$

其中，p(i)表示第p(i)行，q(i)表示第q(i)列，(·)_p(i)q(i)表示矩阵的第p(i)行第q(i)列的元素，|·|表示取绝对值；||·||_∞表示求无穷范数，矩阵的求取过程为：在Φ_i中，筛选出与Ω_i，j-1中元素对应的列向量，然后将Φ_i中筛选出的列向量按照列顺序组成矩阵H表示矩阵的共轭转置；

S423：通过以下迭代关系式得出Ω_i，j和R_i，j：

${\mathrm{\&Omega;}}_{i,j}={\mathrm{\&Omega;}}_{i,j-1}\&cup;{\mathrm{\&beta;}}_i,R_{i,j}=X-{\mathrm{\&Phi;}}_{{\mathrm{\&Omega;}}_{i,j}}{\left({\mathrm{\&Phi;}}_{Q_{i,j}}^H{\mathrm{\&Phi;}}_{{\mathrm{\&Omega;}}_{i,j}}\right)}^{-1}{\mathrm{\&Phi;}}_{{\mathrm{\&Omega;}}_{i,j}}^{H}X$

其中，∪表示集合并运算，矩阵的求取过程为：在Φ_i中，筛选出与Ω_i，j中元素对应的列向量，然后将Φ_i中筛选出的列向量按照列顺序组成矩阵上标-1表示求矩阵的逆；

S424：令j的值自增1，判断j是否等于K+1，K为回波信号Y对应的相干信号源的个数；如果j＝K+1，则求解出Ω_i：Ω_i＝Ω_i，K；然后，执行步骤S43；否则，则返回至步骤S422；

S43：将矩阵定义如下：在矩阵S(i)中筛选出不全为零的列，将矩阵S(i)中筛选出的列按照列顺序组合成矩阵 的计算公式为：

$S{\left(i\right)}_{{\mathrm{\&Omega;}}_i}={\mathrm{\&Phi;}}_{i,{\mathrm{\&Omega;}}_i}{\left({\mathrm{\&Phi;}}_{i,{\mathrm{\&Omega;}}_i}^H{\mathrm{\&Phi;}}_{i,{\mathrm{\&Omega;}}_i}\right)}^{-1}{\mathrm{\&Phi;}}_{i,{\mathrm{\&Omega;}}_i}^{H}X$

然后，将矩阵的各个元素还原至S(i)中，在矩阵S(i)的其余位置补零，从而求解出矩阵S(i)；

S44：根据以下公式得出方向矩阵D的第i次迭代取值D(i)：

s.t.||S||_∞，0＝K

其中，||·||_F表示求F范数；

S45：对角矩阵D(i-1)中主对角线上的元素按顺序排列形成的列向量记为d(i-1)，对角矩阵D(i)中主对角线上的元素按顺序排列形成的列向量记为d(i)；判断||d(i)-d(i-1)||₂是否小于ε，其中，||·||₂表示求2-范数，ε为设定的迭代终止值且ε＞0；如果||d(i)-d(i-1)||₂小于ε，则得出方向矩阵D和信号稀疏矩阵S，D＝D(i)，S＝S(i)；否则，令i的值自增1，返回至步骤S42。

4.如权利要求3所述的相干信号源亚分辨率超分辨到达角估计方法，其特征在于，在步骤S44中，根据以下公式得出方向矩阵D的第i次迭代取值D(i)：

D(i)＝diag(d(i))

其中，||·||₂表示求2-范数，vec(·)表示矩阵的列向量化处理，T表示向量或矩阵的转置，表示K-R积，在得出的D(i)中，d(i)中的元素按照顺序排列在D(i)的主对角线上。

5.如权利要求3所述的相干信号源亚分辨率超分辨到达角估计方法，其特征在于，在步骤S44中，根据以下公式得出方向矩阵D的第i次迭代取值D(i)：

D(i)＝diag(d(i))

其中，||·||₁表示求1-范数，vec(·)表示矩阵的列向量化处理，T表示向量或矩阵的转置，表示K-R积， 表示Kronecker积，I_K为K阶单位矩阵，在得出的D(i)中，d(i)中的元素按照顺序排列在D(i)的主对角线上。

6.如权利要求1所述的相干信号源亚分辨率超分辨到达角估计方法，其特征在于，所述步骤S5具体包括以下步骤：方向矩阵D中主对角线上的非零元素共有K个，将方向矩阵D中主对角线上的非零元素筛选出来，方向矩阵D中主对角线上的第k个非零元素表示为D_k，k取1至K；D_k是方向矩阵D中主对角线上的第k′个元素，1≤k′≤N；则根据以下公式求解得出第k个相干信号源的到达角θ_k：