基于高阶奇异值分解的天波雷达海杂波抑制算法

薄 超 顾 红 苏卫民 吕 婧

（南京理工大学电子工程技术研究中心，江苏 南京210094）

引 言

天波超视距雷达（Over The Horizon Radar，OTHR）是利用电离层对高频段信号的反射作用而进行超视距目标探测的雷达系统，其探测距离范围可以达到800km～3 500km，因此它的一个主要任务是提供大范围的早期预警［1-3］.OTHR的探测目标主要有两类：飞机和舰船，其中飞机的径向运动速度较快，产生的多普勒谱峰和强大海杂波的Bragg峰相距较远，其检测不易受到海杂波影响，而舰船的径向运动速度较慢，产生的多普勒谱峰和强大海杂波的Bragg峰相距较近，易淹没于海杂波中，致使舰船很难被检测［1］.为了增强OTHR的预警能力，实现海杂波背景下的舰船检测尤为重要.

短相干积累时间（Coherent Integration Time，CIT）条件下的舰船检测方法主要有海杂波对消方法［4-5］和现代谱估计方法［6-11］.海杂波对消方法通过估计强海杂波及其谐波分量对应的频率、幅值和初相3个参数，逐次地对消强海杂波及其较强的谐波分量，凸显舰船目标.然而此类方法需要一个判决门限决定其迭代次数，而该判决门限仅能通过经验得出，不恰当的判决门限将导致目标回波的对消，此种情况在多目标检测情况下尤为突出.现代谱估计方法主要有基于线性预测的海杂波抑制方法［6］和基于子空间估计的海杂波抑制方法［7-11］，前者采用线性预测方法对回波数据进行预处理，延长回波数据长度，提高多普勒分辨率，从而将舰船从海杂波中分开，然而此种方法的线性预测模型参数较难选择，影响舰船检测性能；后者根据杂波和舰船在子空间的聚集特性来实现海杂波抑制，文献［7］中方法的子空间估计精度不高，导致海杂波向目标子空间泄露，降低了信号与杂波加噪声之比（Signal to Clutter Plus Noise Ratio，SCNR），且在海态不平稳时，相邻距离单元和方位单元的相关性会降低，文献［8-11］中方法的抑制海杂波性能将下降，上述问题严重影响目标检测和参数估计.

针对上述问题，引入基于高阶奇异值分解（Higher-Order Singular Value Decomposition，HOSVD）子空间估计方法［12］，提高子空间估计精度.与现有子空间类的海杂波抑制方法相比，提出方法的优势主要体现在以下两点：1）子空间估计时仅需要一个距离单元数据，在复杂海态（相邻距离单元和方位单元相关性较弱）下亦能很好抑制海杂波；2）三维正交投影，更加有效地抑制海杂波，从而提高了SCNR和峰值旁瓣电平比（Peak Sidelobe Level Ratio，PSLR）.

1 海杂波产生机理

海浪的运动状态是一种复杂和混沌的随机过程［13］，由不同频率、不同波高和不同传播方向的海浪叠加而成，且这些海浪近似呈现正弦波动，如图1所示.

图1 Bragg峰谐振散射机理

设海浪中某一正弦波动的波长为L，雷达发射信号的波长为λ.当雷达发射源到海浪各反射点的路程差均为半波长λ／2时，各反射回波之间的波程差均为一个波长λ，海浪与高频雷达发射信号产生谐振，回波信号同相相加而加强，形成海杂波的一阶散射Bragg峰.相反，如果雷达发射源到海浪各反射点的路程差不等于半波长λ／2时，回波信号非同相相加，信号功率弱于同相相加回波，产生海杂波的高阶散射.

基于上述散射机理可得

式中，α表示雷达波束的擦地角.根据深水区动力学原理，深水重力波的流速为其中g表示重力加速度.与雷达相比，水流的径向速率为

由式（1）和（2）可得Bragg峰的多普勒频率为

