基于海杂波抑制的双特征海面漂浮小目标检测方法与流程

本发明属于信号处理技术领域，具体涉及一种目标检测方法，可用于对海面漂浮、慢速小目标的识别与跟踪。

背景技术：

海杂波是雷达接收到的从海表面反射回来的雷达回波，海面搜索雷达在对海探测时，特别是对海面上诸如小船、浮冰、漂浮物等漂浮小目标进行检测时，难免会受到海杂波的影响。海杂波的强度会随雷达参数、雷达照射方向、海况等的不同而发生改变。在高分辨海杂波背景下，杂波呈现出较强的非高斯特性，海尖峰的存在导致利用时域能量累积的目标检测方法出现大量虚警，因此，漂浮小目标很难通过传统能量累积的方法进行检测。

为解决这一难题，很多学者对此做出了大量的研究。随着海杂波统计模型的不断完善，许多自适应检测方法被提出，该类方法将海杂波建模为复合高斯模型，由于海杂波的空时非平稳特性，该类方法需要先对海杂波进行抑制，具有一定的局限性。文献“Hu,J.,Tung,W.W.and Gao,J.B.:Detection of low observable targets within sea clutter by structure function based multifractal analysis,IEEE Trans.Antennas Propag.,54(1):136-143,2006.”中提出基于海面分形特征的检测方法，可在观测时间较长时有效的检测目标，然而雷达通常无法对单一波位进行长时间的驻留观测，因此基于分形特征的检测器很难推广到实际应用中。

对海面漂浮小目标的检测，很多方法以海杂波满足某种统计模型为假设前提，然而现有统计模型很难描述海杂波的复杂特性，这导致检测结果具有一定的局限性；自适应检测方法在海况较复杂时，即当目标与杂波在多普勒域无法区分时，无法对海面漂浮或低速小目标进行检测；基于分形的目标检测方法在观测时间较长时能达到很好的检测结果，当观测时间缩短时，检测性能会有明显下降，无法满足对海搜索雷达的要求。

技术实现要素：

本发明的目的在于针对上述已有技术的不足，提出一种基于海杂波抑制的双特征海面漂浮小目标检测方法，以提高对海面漂浮小目标的检测性能，满足雷达对海搜索的要求。

为实现上述目的，本发明的技术方案包括如下：

(1)从回波数据中选取训练单元、参考单元和待检测单元：

利用雷达发射机对海面发送信号，利用雷达接收机接收由海面反射回来的回波数据，该回波数据分为纯杂波数据和包含目标的回波数据；

从纯杂波数据中选取部分距离单元作为一组训练单元，该训练单元时间序列z为：z＝[z(1),z(2),…,z(N)]，从训练单元周围选取Q个临近单元作为参考单元，该参考单元时间序列z_p为：z_p＝[z_p(1),z_p(2),…,z_p(N)]，p＝1,2,…,Q，Q为参考单元数，N为时间序列的长度；

从包含目标的回波数据中选取部分距离单元作为待检测单元T；

(2)对纯杂波回波数据的训练单元时间序列z和参考单元时间序z_p进行块白化，得到白化后的训练单元时间序列和白化后的参考单元时间序列

(3)利用白化后的训练单元时间序列和白化后的参考单元时间序列提取训练单元的相对时频分布的脊能量η₁(z)和训练单元的相对时频分布的脊全变差η₂(z)这两种特征，构造纯杂波数据的相对时频双特征向量η：η＝[η1(z),η₂(z)]^T，[·]^T表示对矩阵进行转置；

(4)利用纯杂波数据的相对时频双特征向量η，在二维特征空间中获得二维凸包并在给定虚警概率P_F下，利用贪婪凸包学习算法对二维凸包进行收缩，并将收缩后的凸包作为检测判决区域Ω；

(5)对待检测单元T提取相对时频分布的脊能量η₁(T)和相对时频分布的脊全变差η₂(T)这两种特征，构造待检测单元的相对时频双特征向量：η_T＝[η₁(T),η₂(T)]^T；

(6)根据构成检测判决区域Ω的凸包的特性和待检测单元的相对时频双特征向量η_T＝[η₁(T),η₂(T)]^T，计算检测统计量ω：

其中，min{·}表示取最小值，det(·)表示求矩阵行列式，r为组成凸包的顶点的个数，x_j表示第j个凸包顶点的相对时频分布的脊能量，y_j表示第j个凸包顶点的相对时频分布的脊全变差，j＝1,2,...,r；

(7)根据检测统计量ω的大小判断目标是否存在：如果检测统计量ω大于零，表明待检测单元的相对时频双特征向量η_T在检测判决区域Ω外，则判定目标存在，反之，判定目标不存在。

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

1)本发明从海杂波序列中提取出两种不同的时频特征，并联合利用这两种时频特征对纯海杂波数据与含有目标的回波数据的区分能力完成对海面目标的检测，相比利用单一特征的传统检测方法，可以在较短观测时间内获得更好的检测效果；

2)本发明利用了贪婪凸包学习算法，由于计算量小，可以迅速的得到检测判决区域，提高了检测器的训练速度，更适合运用于实际应用中。

