CN105116388A - 基于鲁棒主成分分析的天波超视距雷达瞬态干扰抑制方法 - Google Patents

Abstract

本发明公开了一种基于鲁棒主成分分析的天波超视距雷达瞬态干扰抑制方法，主要解决现有技术需先定位瞬态干扰位置导致干扰抑制效果差的问题。其实现步骤是：(1)用距离单元回波信号构造复数形式的Hankel矩阵，并将其分解为实部矩阵和虚部矩阵；(2)计算实部矩阵和虚部矩阵对应的稀疏矩阵和低秩矩阵；(3)根据实部的低秩矩阵与虚部的低秩矩阵获得对瞬态干扰和噪声抑制的回波信号；(4)对每一个距离单元重复步骤1-3，并进行多普勒处理，得到抑制瞬态干扰和噪声后的回波距离多普勒信息。本发明复杂度低，能明显的抑制噪声和多个同时存在的瞬态干扰，可用于对天波或地波超视距雷达中闪电冲击干扰、流星余迹回波干扰和人为冲击干扰的抑制。

Description

基于鲁棒主成分分析的天波超视距雷达瞬态干扰抑制方法

技术领域

本发明属于雷达技术领域，特别是涉及一种天波超视距雷达的瞬态干扰抑制方法，可用于对天波或地波超视距雷达中闪电冲击干扰、流星余迹回波干扰和人为冲击干扰的抑制。

背景技术

天波超视距雷达工作在3-30MHz的高频波段，利用电离层对大功率电波折射弯曲和下视传播而增大探测距离，其探测距离范围为800km～3500km，因此主要用于远程预警，探测的主要目标为飞机、舰船和弹道等。天波超视距雷达容易受到外界环境的干扰，其中瞬态干扰影响最大，其特点是持续时间短、强度大、带宽宽，比接收机的内噪声高20-40dB，因此容易淹没目标信号，极大的影响了天波超视距雷达检测目标的性能。目前可以通过自适应波束形成技术对接收机外部的干扰和噪声进行抑制，但瞬态干扰信号如闪电冲击干扰、流星余迹回波干扰以及人为的冲击干扰等，它们具有很强的能量，因此能够通过天线的副瓣进入接收机，导致了在目标的距离多普勒频谱上出现幅度较高、频带较宽的干扰能量。瞬态干扰通常持续时间较短，例如闪电冲击干扰持续时间为几百毫秒，流星余迹回波干扰的持续时间为0.1s到几秒，因此瞬态干扰只存在于目标回波的1～2个距离单元内。天波超视距雷达的回波信号经过匹配滤波和波束形成之后，一个距离单元的回波信号可以表示为r(m)＝c(m)+i(m)+s(m)+n(m),m＝1,2,…M，其中m表示慢时间序号，M表示一个相干脉冲处理间隔CPI时间内的脉冲数，r(m)表示回波信号，c(m)表示杂波分量，i(m)表示瞬态干扰，s(m)为目标回波分量，n(m)表示加性噪声。研究表明，海面被高频信号照射时，其海面杂波可以表示为两个Bragg分量，即：式中A_I与A_Q分别为两个Bragg分量的复幅度，为Bragg频率，f_c为雷达工作频率，单位为MHz，T为扫频周期。

目前，瞬态干扰抑制的主要方法是在慢时间域对干扰进行抑制，即在匹配滤波和波束形成之后，多普勒处理之前，对各个距离单元的采样数据分别进行瞬态干扰。主要分成两大类，一类是先对干扰位置进行检测，然后将干扰数据剔除掉，最后利用自回归AR模型、神经网络或压缩感知等数据插值的方法对信号进行恢复，例如，邢孟道,保铮,强勇.天波超视距雷达瞬态干扰抑制.电子学报,2002,30(6),pp.823-826，这种方法在检测干扰位置前，需要先进行杂波抑制，由于对杂波抑制的同时会抑制掉部分干扰能量，因此会影响到后续瞬态干扰的位置检测；另外，当定位到瞬态干扰位置后，将瞬态干扰置零，然后利用AR模型对中断数据进行恢复，当存在多个干扰位置的时候，数据的恢复效果并不理想。另一类利用信号分解理论，将回波信号分解为多个分量，根据瞬态干扰的特点，从多个分量中剔除干扰分量，最后利用剩余分量重新合成回波信号，例如周忠根,水鹏朗.基于复数据经验模式分解的天波超视距雷达瞬态干扰抑制.电子与信息学报,2011,33(12),pp.2831-2836，这种方法通过将时域数据分解为多个固有模态函数IMF，该方法无需AR模型进行数据恢复，但是由于目标回波存在于各个IMF分量，将含有干扰的IMF剔除，会导致信杂比降低，且抑制后仍然有较多的干扰残余分量，这种方法的另外一个缺点是算法复杂度较高，工程实现比较困难。

