張 鵬， 黃長強(qiáng)， 魏政磊， 周 歡， 王永乾

(1.空軍工程大學(xué)航空工程學(xué)院， 西安， 710038； 2.95596部隊，河南商丘，476000)

當(dāng)前軍事裝備正向著無人化、智能化的方向發(fā)展[1]，制約無人作戰(zhàn)飛機(jī)作戰(zhàn)模式由對地偵查、打擊向?qū)兆鲬?zhàn)轉(zhuǎn)變的主要因素在于自主決策能力不足，過分依賴地面站遠(yuǎn)程操控，無法滿足高動態(tài)空戰(zhàn)的實時性需求?；诋?dāng)前空戰(zhàn)數(shù)據(jù)實時提取敵我態(tài)勢信息來指導(dǎo)無人平臺做出自主決策是無人自主空戰(zhàn)的關(guān)鍵技術(shù)?？諔?zhàn)態(tài)勢知識提取可以轉(zhuǎn)換為近距空戰(zhàn)態(tài)勢模式的發(fā)現(xiàn)，其本質(zhì)是多元時間序列的分割聚類問題。多元時間序列的聚類算法包括兩種思路：一種是通過對多元時間序列進(jìn)行壓縮降維，再利用傳統(tǒng)的點聚類方法進(jìn)行聚類[2]；另外一種是結(jié)合多元時間序列相似性度量方法改進(jìn)傳統(tǒng)的點聚類方法[3-4]。前者又分為基于特征的聚類和基于模型的聚類，后者為基于初始序列的聚類。基于特征和模型的聚類會忽略原始時間序列的局部特征。基于原始時間序列的聚類一般不對時間序列進(jìn)行壓縮，而是直接進(jìn)行聚類，該方法隨著時間序列增長，其計算代價越大[5]。同時目前的空戰(zhàn)態(tài)勢知識提取文獻(xiàn)大多數(shù)以空戰(zhàn)態(tài)勢信息到固定態(tài)勢類別的映射居多，不能體現(xiàn)聚類中心的客觀性和態(tài)勢變化的連續(xù)性。

針對以上問題，本文提出了一種基于拉普拉斯中心性和Kshape的分層聚類時序分析方法，解決了多維參數(shù)下的態(tài)勢信息提取問題。

1 基于拉普拉斯中心性的聚類數(shù)目確定

目前的分割方法主要是基于特征點的分割和基于監(jiān)督學(xué)習(xí)的分割點檢測[6-7]。然而實際的時序單元特征無法具體描述或者表示，因此采用基于無監(jiān)督學(xué)習(xí)的分割方法是非常有必要的[8]?？紤]到傳統(tǒng)聚類算法以及基于峰值密度算法(density peak clustering,DPC)[9]的聚類個數(shù)無法確定的缺點，本文采用拉普拉斯中心峰值聚類方法(Laplacian peaks clustering algorithm,LPC)[10]，將原始數(shù)據(jù)轉(zhuǎn)換到加權(quán)完全圖中，通過構(gòu)造決策圖可以選取聚類個數(shù)以及聚類中心；LPC在整個聚類過程中不需要任何經(jīng)驗參數(shù)，因此LPC也是一種無參數(shù)聚類算法。

1.1 相關(guān)參數(shù)定義

敵我雙方的態(tài)勢信息主要取決于空戰(zhàn)中雙方的相對運動關(guān)系，如表1所示。

表1 態(tài)勢描述參數(shù)

假設(shè)以我機(jī)質(zhì)心為坐標(biāo)原點O，X軸沿UCAV縱軸指向機(jī)頭方向，Y軸垂直于UCAV對稱面并值向右，Z軸位于UCAV對稱面并垂直于機(jī)體縱軸向下，我方UCAV與敵機(jī)位置坐標(biāo)為Pu=[xu,yu,zu]和Pe=[xe,ye,ze]，則雙機(jī)相對態(tài)勢定義見圖1。

圖1 空戰(zhàn)幾何關(guān)系

