一种基于奇异谱分解的海杂波压制及目标检测方法与流程

本发明属于目标检测技术领域，涉及海杂波背景下的目标检测，为一种基于奇异谱分解的海杂波压制及目标检测方法。

背景技术：

海杂波是指对海观测雷达波束照射到海面时所产生的回波，海杂波的特性与雷达波束覆盖的区域面积、海浪高度(即海况)、雷达工作主频及带宽等因素密切相关，通常海杂波反映了由重力引起的排浪的运动特性和海面张力引起的管状细碎波的微动特性，即海杂波包含2个主分量，海杂波的幅度、带宽等参数随海况的变化呈现出非平稳、非线性特性，如附图1中(a)和(b)所示，为利用加拿大mcmaster大学ipix雷达测量的海杂波数据(http://soma.ece.mcmaster.ca/ipix/dartmouth/datasets.html)的时域幅度图，存在明显的尖峰，且尖峰的强度与数量随海况的增高而增加，附图2为上述海杂波数据的频谱图，如其中(a)和(b)所示，可见高海况时由于海面张力波引起的微动特性加剧而导致其带宽显著展宽。对海杂波背景下弱小目标检测方法的研究一直是雷达信号处理领域的热点，在民用领域内被用于检测海面漂浮物(浮冰、小型船舶和飞行器残骸等)，为民用船舶导航和海上搜救等提供技术支持。在军事领域中，该方法被用于探测海上具有隐身特性的舰船、潜艇潜望镜和贴近海面飞行的飞行器等，传统的方法是将海杂波视为一种复杂的平稳随机过程，通过大量的观测数据建立统计学意义上的概率模型，如weibull分布，log-normal分布、复合k分布模型和pareto分布等，进而利用成熟的检测方法实现目标检测。然而，对海杂波平稳性的假设是建立在短观测时间基础上的，为提高弱小目标的检测概率，雷达不得不增加对目标区域的观测时间，以提高经过相干积累后目标回波能量，但是随观测时间的增长海杂波将不再是平稳随机过程，如附图3中(a)和(b)所示，为南非csirtfc15_023.01海杂波数据集的pareto分布参数，其形状参数k和尺度参数sigma均随距离单元的不同而随机波动，若采用固定模型参数必然会造成较大的测量误差。

综上，对海观测雷达为探测海上弱小目标，通常以增加对海面目标区域的观测时间来增强目标的回波能量，然而在大相干处理时间内所接收的海面回波(海杂波)强度也相应增强，且呈现出非线性、非高斯和非平稳特性，传统的海杂波背景下的弱小目标检测方法均是将海杂波视为服从某种概率密度分布，以有无目标时的概率密度分布特征之间的差异检测目标，传统方法忽略了海杂波的非平稳特性，在高海况时虚警概率过高，不适用于大相干处理时间的目标检测。

为解决上述问题，雷达工程领域内的专家学者进行了大量的研究工作，已有一定数量的研究成果。经典的解决方案之一是对雷达回波的统计学特性进行动态跟踪，按特性相近的原则对其进行分段处理，在工程层面上解决了部分问题，但随着海况的提高，雷达回波数据的分段越来越短，最终与短时观测的情形无本质的差异，致使对弱小目标探测时仍存在较高的虚警概率、较低的发现概率，因此，基于统计学特性的信号处理方法在高海况条件下已不满足工程的需要，必须寻找不依赖统计学模型的海杂波抑制方法，进而解决这一难题。

技术实现要素：

为了克服上述现有技术的缺点，本发明的目的在于提供一种基于奇异谱分解的海杂波压制及目标检测方法，首先计算回波数据的自相关函数，以自相关函数的第一过零点位置作为重构雷达回波数据toeplitz矩阵和轨迹矩阵的维数，进而求得主分量(principalcomponent，pca)数量，依据pca数量判断回波中是否包含目标；而后利用奇异谱分解对雷达回波进行自适应滤波，将处理后的雷达回波进行恒虚警(cfar)目标检测。

为了实现上述目的，本发明采用的技术方案是：

一种基于奇异谱分解的海杂波压制及目标检测方法，包括如下步骤：

s1)对于相干雷达回波，按顺序执行s2)，对于非相干雷达回波，转s3)；

s2)求取雷达回波i、q数据序列的绝对值；

s3)计算回波数据的自相关函数值；

s4)计算回波数据自相关函数的第一过零点位置，估计回波数据的相关时间，作为构造回波数据toeplitz矩阵c的维数和后续奇异谱分解的窗口长度l；

s5)按窗口长度l将雷达回波数据转换成toeplitz矩阵c；

