基于时频脊-Radon变换的海面小目标检测方法

伍僖杰 丁 昊, 2 刘宁波 关 键

(1. 海军航空大学， 山东烟台 264001; 2. 军事科学院国防科技创新研究院， 北京 100071)

1 引言

随着雷达探测环境的复杂化，目标形态的多样化，海面慢速小目标速度低、散射截面积(Radar Cross Section, RCS)小的特性[1]日益凸显，导致雷达对弱小目标的检测愈发成为一个重难点问题。为了更有效地检测海面弱小目标，研究者们提出了很多方法。文献[2]从9个方面出发，对这些方法做了分类综述。其中，从检测机理层面可以分为能量检测方法和特征检测方法。能量检测方法主要依据海杂波的局部幅度或功率水平信息构造似然比，根据门限因子形成检测门限并对目标存在与否做出判决[3-4]。该方法对目标的信杂比(Signal-to-Clutter Ratio, SCR)有较高要求，且在海尖峰密集场景下易引起大量虚警，导致实际检测性能难以达到预期。

特征检测方法通过充分挖掘海杂波与目标回波特征差异性，将二者从高重叠的观测空间转换到低重叠的特征空间，在特征空间内实现目标单元与杂波单元的甄别。与基于严格统计模型的能量检测方法不同，特征提取和融合往往是直观、经验、启发式的[5]，可以根据实际海面环境和雷达系统参数，从雷达回波的幅度、多普勒谱、时频图、极化信息等多个维度提取[6]，由于特征信息不再受低阶矩或统计模型限制，因此具有更大的自由度，在改善弱小目标检测性能方面具有更好的潜力。

目前已提出的特征检测方法有很多种，如时/频域分形特征检测方法[7-9]、多普勒谱非广延熵特征检测方法[10]、三特征和多特征融合检测方法[11-15]等。近年来，基于时频分析的特征检测方法同样受到国内外许多专家学者的关注[16-20]，该方法能够融合应用雷达回波的时域和频域信息，且可以实现目标能量的相干积累。常用的时频分析方法包括短时傅里叶变换(Short-Time Fourier Transform, STFT)、分数阶傅里叶变换(Fractional Fourier Transform, FRFT)、Wigner-Ville分布(Wigner-Ville Distribution, WVD)、平滑伪Wigner-Ville分布(Smoothed Pesudo-Wigner-Ville Distribution, SPWVD)、小波变换(Wavelet Transform, WT)等，通过参数或非参数化方法在二维平面上提取与目标对应的特征来实现海杂波背景下的目标检测[6]。例如，文献[5]提出的特征检测方法融合应用了归一化SPWVD中的3个时频特征，分别为时频脊累积量(Ridge Integration, RI)、连通区域个数(Number of connected Regions, NR)和最大连通区域尺寸(Maximal Size of connected regions, MS)，经实测数据验证，将上述特征应用到海面漂浮小目标检测中，是提高检测性能的一种有效途径。

当雷达驻留时间在秒量级或百毫秒量级时，上述特征检测器大多可在低SCR条件下达到较高的检测概率。然而，在扫描观测模式下，目标驻留时间通常难以达到上述量级，随着积累时间的降低，特征检测方法性能下滑严重。同样以文献[5]提出的方法为例，在虚警率为10-3情况下，当积累时间由1.024 s降为0.128 s时，HH极化数据的检测概率由0.821降为0.629，而VV极化时则由0.789降为0.599。实际上，雷达采用的相干脉冲数通常较少(大多情况下为64个或更少)，可用积累时间基本在几十毫秒量级，在该时间尺度上，易导致特征不稳定，进而影响海杂波与目标可分性。为提升该情形下雷达的检测性能，本文提出一种基于时频脊-Radon变换(Ridges-Radon Transform, RRT)的帧平滑双特征检测方法，在常规杂波抑制和时频域特征提取的基础上，通过多帧扫描历史数据和当前帧数据的综合应用，对特征做帧平滑处理以增强可分性，该方法基本思路为：首先采用分块白化滤波方法对海杂波进行抑制，在获得待检测单元及杂波单元时频分布的基础上，通过时频脊-Radon变换方法提取出时频脊峰值和频宽两个特征，并沿帧时间维进行特征平滑处理，在此基础上进行特征融合，在双特征平面上利用凸包算法进行目标检测。最后，分别采用2级和4级海况实测数据对本文方法的检测性能进行了分析。