假設(shè)原始未分類的空戰(zhàn)態(tài)勢描述參數(shù)數(shù)據(jù)集為Xs={Xs,1,Xs,2,…,Xs,tn}，其中Xs,i=[Ri,Δhi,vu,i,φu,i,qu,i,ve,i,φe,i,qe,i]，i=1,2,…,tn，tn為獲取的空戰(zhàn)數(shù)據(jù)末端時間點。為了挖掘聚類中心，定義2個主要指數(shù)：第i個數(shù)據(jù)點的拉普拉斯中心性ci和最小距離值δi。其次，將數(shù)據(jù)集轉(zhuǎn)換為加權(quán)完全圖G=(Xs,E,W)，將每個數(shù)據(jù)點看作加權(quán)網(wǎng)絡(luò)E的節(jié)點，W表示邊權(quán)值的集合。具體定義如下：

定義1相異度矩陣W。wi,j是數(shù)據(jù)點Xs,i與Xs,j之間的邊權(quán)值，用于衡量多元時間序列相似性。由于原始態(tài)勢參數(shù)序列不等長且具有高動態(tài)性，本文采用動態(tài)時間彎曲距離法(dynamic time warping，DTW)[11]進(jìn)行計算。

(1)

定義2和對角矩陣Z：

Z(G)=diag(z1,z2,…,zn)

(2)

定義3拉普拉斯矩陣L：

L(G)=Z(G)-W(G)

(3)

定義4拉普拉斯能量EL(G)：

(4)

式中：λ1,λ2,…,λtn為拉普拉斯矩陣L(G)的特征值。

定義5拉普拉斯中心性ci。在加權(quán)完全圖G中，刪除節(jié)點Xs,i以及與其相連的邊，得到的加權(quán)完全圖為Gi。節(jié)點Xs,i的拉普拉斯中心性ci可以定義為：

(5)

式中：EL(Gi)表示刪除節(jié)點Xs,i之后的加權(quán)完全圖的拉普拉斯能量，拉普拉斯中心性的意義在于刪除當(dāng)前節(jié)點后的拉普拉斯能量相對下降。當(dāng)前節(jié)點拉普拉斯中心性越高，則該節(jié)點重要程度越高[11]。

定義6最小距離值δi。對于當(dāng)前節(jié)點Xs,i，通過比較拉普拉斯中心性ci，可以得到最小距離值δi，其具體計算如下：

δi=

(6)

1.2 基于正態(tài)分布“3σ”準(zhǔn)則的聚類個數(shù)提取策略

根據(jù)上述定義的分析，具有較高最小距離值和較高拉普拉斯中心性的數(shù)據(jù)節(jié)點可以作為聚類中心，由此可以確定聚類個數(shù)。利用最小距離值δ和拉普拉斯中心性c構(gòu)造決策圖[12]。為了自動完成這一過程，本文采用了基于正態(tài)分布“3σ”準(zhǔn)則的聚類個數(shù)提取策略。為了從數(shù)據(jù)節(jié)點中自動提取聚類個數(shù)，引出另外一個決策量γ，其定義如下：

γi=ci·δi

(7)

為了方便提取聚類中心與聚類個數(shù)，對γ值進(jìn)行放大：

(8)

式中：Eγ,i為節(jié)點Xs,i由決策量γi放大后的值。

假設(shè)數(shù)據(jù)節(jié)點的擴(kuò)大后決策值Eγ服從正態(tài)分布，根據(jù)正態(tài)分布的“3σ”準(zhǔn)則，將(μ-3σ,μ+3σ)區(qū)域看作該隨機(jī)變量實際可能的取值區(qū)間。由于潛在的聚類中心的決策值Eγ遠(yuǎn)大于非聚類中心的Eγ值，可以得到Eγ值大于μ+3σ的數(shù)據(jù)節(jié)點。圖2給出基于拉普拉斯中心性的決策圖，其中紅色的數(shù)據(jù)節(jié)點為Eγ值大于μ+3σ的數(shù)據(jù)節(jié)點，即聚類中心，因此聚類個數(shù)為3。

圖2 某空戰(zhàn)態(tài)勢描述參數(shù)決策圖