s6)对toeplitz矩阵c进行特征值分解，提取特征值σi和对应的特征向量vi，i＝1,2,…l，计算奇异谱依次累加ri，当累加值大于或等于0.9时对应的ri的个数np为主分量σip的数量；定义σip为雷达回波数据的主分量，vip为其特征向量，ip＝1,2,…,np；定义σjm为雷达回波数据的次分量，vjm为其特征向量，jm＝1,2,…,l-np；

在主分量σip、次分量σjm确定后计算：

s6.1:信杂比scr(signaltoclutterratio)

s6.2:信噪比snr(signaltonoiseratio)

s6.3:若回波中仅包含海杂波，则计算杂噪比cnr(cluttertonoiseratio)

s7)对雷达回波实施奇异谱分解，包括以下步骤：

s7.1：将s6)所提取的特征向量按所对应特征值由大到小重新排列v＝(vip,vjm)；根据奇异谱的累加结果，区分了主分量(principalcomponentanalysis,pca)和次分量(minorcomponentanalysis,mca)后，其顺序分别都是由大到小排列的，此处的操作是要相应地调整特征向量的排序，即按ip和jm的顺序排列，以对应pca和mca的位置；

s7.2:将雷达回波x＝{xk}按延迟坐标相空间重构(delaycoordinatephasespacereconstruction)方法转换为2维轨迹矩阵h，k＝1,2,…n，n为雷达回波数据长度；相应的嵌入维(embeddingdimension)m＝l，延迟时间(timedelay)τ＝1；

其中矩阵h的行数m＝n-(m-1)τ。

s7.3:计算2维轨迹矩阵h在特征向量张成的特征空间上的投影w＝hv，按照主、次分量对应的特征向量，w分解成主分量对应的np个子矩阵wip和次分量对应的(l-np)个子矩阵wjm二组矩阵；

s7.4:分别对np个子矩阵wip和(l-np)个子矩阵wjm实施逆对角线元素平均处理(inversehankelizationprocess)将其转换成一维时间序列xp＝{xip}，和xm＝{xjm}，获得主、次分量对应的时间序列；

s8)若雷达回波的主分量个数np>et，保留最大主分量σip对应的时间序列xip，此时ip＝1，舍弃其他分量对应的时间序列；若np≤et，则保留xm＝{xjm}中最后一个时间序列，即最小次分量对应的时间序列，舍弃其他分量对应的时间序列；其中et为判断雷达回波是否包含目标的门限；

优选地，所述et＝2，为判断雷达回波是否包含目标的门限，即海杂波所包含的主分量个数不大于2，2个主分量分别反映海浪的重力波和海面张力波强度。

s9)依次对相干处理时间内所收到的雷达回波进行s2)～s8)处理，对所保留的时间序列按照回波脉冲到达时间排列成2维回波矩阵，而后利用恒虚警检测(cfar)方法进行目标检测；整个雷达回波处理过程定义为ssa-cfar检测。

与现有技术相比，本发明的有益效果是：

本发明是一种不依赖海杂波统计学模型的海杂波压制方法，在给定海杂波背景下，可以检测到传统方法不能发现的弱小目标，在检测精度上有很明显的提高。

附图说明

图1为ipix雷达实测海杂波，其中，(a)为低海况海杂波，(b)为高海况海杂波。

图2为ipix雷达实测海杂波频谱，其中，(a)为低海况海杂波频谱，(b)为高海况海杂波频谱。

图3为南非csirtfc15_023.01海杂波数据集的pareto分布参数，其中，(a)为形状参数k，(b)为尺度参数sigma。

图4为本发明流程图。

图5为经ssa处理后距离单元11内的目标船回波。

图6为经ssa处理后距离单元11内的海杂波，其中，(a)为主分量σ2对应的时间序列，(b)为σ3对应的时间序列。

具体实施方式

下面结合附图和实施例详细说明本发明的实施方式。

本发明为一种基于奇异谱分解(singularspectrumanalysis,ssa)的海杂波压制及目标检测方法，首先计算回波数据的自相关函数，以自相关函数的第一过零点位置作为重构雷达回波数据轨迹矩阵的维数，构建回波数据的toeplitz矩阵、求得主分量(principalcomponentanalysis,pca)数量,依据pca数量判断回波中是否包含目标；而后对回波轨迹矩阵进行奇异谱分解，重构主分量和次分量(minorcomponentanalysis,mca)对应的时间序列，当回波中包含目标时，取最大主分量对应的时间序列、舍去其他分量对应的时间序列；反之，则选择最小次分量对应的时间序列而舍去其他分量对应的时间序列，这一处理过程可被视为不依赖回波统计学模型的自适应滤波，经处理后的雷达回波数据中海杂波强度被有效压制，有利于提高后续弱小目标的检测概率。

本发明的具体步骤如图4所示，包括：