式中：fB表示Bragg峰的多普勒频率，Hz；±表示谐振海浪的运动方向指向及背离雷达波束；f0表示雷达载频，MHz.通常雷达波束的擦地角α较小，则式（3）可近似为

式中fB表示理想Bragg峰的多普勒频率，在实际环境中海浪还受到海风和潮汐等作用，产生多普勒频率为fC的洋流，从而实测数据中海杂波Bragg峰多普勒频率为fB＋fC.

2 基于HOSVD子空间估计的海杂波抑制算法

将舰船所在距离单元的数据s＝［s1，…，sN］按照三阶Hankel数据张量结构排列，应用HOSVD方法求解Hankel数据张量的目标子空间投影矩阵，利用目标子空间投影矩阵将Hankel数据张量投影到目标子空间，求得海杂波抑制后数据张量，采用重构方程进行海杂波抑制后数据张量的矢量化，得到海杂波抑制后的一维数组，对一维数组进行傅里叶变换，求得海杂波抑制后的频谱，此时，频谱上的明显尖峰便是舰船目标.

2.1 一维数据的Hankel张量表示

三阶Hankel张量H∈CⅠ1×Ⅰ2×Ⅰ3内各元素可用hi1i2i3表示，其中Ⅰ1，Ⅰ2和Ⅰ3分别表示张量的长，宽和高，1≤i1≤Ⅰ1，1≤i2≤Ⅰ2和1≤i3≤Ⅰ3表示张量的元素坐标［14］.设舰船所在距离单元的回波信号为s＝［s1，…，sN］，其前P（P＜N）个元素［s1，…，sP］构造Hankel张量的第一个Ⅰ1×Ⅰ2维Hankel矩阵，Hankel张量的第二个Ⅰ1×Ⅰ2维Hankel矩阵由元素［s2，…，sP＋1］构成，以此类推构造Hankel张量的其他Hankel矩阵，直至一维数组s的最后一个元素，这些Hankel矩阵从前至后依次排列得到三阶Hankel张量H∈CⅠ1×Ⅰ2×Ⅰ3，具体排列方法如图2所示.

图2 Hankel张量的排列结构

图2中三阶Hankel张量的每个元素hi1i2i3满足

元素坐标满足1≤i1≤Ⅰ1，1≤i2≤Ⅰ2和1≤i3≤Ⅰ3，张量长宽高之和满足Ⅰ1＋Ⅰ2＋Ⅰ3＝N＋2.如图2所示，分别沿Hankel张量的长i1、宽i2和高i3进行切割，并将各自的切片矩阵向量化表示，可以得到以下三个矩阵，称其为三阶张量的等效模式展开矩阵［15］：

式中H（i）（i＝1，2，3）表示第i等效模式展开矩阵；Hi1，·，·，H·，i2，·和H·，·，i3分别表示沿长i1、宽i2和高i3切割的切片矩阵；vec（·）表示矩阵按列展开；［·］T表示矩阵或向量的转置运算.其中3个等效模式展开矩阵等价，当任一等效模式展开矩阵已知时可求出相应的三阶张量.

2.2 Hankel张量子空间计算和投影

根据文献［15］中的HOSVD定义，三阶Hankel张量H∈CⅠ1×Ⅰ2×Ⅰ3进行HOSVD得

式中：C∈CⅠ1×Ⅰ2×Ⅰ3表示核张量；Uk表示H的第k等效模式展开矩阵H（k）的左奇异矩阵，k＝1，2，3；×k表示张量与矩阵的模k乘积，即H′＝C×kUk⇔H′（k）＝UkC（k），k＝1，2，3.通常，舰船目标回波信号能量远小于海杂波，所以，Uk中大奇异值对应的列向量张成海杂波子空间，小奇异值对应的列向量张成目标子空间.与此同时，核张量可表示为

［·］H表示矩阵或向量的共轭转置运算.根据文献［15］中张量性质，将式（8）代入式（7）可得

