嚴(yán)軍，劉奇濤

(1.南京電子技術(shù)研究所，江蘇 南京 210000；2.西北機(jī)電工程研究所，陜西 咸陽 712099)

炮位偵察校射雷達(dá)的主要作戰(zhàn)使命是對敵方炮位的偵察和己方火炮的校射，通過快速有效地完成火力偵察支撐快速打擊任務(wù)，在現(xiàn)代戰(zhàn)爭中發(fā)揮著重要的作用[1-5]。炮位偵校雷達(dá)通過搜索發(fā)現(xiàn)彈丸目標(biāo)進(jìn)行跟蹤，然后進(jìn)行彈道處理外推炮位和落點(diǎn)。彈道外推的精度決定著炮位偵校雷達(dá)的性能，而運(yùn)動(dòng)模型直接影響彈道外推的精度。國內(nèi)外傳統(tǒng)的炮位偵校雷達(dá)數(shù)據(jù)處理方法通常基于重力平面模型，忽略了其他因素對彈道的影響，因此僅適用于中、近程火炮的偵察校射。近年來，隨著各種遠(yuǎn)程火箭彈、增程彈技術(shù)的突飛猛進(jìn)[6]，平面重力模型已不能很好地滿足精度需求，迫切需要采用高精度彈道模型。

1 彈道模型研究

1.1 理論模型

箭彈的運(yùn)動(dòng)受力非常復(fù)雜，精確的計(jì)算需要同時(shí)考慮質(zhì)心和姿態(tài)的運(yùn)動(dòng)[7]，這需要知道很多彈體相關(guān)的先驗(yàn)信息，對于箭彈的外彈道設(shè)計(jì)來說可以進(jìn)行詳細(xì)的考慮，但是對于炮位雷達(dá)非合作目標(biāo)的彈道外推來說，雷達(dá)觀測條件的限制決定了只能對質(zhì)心運(yùn)動(dòng)進(jìn)行近似的考慮，因此炮位偵校雷達(dá)中通常廣泛地采用簡化的彈道方程描述彈丸質(zhì)心的運(yùn)動(dòng)。

1.1.1 平面運(yùn)動(dòng)模型

傳統(tǒng)炮位偵校雷達(dá)的數(shù)據(jù)處理為了簡化計(jì)算，不考慮地球曲面的影響，采用了簡單的平面模型，并在平面模型中將地球引力取為常值，同時(shí)考慮大氣阻尼和風(fēng)的影響。其運(yùn)動(dòng)方程如下[8-9]：

(1)

1.1.2 考慮地球曲率的運(yùn)動(dòng)模型

彈丸是在地球表面的空間運(yùn)動(dòng)的，彈丸在運(yùn)動(dòng)過程中的地心引力是指向地心的，由于地球表面是復(fù)雜的曲面，因此作用在彈丸上的地心引力指向隨著彈丸運(yùn)動(dòng)不斷變化。為了更好地描述彈丸的運(yùn)動(dòng)，應(yīng)當(dāng)選擇在地心直角坐標(biāo)系中描述彈丸的運(yùn)動(dòng)，一般可以選擇常用的WGS84坐標(biāo)系，相應(yīng)的運(yùn)動(dòng)模型如下[10-11]：

(2)

1.1.3 考慮地球自轉(zhuǎn)的運(yùn)動(dòng)模型

平面模型和考慮地球曲率的模型忽略了地球的自轉(zhuǎn)，實(shí)際情況下地球是存在自轉(zhuǎn)的，地球自轉(zhuǎn)會產(chǎn)生離心力和柯氏力?？紤]到這些力的影響后，在WGS84坐標(biāo)系下，相應(yīng)運(yùn)動(dòng)模型如下：

(3)

式中，we為地球自轉(zhuǎn)角速度。

1.1.4 考慮地球J2引力的運(yùn)動(dòng)模型

地球不是完美的球體，更準(zhǔn)確地描述應(yīng)該是一個(gè)定扁率的橢球體，這導(dǎo)致地球引力對彈丸的影響在質(zhì)心引力外存在附加引力，其中最大的一項(xiàng)為J2項(xiàng)非球形引力。對于彈丸的運(yùn)動(dòng)來說，只需要考慮J2項(xiàng)非球形引力就足夠了，相應(yīng)的運(yùn)動(dòng)模型如下：

