基于L-Kshape-HACA的空戰態勢分割聚類

張 鵬， 黃長強， 魏政磊， 周 歡， 王永乾

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

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

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

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

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

1.1 相關參數定義

敵我雙方的態勢信息主要取決于空戰中雙方的相對運動關系，如表1所示。

表1 態勢描述參數

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

圖1 空戰幾何關系

假設原始未分類的空戰態勢描述參數數據集為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為獲取的空戰數據末端時間點。為了挖掘聚類中心，定義2個主要指數：第i個數據點的拉普拉斯中心性ci和最小距離值δi。其次，將數據集轉換為加權完全圖G=(Xs,E,W)，將每個數據點看作加權網絡E的節點，W表示邊權值的集合。具體定義如下：

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

(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。在加權完全圖G中，刪除節點Xs,i以及與其相連的邊，得到的加權完全圖為Gi。節點Xs,i的拉普拉斯中心性ci可以定義為：

(5)

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

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

δi=

(6)

1.2 基于正態分布“3σ”準則的聚類個數提取策略

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

γi=ci·δi

(7)

為了方便提取聚類中心與聚類個數，對γ值進行放大：

(8)

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

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

圖2 某空戰態勢描述參數決策圖

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

時間序列的聚類算法的提出主要來源于2種途徑[13]：第1種是替代相似性度量方法，屬于基于原始數據的聚類方法；第2種是將時間序列轉換為靜態數據，包括基于特征和基于模型的聚類方法。基于特征和基于模型的聚類方法是對原始時間序列進行特征提取或者模型表示，不能體現出時序信息。本文選擇基于原始數據的時間序列聚類方法，其主要依賴于矩陣的分解。

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

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

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

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

(9)

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

圖3 移位時間序列

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

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

(10)

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

(11)

通過計算式(10)中的CCω(Xs,Ys)中最大的互相關序列值，得到對應的ω。基于ω值，可以得到κ=ω-n，對應的移位序列Xs(κ)是最優的匹配序列。由于問題的不同，則采用的標準化方法不同。本文選取了文獻[15]中的系數標準化方法，即NCCc，具體定義如下：

(12)

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

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

SBD(Xs,Ys)=

(13)

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

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

(14)

(15)

(16)

2.2 聚類中心提取策略

(17)

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

(18)

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

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

3 分層聚類時序分析

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

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

3.1 聚類時序分析方法

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

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

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

(19)

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

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

(20)

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

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

(21)

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

J(v)=

(22)

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

J(v)=

(23)

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

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

3.2 分層策略

圖4 分層聚類時序分析結果對比

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

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

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

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

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

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

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

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

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

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

5 結語

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