一种基于时间序列分割的多功能雷达状态序列聚类方法

本发明涉及雷达电子侦察技术领域，具体涉及一种基于时间序列分割的多功能雷达状态序列聚类方法。

背景技术：

多功能雷达是具有多种动态变化工作模式的复杂传感器，被广泛应用于监视,目标跟踪与识别等领域。多功能雷达可以同时执行多个不同的雷达工作模式，具有捷变的波束调度能力、复杂的信号调制样式、程控的工作模式编排。多功能雷达在时间线上编排多个不同的雷达工作模式，对每一个编排的雷达工作模式，可以基于对当前环境和目标的感知来选择或者优化雷达的控制参数(这些控制参数在侦察方又称为脉冲描述字，pulsedescriptivewords,pdw)，如脉冲重复间隔(pulserepetitioninterval,pri)，射频(radiofrequency,rf)和脉宽(pulsewidth,pw)，来充分挖掘多功能雷达的性能潜力。这些特点使得电子侦察系统接收的雷达脉冲流信号序列中，1)脉冲序列以数据流的形式输入至电子侦察系统；2)输入脉冲流中，包含的雷达工作模式数量未知；3)每个模式的脉内，脉间调制样式和对应调制参数可以不同；4)每个工作模式的持续时间也可能不同。多功能雷达的这些动态特性给传统的电子侦察与对抗系统带来了极大的挑战。基于对接收脉冲流数据的分析，实现对非合作多功能雷达工作模式的在线准确识别，是电子侦察研究领域中的一个热点和难点问题。

多功能雷达具有层次化结构，对多功能雷达工作模式的识别研究，可以从对其层次化结构的建模开始。现有的对多功能雷达工作模式识别方法首先是利用自然语言处理中的随机上下文无关文法来建模多功能雷达层次化结构，从句法分析理论的角度建模了多功能雷达的雷达资源管理功能，从而实现对信号产生流程的模拟。具体实现过程包括两个步骤。首先根据接收到的雷达脉冲序列进行雷达字提取，利用如隐马尔科夫模型或者模板匹配方法进行识别，得到脉冲序列对应的雷达字序列。然后针对雷达字序列，利用预测状态表示算法或者自动机将雷达字序列转换成对应的雷达工作模式。上述实现方法假设每个脉冲序列样本中仅存在一个雷达字类别或者一个雷达工作模式类别。

实际系统中，侦察接收机接收到的是连续不断的多功能雷达脉冲信号序列。层次化模型及对应的工作模式识别方法存在两个限制：1)层次化模型对多功能雷达的建模，需要层次结构中所有元素的先验信息以及这些元素之间的状态转移规则。这些信息在实际情况中通常是很难获取的。2)基于层次化模型的识别方法利用所有可用的先验信息得到各个工作模式的固定模板，然后基于这些固定的模板对各个雷达工作模式进行识别。这些固定模板无法反映当前的多功能雷达信号形式多样，灵活可变的特点。在雷达脉冲多维参数上不同的脉间调制类型，往往直接反映了雷达的动态特征和潜在功能意图。如，参差脉冲重复间隔往往在动目标检测雷达中用于强调盲速的影响。从而，一些研究开始着手于研究对多功能雷达脉冲序列脉间调制样式的识别。

早期对脉间调制识别的研究使用手工设计的特征和人工神经网络，随着深度学习网络如卷积神经网络，循环神经网络的发展，其自动提取特征的特点使得脉间调制识别的性能得到了很大的提升。然而，这些方法会遇到如下四个方面的困难：1)这些识别方法都属于监督学习，需要预先获取训练数据，然后在给定的类别中进行识别。在面临先进多功能雷达如软件自定义雷达或者认知雷达时，一个预先训练好的网络总是会遇到新的工作模式类别；2)对于一些脉间调制类型相同，但脉间调制参数存在区别的雷达工作模式，一个训练好的用于脉间调制类型识别的分类器往往不会进行如此细粒度的区分。因此，亟需研究无监督时间序列聚类方法来实现对多功能雷达工作模式(状态)序列的聚类。

技术实现要素：

本发明提出一种基于时间序列分割的多功能雷达工作模式序列聚类方法，定义各个雷达工作模式为不同pdw参数的不同脉间调制类型或者相同调制类型但调制参数不同，对截获的一段长度输入脉冲数据进行分割聚类，得到状态切换点，并给出每个状态各个脉冲的各个控制参数对应的脉间调制类型以及对应的脉间调制参数。