s1)对于相干雷达回波{yi＝vi+jui}，按顺序执行s2)，对于非相干雷达回波{xi}，转s3)，i＝1,2,…n，vi、ui分别为雷达回波数据的实部i序列和虚部q序列；

s2)求取雷达回波i、q数据序列的绝对值；

s3)计算回波数据的自相关函数值:

s4)计算回波数据自相关函数的第一过零点位置，作为构造回波数据toeplitz矩阵c的维数和后续奇异谱分解的窗口长度l；

s5)按窗口长度l将雷达回波数据转换成toeplitz矩阵c；

s6)对toeplitz矩阵c进行特征值分解，即[v,d]＝eig(c)(matlab函数)，从d的对角线提取特征值σi和对应的特征向量vi，i＝1,2,…l，计算奇异谱依次累加ri，当累加值大于或等于0.9时对应的ri的个数np为主分量σip的数量；定义σip为雷达回波数据的主分量，vip为其特征向量，ip＝1,2,…,np；定义σjm为雷达回波数据的次分量，vjm为其特征向量，jm＝1,2,…,l-np；依据海杂波的主分量分析得到的结论，令et＝2为判断雷达回波是否包含目标的门限，即海杂波所包含的主分量个数不大于2，2个主分量分别反映海浪的重力波和海面张力波强度。在主分量σip、次分量σjm确定后可计算：

s6.1:信杂比scr(signaltoclutterratio)

s6.2:信噪比snr(signaltonoiseratio)

s6.3:若回波中仅包含海杂波，则计算杂噪比cnr(cluttertonoiseratio)

s7)对雷达回波实施奇异谱分解，包括以下步骤：

s7.1：将s6)所提取的特征向量按所对应特征值由大到小重新排列v＝(vip,vjm)；根据奇异谱的累加结果，区分了主分量(principalcomponentanalysis,pca)和次分量(minorcomponentanalysis,mca)后，其顺序分别都是由大到小排列的，此处的操作是要相应地调整特征向量的排序，即按ip和jm的顺序排列，以对应pca和mca的位置；

s7.2:将雷达回波x＝{xk}按延迟坐标相空间重构(delaycoordinatephasespacereconstruction)方法转换为2维轨迹矩阵h，k＝1,2,…n，n为雷达回波数据长度；相应的嵌入维(embeddingdimension)m＝l，延迟时间(timedelay)τ＝1；

其中m表示矩阵h的行数，m＝n-(m-1)τ。

s7.3:计算2维轨迹矩阵h在特征向量张成的特征空间上的投影w＝hv，按照主、次分量对应的特征向量，w分解成主分量对应的np个子矩阵wip和次分量对应的(l-np)个子矩阵wjm二组矩阵；

s7.4:分别对np个子矩阵wip和(l-np)个子矩阵wjm实施逆对角线元素平均处理(inversehankelizationprocess)将其转换成一维时间序列xp＝{xip}，和xm＝{xjm}，获得主、次分量对应的时间序列；

s8)若雷达回波的主分量个数np>et，保留最大主分量σip对应的时间序列xip，此时ip＝1，舍弃其他分量对应的时间序列；若np≤et，则保留xm＝{xjm}中最后一个时间序列，即最小次分量对应的时间序列，舍弃其他分量对应的时间序列；其中et为判断雷达回波是否包含目标的门限；

如附图5所示，为南非csirtfc15_023.01海杂波数据集中距离单元11经ssa处理后内的目标船回波；附图6中为该距离单元分离出的海杂波，其中，(a)为主分量σ2对应的时间序列，(b)为σ3对应的时间序列。

s9)依次对相干处理时间内所收到的雷达回波进行s2)～s8)处理，对所保留的时间序列按照回波脉冲到达时间排列成2维回波矩阵，而后利用恒虚警检测(cfar)方法进行目标检测；整个雷达回波处理过程定义为ssa-cfar检测。

以南非csirtfc15_023.01海杂波数据集为处理对象，该数据集包含31个距离单元的回波i-q数据,数据长度均为33001，距离单元11包含目标船回波，距离单元9～10、12～13为目标船与海面形成的干涉区，经s2)～s8)处理后，分别采用ca-cfar和os-cfar检测目标，训练距离单元数设为16,保护距离单元数设为2,目标虚警概率设为pfa_goal＝10^-3，目标检测扫描距离单元9～13，ca-cfar的发现概率pd＝0.9876,虚警概率pfa＝3.0302×10^-5；os-cfar的发现概率pd＝0.9907,虚警概率pfa＝1.2121×10^-4。

为验证本发明的先进性，直接对上述数据集进行目标检测，ca-cfar的发现概率pd＝0.1857,虚警概率pfa＝6.0604×10^-5；os-cfar的发现概率pd＝0.3268,虚警概率pfa＝1.2121×10^-4。