发明内容

本发明的目的在于针对上述已有技术的不足，提出一种基于鲁棒主成分分析的天波超视距雷达瞬态干扰抑制，以在无需检测瞬态干扰位置的条件下提高抑制效果。

本发明的技术思路是：将任意一距离单元回波信号r(t)构造为Hankel矩阵，根据杂波和信号的特性，及在Hankel矩阵中，杂波信号分量与目标信号分量具有低秩性，而瞬态干扰和加性噪声则表现出稀疏性的特点，采用鲁棒主成分分析RPCA的方法，对Hankel矩阵进行分解得到低秩矩阵L和稀疏矩阵S，从低秩矩阵中恢复出回波信号，完成天波超视距雷达瞬态干扰的抑制。其实现步骤包括如下：

(1)将第n个距离单元回波信号r(m)构造成复数形式的Hankel矩阵H(r)，其中n＝1,2,…,N，N表示距离单元数；

(2)将复数矩阵H(r)分解为实部矩阵H_R(r)和虚部矩阵H_I(r)；

(3)根据鲁棒主成分分析RPCA理论，求解实部矩阵H_R(r)的低秩矩阵L_R(r)和稀疏矩阵S_R(r)以及虚部矩阵H_I(r)的低秩矩阵L_I(r)和稀疏矩阵S_I(r)。

(4)将实部的低秩矩阵L_R(r)和虚部的低秩矩阵L_I(r)组合为复矩阵H′(r)，对复矩阵H′(r)进行求解，得到对瞬态干扰和噪声抑制后的回波信号r′(m)；

(5)对每一个距离单元重复步骤1-4，得到每个距离单元的回波信号；

(6)对所有距离单元的回波信号进行多普勒处理，获得回波的距离多普勒信息。

与现有技术相比，本发明具有以下优点：

第一，本发明针对天波超视距雷达进行瞬态干扰抑制时，无需事先对杂波进行抑制，避免了杂波剔除时引起瞬态干扰能量丢失导致的瞬态干扰位置检测困难，因此本发明能够更好的对瞬态干扰进行抑制。

第二，本发明通过将回波信号构造为Hankel矩阵，利用杂波信号与目标信号的Hankel矩阵具有低秩性，而瞬态干扰和加性噪声的Hankel矩阵具有稀疏性特点，采用鲁棒主成分分析算法对矩阵进行求解，无需定位瞬态干扰的位置，克服了现有技术在检测弱干扰信号时性能较差的缺点。

第三，本发明在对瞬态干扰信号抑制的同时，还可以对噪声进行抑制，提高了信噪比，与现有技术相比，可以更好的用于对较小的动目标进行检测。

第四，本发明无需利用自回归AR模型或其他方法对信号进行插值恢复，克服现有技术在抑制多个瞬态干扰时，中断数据恢复效果较差的缺点。

附图说明

图1为本发明的实现流程图；

图2为本发明仿真实验中第21个距离单元未添加瞬态干扰分量时的距离多普勒图；

图3为本发明仿真试验中第21个距离单元添加了瞬态干扰分量时的距离多普勒图；

图4为未经过本发明处理的目标回波距离多普勒图；

图5为经过本发明处理后的目标回波距离多普勒图；

图6为经过本发明处理后的第21个距离单元的回波信号多普勒图；

图7为经过本发明处理前后第21个距离单元的回波信号多普勒结果对比图。

具体实施方式

以下参照附图对本发明作进一步详细的描述。

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

步骤1：将第n个距离单元回波信号r(m)构造成复数形式的Hankel矩阵H(r)，其中n＝1,2,…,N，N表示距离单元数。

(1.1)将第n个距离单元回波信号r(m)表示为：

r(m)＝c(m)+i(m)+s(m)+n(m),m＝1,2,…M

其中m表示慢时间序号，M表示一个相干脉冲处理间隔CPI时间内的脉冲数，c(m)表示杂波分量，i(m)表示瞬态干扰，s(m)为目标回波分量，n(m)表示加性噪声；