一种基于时间序列分割的多功能雷达状态序列聚类方法，包括如下步骤：

s1、接收到长度为t的脉冲序列p＝(p1,p2,…,pt)，每个脉冲由某一个pdw参数表征；

s2、设定算法参数窗长为w，分割点集合splitpoint＝[1,t]，以及设定阈值为ε，令k＝1，k为已有的片段数，然后令i表示已有片段的第i个片段，则，k＝1时，i＝1；k＝2时，i＝1，2，以此类推；另设i的初始值为i＝1；ts＝w；ts表示第i个片段中的索引t的初始值；

s3、如果k≤k,继续执行s4；k为设定的脉冲序列p的最大分割数；否则，执行s10；

s4、如果i＞k，执行s8；

如果i≤k，继续执行s5，则第i个片段为：

seqtosplit(i)＝p(splitpoint(i):splitpoint(i+1))

记len为seqtosplit(i)的长度，splitpoint(i)表示第i个片段的起始脉冲索引号；splitpoint(i+1)表示第i个片段的终止脉冲索引号；p(splitpoint(i):splitpoint(i+1))表示p中第splitpoint(i)个脉冲到第splitpoint(i+1)个脉冲；

s5、基于seqtosplit(i)估计参数化模型θall；

设t′＝len-w，如果t＝t′,则执行s7；

如果t＜t′，则执行s6；

s6、基于s5选择的模型θall，计算在seqtosplit(i)的第t个脉冲处进行分割后得到的指标增益improvement(t)，然后在seqtosplit(i)的t+1个脉冲处执行反复执行s6，直到第t′个脉冲处结束，由此得到了h个指标增益improvement(t)，其中，h＝t′-ts+1；

s7、在h个指标增益improvement(t)中，寻找最大增益，定义为segmentimprovement(i)，并记录最大增益对应的脉冲索引值，记为索引maxindex(i)；令i＝i+1，执行s4-s7；直到i＞k，由此得到k个segmentimprovement(i)；

s8、当k个最大增益segmentimprovement(i)中的最大值超过了阈值ε，则对应的索引maxindex(i)被添加为一个新的切割点，记为tnew，更新切割点集合：

splitpoint＝[splitpoint,tnew]

将splitpoint中的切割点按升序排列；

然后执行s9；

当k个最大增益segmentimprovement(i)中的最大值不超过阈值ε，则不添加新的切割点，执行s10；

s9、当更新切割点集合后，设当前已有的片段数为k′，k′＝k+1，反复执行s3-s9，直到k>k时结束循环；得到最终的切割点集合；

s10、利用以上步骤得到的切割点集合，对p进行切割，得到脉冲序列片段；对每个脉冲序列片段，估计得到对应的参数化模型，由此完成了基于时间序列分割的多功能雷达状态序列聚类。

较佳的，所述s5中，估计参数化模型的方法为：

基于所有可能的脉间调制类型，计算每一个类型下脉冲序列seqtosplit对应的模型以及模型参数；最后分别计算该脉冲序列在各模型下的bic值，取bic值最小的模型作为该脉冲序列所属模型；

bic＝-2ln(l)+numparas×ln(w′)

这里，ln(l)为seqtosplit在估计得到的不同模型下的对数似然值；numparas为各个模型的参数数目，w′为seqtosplit脉冲数目。

较佳的，参数化模型包括高斯抖动调制类型下的参数化模型；其中，在估计参数模型时，高斯抖动调制类型的概率密度函数由截断高斯分布描述，即高斯分布在横坐标等于0处截断。

较佳的，参数化模型包括正向滑变调制类型下的参数化模型；其中，在估计参数模型时，当使用的脉冲序列属于多个滑变周期的脉冲时，排除从第二个周期往后的各个周期中初始脉冲的影响，即对使用的脉冲序列，首先计算一阶差分，然后将一阶差分中大于固定阈值的脉冲删除。

较佳的，参数化模型包括参差调制类型下的参数化模型；其中，参差调制对应一系列有序切换的离散pri值，离散pri值使用带高斯发射的隐马尔可夫模型hmm描述。

