基于聚類思想的彈道外推算法研究

李同亮 ， 朱 勇 ， 于 瓊

(1.中國電子科技集團公司第三十八研究所， 安徽合肥 230088；2.中國人民解放軍陸軍工程大學， 江蘇南京 210007)

0 引言

炮位偵校雷達具備全天候、全天時、火炮類目標偵察校射能力，為炮兵火力打擊提供重要保障。彈道外推算法是決定炮位偵校雷達性能的關鍵技術之一，炮位偵校雷達的定位精度性能直接受其影響。

當前的彈道外推算法主要包括兩大類：一是基于先驗樣本的機器學習算法；二是基于彈道模型的卡爾曼濾波估計算法。前者受實際彈道樣本的限制，算法實用性不高，目前在炮位偵察校射雷達彈道外推領域，更多地采用卡爾曼濾波及其衍生算法。然而，基于卡爾曼濾波的外推算法受彈道、氣象條件以及量測點跡質量影響較大，故外推精度會受到影響。針對不同彈道，如何提升外推結果的一致性及外推精度是炮位偵校雷達領域亟待解決的一大難題。

為此，本文將聚類思想引入至彈道外推過程中，將聚類思想與七態UKF濾波算法相結合。對量測數據進行多次濾波外推，獲取多個外推結果，然后采用-均值聚類算法對外推結果進行聚類處理，最后采用綜合多因子方法計算簇品質，選取最優簇對應的聚類中心作為最終的火炮位置進行輸出。經實驗仿真驗證，本算法有效提升了外推結果的一致性及定位精度。

1 彈道運動模型

在炮位偵察模式下，因炮彈種類未知，相關彈道參數均未知，射表中的標準條件與實際條件也不能保證完全一致，并且在炮位偵察過程中炮口是未知的，這時，我們采用前序數個量測點粗略擬合出炮口位置，并以此作為坐標中心進行后續濾波預測及外推處理，故我們以炮口為坐標中心，基于彈丸質心運動微分方程組(非標準條件下)建立彈道運動模型：

(1)

式中：，，為炮彈離開炮口后在射向上的水平分量、垂直分量和側偏分量；，，為對應的水平速度、垂直速度以及側偏速度；為彈道上對應時刻的合成速度；為彈道系數；()為空氣密度函數；()為空氣阻力函數；為標準狀態下的虛溫，對應于不同高程上的虛溫；為重力加速度。

圖1表示炮口坐標系。其中，為炮口中心，對應射擊方向，垂直水平面向上，基于右手螺旋法則，指向右方。彈道方程組的解算基于炮口坐標系完成。射擊方向與坐標北的夾角為射向(順時針為正)。

圖1 炮口坐標系

2 彈道外推

彈道外推的主要任務是接收彈道量測數據，根據彈道目標飛行特征采用相應的濾波模型進行跟蹤濾波處理，完成彈丸類目標軌道的正反向外推，獲取炮彈發點(偵察)或落點(校射)。本文彈道處理算法專注于偵察模式發點推算，具體處理流程如圖2所示。

圖2 彈道外推流程圖

彈道外推軟件實時接收來自彈丸飛行過程中的雷達量測數據，當量測點數積累至一定數量后，開始啟用彈道處理算法。首先，基于彈丸質心運動模型創建UKF濾波器進行彈道濾波處理。然后，采用四階龍格-庫塔算法進行彈道解算并依據炮位高程計算炮位發點。隨雷達觀測點增加個數，上述過程重復次。最后，對前序過程求解出的炮位發點進行-均值聚類處理，并采用綜合多因子算法選取最優簇，以該簇對應的聚類中心作為最終外推結果輸出。

2.1 七態UKF濾波

UKF濾波算法對非線性問題的處理效果顯著且該方法估計精度至少達到2階，計算量與EKF 算法處于同一量級，不需要計算Jacobian 矩陣，可以處理不可導的非線性函數。基于上述優點以及彈道運動方程非線性特性，本文采用七態UKF濾波方法，將彈道系數作為第七維狀態進行濾波，實時估計并更新彈道系數。

基于第1節中彈道運動模型建立狀態方程和量測方程。