由H′＝C×kUk⇔H′（k）＝UkC（k）可知，式中UkUHk是H的第k等效模式展开矩阵H（k）的投影矩阵，k＝1，2，3.依据多重信号分类（Multiple Signal Classification，MUSIC）算法，三阶数据张量中各个模式展开矩阵中海杂波子空间投影矩阵可表示为，目标子空间投影矩阵可表示为，其中Uks由Uk中大特征值对应的列向量组成，k＝1，2，3.将三阶Hankel张量各模式展开矩阵H（k）投影到其目标子空间，k＝1，2，3，从而抑制海杂波后的数据张量为

式中，张量H″∈CⅠ1×Ⅰ2×Ⅰ3是张量H在其目标子空间的投影张量.

2.3 Hankel张量重构

根据重构方程，将海杂波抑制后的Hankel张量H″∈CⅠ1×Ⅰ2×Ⅰ3重构为一维数据，其中重构方程为

式中：s″n表示海杂波抑制后的一维数据元素；表示Hankel张量H″中长、宽和高分别为i1、i2和i3的元素；表示Hankel张量H″中元素系数满足i1＋i2＋i3-2＝n的元素之和；kn表示Hankel张量H″中元素系数满足i1＋i2＋i3-2＝n的元素个数.对s″＝进行傅里叶变换可以得到海杂波抑制后的频谱，此时频谱上的明显尖峰便是舰船目标.

2.4 算法步骤

基于HOSVD子空间估计的海杂波抑制算法基本步骤：

1）根据式（5）中一维数据与Hankel张量之间的转换关系，将目标所在距离单元的一维数据转换为Hankel数据张量.

2）依据式（6）求得Hankel数据张量的各个模式展开矩阵，应用HOSVD计算各个模式展开矩阵的左奇异矩阵Uk，进而得到海杂波子空间投影矩阵，k＝1，2，3.

4）依据式（11）重构方程将Hankel数据张量H″转换为一维数据s″，采用傅里叶变换处理数据s″，获得海杂波抑制后的频谱，此时频谱上的明显尖峰便是舰船目标.

3 理论仿真实验和基于实测数据的仿真实验的处理结果与分析

3.1 理论仿真实验处理结果与分析

通过理论仿真实验对比奇异值分解（Singular Value Decomposition，SVD）［7］算法和HOSVD算法的海杂波抑制性能.载频为f0＝7.5MHz；脉冲重复周期为0.7s；脉冲重复周期数为64；Bragg峰的多普勒频率为fB＝±0.28Hz；Bragg峰信噪比（Signal to Noise Ratio，SNR）分别为20和10dB；加入噪声为复高斯白噪声；Hankel矩阵列数为32；目前，没有详细理论推导说明如何选择Ⅰ1，Ⅰ2和Ⅰ3能得到Hankel张量子空间估计的最优解，而Ⅰ1＝Ⅰ2＝Ⅰ3时可用优化算法降低Hankel张量HOSVD的计算量［16］，理论仿真实验中张量维数为Ⅰ1＝Ⅰ2＝Ⅰ3＝22.

采用理论仿真实验，分别从海杂波抑制后输出的PSLR和SCNR两个方面说明本算法的优越性.

仿真实验1：在仿真海杂波数据中加入三个目标，目标多普勒频率分别为0.48、0.2和-0.17Hz，目标SNR分别为-3、-2和1dB，且SNR均小于海杂波Bragg峰.奇异值分解时大特征值对应能量较大（即频谱中幅值较高）的频率点，而从图3（a）可以看出海杂波Bragg峰的幅值较大，且高于目标幅值，选择特征个数为2进行海杂波抑制.海杂波抑制后的傅里叶频谱如图3（b）所示，图中显示SVD算法和HOSVD算法均能抑制海杂波，使目标凸显出来，然而HOSVD算法抑制海杂波后的PSLR较高，且在海杂波Bragg峰的位置产生较深的零陷.与SVD算法相比，HOSVD算法的PSLR提高了4dB左右.上述现象是由于HOSVD算法分别计算数据张量各个模式展开矩阵的海杂波子空间和目标子空间，并将数据张量投影到目标子空间，减少海杂波子空间向目标子空间扩散.目标与Bragg峰重合时两种方法均失效.