2 分块白化杂波抑制方法

假设对海雷达在一个波束内发射一系列相干脉冲，并在每个距离单元接收到长度为N的回波时间序列，则海杂波中的目标检测可以归结为以下的二元假设检验问题：

(1)

式中，x(n)、xr(n)分别代表检测单元、参考单元信号，s(n)、c(n)分别代表目标信号、海杂波信号，R为参考单元数目。

复合高斯模型是对海杂波进行建模的常用模型，该模型将海杂波分解为两种分量：散斑分量(Speckle Component, SC)和纹理分量(Texture Component, TC)，其中SC为零均值复高斯过程，TC为非负随机过程。当雷达的积累时间较短时，由于纹理分量可近似看作常数，复合高斯模型可以简化为球不变随机向量(Spherically Invariant Random Vector)模型，自适应检测器基于该模型利用待检测单元周围的参考单元数据估计杂波的协方差矩阵(Sample Co-variance Matrix, SCM)，完成杂波白化。然而，在高分辨海杂波情形下，参考单元数往往小于雷达的相干脉冲数，杂波抑制可能出现较大误差，解决该问题的一种有效方法是分块白化[21]。具体流程如下：

首先将待测单元序列x(n)和参考单元序列xr(n)截成长度为L的不重叠向量ui、up,i，如下式：

为保证检测器性能，参考单元数R与截断长度L需满足关系：R≥2L。待检测单元的杂波协方差矩阵可用参考单元的归一化样本协方差矩阵(Normalized SCM, NSCM)估计，即：

(3)

(4)

类似地，对所有距离单元上的数据进行分块白化杂波抑制，则海面小目标检测问题可以转化为杂波白化条件下的目标检测问题：

(5)

3 双特征融合检测方法

3.1 时频脊-Radon变换提取特征

利用Cohen类时频分布方法对回波数据进行分析和特征提取。虽然Wigner-Ville分布具有良好的时频分辨率，但由于其具有非线性特性，在实际应用时会出现严重的交叉项干扰，因此本文采用平滑伪Wigner-Ville分布作为时频分析工具进行处理。分块白化后回波信号的SPWVD离散形式为：

(6)

定义时频脊序列Ridges(n)为DSPWVD(t,f)沿时间维的极大值集合，表示为：

Ridges(n)=max(DSPWVD(t,f)|t=t(n))，n=1,2,…,N

(7)

(8)

取平均信杂比取3 dB，速度取1.5 m/s。图1显示了白化前后目标单元与杂波单元的典型时频分布图像，其中相干脉冲数为64，频率平滑窗和时间平滑窗分别采用长度5点和31点的Hamming窗，其中时频脊序列用绿点标识。

由图1发现，在当前目标运动状态下，目标单元回波在分块白化滤波后，时频脊线向频率轴两端扩散程度较轻，基本呈线型；海杂波单元的时频脊线经分块白化滤波后，在整个平面的分布无明显规律。经大量实测数据统计分析可知，白化滤波后目标单元和杂波单元在时频脊分布特性上差异较大，因此目标检测问题可以转化为时频平面上检测时频脊线性程度的问题，本文采用时频脊-Radon变换方法提取特征以对其线性程度进行定量描述，该变换表示为：

图1 块白化前后杂波单元与目标单元的时频分布图Fig.1 The SPWVD of the sea clutter cell and target cell before and after block-whitening

RRidges(ρ,θ)=∬DRidgesDRidges(t,f)δ(ρ-tcosθ-fsinθ)dtdf

(9)

其中，DRidges(t,f)是从SPWVD中提取的时频脊图像，表示为：

(10)

式中，location(Ridges(n))表示时频脊序列Ridges(n)在时频平面中的(t,f)坐标位置集合。由于目标单元的时频脊具有良好的聚集性，近似呈水平线，故在实际检测中可以通过限制Radon的角度θ以达到减少计算量的目的。