狀態方程：是描述動態系統各狀態變量和輸入之間關系的方程。公式(2)描述了關于七維向量= [,,,,,,]狀態方程:

(2)

式中，為高斯白噪聲，協方差矩陣為。

量測方程：量測方程用來表示狀態向量、測量向量以及輸入之間關系。本文算法中量測向量= []，量測方程描述如下：

=

(3)

(4)

式(3)中為坐標轉移矩陣，將炮口坐標系下的狀態向量轉移至以雷達站址為中心的東北天坐標系。式(4)中代表高斯白噪聲，協方差矩陣為。

考慮如下非線性模型：

(5)

式中，∈R為七維狀態向量，為七維過程噪聲，-1為七維狀態向量函數，為狀態與量測轉換函數，為三維隨機量測噪聲。其中，過程噪聲與量測噪聲為相互獨立，不相關的高斯白噪聲，和分別為對應的協方差陣。具體濾波過程描述如下：

1) 狀態初始化。

(6)

(7)

(8)

(9)

②公式(10)描述了一步提前預測計算過程。

(10)

4) 濾波更新。

(11)

式中，表示濾波器增益矩陣。

各采樣點權值計算：

(12)

2.2 彈道解算

目前彈道解算主要通過數值計算方法求解彈道方程組，進而完成發點或落點的推算。彈道方程組通常由一階變系數聯立方程組表示，一般情況下采用數值方法求解其數值解是唯一手段，在某些特殊情況下對其簡化處理后才能求得近似解析解。對微分方程求解的方法較多，如阿當姆斯預報-校正法和龍格-庫塔法，而四階龍格-庫塔法在工程中得到廣泛應用，故本文算法采用龍格-庫塔法進行彈道解算。

龍格-庫塔法是基于泰勒級數的一種改進算法，四階龍格-庫塔法的具體描述如公式(13)～(15)所示，其中公式(13)表示微分方程組及初值：

0(=1,2,…,)

(13)

方程組在第點處的所有變量的值為 (,1,2,…,)，那么+1點處各變量的四階龍格-庫塔計算公式為

(14)

其中：

(15)

式中，代表步長，步長的選取會影響四階龍格-庫塔法的計算精度。步長越小計算精度越高，然而過小的步長會導致迭代計算過程中的累計誤差變大，同時會消耗更多的計算時間，影響時效性。本文的選擇分為兩步，先進行粗粒度外推，再進行細粒度外推，這樣既保證了外推精度，又不會影響時效性。

2.3 K-均值聚類處理

為了進一步提升彈道外推精度及外推結果的一致性，引入聚類思想，采用-均值聚類算法對多次外推結果進行聚類處理，選取最優簇對應的聚類中心作為外推結果輸出。-均值聚類算法屬于劃分型的動態聚類算法，其計算過程需要給定待聚類的數目，通過聚類處理后數據集被劃分成個不同的類。該算法的核心思想是：給定聚類的個數后，第一次迭代的個中心點被隨機選取。這時，依次計算個中心點與數據集中其他數據之間的距離，通過比較篩選，數據被劃分至距離其最近的類中。然后，計算新生成的數據類的聚類中心，同時調整數據集。倘若新舊類之間的聚類中心沒有變化或者在某個較小的范圍內變化，那么聚類完成。

圖3描述了在彈道外推過程中-均值聚類的具體處理步驟，詳細處理流程如下：

圖3 K-均值聚類算法流程

輸入： 彈道外推發點數據集={,,…,}，待聚類的個數為，本文算法值可以取2～5；

輸出：個聚類{,,…,};

1) 在給定的彈道外推發點數據集={,,…,}中，隨機抽取個不相同的發點樣本用作初始聚類中心點{,,…,};

2) 計算,,…,這個中心點與數據集中其他數據之間的距離，記作(,)，其中，= 1,2,…,，= 1,2,…,。當

(,)=min{(,)},

=1,2,…,

(16)

則表明屬于類。

3) 當所有數據的歸屬類調整之后，利用公式(18)重新生成個類的聚類中心：

(17)

并計算誤差平方和準則函數()，()代表第次迭代計算得出的誤差平方和值。其中公式(18)描述了某一次誤差平方和的計算過程。

(18)