图3 64个脉冲仿真实验的SVD算法和HOSVD算法抑制海杂波性能对比

仿真实验2：在仿真海杂波数据中加入三个目标，目标多普勒频率分别为0.48、0.2和-0.17Hz，仿真实验中加入目标SNR均小于海杂波Bragg峰，选择特征值个数为2，仿真实验中蒙特卡洛实验次数为100.图4是分别运用SVD算法和HOSVD算法抑制海杂波后输出的SCNR变化曲线，其中海杂波抑制后输出SCNR的详细计算方法见附录A.图4中显示无论目标远离还是靠近Bragg峰，HOSVD算法的海杂波抑制性能均优于SVD算法，其数值大小为2.5dB，且具有较好的稳定性.上述现象是由于HOSVD算法能够在海杂波Bragg峰位置产生较深的零陷，减少了海杂波的残余分量，从而提高了海杂波抑制后输出的SCNR.

图4 64个脉冲仿真实验时两种算法抑制海杂波后SCNR变化曲线

3.2 基于实测数据的仿真实验处理结果与分析

通过基于实测数据的仿真实验，对比SVD［7］算法和HOSVD算法的海杂波抑制性能.实验采用武汉大学提供的高频雷达海杂波实测数据对文中算法进行检验，载频为f0＝7.5MHz；脉冲重复周期为0.726 4s；Bragg峰的多普勒频率为fB＝±0.29 Hz.天波雷达短CIT的取值范围为4～10s，舰船检测时信号重复周期为100～200ms甚至更长［1］227，230，短CIT的脉冲数目取值范围应小于100.由于实测数据有限，仅有信号重复周期为0.726 4s的海杂波实测数据，选取脉冲数目为64的数据进行算法验证，其取值小于100，属于短CIT范畴，且实验中Hankel矩阵列数为32，张量维数Ⅰ1＝Ⅰ2＝Ⅰ3＝22.

采用实测海杂波数据验证算法，分别从海杂波抑制后输出的PSLR和SCNR两个方面说明本算法的优越性.

实测实验1：在某一个距离单元内海杂波数据中加入SCNR分别为-34、-40和-36dB的目标，其多普勒频率分别为0.47、0.19和-0.18Hz，实验中大特征值个数从1开始依次加1，当大特征值个数为3时海、地杂波均被抑制，三个目标凸显出来.目标所在距离单元海杂波抑制前后的傅里叶频谱如图5所示，图5（a）中两个较强的谱峰是海杂波的正负Bragg峰，零频附近较强的谱峰是地杂波.分别采用SVD算法和HOSVD算法抑制海、地杂波，抑制后的对比结果如图5（b）所示，图中显示两种算法均能抑制强海杂波凸显目标，即在目标所在位置出现尖峰，然而HOSVD算法抑制海杂波后的PSLR相对较高，且在海杂波和地杂波所在频率点上产生较深的凹口，图中0.47Hz处目标峰值低于-0.18Hz处目标，是由于0.47Hz处高阶海杂波幅值较强所致.与SVD方法相比，HOSVD算法的PSLR提高了5dB左右.上述现象与仿真实验1基本相符，说明HOSVD算法在实测数据实验中亦优于SVD算法.

图5 64个脉冲实测数据的SVD方法和HSVD方法抑制海杂波性能对比

图6 64个脉冲实测数据时两种方法抑制海杂波后SCNR变化曲线

