阳 榴, 朱卫纲, 吕守业, 马 爽
(1. 航天工程大学电子与光学工程系, 北京 101416; 2. 北京遥感信息研究所, 北京 100192)
随着相控阵技术的发展,多功能雷达(multi-function radar,MFR)广泛部署于各军事系统中,战场电磁环境日益复杂。MFR的行为分析与辨识[1-2]成为电子侦察领域近年来的研究热点。其中,波形单元作为MFR信号的基本构成,其提取效果直接影响到后续MFR“雷达字”集构建、工作模式识别、威胁等级判定等环节的准确性,具有重要的研究意义。
当前MFR波形单元提取的方法主要有以下两类。一类是有监督学习框架下的匹配分类技术,Visnevski等[3]在波形库已知的条件下,对MFR不同“雷达字”分别建立隐马尔可夫链模板,通过计算与模板的匹配程度提取波形单元。刘海军[4]利用三级匹配的方法进一步提高了匹配精度。文献[5]将降噪编码器与分类算法结合获取了较好的鲁棒性。文献[6]利用已知标签的训练数据对循环神经网络(recurrent neural network, RNN)模型调参,再通过RNN分类提取波形单元。另一类方法是在缺乏先验信息的条件下,利用不同波形单元间脉冲列的参数变化进行提取。文献[7]通过检测单一脉幅序列的多个统计变化点,实现了波形库未知的MFR波形单元提取。文献[8]用定长的滑窗遍历MFR脉冲描述字序列,寻找滑窗数据的变化来划分脉冲列,但算法中如何选取滑窗长度及步长缺乏相应的参考依据。除此之外,文献[2]通过数值点上的聚类提取MFR波形单元,但难以应用于采用随机调制的雷达。文献[9]结合脉冲到达时间和脉宽,提出一种事件驱动的波形单元提取方法,对虚假脉冲和漏脉冲有较好的适应性。
上述研究取得了一定的成果,但在实际应用中探测目标多为非协作MFR,事先难以获取雷达的波形单元集合。因此,在缺乏匹配标准或标签数据的条件下,基于有监督学习的匹配分类算法将难以适用。另一方面,寻找序列变化点的方法需要多次重复遍历、统计、判决的步骤,计算复杂度随着脉冲个数的增加迅速上升。同时,变化点检测依据的是序列的局部变化规律,没有充分利用整个序列的数据分布特点,加上目前较多算法只对单一参数序列进行统计判决,因此在具有测量误差和脉冲丢失的实际应用场景中提取效果往往不太理想。
针对以上问题,首先构建了一种基于MFR多参数序列的二分类模型,将多个变化点检测问题转化成了一个二分类问题。在此基础上针对实际应用中非协作MFR缺乏先验知识的问题,提出了一种参数自适应的密度聚类算法(density adaptive density-based spatial clustering of applications with noise,DA-DBSCAN)进行分类,能够有效地提取非协作MFR的波形单元。
多功能雷达信号具有功能(function,F)、任务(progress,P)、波形(waveform,W)的层级结构模型[7],如图1所示。
图1 MFR信号的层级模型Fig.1 Hierarchical model of MFR signal
MFR通过执行一系列雷达任务P来实现搜索、目标跟踪、制导等功能F。执行任务时,MFR从波形库[10]中选择合适的波形W在指定方向上发射出去,从而在波束驻留时间内对目标完成相应的动作。MFR波形单元提取则是指从侦收的脉冲列中获取上述波形单元的过程。
由于任务需求和波位指向的差异,不同的波形单元在射频(radio frequency,RF)、脉宽(pulse width,PW)、脉冲重复间隔(pulse repetition interval,PRI)、脉幅(pulse amplitude,PA)等参数上将发生变化[11]。以往的研究大多采用脉幅作为划分雷达波形单元的依据[7]。但实际探测过程中伴随着杂波、干扰和噪声,在功率维度具有较强干扰性,脉幅出现较大畸变,单一的脉幅数据难以准确地反映波形单元的变化情况。
通过分析某探测系统实际获取的MFR信号,发现对于采用随机调制的新体制雷达,同一波形单元内其RF和PRI也存在一定程度的随机波动,具有更高的抗干扰性能。但由于各波形单元在特定波位指向上完成不同的任务,波形单元转换时其参数波动更为剧烈。因此,形成了整体震荡剧烈、局部呈“簇”状的时序特点。如图2所示,根据某探测系统对某二维相控阵雷达的实际获取数据特点,仿真得到其脉冲列时序图。其中,PRI、RF、PW时序图整体震荡剧烈。对范围在500~750 μs 的PRI时序图放大,可以明显看到局部呈“簇”状的时序特点,如图2(b)所示。由于波形单元内部的随机调制,RF在遍历整个参数范围上取值,难以通过数值点上的聚类达到提取波形单元的目的,如图2(c)所示。基于以上分析,本文将利用PRI、RF、PW以及参数间的联合变化规律作为衡量波形单元转换的依据,采用数据驱动的思路,从数据集本身的特点出发构建波形单元提取的二分类模型。
图2 某二维相控阵雷达脉冲列的时序图Fig.2 Pulse sequencediagram of a two dimensional phased array radar
假设某多功能雷达侦收脉冲列P={p1,p2,…,pn}由三元时间序列{PRI,RF,PW}构成,包含n个脉冲和m个波形单元。其中:
{
PRI={x1,x2,…,xn}
RF={y1,y2,…,yn}
PW={z1,z2,…,zn}
(1)
令label={l1,l2,…,ln},l1=1,
(2)
其中,i=2,3,…,n。
若脉冲pi对应的li=1,则该脉冲为某波形单元的起始脉冲;若脉冲pi对应的li=0,则该脉冲为位于某波形单元内部的非起始脉冲。于是,将从脉冲列中寻找m个统计变化点的问题转化成一个二分类问题。
基于不同波形单元转换时参数波动大,整体震荡剧烈,波形单元内部参数波动较小,局部呈“簇”状的数据特点,利用相邻脉冲间参数的差值{ΔPRI,ΔRF,ΔPW}来衡量脉冲列时序变化的波动程度。其中:
(3)
图3 某二维相控阵雷达联合分布情况Fig.3 Distribution map of a two dimensional phased array radar
图3中,红点是label值为1的波形单元的起始脉冲,黑点是label值为0的波形单元内部脉冲。由于波形单元转换时参数波动剧烈,ΔPRI、ΔRF、ΔPW数值较大,因而红点分布分散、密度小;波形单元内参数波动较小,黑点分布集中、密度大。同时可以看到,图3中黑点聚集内部和边界也都分布着些许红点,这是因为波形单元转换时部分起始脉冲的某个或某两个参数变化较小,从而在二维的分布图中难以与黑点分离,因此二分类时需要同时考虑ΔPRI、ΔRF、ΔPW,利用参数间的联合变化规律获取更高的准确率。
针对非协作MFR缺乏匹配标准和标签数据的问题,采用密度聚类的无监督学习算法,将非起始脉冲构成的大密度“簇”从起始脉冲构成的噪声背景中分离出来,实现二分类,从而提取MFR波形单元。
根据邻域半径Eps和邻域点数MinPts的输入参数,DBSCAN聚类算法[12]将高密度数据集合从低密度区域中划分出来。传统的DBSCAN聚类算法需要人为设置两个输入参数,文献[13]采用密度阈值Density综合考量Eps和MinPts,定义Density为以Eps为半径的圆存在MinPts个数据点,即
(4)
Density太大可能导致同簇集合内部被划分为多个簇,太小则会带来不同簇合并的问题,因此选取适应于待聚类数据集的参数尤为重要。
在此基础上,本文提出了DA-DBSCAN聚类算法自适应确定密度阈值Density。由式(4)可知,Density由MinPts和Eps共同确定。对任意固定的MinPts值,可以通过寻找最优的Eps值获取合适的密度阈值Density。其中,MinPts一般不小于待聚类数据点的维度,综合考虑计算复杂度,现固定MinPts=3寻找最优的Eps,从而获取适应于待聚类数据集密度阈值Density。
(5)
3-dist值是指对象q∈D与数据集中到q第3近的对象p之间的距离值。3-dist排序图是所有对象的3-dist值升序排列构成的有序图[12]。数据集D的3-dist排序图如图4(a)所示,“簇”中数据点密度大3-dist值小,而噪声点密度小3-dist值大,曲线由平缓到突增存在着明显突变,因此突变区域中存在着某个阈值点将数据集的对象分为密度差异较大的两部份,该阈值便是最优的Eps值。
图4 3-dist值及3-dist差值排序曲线Fig.4 3-dist value and 3-dist variation sorting curve
(6)
如图4(b)所示,初期ΔR小范围内波动,随着序号增加,ΔR开始剧烈波动。ΔR依次表征了该簇的核心点、簇边界点以及噪声点3-dist值递增的速度。不规则的簇边界导致边界点对应的ΔR差异较大,因此当横轴序号遍历到边界点时,ΔR开始剧烈波动。
计算方差Var刻画ΔR的波动程度。令Var={v1,v2,…,vn},其中:
(7)
方差Var的分布如图5(a)所示,曲线突变点以前对应着簇的核心点和边界点,突变点以后对应着噪声点,因此该突变点即为所求阈值点。现引入熵E={e2,e3,…,en}来求解该突增点,设数据集的集合M={D2,D3,…,Dn-1},其中Di={v2,v3,…,vi},i=2,3,…,n。令
图5 方差和熵的分布图Fig.5 Diagram of variance and entropy
(8)
(9)
熵E的分布如图5(b)所示,假设突变点序号为m,由于{v2,v3,…,vm-1}在0附近小范围波动,∀j∈{2,3,…,i},Pij(i={2,3,…,m-1})分布较均匀,因此随着序号i的增大,Pij的个数增多,熵ei增大。
当i=m时,vm≫vj(j={2,3,…,m-1}),则Pmm≫Pmj(j={2,3,…,m-1}),熵ei开始减小,于是有m=arg max(E)。E取最大值对应的序号即为所求阈值点的序号,因此最优的Eps值为
Epsbest=rarg max(E)
(10)
至此DA-DBSCAN算法的两个输入参数已确定,聚类后得到簇Ci和噪声点N。对于二分类模型,只有波形单元非起始脉冲构成的一个簇,因此将簇Ci都合并为一个簇C。序号属于簇C对应非起始脉冲,而序号属于噪声点集合N则为起始脉冲。
对于二分类问题,真正例(true positive,TP)是指正例的数据点被标记为正例;假正例(false positive,FP)是指反例的数据点被标记为正例;真反例(true negative,TN)是指反例的数据点被标记为反例;假反例(false negative,FN)是指正例的数据点被标记为反例。
通过真正率(true positive rate,TPR)、假正率(false positive rate,FPR)、F-值(F-score)[14]3个指标来度量波形单元提取模型及算法的性能。相关定义如下:
(11)
(12)
F-score是一种常用于评价聚类结果的指标[15],公式如下:
(13)
TPR衡量了模型将波形单元起始脉冲检测出来的能力,FPR衡量了模型将非起始脉冲误判为起始脉冲的虚警率,而F-score综合评价了聚类的分类结果对实际数据类型标签的逼近能力。
在检测波形单元起始脉冲的问题中,检测率和虚警率难以同时优化,强检测能力对应着高虚警,而降低虚警也会导致检测能力的下降。具体到聚类算法中,第2.1节中选取MinPts=3,Eps=Epsbest作为输入参数,若增大Eps的取值,Density将减小。DA-DBSCAN算法通过密度可达形成簇,此时过渡区域的对象更容易被判定属于簇,模型FPR减小,具有较低的虚警率,同时TPR也减小,检测能力也将下降。因此,可以针对实际应用场景调节Eps获得预期的用户偏向需求。
设L为熵的最大序号m′与最大值点em的序号差m′-m,如图5(b)中红色部分的序号长度。输入参数λ∈(-1,1),则
Eps=r(m+λ L)
(14)
当λ=0时,模型不体现额外的需求偏向;当λ为负时,模型偏向具有更高的TPR和检测能力;当λ为正时,模型偏向具有更低的FPR和虚警率。且有
(15)
MFR波形单元提取算法的实现步骤如表1所示。
算法 1 MFR波形单元提取算法输入 多功能雷达脉冲列P,参数λ(默认为0)输出 波形单元的起始脉冲和非起始脉冲1.计算多功能雷达脉冲列P的{ΔPRI,ΔRF,ΔPW},建立数据集D;2.计算数据集D中所有对象的3-dist值,由小到大排序后得到R;3.对R做差得到ΔR;4.计算ΔR的区间方差Var;5.根据Var计算熵E,并寻找最大熵值对应的序号m;若λ≠0,计算熵E最大序号m′与最大熵值点的序号差L=m′-m;6.确定聚类算法的输入参数Eps=r(m+λL)和MinPts=3;7.对数据集D进行密度聚类,得到簇Ci和噪声点集合N;8.将簇Ci合并为一个簇C,令序号属于簇C的脉冲label值为0,即为非起始脉冲;令序号属于噪声点集合N的脉冲label值为1,即为起始脉冲。
MFR波形单元提取算法采用R语言实现,为了验证其有效性、适应性,现通过以下实验分别研究脉冲缺失、测量误差和参数选择对提取效果的影响。上述二维相控阵雷达的仿真波形库如表1所示,共包含5个类型40种波形单元。从波形库中随机均匀地选取100个波形单元构成脉冲列进行实验。
表1 某二维相控阵雷达的仿真波形库Table 1 Simulation waveform library of a two dimensional phased array radar
为验证本文构建的二分类模型和DA-DBSCAN聚类算法的性能,现进行以下对比实验。一方面,将DA-DBSCAN聚类算法与传统DBSCAN、层次聚类算法进行性能对比,从而验证DA-DBSCAN聚类算法的有效性。另一方面,对同一聚类算法分别将D′={PRI,RF,PW}和D={ΔPRI,ΔRF,ΔPW}作为输入数据集,验证所提二分类模型的科学性。其中,传统DBSCAN算法的参数MinPts设置为6时,参数Eps的取值范围大致位于区间[5,9](数据集为D′)和[3,8](数据集为D),随机选取该区间内3组距离值进行实验。提取效果的相关指标如表2所示。
表2 不同聚类算法提取性能对比Table 2 Performance comparison of different clustering algorithms
由表2可知,同一聚类算法采用二分类模型后,其F-score都得到了明显提高。二分类模型将波形单元内部和转换时的波动差异作为先验知识引入到无监督聚类算法中,充分利用了数据集的整体分布信息,从而大幅优化了提取效果。另一方面,传统DBSCAN聚类算法的输入参数依赖人为设置,参数设置质量直接影响了聚类精度,难以应用于工程实践。DA-DBSCAN聚类算法能够根据数据集的分布特点,自适应确定输入参数,实现了波形单元的全自动提取。同时,与传统DBSCAN和层级聚类算法相比,F-score达到0.979,具有最高的聚类精度。
脉冲丢失率(ratio of dropped pulses,RDP)[4]和误差偏离水平(error deviation level,EDL)[4]的定义如下:
(16)
(17)
RDP和EDL取值为0%~50%,MinPts=3,λ=0,Eps=Epsbest,重复1 000次实验的结果如表3和图6所示。
表3 不同RDP/EDL下的提取结果Table 3 Extraction results under different RDP/EDL
图6 不同RDP/EDL下的提取效果Fig.6 Extraction effect under different RDP/EDL
总体上随着RDP和EDL的增加,算法检测能力下降,虚警升高,聚类结果的F-score下降。但脉冲缺失时,模型将采用缺失脉冲前后脉冲的参数差值表征波形单元,依然能够较好地刻画波形单元内部和转换时的差异。因此,当RDP取到50%,TPR仍在91.9%以上,FPR低于2.2%,F-score大于0.97,算法在高脉冲缺失率的条件下仍有较好的表现。相似地,算法对测量误差也有良好的适应性。
本实验通过研究不同参数下的提取效果来考察算法分类结果对实际数据标签的逼近能力,从而验证算法自适应确定参数的合理性。EDL为10%,RDP为20%,Eps取值为(r(m-L),r(m+L))的实验结果如图7所示。
图7 不同Eps值下的提取效果Fig.7 Extraction results under different Eps values
由于检测能力和虚警无法同时优化,随着Eps增大,TPR和FPR都下降,F-score综合评价了算法性能,呈现先增大后减小的趋势。图7中,红线相交处为Eps=Epsbest时对应的评价指标,此时F-Score临近最大值,验证了算法自适应确定Eps的合理性。此次实验Eps=Epsbest时聚类算法的结果如图8所示,算法有效地将非起始脉冲构成的大密度“簇”(黑色)从起始脉冲构成的噪声背景(红色)中分离出来,实现了二分类从而提取非协作MFR的波形单元。
图8 Eps=Epsbest的聚类结果Fig.8 Clustering result of Eps=Epsbest
另外,在第2.3节中提出通过输入参数λ调节Eps值获得预期的性能偏向。现以1/10的步长遍历参数区间λ∈[-1,1],令EDL为10%,RDP为20%,1 000次实验的结果如图9和表4所示。
图9 不同参数λ下的提取效果Fig.9 Extraction results under different λ values
表4 不同参数λ下的提取效果
当λ=0(Eps=Epsbest)时,起始脉冲的平均检测能力达到88.59%,平均虚警率为3.36%,算法在具有测量误差和丢失脉冲的条件下仍然表现出较好的提取性能;当λ∈[-1,0)时,TPR和FPR都增大,模型偏向牺牲虚警率获取更高检测率;当λ∈(0,1]时,TPR和FPR都减小,模型偏向牺牲检测率获取更低虚警率,且偏向程度与|λ|呈正比。因此在实际应用中,用户可以根据当前提取效果和表中评价指标的变化,选取合适的参数λ获得预期的算法性能偏向。
本文对MFR波形单元提取的问题展开了研究,首先构建了基于PRI、RF、PW序列的二分类模型,该模型充分利用了数据集的整体分布信息和多参数间的联合变化,能够较好地适应实际探测过程中的脉冲丢失和测量误差。然后,针对该模型提出了一种参数自适应确定的DA-DBSCAN聚类算法进行分类,从而对非协作MFR实现了波形单元提取。最后,通过引入输入参数λ可针对不同的用户需求调整算法性能。
仿真实验表明,所提方法有效提取了非协作MFR的波形单元,较好地应对了实际工程应用中缺乏波形库先验知识的问题。同时,能够适应测量误差和脉冲丢失带来的干扰,具有良好的准确性和鲁棒性,因而对实际探测获取的数据有较好地适用基础,具有重要的现实意义。另外,由于实际获取的MFR脉冲列还具有不同程度的虚假脉冲,在对实际获取数据进行波形单元提取前可以通过如TTP变换等参数过滤的预处理剔除虚假脉冲,降低算法的虚警率。