(1.2)用回波信号r(m)构成复数形式的Hankel矩阵H(r)，表示为：

其中，表示复数域，P为Hankel矩阵的行数，Q为Hankel矩阵的列数，且M＝P+Q-1；

(1.3)根据回波信号r(m)＝c(m)+i(m)+s(m)+n(m)，将复数形式的矩阵H(r)写为：

H(r)＝H(c)+H(i)+H(s)+H(n)，

式中，为回波信号r(m)的Hankel矩阵，为杂波信号分量c(m)的Hankel矩阵，为瞬态干扰i(m)的Hankel矩阵，为目标信号分量s(m)的Hankel矩阵，为加性噪声n(m)的Hankel矩阵，rank(H(c))＜＜min(P,Q)，rank(H(s))＜＜min(P,Q)，rank表示求矩阵的秩，min表示取最小值。

步骤2：将复数形式的矩阵H(r)分解为实部矩阵H_R(r)和虚部矩阵H_I(r)：

H_R(r)＝Re(H(r))，H_I(r)＝Im(H(r))

其中，Re表示求矩阵H(r)的实部，Im表示求矩阵H(r)的虚部。 表示实数域。

步骤3：根据鲁棒主成分分析RPCA理论，求解实部矩阵H_R(r)的低秩矩阵L_R(r)和稀疏矩阵S_R(r)以及虚部矩阵H_I(r)的低秩矩阵L_I(r)和稀疏矩阵S_I(r)。

鲁棒主成分分析RPCA理论最早由Candès等人提出，用于从较大但存在稀疏误差的矩阵中恢复出本质上低秩的数据矩阵。其中恢复低秩矩阵是一个双目标优化问题，需要求解低秩矩阵核范数和稀疏矩阵的l₁范数联合最小化。经典的RPCA求解算法包括：奇异值阈值法SVT、近似梯度算法APG、交替方向算法ADM、精确拉格朗日乘子算法ALM以及不精确拉格朗日乘子算法IALM。在本实例中采用不精确拉格朗日乘子算法IALM，不精确拉格朗日乘子算法IALM通过构造增广拉格朗日函数对双目标优化问题进行求解。

增广拉格朗日函数如下：

$L(A,E,Y,\mathrm{\&mu;})=\left|\right|A|{|}_{*}+\mathrm{\&lambda;}|\left|E\right|{|}_1+<Y,D-A-E>+\frac{\mathrm{\&mu;}}{2}\left|\right|D-A-E|{|}_F^2$

式中，D表示观测数据矩阵，A表示低秩矩阵，E表示稀疏矩阵，Y为拉格朗日乘子矩阵，μ>0为惩罚算子，λ为正的实数，||·||_*表示矩阵的核范数，||·||₁示矩阵的l₁范数，<·,·>表示求矩阵的内积运算，||·||_F表示矩阵的F范数。

本步骤的实现如下：

(3.1)通过不精确拉格朗日乘子算法IALM对下述凸优化问题的求解,获得实部矩阵H_R(r)的低秩矩阵L_R(r)和稀疏矩阵S_R(r)：

$\underset{L_R\left(r\right),S_R\left(r\right)}{\min }\left|\right|L_R\left(r\right)|{|}_{*}+\mathrm{\&lambda;}|\left|S_R\left(r\right)\right|{|}_1,s.t.H_R\left(r\right)=L_R\left(r\right)+S_R\left(r\right)$

其中，||L_R(r)||_*表示矩阵L_R(r)的核范数，即矩阵的奇异值之和；||S_R(r)||₁表示矩阵S_R(r)的l₁范数；λ为正的实数，取值为max表示求最大值操作；s.t.表示约束条件；

(3.2)通过不精确拉格朗日乘子算法IALM对下述凸优化问题的求解，获得虚部矩阵H_I(r)的低秩矩阵L_I(r)和稀疏矩阵S_I(r)：

$\underset{L_I\left(r\right),S_I\left(r\right)}{\min }\left|\right|L_I\left(r\right)|{|}_{*}+\mathrm{\&lambda;}|\left|S_I\left(r\right)\right|{|}_1,s.t.H_I\left(r\right)=L_I\left(r\right)+S_I\left(r\right),$