实测实验2：在某一个距离单元内海杂波数据中加入目标，多普勒频率分别为0.47、0.19和-0.18Hz，实验中海杂波抑制前输入目标SCNR区间为-40至-30dB，选择特征值数为3时均能凸显目标抑制海、地杂波，实验中选择特征值个数为3.图6是分别运用SVD算法和HOSVD算法抑制海杂波后的SCNR变化曲线.图6中显示无论目标远离还是靠近Bragg峰，HOSVD算法的海杂波抑制性能均优于SVD算法，海杂波抑制后输出SCNR提高2dB，与仿真实验2基本相符.当目标距海杂波较近时SVD算法输出的SCNR下降较大，且SVD算法海杂波抑制后的SCNR输出随目标所在位置变化较大，即SVD算法对目标所在位置较敏感.上述现象是由于实测海杂波数据中回波信号除了Bragg峰外还存在高阶海杂波分量，且噪声形式较为复杂，不一定为高斯白噪声，从而影响了算法的抑制性能.

4 结 论

在子空间算法的基础上，提出了基于HOSVD子空间估计的海杂波抑制算法，对该算法进行理论分析和推导，并应用理论仿真和海杂波实测数据验证了算法的有效性.理论仿真实验结果和基于海杂波实测数据的实验结果均表明，与SVD海杂波抑制方法相比，该方法有利于凸显目标和提高输出SCNR，但本文算法需要三次奇异值分解，计算量较大，因此降低计算复杂度是下一步工作重点.

致谢：衷心感谢武汉大学提供高频雷达实测数据，使得本文算法得以实测验证.

［1］周文瑜，焦培楠.超视距雷达技术［M］.北京：电子工业出版社，2008：224.

［2］王子曦，胡进峰，肖赛军，等.天波雷达在不规则地形中的接收阵列天线综合［J］.电波科学学报，2012，27（4）：672-679.WANG Zixi，HU Jinfeng，XIAO Saijun，et al.Synthesis of OTHR receive antenna arrays in irregular terrain［J］.Chinese Journal of Radio Science，2012，27（4）：672-679.（in Chinese）

［3］游 伟，何子述，陈绪元，等.基于三次相位建模的天波雷达污染 校 正［J］.电 波 科 学 学 报，2012，27（5）：875-880.YOU Wei，HE Zishu，CHEN Xuyuan，et al.Skywave radar decontamination based on the cubic phase model［J］.Chinese Journal of Radio Science，2012，27（5）：875-880.（in Chinese）

［4］郭 欣，倪晋麟，刘国岁.短相干积累条件下天波超视距雷达的舰船检测［J］.电子与信息学报，2004，26（4）：613-618.GUO Xin，NI Jinlin，LIU Guosui.The ship detection of sky wave over-the-horizon radar with short coherent interaction time［J］.Journal of Electronics ＆Information Technology，2004，26（4）：613-618.（in Chinese）

［5］仇永斌，张 宁，张树春.双基地高频地波雷达海杂波抑制［J］.哈尔滨工业大学学报，2012，44（1）：71-77.CHOU Yongbin，ZHANG Ning，ZHANG Shuchun.Ocean clutter suppression for a bistatic HF ground wave radar［J］.Journal of Harbin Institute of Technology，2012，44（1）：71-77.（in Chinese）

［6］CHEN J W，GAO S.Detection of ship for OTHR based on AR-MUSIC algorithm［C］／／Proc of the International Conference on Wireless Communications ＆Signal Processing.Nanjing，2009：1-4.

［7］LU K，LIU X Z，LIU Y T.Ionospheric decontamination and sea clutter suppression for HF skywave radars［J］.IEEE Journal of Oceanic Engineering，2005，30（2）：455-462.

［8］WANG G，XIA X G，ROOT B T，et al.Maneuvering target detection in over-the-horizon radar using adaptive clutter rejection and adaptive chirplet transform［J］.IEE Proceedings-Radar，Sonar and Navigation，2003，150（4）：292-298.

［9］ZHAO Z G，CHEN J W，BAO Z.A method to estimate subspace via Doppler for ocean clutter suppression in skywave radars［C］／／Proc of the IEEE CIE International Conference on Radar.Beijing，2011：145-148.

［10］赵志国，陈建文，鲍 拯.一种改进的OTHR自适应海杂波抑制算法［J］.系统工程与电子技术，2012，34（5）：909-914.ZHAO Zhiguo，CHEN Jianwen，BAO Zheng.Modified adaptive ocean clutter suppression approach in OTHR［J］.Systems Engineering and Electronics，2012，34（5）：909-914.（in Chinese）