4) 倘若第和+1次迭代，誤差平方和沒有發生較大的變化，則表明誤差平方和已經收斂，可以跳出迭代過程。反之，進入步驟2)繼續迭代。

5) 當聚類收斂后采用綜合多因子方法選取最優簇。為此定義了簇品質以及三類因子：

簇品質_：用于表示簇優劣程度的量稱為簇品質。

簇體積因子_：假設簇成員體積大小相同，均為1 m，則此時，可用簇成員個數表示簇體積。

簇質量因子_：用簇成員質量平均值表示簇質量。簇成員質量定義為參與濾波的量測點跡的數目，簇成員中參與濾波的量測點跡個數越多，表示該成員質量越大。

簇一致性因子_：用該簇外推發點位置的方差表示簇一致性。方差越小，一致性越好。

對上述三類因子分別做歸一化處理，得到歸一化簇體積因子_=[_，_,…,_res]，歸一化簇質量因子_=[_,_,…,_res]和歸一化簇一致性因子=[__,…,_res]，然后利用公式(19)計算簇品質：

_=[___][]

(19)

式中，_=[_,_,…,_]為簇品質向量，,,分別表示歸一化簇體積因子、歸一化簇質量因子及歸一化簇一致性因子對應的權值。最優簇的品質即為

__=max{_,_,…,_},

_∈[1,]

(20)

式中，_對應的簇即為最優簇。

3 算法仿真與分析

本文以155榴彈炮為例，利用彈道仿真軟件生成雷達量測數據，進行數值仿真計算，仿真平臺為Intel Core i7-5600U、主頻2.6 GHz、四核CPU 計算機，仿真軟件為MATLAB。具體仿真條件如下：

1) 155榴彈發射條件: 初速設置為900 m/s，射角設置為35°；

2) 取雷達隨機測量誤差:= 76 m，= 17 mil，= 1.6 mil;

3) 雷達距離炮48 km。

基于上述數據對六態_EKF(原算法)、六態_UKF、七態_UKF算法、聚類_七態_UKF算法分別進行10 000次蒙特卡洛仿真實驗，并采用炮位偵察定位精度計算方法圓中間誤差(ECP)統計外推精度。

圖4在ENU坐標系下分析了六態_EKF(原算法)與聚類_七態_UKF彈道濾波情況，可以看出相比六態_EKF算法，聚類_七態_UKF算法濾波結果更接近真實值，效果更優。

圖4 原算法與本文算法東北天坐標系下濾波結果

圖5從濾波后位置誤差角度分析了二者濾波效果。結果表明本文算法在3個坐標軸上的濾波誤差明顯小于原算法，進一步驗證本文算法濾波效果優于原算法。

圖5 原算法與本文算法濾波后位置誤差比對

圖6和圖7給出了七態_UKF算法與聚類_七態_UKF算法目標外推發點位置散布情況，仿真結果可以看出，聚類_七態_UKF算法發點外推結果一致性明顯優于原算法，驗證了聚類處理的有效性。

圖6 七態_UKF算法外推發點散布情況

圖7 聚類_七態_UKF外推發點散布情況

表1 給出了本文算法與其他算法的比較結果，幾種算法相比，六態_EKF算法圓中間誤差最大，聚類_七態_UKF圓中間誤差最小，定位精度相比原算法提升42.62%，精度最佳，六態_UKF和七態_UKF算法較原算法定位精度也有明顯提升。然而算法耗時角度分析，六態_EKF算法耗時最短，六態_UKF和七態_UKF次之，聚類_七態_UKF耗時最長。

表1 外推算法仿真結果對照表

4 結束語

本文針對一定條件下炮位偵校雷達定位結果一致性較差及定位精度偏低的問題，提出一種基于-均值聚類的彈道外推算法。該算法對單發炮彈軌跡進行多次反向UKF濾波和外推處理，獲取多個發點。然后將聚類思想引入其中，采用-均值聚類及綜合多因子算法獲取最優結果，起到剔除奇異值以及較差發點位置的效果，同時可以進一步消除隨機誤差。仿真結果表明，本文算法顯著提升了彈道外推的定位精度及一致性。本文所采用的外推算法雖然具有較高的精度和一致性，但是對于當前各種新型炮彈的彈道濾波及外推處理技術仍需作進一步研究。