李同亮 ， 朱 勇 ， 于 瓊

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

0 引言

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

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

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

1 彈道運(yùn)動(dòng)模型

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

(1)

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

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

圖1 炮口坐標(biāo)系

2 彈道外推

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

圖2 彈道外推流程圖

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

2.1 七態(tài)UKF濾波

UKF濾波算法對(duì)非線性問(wèn)題的處理效果顯著且該方法估計(jì)精度至少達(dá)到2階，計(jì)算量與EKF 算法處于同一量級(jí)，不需要計(jì)算Jacobian 矩陣，可以處理不可導(dǎo)的非線性函數(shù)。基于上述優(yōu)點(diǎn)以及彈道運(yùn)動(dòng)方程非線性特性，本文采用七態(tài)UKF濾波方法，將彈道系數(shù)作為第七維狀態(tài)進(jìn)行濾波，實(shí)時(shí)估計(jì)并更新彈道系數(shù)。

基于第1節(jié)中彈道運(yùn)動(dòng)模型建立狀態(tài)方程和量測(cè)方程。

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

(2)

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

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

=

(3)

(4)

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

考慮如下非線性模型：

(5)

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

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

(6)

(7)

(8)

(9)

②公式(10)描述了一步提前預(yù)測(cè)計(jì)算過(guò)程。

(10)

4) 濾波更新。

(11)

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

各采樣點(diǎn)權(quán)值計(jì)算：

(12)

2.2 彈道解算

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

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

0(=1,2,…,)

(13)

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

(14)

其中：

(15)

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

2.3 K-均值聚類處理

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

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

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

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

輸出：個(gè)聚類{,,…,};

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

2) 計(jì)算,,…,這個(gè)中心點(diǎn)與數(shù)據(jù)集中其他數(shù)據(jù)之間的距離，記作(,)，其中，= 1,2,…,，= 1,2,…,。當(dāng)

(,)=min{(,)},

=1,2,…,

(16)

則表明屬于類。

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

(17)

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

(18)

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

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

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

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

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

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

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

_=[___][]

(19)

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

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

_∈[1,]

(20)

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

3 算法仿真與分析

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

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

2) 取雷達(dá)隨機(jī)測(cè)量誤差:= 76 m，= 17 mil，= 1.6 mil;

3) 雷達(dá)距離炮48 km。

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

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

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

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

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

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

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

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

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

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

4 結(jié)束語(yǔ)

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