［11］邢孟道，保 铮，强 勇.天波超视距雷达瞬态干扰抑制［J］.电子学报，2002，30（6）：823-826.XING Mengdao，BAO Zheng，QIANG Yong.Transient interference excision in OTHR［J］.Acta Electronica Sinica，2002，30（6）：823-826.（in Chinese）

［12］HAARDT M，ROEMER F，GALDO G D.Higherorder SVD-based subspace estimation to improve the parameter estimation accuracy in multidimensional harmonic retrieval problems［J］.IEEE Transactions on signal processing.2008，56（7）：3198-3213.

［13］盛 文，任 吉.天波超视距雷达海杂波的混沌动态特性分析［J］.电波科学学报，2012，27（2）：350-357.SHENG Wen，REN Ji.Analysis of chaotic dynamics of skywave over-the-horizon radar sea clutter［J］.Chinese Journal of Radio Science，2012，27（2）：350-357.（in Chinese）

［14］PAPY J M，LATHAUWER L D，HUFFEL S V.Exponential data fitting using multilinear algebra：The single-channel and multi-channel case［J］.Numerical Linear Algebra with Applications，2005，12（8）：809-826.

［15］LIEVEN D L，BART D M，JOOS V.A multilinear singular value decomposition［J］.SIAM Journal on Matrix Analysis and Applications，2000，21（4）：1253-1278.

［16］BADEAU R，BOYER R.Fast Multilinear singular value decomposition for structured tensors［J］.SIAM J Matrix Anal Appl，2008，30（3）：1008-1021.

附录A

输出SCNR计算方法

某一距离单元的回波信号可表示为

式中：st表示目标回波信号矢量；sc表示海、地杂波回波信号矢量；n表示噪声矢量.

SVD算法：

目标回信号是人为加入，因此可将目标回波信号矢量st和总回波信号矢量s分别应用Hankle矩阵表示为：

式中：Ht表示目标回波信号Hankel矩阵；H表示总回波信号Hankel矩阵；C为Hankel矩阵的列数.对H进行奇异值分解可得

式中：U表示左奇异矩阵；Σ表示对角矩阵；V表示右奇异矩阵.假设海、地杂波对应大奇异个数为K，则取左奇异矩阵的前K列为

根据U′求得目标子空间的投影矩阵为I-U′U′H，将H和Ht分别投影到目标子空间可得：

式中：H′表示海杂波抑制后的总回波信号矩阵；Ht′表示海杂波抑制后的目标信号矩阵.应用Hankel矩阵重构方程将上述两个矩阵重构，得到一维数据s′和s′t，其中s′包含目标信号、杂波残余分量和噪声；s′t包 含目标信号，将s′中的目标信号减去便可求得杂波残余分量和噪声，如下式所示：

式中：s′c表示杂波残余分量矢量；n′表示投影后噪声矢量.求得目标信号和杂波残余分量加噪声信号后通过下式计算海杂波抑制后输出SCNR：

HOSVD算法：

目标回信号是人为加入，因此可将目标回波信号矢量st和总回波信号矢量s分别应用Hankle张量表示为Ht和H.计算张量H各模式展开矩阵的左奇异矩阵.假设海、地杂波对应大奇异个数为K，则取各模式展开矩阵左奇异矩阵的前K列为

根据Uks求得目标子空间的投影矩阵为I-将张量H和Ht分别投影到目标子空间可得：

式中：H″表示海杂波抑制后的总回波信号张量；Ht″表示海杂波抑制后的目标信号张量.应用Hankel张量重构方程将上述两个张量重构，得到一维数据s″和s″t，其中s″包含目标信号、杂波残余分量和噪声；s″t包含目标信号，将s″中的目标信号减去便可求得杂波残余分量和噪声，如下式所示为

求得目标信号和杂波残余分量加噪声信号后通过下式计算海杂波抑制后输出SCNR.