2 基于拉普拉斯中心性的聚類中心確定

時間序列的聚類算法的提出主要來源于2種途徑[13]：第1種是替代相似性度量方法，屬于基于原始數(shù)據(jù)的聚類方法；第2種是將時間序列轉(zhuǎn)換為靜態(tài)數(shù)據(jù)，包括基于特征和基于模型的聚類方法?；谔卣骱突谀Ｐ偷木垲惙椒ㄊ菍υ紩r間序列進(jìn)行特征提取或者模型表示，不能體現(xiàn)出時序信息。本文選擇基于原始數(shù)據(jù)的時間序列聚類方法，其主要依賴于矩陣的分解。

針對以上問題的分析，本文提出了一種基于原始數(shù)據(jù)的多元時間序列聚類方法，即Kshape。首先，提出一種基于互相關(guān)度量的形狀距離度量方法；其次，基于形狀距離度量方法，設(shè)計一種計算時間序列聚類中心的計算方法；最后，根據(jù)相似性度量方法與聚類中心計算方法，給出了基于Kshape的多元時間序列聚類算法具體步驟。

2.1 時間形狀相似性度量方法

為了提高時間序列相似性度量計算效率，本文設(shè)計了一種互相關(guān)度量方法。

對于序列Xs={Xs,1,Xs,2,…,Xs,n}與Ys={Ys,1,Ys,2,…,Ys,m}，為了實現(xiàn)位移不變性，固定時間序列Ys不動，滑動Xs，互相關(guān)度量就是計算滑動過程中對齊的元素之間的內(nèi)積。Xs時間序列的移位序列計算表達(dá)式為：

(9)

式中：Xs(κ)表示移位κ之后的時間序列。一般情況下，κ∈(-n,n)或者κ∈[-n+1,n-1]。假設(shè)n=3，那么其移位時間序列的個數(shù)為2n-1=5，為了方便理解，給出示意圖見圖3。移位序列中，Xs(0)就是原始時間序列，其他移位序列的某些元素需要填充為0。

圖3 移位時間序列

根據(jù)移位序列Xs(k)，可以得到互相關(guān)序列CCω(Xs,Ys)=(c1,c2,…,cω)，其長度為2n-1，與移位序列的個數(shù)相等。互相關(guān)序列可以定義如下：

CCω(Xs,Ys)=Rω-n(Xs,Ys)

(10)

式中：ω∈{1,2,…,2n-1}。同時，Rω-n(Xs,Ys)的計算公式如下：

(11)

通過計算式(10)中的CCω(Xs,Ys)中最大的互相關(guān)序列值，得到對應(yīng)的ω?；讦刂?，可以得到κ=ω-n，對應(yīng)的移位序列Xs(κ)是最優(yōu)的匹配序列。由于問題的不同，則采用的標(biāo)準(zhǔn)化方法不同。本文選取了文獻(xiàn)[15]中的系數(shù)標(biāo)準(zhǔn)化方法，即NCCc，具體定義如下：

(12)

式中：NCCc(Xs,Ys)的值域為[-1,1]。NCCc(Xs,Ys)=1說明兩個時間序列相似，NCCc(Xs,Ys)=-1說明兩個時間序列負(fù)相關(guān)。

為了設(shè)計一種基于形狀距離的度量方法，將NCCc(Xs,Ys)取值映射到[0,2]，具體計算公式如下：

SBD(Xs,Ys)=

(13)

為了提高基于形狀距離(shape-based distance,SBD)的相似性度量計算過程，采用Xs與Ys離散傅里葉變換內(nèi)積的離散傅里葉逆變換計算[14]。具體的計算公式如下：

CC(Xs,Ys)=F-1{(Xs)(Ys)}

(14)

(15)

(16)

2.2 聚類中心提取策略

(17)

代入到式(17)中，得到：

(18)

2.3 基于Kshape的多元時間序列聚類

