李同亮 ， 朱 勇 ， 于 瓊

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

0 引言

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

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

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

1 彈道運動模型

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

(1)

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

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

圖1 炮口坐標系

2 彈道外推

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

圖2 彈道外推流程圖

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

2.1 七態(tài)UKF濾波

UKF濾波算法對非線性問題的處理效果顯著且該方法估計精度至少達到2階，計算量與EKF 算法處于同一量級，不需要計算Jacobian 矩陣，可以處理不可導的非線性函數(shù)?；谏鲜鰞?yōu)點以及彈道運動方程非線性特性，本文采用七態(tài)UKF濾波方法，將彈道系數(shù)作為第七維狀態(tài)進行濾波，實時估計并更新彈道系數(shù)。

基于第1節(jié)中彈道運動模型建立狀態(tài)方程和量測方程。

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

(2)

式中，為高斯白噪聲，協(xié)方差矩陣為。

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

=

(3)

(4)

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

考慮如下非線性模型：

(5)

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

1) 狀態(tài)初始化。

(6)

(7)

(8)

(9)

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

(10)

4) 濾波更新。

(11)

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

各采樣點權(quán)值計算：

(12)

2.2 彈道解算

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

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

0(=1,2,…,)

(13)

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

(14)

其中：

(15)

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

2.3 K-均值聚類處理

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

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

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

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

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

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

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

(,)=min{(,)},

=1,2,…,

(16)

則表明屬于類。

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

(17)

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

(18)

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

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

簇品質(zhì)_：用于表示簇優(yōu)劣程度的量稱為簇品質(zhì)。

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

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

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

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

_=[___][]

(19)

式中，_=[_,_,…,_]為簇品質(zhì)向量，,,分別表示歸一化簇體積因子、歸一化簇質(zhì)量因子及歸一化簇一致性因子對應的權(quán)值。最優(yōu)簇的品質(zhì)即為

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

_∈[1,]

(20)

式中，_對應的簇即為最優(yōu)簇。

3 算法仿真與分析

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

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

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

3) 雷達距離炮48 km。

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

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

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

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

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

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

圖6 七態(tài)_UKF算法外推發(fā)點散布情況

圖7 聚類_七態(tài)_UKF外推發(fā)點散布情況

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

表1 外推算法仿真結(jié)果對照表

4 結(jié)束語

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