其中，||L_I(r)||_*表示矩阵L_I(r)的核范数，即矩阵的奇异值的和；||S_I(r)||₁表示矩阵S_I(r)的l₁范数； $\mathrm{\&lambda;}=1/\sqrt{max(P,Q)}.$

步骤4：根据实部的低秩矩阵L_R(r)和虚部的低秩矩阵L_I(r)，得到对瞬态干扰和噪声抑制后的回波信号r′(m)。

(4.1)将实部的低秩矩阵L_R(r)和虚部的低秩矩阵L_I(r)组合为复矩阵H′(r)：

H′(r)＝L_R(r)+j*L_I(r)

式中， $j=\sqrt{-1};$

(4.2)计算复矩阵H′(r)中对应辅对角线上J个元素的平均值，得到恢复抑制瞬态干扰和噪声后的回波信号：

式中， $H=\{\begin{array}{cc}m,& 1\&le;m<Q\\ M-m+1,& Q\&le;m\&le;M\end{array},$ m＝1,2,…M，p和q分别为复矩阵H′(r)的行下标和列下标，且满足p+q-1＝m。

步骤5：重复步骤1-4，对每一个距离单元回波信号进行处理，得到抑制瞬态干扰和噪声后的回波信号。

步骤6：对所有距离单元抑制瞬态干扰和噪声后的回波信号进行多普勒处理，即对所有距离单元抑制瞬态干扰和噪声后的回波信号做FFT运算，获取回波的距离多普勒信息。

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

1)仿真条件与参数

所有仿真实验均在Windows8.1操作***下采用Matlab2012b软件实现。

雷达工作频率为20MHz，扫频周期T为0.025s，每一个相干脉冲处理间隔CPI内的脉冲数M＝512，距离单元数N＝32。信杂比为-11dB，杂噪比为-13dB。仿真目标出现在第21个距离单元，其多普勒频率为10Hz。在第21个距离单元内添加两个瞬态干扰分量，其表达式i₁(m)和i₂(m)由下式给出：