(4)

式中：J2為地球非球形引力系數(shù)；ae為地球赤道半徑。

1.2 仿真分析

為了比較上述不同模型的差異，選取一射程為35 km的彈丸進(jìn)行仿真，仿真中采用相同的初始條件，比較不同模型預(yù)報(bào)的軌道差異，仿真參數(shù)如下：

炮位位置：經(jīng)度45°，緯度45°；

射角：35°；

射向：45°；

初速：1 000 m/s。

1.2.1 地球曲率的影響分析

圖1給出的是上述仿真條件下，采用平面模型和考慮地球曲率的運(yùn)動(dòng)模型，軌道預(yù)報(bào)結(jié)果在地心直角坐標(biāo)系中3個(gè)坐標(biāo)分量的差異。

可以看到，在給定的仿真條件下，采用平面模型和曲面模型，位置預(yù)報(bào)差異可達(dá)50 m左右。

1.2.2 地球自轉(zhuǎn)的影響分析

圖2給出的是上述仿真條件下，采用球面模型和同時(shí)考慮地球自轉(zhuǎn)的運(yùn)動(dòng)模型，軌道預(yù)報(bào)結(jié)果在地心直角坐標(biāo)系中3個(gè)坐標(biāo)分量的差異。

可以看到，考慮地球自轉(zhuǎn)和不考慮地球自轉(zhuǎn)，預(yù)報(bào)結(jié)果差異在150 m量級。

1.2.3J2項(xiàng)引力的影響分析

圖3給出的是上述仿真條件下，采用考慮地球自轉(zhuǎn)的運(yùn)動(dòng)模型和同時(shí)考慮地球J2項(xiàng)引力運(yùn)動(dòng)模型，軌道預(yù)報(bào)結(jié)果在地心直角坐標(biāo)系中3個(gè)坐標(biāo)分量的差異。

可以看到，考慮地球J2項(xiàng)引力和不考慮地球J2項(xiàng)引力，預(yù)報(bào)結(jié)果差異在30 m量級。

2 參數(shù)估計(jì)方法

建立彈道運(yùn)動(dòng)方程后，為了進(jìn)行彈道外推，首先要進(jìn)行彈道參數(shù)估計(jì)。由前面的介紹可知，目標(biāo)運(yùn)動(dòng)方程是非線性的，觀測方程一般也是非線性的，因此是一個(gè)非線性參數(shù)估計(jì)問題。炮位雷達(dá)探測目標(biāo)時(shí)獲得的觀測數(shù)據(jù)帶有誤差，常用的非線性參數(shù)估計(jì)方法有最小二乘估計(jì)、擴(kuò)展卡爾曼濾波(EKF)、無跡卡爾曼濾波(UKF)和粒子濾波(PF)等[12-14]。通過參數(shù)估計(jì)得到彈道目標(biāo)的運(yùn)動(dòng)參數(shù)，利用彈道運(yùn)動(dòng)方程采用成熟的數(shù)值積分方法即可進(jìn)行彈道外推，獲得炮位和落點(diǎn)等參數(shù)。

2.1 最小二乘擬合

最小二乘估計(jì)是一種常用的參數(shù)估計(jì)方法，以預(yù)測數(shù)據(jù)與觀測數(shù)據(jù)的殘差平方和最小為目標(biāo)函數(shù)來進(jìn)行未知參數(shù)估計(jì)。假設(shè)測量值Z與估計(jì)參數(shù)X滿足關(guān)系：

Z(j)=h(j,X)+V(j),j=1,2,…，

(5)

式中，V(j)為測量噪聲。

參數(shù)估計(jì)結(jié)果X(k)為使得k次觀測擬合殘差平方和達(dá)到最小的解，即

(6)

炮彈飛行軌跡可以近似為一拋物線，首先利用雷達(dá)測量數(shù)據(jù)對拋物線的參數(shù)進(jìn)行最小二乘估計(jì)，然后用估計(jì)得到的參數(shù)曲線即可外推算出彈丸發(fā)射點(diǎn)或落點(diǎn)。