3)本发明利用块白化来抑制杂波，同时增强目标回波，利用平均散斑协方差矩阵来代替直接估计得到的协方差矩阵，可以有效减小白化过程中目标能量的分块效应。

附图说明

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

图2是在观测时间为512ms时，用本发明与现有两种检测方法在四种极化下的检测性能比较图。

图3是在观测时间为1024ms时，用本发明与现有两种检测方法在四种极化下的检测性能比较图。

具体实施方式

参照图1，本发明包括训练和检测两部分，其具体步骤如下：

步骤1，获取回波数据，并从回波数据中选取训练单元、参考单元和待检测单元。

利用雷达发射机对海面发送信号，利用雷达接收机接收由海面反射回来的回波数据，该回波数据分为纯杂波数据和包含目标的回波数据；

从纯杂波数据中选取部分距离单元作为一组训练单元，该训练单元时间序列z为：z＝[z(1),z(2),…,z(N)]，从训练单元周围选取Q个临近单元作为参考单元，该参考单元时间序列z_p为：z_p＝[z_p(1),z_p(2),…,z_p(N)]，p＝1,2,…,Q，Q为参考单元数，N为时间序列的长度；

从包含目标的回波数据中选取部分距离单元作为待检测单元T。

步骤2，对纯杂波回波数据的训练单元时间序列z和参考单元时间序z_p进行块白化。

块白化是指：将时间序列分为互不重叠的短向量块，利用参考单元时间序列估计出的每一短向量块的散斑协方差矩阵，对每个短向量块进行白化，其步骤如下：

(2a)将训练单元时间序列z和参考单元时间序列z_p分别平均分成长度为M的互不重叠的短向量，即：

z＝[z₁,z₂,…,z_m,…,z_N/M]，

z_p＝[z_p,1,z_p,2,…,z_p,m,…,z_p,N/M]，p＝1,2,…,Q

其中，z_m表示训练单元时间序列的第m个短向量，z_p,m表示参考单元时间序列的第m个短向量，m＝1,2,…,N/M；

(2b)利用上述每个短向量，对训练单元时间序列z和参考单元时间序z_p进行块白化，得到白化后的训练单元时间序列和白化后的参考单元时间序列

其中，表示白化后的训练单元时间序列的第m个短向量，表示白化后的参考单元时间序列的第m个短向量，表示参考单元时间序列的第m个短向量z_p,m的散斑协方差矩阵估计；

所述的对参考单元时间序列的第m个短向量z_p,m的散斑协方差矩阵进行估计，可通过现有多种方法实现，例如最大似然估计法、归一化样本协方差矩阵估计法和近似最大似然估计法，本实例采用近似最大似然估计法得到z_p,m的散斑协方差矩阵估计

步骤3，利用上述白化后的训练单元时间序列和白化后的参考单元时间序列提取训练单元的相对时频分布的脊能量η₁(z)和训练单元的相对时频分布的脊全变差η₂(z)这两种特征，构造纯杂波数据的相对时频双特征向量η。

(3a)计算训练单元的平滑维格纳-威利分布和参考单元的平滑维格纳-威利分布

其中，上标*表示共轭，g(m)为时间平滑窗，h(k)为频率平滑窗，E表示时间平滑窗长的一半，F表示频率平滑窗长的一半，Δf_d为归一化多普勒频率的采样间隔，n＝1,2,…,N，l＝1,2,…,N；

这里的时间平滑窗和频率平滑窗采用但不限于汉宁窗，海明窗，布莱克曼窗，凯瑟窗，本实例采用长度为31的凯瑟窗作为时间平滑窗，长度为63的凯瑟窗作为频率平滑窗；

(3b)利用参考单元的平滑维格纳-威利分布估计训练单元的平均时频分布A(n,l)：

其中，n′的取值为-1，0，1，l′的取值为-1，0，1；

(3c)利用训练单元的平滑维格纳-威利分布和训练单元的平均时频分布A(n,l)，计算训练单元的相对时频分布

(3d)计算训练单元的相对时频分布的时频脊R(n)：

其中，表示取最大值时对应l的取值；

(3e)利用训练单元的相对时频分布和训练单元的相对时频分布的时频脊R(n)，计算训练单元的相对时频分布的脊能量η₁(z)和训练单元的相对时频分布的脊全变差η₂(z)：

(3f)利用训练单元的相对时频分布脊能量η₁(z)和训练单元的相对时频分布脊全变差η₂(z)，构造纯杂波数据的相对时频双特征向量η：

η＝[η₁(z),η₂(z)]^T，

其中，[·]^T表示对矩阵进行转置。

步骤4，利用纯杂波数据的相对时频双特征向量η，通过凸包学习算法确定检测判决区域Ω。

现有的凸包学习算法有：快速凸包学习算法，贪婪凸包学习算法，Graham扫描法，礼品包裹法，本实例采用贪婪凸包学习算法确定检测判决区域Ω，其步骤如下：

(4a)利用纯杂波数据的相对时频双特征向量η，组成训练样本集合S：

其中，ηi表示第i个训练样本的相对时频双特征向量，表示第i个训练样本的相对时频分布的脊能量，表示第i个训练样本的相对时频分布的脊全变差，I为训练样本数目，本实例取训练样本数为20000个；