基于SBD度量方法與聚類中心提取策略的Kshape聚類算法類似于k-means，可以更新得到聚類劃分較好的聚類中心，且可以有效地處理時間序列的縮放、平移、轉(zhuǎn)移不變性。Kshape聚類過程是一個迭代尋優(yōu)的過程，主要的步驟為：①時間序列樣本分配過程：根據(jù)聚類中心，分配時間序列Xs到距離最近聚類中心；②聚類中心更新過程。依據(jù)上一迭代中聚類關(guān)系的變化，重新計算聚類中心。Kshape算法重復(fù)步驟①與步驟②，直到聚類關(guān)系不再發(fā)生變化(各類別樣本到聚類中心的SBD距離平方和不再發(fā)生變化)或者達(dá)到最大迭代次數(shù)。

3 分層聚類時序分析

從一組空戰(zhàn)數(shù)據(jù)中尋找態(tài)勢時序單元(未知)是一個多元時間序列分割聚類問題。時間序列的分割聚類是一種非監(jiān)督學(xué)習(xí)問題，主要應(yīng)用于計算機(jī)視覺與計算機(jī)圖形學(xué)[15]?？諔?zhàn)態(tài)勢描述參數(shù)具有時序性，也具有空間性，可以分解為多種態(tài)勢時序單元，目前面臨的問題主要體現(xiàn)在：①態(tài)勢時序單元組合呈指數(shù)增長；②每個相同場景的態(tài)勢描述參數(shù)時間序列不一樣，具有時間的伸縮性；③態(tài)勢時序單元劃分粒度可變性。

針對以上問題，本文提出了一種分層的聚類時序分析方法(HACA)。該方法結(jié)合聚類算法與分割算法，提出坐標(biāo)下降法在計算分割段以及其聚類過程中選擇求解。

3.1 聚類時序分析方法

假設(shè)給定多元時間序列Xs={Xs,1,Xs,2,…,Xs,n}，其中Xs∈Rd×n，n為時間序列的時間長度，d=8 表示態(tài)勢描述參數(shù)的維度。將Xs分割為mc個態(tài)勢時序單元，其中每個單元對應(yīng)k聚類類別中的某一類。假設(shè)第i個分割態(tài)勢時序單元Yi=Xs,[si,si+1)=[xsi，…,xsi+1-1∈Rd×n]，該態(tài)勢時序單元的開始時刻為si，結(jié)束時刻為si+1-1，其長度ni=si+1-si必須滿足ni≤nmax，其中nmax為態(tài)勢時序最大長度且控制時間序列分割粒度。為了方便表示與計算每個態(tài)勢時序單元的所屬類別，則定義聚類指示矩陣G∈{0,1}k×mc，其中聚類指示矩陣元素gci=1表示Yi屬于類別c，否則gci=0。結(jié)合第2節(jié)的時間序列聚類算法Kshape，基于聚類時序分析的聚類分割方法目標(biāo)函數(shù)為：

SBD([Y1,…,Ymc],UcG)

s.t.GT1k=1mandsi+1-si∈[1,nmax]

(19)

式中：G∈{0,1}k×mc為聚類類別指示器，s∈Rmc+1為態(tài)勢時序單元的分割起始點與末端點時刻，Yi=Xs,[si,si+1]表示第i個分割段(態(tài)勢時序單元)，Uc={u1,u2,…,uk}為聚類中心集合。本文采用的多元時間序列相似性度量方法為基于形狀距離(SBD)的度量方法。

式(19)是一個關(guān)于G和s整數(shù)規(guī)劃問題。根據(jù)上述的定義，G表示態(tài)勢時序單元與聚類中心的指示矩陣，s表示態(tài)勢描述參數(shù)點與態(tài)勢時序單元之間的對應(yīng)關(guān)系，優(yōu)化關(guān)于G和s的聚類分割方法目標(biāo)函數(shù)是一個非線性多項式NP難問題。為了解決該問題，本文采用坐標(biāo)下降法框架，利用動態(tài)規(guī)劃策略來優(yōu)化計算s，利用winner-take-all策略來優(yōu)化計算G。每一迭代次數(shù)具體優(yōu)化如下表示：

(20)