图2给出了目标单元与海杂波单元的时频脊-Radon变换图像，其中，θ∈[80°,100°]。可以看出，目标单元的RRT峰值明显高于海杂波单元，这主要与目标单元时频聚集性较强有关。相比之下，海杂波单元的峰值较小。经大量实测数据验证，这种规律普遍存在。因此，本文将RRT峰值作为检测器的特征之一，定义为：

图2 杂波单元与目标单元的时频脊-Radon变换图像Fig.2 Images of sea clutter and target cell after RRT

ξMAX=maxRRidges(ρ,θ)}

(11)

式中，ξMAX代表峰值特征，RRidges(ρ,θ)代表时频脊-Radon变换后的特征空间。

显然，峰值特征体现了时频脊能量的集中性，它意味着时频脊序列Ridges(n)最多有ξMAX个脊点在同一条直线上。但该特征并不足以用于精确描述时频脊的线性程度，因为它的产生只利用了部分脊点，而未考虑其他脊点的分布特性。下面保持峰值点对应的θMAX角不变，以径向坐标ρMAX为中心，向两端扩散，记录每个径向坐标上存在的脊点数，可以得到脊点数随径向坐标ρ变化的散布图，如图3。可以看出，目标单元的时频脊散布性较杂波单元要弱，这从另一角度描述了时频脊能量的集中性。因此为了补足单峰值特征的缺陷，本文定义频宽特征ξBW用于衡量时频脊序列在径向坐标轴ρ上的散布性，如式(12)所示：

图3 杂波单元与目标单元的脊点散布图Fig.3 Ridge points distribution of sea clutter and target cell

(12)

式中，COEFF是为了避免个别脊点分布的偶然性而引入的控制量，称作频宽系数；P为时频脊点的总数目(数值上与每一帧的信号长度相等)；ρMAX和θMAX分别是峰值特征ξMAX在Radon空间中的坐标。

显然，与峰值特征ξMAX不同，频宽特征ξBW是从宏观角度对时频脊的集中性进行描述，它考察了所有脊点在径向坐标轴上的分布特性，这是只分析某一直线上脊点数的峰值特征无法做到的，而且与峰值特征ξMAX的表现相反，由于时频脊能量聚集性的差异，目标单元的ξBW特征值较杂波单元小。

3.2 帧间特征平滑方法

图4 待检测单元的时频特征直方图对比Fig.4 Time-frequency features’ histogram comparison of detection cell

为使特征检测器在短脉冲数条件下仍能具备较好的目标区分能力，这里通过多帧扫描历史数据和当前帧数据的综合应用，对RRT特征做帧平滑处理以增强可分性，即保留多帧扫描历史数据得到的特征值，并与当前帧特征值进行平滑处理，最终得到的特征提取结果表示为：

(13)

图5 帧平滑后待检测单元的时频特征直方图对比Fig.5 Time-frequency features’ histogram comparison of detection cell after frame smoothing

3.3 凸包算法检测分类

构造一个双特征平面，横轴代表RRT峰值特征，纵轴代表RRT频宽特征。由前文分析可知，当待测单元回波中包含目标时，由于在时频平面上良好的能量聚集性，其峰值特征较大，频宽特征较小，理应分布在特征平面的右下方；相对的，杂波单元的峰值特征较小，频宽特征较大，应分布在特征平面的左上方。基于此，本文根据设定的虚警率利用凸包学习算法[5,23]对特征点进行检测分类，流程如下：

1) 初始化：设置海杂波数据的特征训练集合为H，个数为W，计算虚警数为L=W·pF，其中pF为设定的虚警率。令l=0。

2) 计算当前数据点的凸包CH(H)，其中CH(H)的顶点为{v1,v2,…,vr}。然后统计落入凸包CH(H)中的特征点数量，设为nall。

3) 设置一个循环变量q，令其从1到r，计算凸包CH(H-{vq})，即从集合H中去除顶点vq，然后计算处于新凸包CH(H-{vq})中的特征点数量，设为nq。

4) 比较集合{nall-n1,…,nall-nr}，去掉该集合中最大数值所对应的顶点vi。

5) 令H-{vi}=H,l+1=l。

6) 若l