较佳的，参数化模型包括正弦调制类型下的参数化模型；其中，在估计参数模型时，用f来产生pri序列，其中：fc为载频，fs为对应的采样频率。

进一步的，还包括s11，具体为：

对s10中t个脉冲得到的各个相同调制类型的模型，基于估计得到的模型参数再次进行聚类，将聚成一类的模型对应的脉冲标记为同一个标签，对存在的所有调制类型的模型进行聚类后，就得到了最终的脉冲标签序列然后将这些聚成一类的多个模型合并，利用对应的脉冲重新估计模型参数

较佳的，指标增益improvement(t)的计算方式如下：

将t作为分割点，将seqtosplit划分为两个序列，从第1到t个脉冲的脉冲序列为seqleft,t+1到最后一个脉冲的脉冲序列为seqright；分别基于seqleft和seqright，根据s5中的模型估计方法，估计两个序列对应的参数化模型θleft和θright，然后，根据θleft和θright计算seqleft和seqright在这两个模型下的bic指标值bicleft和bicright，得到improvement(t)＝bicall-(bicleft+bicright)。

较佳的，对处理的t个脉冲，根据这t个脉冲基于一定的评价指标优化阈值，并将优化的阈值再次用于这t个脉冲的分割聚类，阈值的优化具体为：

当score为bic时，

bic＝-2ln(l)+numparas×ln(t)

其中，ln(l)为p在估计得到的模型θ和标签序列d下的对数似然值,k表示估计得到的模型数目，numparas表示模型θ包含的参数数目。

本发明的有益效果：

本发明提出了一种基于时间序列分割的多功能雷达状态序列聚类方法，可以对输入所含模式个数可变、模式长度可变、模式调制参数可变的复杂多功能雷达工作模式脉冲序列，实现各个pdw参数上，脉冲级别的无监督时间序列分割聚类。本方法利用各个pdw参数的脉间调制类型，构建相应的时间序列参数化模型并给出对应的模型参数估计方法。各个参数化模型能够表示不同脉间调制类型的时间序列特征，对应的模型参数估计方法也是目前估计方法中效果最优或者比较好的。本方法的聚类结果标注到脉冲级别，一方面能够将长脉冲序列中属于不同雷达工作模式的各个脉冲进行聚类，另一方面，聚类结果能够给出各个参数化模型的模型参数，给后续脉冲序列更精细地分析和处理提供了条件。综上，本发明给出的多功能雷达工作模式分割聚类方法，可以为后续的多功能雷达系统行为的辨识和推理提供技术手段支持。

附图说明

图1为本发明实现雷达工作模式分割聚类的功能图。

图2为实施例中得到的雷达工作模式分割聚类结果图。

具体实施方式

本发明提供了一种对多功能雷达工作模式的在线聚类方法。

一种对多功能雷达工作模式的在线聚类方法，包括如下步骤：

s1、接收到长度为t的脉冲序列p＝(p1，p2,…,pt)，每个脉冲由某一个pdw参数表征；

s2、设定算法参数窗长为w，分割点集合splitpoint＝[1，t]，以及设定阈值为ε，令k＝1，k为已有的片段数，然后令i表示已有片段的第i个片段，则，k＝1时，i＝1；k＝2时，i＝1，2，以此类推；另设i的初始值为i＝1；ts＝w；ts表示第i个片段中的索引t的初始值；

s3、如果k≤k,继续执行s4；k为设定的脉冲序列p的最大分割数；否则，执行s10；

s4、如果i＞k，执行s8；

如果i≤k，继续执行s5，则第i个片段为：

seqtosplit(i)＝p(splitpoint(i):splitpoint(i+1))

记len为seqtosplit(i)的长度，splitpoint(i)表示第i个片段的起始脉冲索引号；splitpoint(i+1)表示第i个片段的终止脉冲索引号；p(splitpoint(i):splitpoint(i+1))表示p中第splitpoint(i)个脉冲到第splitpoint(i+1)个脉冲；

s5、基于seqtosplit(i)估计参数化模型θall；

设t′＝len-w，如果t＝t′,则执行s7；

如果t＜t′，则执行s6；