式中：uc是類別為c的聚類中心，該聚類中心可以由基于SBD的聚類中心提取策略，由上一迭代次數(shù)的解(G′,s′)得到。式(20)的計算復(fù)雜度隨著原始時間序列的增長而增加，因此，本文采用了一種基于動態(tài)規(guī)劃的算法去檢測可能的態(tài)勢時序單元。為了進(jìn)一步提高求解(G,s)的效率，引入一種輔助函數(shù)：

計算機(jī)作為高新技術(shù)應(yīng)用，對于提高學(xué)生社會實踐能力及計算機(jī)技術(shù)水平具有良好的推動作用。尤其是極域多媒體電子教室軟件應(yīng)用，對解決計算機(jī)基礎(chǔ)性教學(xué)的諸多難題有一定的幫助，使學(xué)生在計算機(jī)知識學(xué)習(xí)方面，不再受到基礎(chǔ)學(xué)習(xí)環(huán)境及知識學(xué)習(xí)能力的限制，實現(xiàn)對學(xué)生計算機(jī)知識的立體化教學(xué)。通過理論知識講解與教學(xué)實現(xiàn)對學(xué)生計算機(jī)知識的掌握能力及應(yīng)用能力進(jìn)行深度提高。

(21)

式中：Xs,[i,v]=[x1,x2,…,xv]表示末端時刻為v的態(tài)勢時序單元，通過優(yōu)化求解v來替代優(yōu)化求解s。利用動態(tài)規(guī)劃的思想，可以將式(21)改寫為：

J(v)=

(22)

式中：時間序列Xs被劃分為Xs,[1,i-1]與Xs,[i,v]，優(yōu)化公式(22)使其最小化，得到最優(yōu)的子序列Xs,[1,i-1]與Xs,[i,v]。盡管分割時間序列Xs的可能情況隨著時間序列的時間長度增加n，但動態(tài)規(guī)劃方法利用Bellman方程解決最優(yōu)問題仍然是一個有效的方法。

J(v)=

(23)

步驟2后向過程。根據(jù)步驟1中得到的最優(yōu)分割點位置與聚類類別，從右向左依次劃分，通過計算得到s與G。

重復(fù)步驟1與步驟2直到J(n)收斂不變或者迭代次數(shù)達(dá)到最大。

3.2 分層策略

圖4 分層聚類時序分析結(jié)果對比

4 基于L-Kshape-HACA的分層聚類時序分析仿真

結(jié)合基于拉普拉斯中心性的聚類個數(shù)確定方法與基于Kshape的多元時間序列聚類算法，本文提出了一種基于L-Kshape的分層聚類時序分析方法(L-Kshape-HACA)，用于解決長時域態(tài)勢描述參數(shù)序列的分割聚類問題。L-Kshape-HACA分割聚類方法流程圖如圖5所示，其中迭代停止標(biāo)準(zhǔn)均為J(n)不再發(fā)生變化或者達(dá)到最大迭代次數(shù)。

圖5 基于L-Kshape-HACA的分層聚類流程圖