(4b)利用训练样本集合S，得到一个属于特征平面的集合Ψ：

其中，表示训练样本的相对时频分布的脊能量的平均值，表示训练样本的相对时频分布的脊全变差的平均值；

(4c)计算训练样本集合S组成的凸包

其中，polygon{·}表示由·组成的凸多边形，v_j为组成凸包的第j个顶点，j＝1,2,…,r,r为组成凸包的顶点的个数；

(4d)计算在凸包中去掉一个凸包顶点v_j后的凸包面积减小量Δ(v_j)：

其中，area表示凸包的面积，∩表示集合的取交运算，表示从凸包中删除凸包顶点v_j后的集合；

(4e)找出使凸包面积减少量Δ(v_j)最大的顶点v_*：

其中表示Δ(v_j)取最大值时对应j的取值；

(4f)将上述顶点v_*从训练样本集合S中去除，得到新的训练样本集合S'＝S-{v_*}；

(4g)重复步骤4a)至4f)，总共去除C＝[I×P_F]个相应的凸包顶点，最终得到收缩后的凸包即为满足虚警概率P_F的检测判决区域Ω，其中，[I×P_F]表示取I×P_F的整数部分。

步骤5，提取待检测单元T的时频双特征向量η_T。

在训练环节中，利用贪婪凸包学习算法得到判决区域Ω后，即可进行检测，检测时需先对待检测单元T提取相对时频分布的脊能量η₁(T)和相对时频分布的脊全变差η₂(T)这两种特征，构造待检测单元的相对时频双特征向量η_T：

5a)待检测单元T的时间序列为：T＝[T(1),T(2),…,T(N)]，在待检测单元T的周围选取Q个临近单元作为参考单元，该参考单元时间序列为：T_p＝[T_p(1),T_p(2),…,T_p(N)]，p＝1,2,…,Q；N为时间序列的长度；

5b)对待检测单元T的时间序列和参考单元时间序T_p进行块白化，得到白化后的待检测单元时间序列和白化后的参考单元时间序列；

5c)利用白化后的待检测单元时间序列和白化后的参考单元时间序列，提取待检测单元的相对时频分布的脊能量η₁(T)和待检测单元的相对时频分布的脊全变差η₂(T)这两种特征，得到待检测单元T的相对时频双特征向量：η_T＝[η₁(T),η₂(T)]^T。

步骤6，根据构成检测判决区域Ω的凸包的特性和待检测单元的相对时频双特征向量η_T计算检测统计量ω。

(6a)利用待检测单元的相对时频双特征向量η_T和凸包顶点计算矩阵行列式：

其中，det(·)表示求矩阵行列式，x_j表示第j个凸包顶点的相对时频分布的脊能量，y_j表示第j个凸包顶点的相对时频分布的脊全变差，j＝1,2,...,r；

(6b)取上述矩阵行列式的最小值，得到检测统计量ω：

其中，min{·}表示取最小值。

步骤7，根据检测统计量ω的大小判断目标是否存在：如果检测统计量ω大于零，表明待检测单元的时频双特征向量η_T在检测判决区域Ω外，则判定目标存在，反之，判定目标不存在。

基于步骤1到步骤7，实现了基于杂波抑制的双特征海面漂浮小目标的检测。

一.实验数据

本实例所用数据为12组IPIX雷达获取的实测海杂波数据，雷达架设高度为30m，脉冲重复频率为1000Hz，距离分辨率为30m；每组数据包含四种极化数据，其中两种为同极化数据HH、VV，两种为交叉极化数据HV、VH。其中有10组数据是93年采集的海杂波数据，每种极化数据包括14个距离单元，数据长度为2¹⁷，目标是直径为1米的圆球，表面用金属丝包裹；剩余2组数据是98年采集的海杂波数据，每种极化数据包括28个距离单元，数据长度为60000，目标是一艘漂浮的小船。

二.仿真实验

仿真1，在观测时间为512ms时，利用本发明与基于分形的检测方法和基于三特征的检测方法，在四种极化数据下对雷达检测性能进行仿真对比，结果如图2所示。其中图2(a)为同向HH极化数据下的雷达检测性能比较图；图2(b)为同向VV极化数据下的雷达检测性能比较图；图2(c)为异向HV极化数据下的雷达检测性能比较图；图2(d)为异向VH极化数据下的雷达检测性能比较图；

从图2可以看出，本发明对海面漂浮小目标的检测性能优于现有两种检测方法的检测性能。

仿真2，在观测时间为1024ms时，利用本发明与基于分形的检测方法和基于三特征的检测方法，在四种极化数据下对雷达检测性能进行仿真对比，结果如图3所示。其中图3(a)为同向HH极化数据下的雷达检测性能比较图；图3(b)为同向VV极化数据下的雷达检测性能比较图；图3(c)为异向HV极化数据下的雷达检测性能比较图；图3(d)为异向VH极化数据下的雷达检测性能比较图；

从图3可以看出，本发明与现有两种检测方法相比，对海面漂浮小目标的检测性能更好。