综合以上流程，本文提出一种海面小目标的双特征检测器，检测器原理框图如图6所示。在检测时，先对当前环境下的海杂波数据进行仿真分析，一方面，通过分析形成给定虚警概率条件下的判决空间，另一方面，确定出频宽系数、平滑因子的适宜取值。然后，分别对待检测单元数据做分块白化滤波、时频脊-Radon变换和特征提取、帧间特征平滑处理，形成待检测单元的双特征空间，若落入判决空间内，判定为杂波单元，假设成立；否则为目标单元，假设成立。需要说明的是，当海杂波特征因海况或雷达视角的改变而发生明显变化时，需要重新采集海杂波数据以确定新的凸包判决区域。本文提出的检测器虽然历经的步骤较多，但主要的时间消耗依然在特征提取阶段，而由于雷达的相参积累时间较短，极大削减了算法运行所需的时间，因此该检测器可以满足雷达实时检测的要求，具备一定的实用性。

图6 双特征融合检测器框图Fig.6 Block diagram of double feature fusion detector

4 实测数据分析结果

4.1 低海况数据分析

该部分采用低海况数据对本文检测器性能进行分析，数据来源于某X波段对海雷达试验数据，相参体制，脉冲重复频率(Pulse Repetition Frequency, PRF)4 kHz，带宽10 MHz，HH极化，海况约为2级。合作目标为渔船，匀速运动且运动速度较慢，图7是回波数据的距离-时间图像。经估算，目标SCR约为5～6 dB。在后续分析时，假定每帧数据均包含64个脉冲，对应的相干处理时间为16 ms。

图7 某型X波段雷达雷达接收数据的距离-时间图像Fig.7 The range-time image of data received by an X-band radar

4.1.1检测器性能的影响因子分析

首先抽取该数据集中的部分海杂波数据添加仿真目标(以1.5 m/s的速度匀速运动)，分析虚警率为10-3时频宽系数对本文检测器性能的影响，如图8(a)。显然，频宽系数在0.75～0.85范围内时，检测器性能较好，图8(b)继续对该区间内的检测器性能作了详细分析，发现在该海况条件下较为合理的频宽系数区间为[0.79, 0.81]，故本文取COEFF=0.8。

事实上，经大量数据集验证，该频宽系数的取值具有较强的普适性。经理论分析可知，频宽特征是从宏观上对时频脊的集中性进行描述，其准确性依赖于所使用的脊点数，当COEFF值过小时，频宽特征的形成只使用了少部分脊点，而未考虑其他可能有效的脊点，必然导致其准确性不足，使目标单元与海杂波单元之间的可分性不够，检测器性能下滑，如图8(a)中COEFF∈[0.2,0.6]区间；反之，若COEFF值过大，目标单元的时频谱中也可能存在个别脊点严重偏离直线(θMAX,ρMAX)，同样会引起检测器性能下滑，如图8(a)中COEFF∈[0.9,1]。因此后文不再对频宽系数作分析，均取定COEFF=0.8。

图8 双特征检测器性能随频宽系数的变化曲线Fig.8 The performance curve of the double feature detector varies with the bandwidth coefficient

进一步对平滑因子的取值影响作分析。分别取仿真目标的速度为1.5 m/s、0 m/s，图9给出了平滑因子对检测性能的影响。显然，随着M值增大，检测器对海上目标的检测性能呈提升趋势，且在[5,30]区间内改观尤为明显。当M∈[30,50]时，检测器性能的提升程度非常小，说明帧间特征平滑方法对检测器性能的改善是有限的。考虑到在实际场景中目标还有可能出现跨距离单元走动现象，进而影响能量累积和特征可分性，因此在兼顾检测性能和处理复杂度的情况下，取平滑因子M=30。

图9 平滑因子对检测器性能的影响Fig.9 The effect of smoothing factor on detector performance

4.1.2检测器的检测结果及对比分析

确定好频宽系数、平滑因子的取值之后，利用前文所述方法提取出雷达数据的RRT峰值和频宽特征，并在双特征平面上利用凸包算法进行检测分类。图10分别展示了本文提出的检测器、时频三特征检测器[5]、频域CFAR检测器、分形检测器在虚警率为10-3时的检测结果图，可以看出，本文方法得到的检测结果较为连续、稳定，且虚警较少，在性能上优于已有方法。

图10 四类检测器的检测结果Fig.10 The detection results of four detectors

