嚴軍，劉奇濤

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

炮位偵察校射雷達的主要作戰(zhàn)使命是對敵方炮位的偵察和己方火炮的校射，通過快速有效地完成火力偵察支撐快速打擊任務，在現(xiàn)代戰(zhàn)爭中發(fā)揮著重要的作用[1-5]。炮位偵校雷達通過搜索發(fā)現(xiàn)彈丸目標進行跟蹤，然后進行彈道處理外推炮位和落點。彈道外推的精度決定著炮位偵校雷達的性能，而運動模型直接影響彈道外推的精度。國內(nèi)外傳統(tǒng)的炮位偵校雷達數(shù)據(jù)處理方法通?；谥亓ζ矫婺Ｐ停雎粤似渌蛩貙椀赖挠绊?，因此僅適用于中、近程火炮的偵察校射。近年來，隨著各種遠程火箭彈、增程彈技術的突飛猛進[6]，平面重力模型已不能很好地滿足精度需求，迫切需要采用高精度彈道模型。

1 彈道模型研究

1.1 理論模型

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

1.1.1 平面運動模型

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

(1)

1.1.2 考慮地球曲率的運動模型

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

(2)

1.1.3 考慮地球自轉(zhuǎn)的運動模型

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

(3)

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

1.1.4 考慮地球J2引力的運動模型

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

(4)

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

1.2 仿真分析

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

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

射角：35°；

射向：45°；

初速：1 000 m/s。

1.2.1 地球曲率的影響分析

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

可以看到，在給定的仿真條件下，采用平面模型和曲面模型，位置預報差異可達50 m左右。

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

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

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

1.2.3J2項引力的影響分析

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

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

2 參數(shù)估計方法

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

2.1 最小二乘擬合

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

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

(5)

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

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

(6)

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

2.2 擴展卡爾曼濾波

2.2.1 狀態(tài)方程

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

(7)

(8)

f(θ)的具體表達式由選擇的運動模型決定，該非線性方程只是目標真實運動模型的近似，有一定誤差，為了補償這個誤差，引入一個模型噪聲W，則變?yōu)?/p>

(9)

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

2.2.2 觀測方程

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

(10)

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

Z=h(θ)+V，

(11)

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

(12)

2.2.3 擴展的卡爾曼濾波器

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

E(Vk)=0，

(13)

(14)

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

E(Wk)=0，

(15)

(16)

擴充的卡爾曼濾波方程：

1)預測方程為

(17)

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

(18)

3)觀測量預測方程為

(19)

4)濾波方程為

(20)

式中，Wk+1為加權矩陣，

(21)

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

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

(22)

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

2.3 各種方法比較

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

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

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

3 結(jié)束語

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