2.2 擴(kuò)展卡爾曼濾波

2.2.1 狀態(tài)方程

取x，y，z，vx，vy，vz，c作為卡爾曼濾波的狀態(tài)變量。

(7)

(8)

f(θ)的具體表達(dá)式由選擇的運(yùn)動(dòng)模型決定，該非線性方程只是目標(biāo)真實(shí)運(yùn)動(dòng)模型的近似，有一定誤差，為了補(bǔ)償這個(gè)誤差，引入一個(gè)模型噪聲W，則變?yōu)?/p>

(9)

式中，W為零均值高斯白噪聲，W～N(0,Q)。

2.2.2 觀測方程

設(shè)雷達(dá)量測值為斜距R，方位角β，高低角ε，雷達(dá)球坐標(biāo)與直角坐標(biāo)之間的轉(zhuǎn)換方程為

(10)

令觀測矢量Z=(R,β,ε)T，則有

Z=h(θ)+V，

(11)

式中：V為雷達(dá)測量噪聲，假定為零均值高斯白噪聲，V～N(0,R)；h(θ)為三維矢量函數(shù)：

(12)

2.2.3 擴(kuò)展的卡爾曼濾波器

設(shè)離散化的量測噪聲Vk有如下均值和方差：

E(Vk)=0，

(13)

(14)

離散形式的狀態(tài)擾動(dòng)Wk有如下均值和方差：

E(Wk)=0，

(15)

(16)

擴(kuò)充的卡爾曼濾波方程：

1)預(yù)測方程為

(17)

2)預(yù)測協(xié)方差為

(18)

3)觀測量預(yù)測方程為

(19)

4)濾波方程為

(20)

式中，Wk+1為加權(quán)矩陣，

(21)

5)濾波協(xié)方差：

Pk+1=[I-Wk+1·Hk+1]·Pk+1|k.

(22)

通常EKF是逐漸收斂的過程，濾波終點(diǎn)的參數(shù)估計(jì)精度高于中間點(diǎn)，一般采用濾波終點(diǎn)的參數(shù)進(jìn)行彈道外推，為了使外推距離更短，可以從數(shù)據(jù)末點(diǎn)開始進(jìn)行反向?yàn)V波，有利于提高精度。

2.3 各種方法比較

參數(shù)估計(jì)方法中，最小二乘估計(jì)是一種常用的傳統(tǒng)方法，由于無法估計(jì)彈道系數(shù)等其他彈道參數(shù)，外推精度相對較差，它的優(yōu)勢是在數(shù)據(jù)點(diǎn)很少的情況下依然可以得到外推結(jié)果。

EKF算法計(jì)算量較小，但是需要計(jì)算雅可比矩陣，對于復(fù)雜的非線性方程，雅可比矩陣的計(jì)算比較復(fù)雜。EKF算法適用于弱非線性問題,在非線性較大的場景，可能會出現(xiàn)濾波發(fā)散的情況。

除了上面的方法，無跡卡爾曼濾波(UKF)算法和粒子濾波(PF)也被應(yīng)用到彈道參數(shù)估計(jì)中。UKF算法適用于各種非線性高斯噪聲問題，計(jì)算量比EKF大，對非線性問題的適應(yīng)性更好，其濾波精度和外推精度都比較好。PF算法由于需要采樣大量的粒子，計(jì)算量明顯較大。對于高度非線性問題，這兩者較EKF方法有一定的優(yōu)勢，但對于火箭彈的彈道外推來說，非線性程度還無法體現(xiàn)兩者的優(yōu)勢，采用EKF、UKF和PF沒有本質(zhì)區(qū)別。

3 結(jié)束語

彈丸在運(yùn)動(dòng)過程中受到多種力的影響，通過建立考慮不同因素的運(yùn)動(dòng)模型，采用相同初值進(jìn)行軌道預(yù)測并比較。仿真表明，在給定仿真條件下，地球曲率的影響在50 m左右、地球自轉(zhuǎn)的影響在150 m左右、J2項(xiàng)引力的影響都在30 m左右。對于射程更遠(yuǎn)的彈丸，這些力的影響可能更大，隨著裝備的發(fā)展，這些力的影響不可忽略，需要進(jìn)行考慮。