估计参数化模型θall时，当脉冲序列所属的脉间调制样式未知，在创造模型并估计模型参数时，需要使用模型选择方法。以这里的四种脉间调制类型为例说明，对脉冲序列seqtosplit，需要构建四个相应的模型，然后利用seqtosplit估计模型参数，最后分别计算seqtosplit在四个估计得到的模型下的bic值，取bic值最小的模型作为创造的新模型。优化问题如下:

这里，ln(l)为seqtosplit在估计得到的不同模型下的对数似然值。numparas为各个模型的参数数目，w为seqtosplit脉冲数目。

以pri参数为例，给出四种典型的脉间调制类型对应的参数化模型及模型参数估计方法。

1)高斯抖动调制

高斯抖动调制基于均值μ和方差σ²产生对应的pri序列。由于pri序列中的每个pri取值都为正数，所以高斯抖动调制类型的概率密度函数(probabilitydensityfunction,pdf)可以由一个截断高斯分布描述。这个高斯分布在横坐标等于0处截断。从而，对一个给定的pri的值pt，对应的pdf为：

其中，φ是正态累积分布函数，当μ＞＞0且μ＞＞σ时，因此，当输入脉冲数据为pwi＝(pi,pi+1,…,pi+w-1)时，pwi对应的对数似然函数为：

通过该对数似然函数，可以求得给定pwi情况下参数μ和σ²的最大似然估计。如下：

2)正向滑变调制

在一个滑变周期内，脉冲pt可以表示为

pt＝α+pt-1+ωt

其中，α>0是滑变步长，是零均值，方差为σ²的高斯白噪声。其他滑变类型如负向滑变调制，三角滑变调制(如正向滑变调制后接负向滑变调制)等可以由正向滑变调制的参数化模型变换得到。

我们首先考虑，输入的pwi中不含属于多个滑变周期的脉冲，则给定pwi时，脉冲pt的条件概率分布为：

则pwi的联合概率密度分布为：

我们应该忽略f(p1；σ²)的影响，因此p1条件下pwi的对数似然函数为：

可以得到对应的滑变步长α和噪声方差σ²的估计值：

当pwi中包含属于多个滑变周期的脉冲时，我们还需要排除从第二个周期往后的各个周期初始值的影响。如p7是第二个周期的起始脉冲，p6是第一个周期最后一个值。在估计值公式中，我们需要排除p7-p6的影响，因为这一项远大于其他项，且和步进值没有直接关系。对pwi，我们可以首先计算pwi的一阶差分，然后在一阶差分中，利用固定阈值(在正向滑变中，可以用0，或者均值)排除这些项的影响。由于我们这里是在线聚类，窗长w不会太长，一个窗内包含来自两个周期脉冲的情况仅出现在窗位于上一个周期结束和下一个周期开始的部分。且窗内始终只有一个周期交替点。

3)参差调制

从雷达系统的角度，参差调制对应一系列有序切换的离散pri值，可以用隐马尔可夫模型(hiddenmarkovmodel,hmm)描述。但是从电子侦察系统的角度，接收到的pri序列是被噪声污染的，因此也就不能用hmm中的离散观测矩阵描述。我们这里使用带高斯发射的hmm描述，记为gaussianhmm。gaussianhmm可以由三元组描述：

θ＝<a,b,π>

其中a＝[aij]m×m为状态转移矩阵，m为状态数目。π＝π1,π2,…,πm为初始状态分布。记q＝(q1,q2,…,qm)为状态集合，每个状态对应一个均值为μm，方差为的高斯模型。这些高斯模型为b＝(φ1,φ2,…,φm)，其中

这些状态就对应了参差调制中各个参差点上的pri取值。则，对输入数据pwi而言，其对应的对数似然函数为：

其中

对于t＝1,2,…,n-1

通过期望最大算法进行模型参数的估计。在期望最大算法每轮迭代中θ的估计值为：

其中：

γt(i)＝f(it＝qi|pwi,θ)

ξt(i,j)＝f(it＝qi,it+1＝qj|pwi；θ)

4)正弦调制

正弦调制的pri序列可以表示为

pt＝asin(2πft+φ)+c+ωt

其中，从雷达的角度，雷达选择正弦调制的载频fc和对应的采样频率fs来产生正弦pri序列。然而，从侦察方的角度，fs是未知的，因此我们用f来表征fc和fs，用以产生pri序列。a为幅度，c为常数项，φ为相位。同时，为了满足pri为正值的需求，c>a。