进一步对检测性能进行量化分析，结果如图11所示，表1展示了其中三种虚警率条件下的检测概率，虚警率最高者用粗体标识。显然，四类检测器性能均随着虚警率增加逐步提升；同时经过对比发现，本文提出的检测器性能在虚警率pF∈[10-3,0]时要优于其他三类检测器，且由表1可知，性能差距超过20%，这是本文综合运用雷达历史扫描数据的做法体现出的优势；时频三特征检测器与频域CFAR检测器性能在不同虚警率时各有优劣，而分形检测器在低虚警率检测概率接近于零，无法满足目标检测的需要。

表1 四类检测器对低海况数据的检测概率

图11 四类检测器的性能对比Fig.11 Performance comparison of four detectors

4.2 高海况数据分析

本节继续采用高海况数据对双特征检测器性能进行分析，数据来源于CSIR数据库，该数据库由Fynmeet雷达于2006年在南非西南海岸采集，该雷达为VV极化，相参体制。本文采用文件名为TFA17-007的数据进行分析，对应的雷达工作参数和环境参数如表2所示[24-26]，海况约为4级。图12给出了该组数据的距离-时间图像，其中目标在第24距离单元匀速运动。经估算，目标SCR约为3 dB。在后续分析时，同样假定每帧数据均包含64个脉冲，对应的相干处理时间为12.8 ms，虚警率为10-3。

图12 TFA17-007数据的距离-时间图像Fig.12 The range-time image of TFA17-007 dataset

表2 雷达性能及环境参数

与4.1.1节类似，首先抽取检测当时的海杂波

数据，添加匀速运动(1.5 m/s)的仿真目标进行分析，结果如表3、图13所示。显然，此时由于海况等级较高、海浪起伏周期长的影响，平滑因子对检测器性能的改善程度并不如图9明显，且存在平滑因子、信杂比增大时检测概率反而下降的异常现象。进一步分析可知，平滑因子M=50、60、70时，本文检测器的性能在信杂比scr∈[0,6]区间内相近，由此推断，当M>50时，平滑因子对检测器性能的改善是微弱的。而当0

表3 本文检测器的检测概率

图13 平滑因子对检测器性能的影响Fig.13 The effect of smoothing factor on detector performance

在图14中，依旧通过量化分析给出了四类检测器的ROC曲线，且用表4展示三种虚警率条件下的检测概率。可以看出，在高海况条件下，时频三特征检测器性能不佳，说明时频特征容易受到海况条件的影响，且在低虚警率范围内该检测器与本文检测器的检测概率均无明显变化，这是因为部分属于目标单元的特征点落入凸包判决空间的位置较深，只有当虚警率较大时，凸包学习算法形成的判决空间较小，才能将它们正确分类。而频域CAFR检测器、分形检测器分别因为检测机理及判决方式的不同，其检测性能均随虚警率增加稳步上升，未出现此种现象，但在虚警率较低时，两者检测概率不足10%，无法用于检测目标。进一步分析表4可知，由于平滑因子的作用，本文检测器相对于其他三类检测器的性能改善程度在pF=10-3时将近60%，能够较好满足目标检测的需求。

图14 四类检测器的性能对比Fig.14 Performance comparison of four detectors

表4 四类检测器的检测概率

5 结论

针对相干脉冲数较少条件下特征检测器性能损失严重的问题，本文提出一种基于时频脊-Radon变换的帧平滑双特征检测方法。该方法采用分块白化滤波实现初步的海杂波抑制和SCR提升，以SPWVD为时频分析手段，对时频脊做Radon变换并提取出时频脊峰值和频宽两个特征。为使特征检测器在短脉冲数条件下仍能具备较好的目标区分能力，本文通过多帧扫描历史数据和当前帧数据的综合应用，对RRT特征做帧平滑处理以增强可分性。最后，在双特征平面上利用凸包算法形成判决空间并实现目标检测。经2级、4级海况实测数据验证，当相干脉冲数仅有64个的条件下，本文检测器的检测性能明显优于频域CFAR检测器以及已有的两类特征检测器(分别是时频三特征检测器和分形检测器)。后续还需进一步深入研究不同环境因素和不同目标运动状态下双特征检测器参数优化选取问题。