態(tài)勢聚類的有效性是實現(xiàn)正確的近距空戰(zhàn)決策的數(shù)據(jù)基礎(chǔ)。為了驗證L-Kshape-HACA分割聚類方法解決近距空戰(zhàn)態(tài)勢描述參數(shù)態(tài)勢時序單元提取問題的可行性與有效性，利用基于貝葉斯推理的態(tài)勢評估方法[16]、基于多維點數(shù)據(jù)的聚類算法(聚類數(shù)據(jù)由多維空間中點數(shù)據(jù)組成)、基于模型表示的時間序列聚類算法和于多元時間序列數(shù)據(jù)的分割聚類算法與L-Kshape-HACA分割聚類算法進(jìn)行比較。為了對比基于數(shù)據(jù)驅(qū)動的態(tài)勢聚類分析與基于推理的態(tài)勢評估，本仿真實驗的對比算法包括文獻(xiàn)[16]中基于貝葉斯推理的態(tài)勢評估方法；為了對比基于時序數(shù)據(jù)的態(tài)勢聚類與基于點數(shù)據(jù)的態(tài)勢聚類，本仿真實驗比較的基于點數(shù)據(jù)態(tài)勢聚類算法還包括自組織特征映射神經(jīng)網(wǎng)絡(luò)(self-organizing feature map,SOM)[17]、層次聚類(hierachical clustering,HC)[18]、模糊C均值的態(tài)勢聚類[19]、高斯混合的態(tài)勢聚類[20]、k-means、DBSCAN以及文獻(xiàn)[9]中基于密度峰值的態(tài)勢聚類；為了比較基于分割聚類的時序態(tài)勢知識提取方法的好壞，本仿真實驗比較的分割聚類方法包括ACA[8]、基于譜聚類的分層ACA(HACA)[16]、基于k-means的分層ACA(k-means-HACA)以及基于密度峰值的分層ACA(DPC-HACA)；除此之外，為了測試基于拉普拉斯中心性聚類個數(shù)的確定方法是否有效，則需要比較基于Kshape的分層ACA(Kshape-HACA)。本文仿真實驗數(shù)據(jù)來源于某近距空戰(zhàn)模擬訓(xùn)練的對抗數(shù)據(jù)，該實驗選取其中的12組近距空戰(zhàn)對抗數(shù)據(jù)，其皆為紅方(我機(jī))獲取勝利的數(shù)據(jù)，實驗仿真環(huán)境為Windows 10，CPU為2.80 GHz,8 GB內(nèi)存，編程語言Matlab。

表2 L-Kshape-HACA不同時序分割最大值設(shè)置

圖6 不同時序分割長度最大值的影響

基于L-Kshape-HACA分割聚類方法與比較算法的仿真測試結(jié)果如圖7和圖8所示。

本文選取了4組空戰(zhàn)對抗數(shù)據(jù)進(jìn)行了對比，圖7給出了基于L-Kshape-HACA分割聚類算法與其他算法的時間序列分割聚類示意圖，不同的顏色表示不同的類別，可以看出，與基于先驗知識貝葉斯推理的態(tài)勢評估結(jié)果相比，基于L-Kshape-HACA的態(tài)勢描述參數(shù)序列分割聚類結(jié)果是最接近的。圖8給出了態(tài)勢描述參數(shù)歸一化之后的特征空間，可以看出，整個態(tài)勢描述空間可以聚類為4類，態(tài)勢描述參數(shù)序列可以分割為9個態(tài)勢時序單元。為了表達(dá)這種抽象的態(tài)勢時序單元，結(jié)合基于貝葉斯推理的態(tài)勢評估結(jié)果，根據(jù)圖8所示的對抗軌跡分割聚類結(jié)果，得到以下態(tài)勢時序單元特征：1)Cluster #1：我機(jī)其他狀態(tài)-互不安全的轉(zhuǎn)換；2)Cluster #2：我機(jī)中立互為安全-守勢-中立相向狀態(tài)的轉(zhuǎn)換；3)Cluster #3：我機(jī)攻勢-互為安全狀態(tài)的轉(zhuǎn)換；4)Cluster #4：我機(jī)互不安全-攻勢狀態(tài)的轉(zhuǎn)換。

圖7 L-Kshape-HACA與比較算法的分割聚類結(jié)果對比

圖8 Case(2)空戰(zhàn)對抗組的態(tài)勢時序單元分析

5 結(jié)語

無人自主空戰(zhàn)是未來智能作戰(zhàn)的關(guān)鍵技術(shù)，本文設(shè)計了一種基于數(shù)據(jù)的無監(jiān)督態(tài)勢信息提取新方法，為無人空戰(zhàn)決策提供了更加客觀可靠的數(shù)據(jù)基礎(chǔ)。本文研究的泛化分析算法會在態(tài)勢提取過程中忽略飛機(jī)獨有的戰(zhàn)術(shù)機(jī)動性能，因此，下一步研究必須要針對特定機(jī)型的空戰(zhàn)數(shù)據(jù)進(jìn)行研究，從而提高戰(zhàn)術(shù)機(jī)動態(tài)勢提取的準(zhǔn)確度。