从而，我们可以计算给定pwi情况下的对数似然函数：

其中，θ＝(a，f,c，σ²，φ)，θ的最大似然估计值为：

s6、基于s5选择的模型θall，计算在seqtosplit(i)的第t个脉冲处进行分割后得到的指标增益improvement(t)，然后在seqtosplit(i)的t+1个脉冲处执行反复执行s6，直到第t′个脉冲处结束，由此得到了h个指标增益improvement(t)，其中，h＝t′-ts+1；

其中，指标增益improvement(t)的计算方式如下：

将t作为分割点，将seqtosplit划分为两个序列，从第1到t个脉冲的脉冲序列为seqleft,t+1到最后一个脉冲的脉冲序列为seqright。分别基于seqleft和seqright，根据s5中的模型估计方法，估计两个序列对应的参数化模型θleft和θright。然后，根据θleft和θright计算seqleft和seqright在这两个模型下的bic指标值bicleft和bicright，得到improvement(t)＝bicall-(bicleft+bicright)。

s7、在h个指标增益improvement(t)中，寻找最大增益，定义为segmentimprovement(i)，并记录最大增益对应的脉冲索引值，记为索引maxindex(i)；令i＝i+1，执行s4-s7；直到i＞k，由此得到k个segmentimprovement(i)；

s8、当k个最大增益segmentimprovement(i)中的最大值超过了阈值ε，则对应的索引maxindex(i)被添加为一个新的切割点，记为tnew，更新切割点集合：

splitpoint＝[splitpoint，tnew]

将splitpoint中的切割点按升序排列；

然后执行s9；

当k个最大增益segmentimprovement(i)中的最大值不超过阈值ε，则不添加新的切割点，执行s10；

s9、当更新切割点集合后，设当前已有的片段数为k′，k′＝k+1，反复执行s3-s9，直到k>k时结束循环；得到最终的切割点集合；

s10、利用以上步骤得到的切割点集合，对p进行切割，得到脉冲序列片段；对每个脉冲序列片段，估计得到对应的参数化模型，由此完成了基于时间序列分割的多功能雷达状态序列聚类。

进一步的，还包括s11：每t个脉冲，对s10中对该段t个脉冲得到的各个相同调制类型的模型，可以基于估计得到的模型参数再次进行聚类(如基于距离的聚类)。将聚成一类的模型对应的脉冲标记为同一个标签。对存在的所有调制类型的模型进行聚类后，就得到了最终的脉冲标签序列然后将这些聚成一类的多个模型合并，利用对应的脉冲重新估计模型参数

实施例：

将pri脉间调制样式不同，或者相同脉间调制样式不同调制参数的脉冲片段定义为不同的雷达工作模式片段。本例中给出脉冲序列中各个工作模式片段均为正弦调制，但各个片段调制参数不同且未知的情况。这种情况中，我们产生一个待测样本，每个样本包含五个正弦片段，每个片段包含120个脉冲。每个待测样本对应各个工作模式片段设置为θp＝{(2，0.25，50，0.5²)，(3,0.2,40,0.5²),(4,0.15，30，0.5²)，(5,0.1,20,0.5²)}，其中第五个片段为第二个片段重复一次。

具体的分割聚类方法如下：

先接收到长度为t的脉冲序列p＝(p1，p2，…，pt)，其中pi为第i个脉冲，每个脉冲由pri参数表征。设置w＝20，高斯抖动的阈值ε＝10^-8,k＝10，t＝600。

然后按照步骤s2-s10执行，估计得到对应的参数化模型，该样本得到的聚类输出结果如图2所示。

综上所述，以上仅为本发明基于所选定的pri调制类型及调制参数定义的雷达工作模式的实施示例而已，并非用于限定本发明的保护范围。基于脉冲序列各个pdw参数和参数的脉间调制类型、调制类型参数定义多功能雷达工作模式、基于参数化模型表征雷达pdw参数的脉间调制类型和对应的调制参数、基于接收脉冲序列对雷达pdw参数的参数化模型进行模型选择及参数估计、基于参数化模型的雷达工作模式分割聚类方法是本发明的核心要点。凡在本发明上述设计原则和实现要点之内，经过相应修改、替换、改进等操作形成的多功能雷达工作模式在线聚类方法，均应包含在本发明的保护范围之内。