$i_1\left(m\right)=30\exp \left(\frac{{(mT-2.5)}^2}{2\&times;{0.01}^2}\right)\exp (-j2\mathrm{\&pi;}(-mT+0.047{\left(mT\right)}^{1.5})$

$i_2\left(m\right)=3\exp \left(\frac{{(mT-1.2)}^2}{2\&times;{0.05}^2}\right)\exp (-j2\mathrm{\&pi;}(-mT+0.06{\left(mT\right)}^{2.3})$

添加的两个瞬态干扰位置出现在第21个距离单元中的1.2s和2.5s的位置。

仿真使用距离多普勒谱图如图2-图4所示，其中：

图2为第21个距离单元未添加瞬态干扰分量时距离多普勒谱图，从图2中可以看出，目标信号出现在多普勒频率为10Hz的地方。

图3给出了第21个距离单元添加了瞬态干扰分量时的距离多普勒谱图，从图3中可以明显看出，瞬态干扰将多普勒谱的噪底抬高，目标信号已经被完全淹没。

图4给出了目标回波的距离多普勒谱图，从图4中也可以看出瞬态干扰已经将目标信号覆盖。

2)仿真内容

仿真1：采用本发明方法对32个距离单元的回波信号按照步骤1-5进行处理，得到处理后的目标回波信号，对处理后的目标回波信号进行距离多普勒处理仿真，得到图5所示距离多普勒图，对比图4与图5可以看出，瞬态干扰已经被抑制。

仿真2：采用本发明方法对处理后的目标回波第21个距离单元信号进行距离多普勒处理仿真，得到图6所示回波信号多普勒图。将图6与图2、图3对比可知，经过本发明方法处理后的第21个距离单元的距离多普勒图上可以很容易看出目标的多普勒频率为10Hz，从而验证了本发明方法可以用于对瞬态干扰的抑制。

仿真3：采用本发明方法对经过本发明处理前后的第21个距离单元的回波信号进行多普勒处理仿真，得到图7所示多普勒结果对比图。从图7可以看出，本发明方法不但对瞬态干扰有很好的抑制效果，而且对噪声具有一定的抑制作用，这对于能量较小的动目标检测具有很大的优势。

Claims (6)

1.一种基于鲁棒主成分分析的天波超视距雷达瞬态干扰抑制方法，包括如下步骤：

(1)将第n个距离单元回波信号r(m)构造成复数形式的Hankel矩阵H(r)，其中n＝1,2,…,N，N表示距离单元数；

(2)将复数矩阵H(r)分解为实部矩阵H_R(r)和虚部矩阵H_I(r)；

(3)根据鲁棒主成分分析RPCA理论，求解实部矩阵H_R(r)的低秩矩阵L_R(r)和稀疏矩阵S_R(r)以及虚部矩阵H_I(r)的低秩矩阵L_I(r)和稀疏矩阵S_I(r)。

(4)用实部的低秩矩阵L_R(r)和虚部的低秩矩阵L_I(r)组合为复矩阵H′(r)，对H′(r)进行求解，得到对瞬态干扰和噪声抑制后的回波信号r′(m)；

(5)对每一个距离单元重复步骤1-4，得到每个距离单元的回波信号；

(6)对所有距离单元的回波信号进行多普勒处理，获得回波的距离多普勒信息。

2.根据权利要求1所述的基于鲁棒主成分分析的天波超视距雷达瞬态干扰抑制方法，其中所述步骤(1)中将一个距离单元回波信号r(m)构造成复数形式的Hankel矩阵H(r)，按如下步骤进行：

(2.1)将第n个距离单元回波信号r(m)表示为：

r(m)＝c(m)+i(m)+s(m)+n(m),m＝1,2,…M

其中m表示慢时间序号，M表示一个相干脉冲处理间隔CPI时间内的脉冲数，r(m)表示回波信号，c(m)表示杂波分量，i(m)表示瞬态干扰，s(m)为目标回波分量，n(m)表示加性噪声。

用回波信号r(m)构成Hankel矩阵H(r)，表示为：

其中，表示复数域，P为Hankel矩阵的行数，Q为Hankel矩阵的列数，M为一个相干脉冲处理间隔CPI内的脉冲数，且M＝P+Q-1；

(2.2)根据回波信号r(m)＝c(m)+i(m)+s(m)+n(m)，将矩阵H(r)写为：

H(r)＝H(c)+H(i)+H(s)+H(n)

式中，为回波信号r(m)的Hankel矩阵，为杂波信号分量c(m)的Hankel矩阵，为瞬态干扰i(m)的Hankel矩阵，为目标信号分量s(m)的Hankel矩阵，为加性噪声n(m)的Hankel矩阵，rank(H(c))＜＜min(P,Q)，rank(H(s))＜＜min(P,Q)，rank表示求矩阵的秩，min表示取最小值。

3.根据权利要求1所述的基于鲁棒主成分分析的天波超视距雷达瞬态干扰抑制方法，其中所述步骤(3)中根据鲁棒主成分分析RPCA理论计算实部矩阵H_R(r)的低秩矩阵L_R(r)和稀疏矩阵S_R(r)，是通过不精确拉格朗日乘子算法IALM对下述凸优化问题的求解获得：

其中，||L_R(r)||_*表示矩阵L_R(r)的核范数，即矩阵的奇异值的和；||S_R(r)||₁表示矩阵S_R(r)的l₁范数；λ为正的实数，取值为max表示求最大值操作；s.t.表示约束条件。

4.根据权利要求1所述的基于鲁棒主成分分析的天波超视距雷达瞬态干扰抑制方法，其中所述步骤(3)中根据鲁棒主成分分析RPCA理论计算虚部矩阵H_I(r)的低秩矩阵L_I(r)和稀疏矩阵S_I(r)，是通过不精确拉格朗日乘子算法IALM对下述凸优化问题的求解获得：

其中，||L_I(r)||_*表示矩阵L_I(r)的核范数，即矩阵的奇异值的和；||S_I(r)||₁表示矩阵S_I(r)的l₁范数；

5.根据权利要求1所述的基于鲁棒主成分分析的天波超视距雷达瞬态干扰抑制方法，其中所述步骤(4)中的复矩阵H′(r)，表示如下：

H′(r)＝L_R(r)+j*L_I(r)

式中，L_R(r)为实部的低秩矩阵，L_I(r)为虚部的低秩矩阵。

6.根据权利要求1所述的基于鲁棒主成分分析的天波超视距雷达瞬态干扰抑制方法，其中所述步骤(4)中对复矩阵H′(r)进行求解，是计算复矩阵H′(r)中对应辅对角线上的J个元素的平均值，表示为：

式中，r′(m)表示回波信号，m＝1,2,…M，p和q分别为复矩阵H′(r)的行下标和列下标，且满足p+q-1＝m。