基于剛體彈道模型的防空火控解算方法

何山，吳盤龍，李星秀，惲鵬，薄煜明

(1.南京理工大學 自動化學院, 江蘇 南京 210094; 2.南京理工大學 理學院, 江蘇 南京 210094)

0 引言

隨著網(wǎng)絡化、信息化技術的快速發(fā)展，防空反導作為現(xiàn)代戰(zhàn)爭的重要任務之一，發(fā)揮著越來越重要地位[1]。高炮以其靈活性高、射速快、機動性強和火力密集等優(yōu)點成為了當前低空防御的主力[2]。

火控系統(tǒng)作為武器系統(tǒng)的神經(jīng)中樞，其解算出的射擊諸元精度直接影響武器系統(tǒng)的命中概率[3]。一般來說，火控系統(tǒng)解算的難點在于目標航跡和彈丸軌跡的準確預測[4]。而隨著末端突防能力的日漸提升，目標強機動性使得在航跡的建模和狀態(tài)估計上變得尤為困難；另一方面，復雜的戰(zhàn)場環(huán)境和彈丸的一致性直接影響彈丸零升阻力、誘導阻力、升力以及馬格努斯力等，從而給彈丸的彈道模型建立帶來了極大的挑戰(zhàn)。

現(xiàn)行防空火控系統(tǒng)中的彈道模型多以射表函數(shù)逼近為主，極大縮短了解算時間[5-6]。但在現(xiàn)今信息化作戰(zhàn)的背景下，這一方法存在著精度有限、通用性差、設計過程繁瑣等缺點，尤其是只能對各個彈種的具體彈道和各個修正量分別進行解算，極大地制約了火控系統(tǒng)的開發(fā)[7]。對此，基于彈道微分方程組的火控外彈道模型在該背景下應運而生，它具有通用性好、精度高等特點，但解算處理復雜，耗費時間較多，尤其是對動目標射擊時，需以逆解法[6]為工具，通過試探法不斷迭代求解出命中點的估計值，因此計算量更大。

目前工程應用中所選用的彈道微分方程組主要是質點彈道方程組和修正的質點彈道方程組[8]。質點彈道方程組是將彈丸作為質點處理，忽略了作用在彈丸上的升力、馬氏力和攻角引起的誘導阻力，使彈道計算結果產生了一定的誤差；修正質點彈道模型在質點彈道模型的基礎上綜合了這些因素，有效地提高了計算精度。然而對于高轉速彈丸，修正質點彈道模型在偏流的計算上存在著一定的誤差。對此有學者將彈丸在空中的運動視為一般剛體的運動，考慮彈丸飛行期間產生的進動、章動和自轉等，極大地提高了對彈丸軌跡的預測精度[9]，但復雜的剛體彈道方程解算很難在實際火控系統(tǒng)的一個工作周期內完成。

鑒于此，本文將剛體彈道方程引入到防空火控系統(tǒng)的彈道模型中，極大地提升了火控系統(tǒng)整體的通用性與精確性。同時為了避免剛體彈道方程計算時給火控帶來的實時性問題，本文從火控解算的方法入手，將無偏轉換得到的噪聲協(xié)方差矩陣做解耦，降低了目標狀態(tài)估計時濾波的復雜度；并在割弦迭代法中，將虛擬脫靶量引入到射角和射向的迭代- 修正過程中，減少了火控解算中對彈道方程計算次數(shù)，加快射擊諸元求解的收斂速度。從而在算法執(zhí)行機制上，有效地克服了基于剛體彈道方程火控解算實時性不足的缺陷。

1 彈道建模

由于影響彈丸在空氣中運動的因素極為錯綜復雜，本文針對某型防空高炮，根據(jù)各個因素對彈丸飛行影響的大小，在六自由度剛體彈道方程的基礎上[10],忽略科式慣性力和赤道阻尼力矩；同時為了增大龍格- 庫塔法求解時的步長，忽略彈軸高低角和方位角的2階導項，建立如下彈道方程：

(1)

2 脫靶量

如圖1所示，以武器回轉中心O為坐標原點構建地面基準坐標系Oxyz，其中：Ox軸為射面與炮口水平面的交線，指向射擊方向為正，Oy軸在射擊面內，垂直O(jiān)x軸，向上為正；彈目偏差觀測坐標系TxEyEzE以目標中心T為原點，沿著觀目矢量OT為彈目偏差觀測坐標系的縱深軸yE，過OT的鉛錘面與迎光面(也叫Q平面)的交線為高低軸zE，向上為正，在迎光面內，垂直于zE軸的規(guī)定為方位軸xE，方向滿足右手法則[11-12]；當彈丸以速度v發(fā)射時，記T*為目標T在Oxz面上的投影，Z*為彈跡穿過迎光面的穿越點，Δβ和Δε分別表示觀目矢量OT與觀炸矢量OZ*之間張角的水平分量和高低分量，βh為觀目矢量在東、北、天地理坐標系下的方位角，εh為觀目矢量在東、北、天地理坐標系下的高低角。則彈目偏差TZ*為

TZ*=(xE,0,yE,0,zE,0)=
(|OT|tan Δβ,0,|OT|tan Δε).

(2)

圖1 脫靶量示意圖Fig.1 Schematic diagram of miss distance

而在火控系統(tǒng)中，彈目偏差的數(shù)據(jù)處理都是建立在東、北、天地理坐標系下的[6]。因此，實際地面基準坐標系下的脫靶量Δβ′、Δε′分別為

(3)

3 快速防空火控解算方法

3.1 解耦無偏轉換卡爾曼濾波

為了準確描述目標狀態(tài)隨著時間變化的過程，將目標位置、速度和加速度作為狀態(tài)量，構建目標的運動方程

Xk=Fk-1Xk-1+Wk-1，

(4)

火控系統(tǒng)中對目標狀態(tài)的估計離不開探測設備對目標的跟蹤，在探測系統(tǒng)中，傳感器的量測信息通常是建立在三維球坐標中，目標與傳感器之間的相對位置如圖2所示。圖2中，rk、ηk、βk分別為徑向距離、高低角、方位角。

圖2 傳感器測量坐標系Fig.2 Sensor measurement coordinate system

傳感器的量測值Zk主要包括徑向距離rk、高低角ηk和方位角βk，其測量方程[13-14]為

(5)

式中：Vrk、Vηk和Vβk是相互獨立、均值為0且恒定方差為σrk、σηk和σβk的高斯白噪聲。

球坐標可通過(6)式轉化到笛卡爾坐標系下的量測：

(6)

因此，根據(jù)(5)式的量測值Zk，對真實均值和協(xié)方差矩陣求數(shù)學期望得到無偏量測偏差μk和協(xié)方差Rk[15]：

(7)

(8)

式中：

(9)

3.2 基于剛體彈道模型的快速射擊諸元求解

3.2.1 命中方程建立

假設觀瞄設備中心O′與武器回轉中心O不在同一位置，記J=OO′為觀炮基線，則實際防空武器的命中多邊形如圖3所示。

圖3 命中多邊形Fig.3 Hit polygon

相應的命中矢量方程為

Dh(t)-D(t)-J-Sh(t)=0，

(10)

式中：Dh(t)=D(t+tf)，tf為彈丸飛行時間。

為了計算方便，假設目標在彈丸飛行時間內做勻速直線運動，則

(11)

通過將(10)式中的矢量在東、北、天地理坐標系上進行投影，得到標量命中方程組

(12)

式中：(xoo′,yoo′,zoo′)和(xTh,yTh,zTh)分別為瞄具中心O和命中點Th的坐標；瞄準矢量D(t)=(T,T,T).

由(1)式可知，當射角φh和彈丸飛行時間tf已知時，利用4階龍格- 庫塔法可求得彈丸軌跡上任意一點坐標，即

(13)

再根據(jù)射向方位角βh，命中點坐標可寫為

(14)

由(12)式、(13)式和(14)式即可求得射擊諸元(φh,βh,tf)。

3.2.2 基于虛擬脫靶量的迭代——修正算法

考慮到火控解算的精度與計算速度等要求，可利用割弦迭代法計算出命中點坐標，在每次割弦過程中，需利用多次修正龍格- 庫塔法求解彈道模型時的初值，才能計算出符合當前預測未來點的射擊諸元。因此，本文引入虛擬脫靶量概念，即目標真實航跡與計算機仿真彈丸軌跡間偏差量，對初始射擊射角和射向進行修正，極大地提高了計算當前預測未來點射擊諸元的收斂速度。具體求解過程如下：

步驟2割弦迭代法初值的選取。假設在彈丸飛行時間內目標處于靜止狀態(tài)，即目標運動時間Tf=0 s，則命中點坐標為(xTh,k,yTh,k,zTh,k)，有

(15)

這時可將火控解算問題轉化為一個外彈道學中的兩點邊值問題，即從上述邊值條件出發(fā)，求出射擊諸元。

步驟3彈道方程初值給定。根據(jù)所給的目標未來點坐標預估初始射角φh，0和射向方位角βh，0:

(16)

步驟4彈道方程求解。利用4階龍格- 庫塔法求解微分方程組(1)式，得到地面基準坐標系下彈道諸元的序列值：

{xi,yi,zi,vx,i,vy,i,vz,i,tf,i},i=0,1,2,…，

(17)

式中：i為龍格- 庫塔法定步長求解時的序列號。

步驟5虛擬脫靶量計算。

1)將步驟4得到的彈道諸元序列通過(18)式轉化到彈目偏差觀測坐標系中，

(18)

2)利用拉格朗日插值，可計算出在Q平面內的二維彈目偏差(xE,0,zE,0)和彈丸飛行時間tf,i,0，利用(3)式將(xE,0,zE,0)轉換到東、北、天地理坐標系下的虛擬脫靶量(Δβ′0,Δε′0)。

步驟6射角與射向迭代修正。判斷(Δβ′0,Δε′0)是否滿足系統(tǒng)所需的精度(一般取0.1 mil)。若滿足，則射擊諸元為φh=φh,0，βh=βh,0，tf=tf,i,0，并轉到步驟7；否則對初始射角和射向進行修正，得到新的射角和射向方位角：

(19)

將修正后的射角和射向取代原來的變量，并返回步驟4.

至此，一個周期內的火控解算完成，相比傳統(tǒng)火控解算的迭代修正算法，步驟5中利用虛擬脫靶量進行射角和射向的迭代修正，減少了火控解算中對彈道方程計算次數(shù)，極大提高了基于剛體彈道方程的防空火控解算方法求解速度。

4 仿真分析

4.1 仿真場景

本文以某型防空高炮為研究對象，構建三維仿真場景，對所提火控解算方法進行仿真驗證。假設目標做勻速運動，起始位置為(3 500 m, 3 500 m, 1 000 m)，初始速度為(-300 m/s, -200 m/s, -25 m/s)；測量傳感器采用雷達，采樣周期為0.02 s，距離誤差為5 m，方位和俯仰角誤差為1 mrad，目標在雷達界面的顯示如圖4所示。

圖4 目標在雷達界面顯示Fig.4 Target on radar interface

彈丸初速v0=1 050 m/s，龍格- 庫塔法的解算步長為0.01 s，在標準氣象條件下，不同射角的彈丸軌跡如圖5所示。

圖5 彈丸軌跡Fig.5 Projectile trajectory

4.2 仿真分析

4.2.1 狀態(tài)濾波

在實際工程中，通過傳感器測量的目標運動參數(shù)包含眾多干擾噪聲，這些噪聲不僅影響系統(tǒng)的計算精度，嚴重時系統(tǒng)將失穩(wěn)。因此對目標進行狀態(tài)濾波是準確預測目標未來點位置的關鍵技術。本文通過無偏轉換將雷達對目標運動狀態(tài)的非線性量測方程圍繞目標運動狀態(tài)進行展開，并將無偏轉換得到的噪聲協(xié)方差矩陣做解耦，降低了濾波算法計算的復雜度，其位置均方根誤差(RMSE)、速度RMSE和濾波處理時間如圖6、圖7和圖8所示。由圖6～圖8可以看出，解耦無偏轉換卡爾曼濾波與無偏轉換卡爾曼濾波算法[15]的誤差差值不大于總誤差的1%，而濾波時間減少了15%. 因此，該方法具有較高的精度和良好的性能。

圖6 位置RMSEFig.6 RMSEs of position

圖7 速度RMSEFig.7 RMSEs of velocity

圖8 濾波的處理時間Fig.8 Processing time of filtering

4.2.2 火控解算誤差

根據(jù)上述仿真場景，以2 s穩(wěn)定建航后，分別用傳統(tǒng)的火控解算方法[6]和本文所提方法對不同時間點進行火控解算，其解算點的射擊諸元對比如表1所示，火控解算的誤差如圖9所示。

從圖9仿真結果可以看出，本文所提方法的有效性和可行性。由表1可以看出：當目標在較遠處，其解算誤差很大，但隨著目標越來越接近，彈丸飛行時間越來越短，誤差呈現(xiàn)減小趨勢，這是由于測量傳感器所得目標數(shù)據(jù)存在著一定的量測誤差，導致解算誤差大部分來源于濾波誤差；當目標距離3 000 m以內時，單軸方向上的解算誤差不大于0.5 m，該方法可以很好地滿足實際防空作戰(zhàn)需求。

表1 射擊諸元對比Tab.1 Comparison of firing data

圖9 火控解算誤差Fig.9 Errors of fire control calculation

4.2.3 解算時間

本文通過對不同距離目標的火控解算進行仿真實驗來有效地說明所提出方法的實時性優(yōu)勢。在Matlab2014環(huán)境下(Windows 7, i5-3210M CPU, 2.5 GHz主頻, 8 GB RAM)將傳統(tǒng)的火控解算方法[6]與本文方法進行比較，解算時間對比如表2所示。

表2 Windows平臺下的解算時間對比Tab.2 Comparison of calculation times based on Windows s

從表2中可以看出，本文所提方法在計算量上遠遠低于傳統(tǒng)的火控解算方法[6]，這是由于火控諸元的解算主要包括狀態(tài)濾波與命中點解算兩大步驟。對于狀態(tài)濾波，通過將無偏轉換卡爾曼濾波的噪聲協(xié)方差做解耦，一定程度上降低了濾波計算的復雜度；對于命中點解算，本文提出基于虛擬脫靶量的迭代——修正算法，利用虛擬脫靶量進行射角和射向的迭代修正，減少了火控解算中對彈道方程計算次數(shù)，大幅度提高火控求解速度。

為了有效說明該方法在工程時是可行的，將代碼改編成C語言并移植到配置為i7-4790CPU、SylixOS實時操作系統(tǒng)的工控機上運行，解算時間對比如表3所示。從表3實驗的數(shù)據(jù)中可以看出，本文所提方法在工程上具有可觀的應用前景。

5 結論

為了提高防空火控解算的精確性與實時性，本文提出了基于剛體彈道方程的防空火控解算方法，仿真結果證明了該方法的有效性和可行性。然而，對于實際的防空系統(tǒng)，目標的強機動特性使得基于傳統(tǒng)運動模型的濾波算法出現(xiàn)精度下降甚至發(fā)散現(xiàn)象，因此后續(xù)的工作是考慮在該類情況下，提高目標狀態(tài)估計與軌跡預測的準確性，并在一定的基礎上降低算法的復雜度，使得該算法可用于未來防空武器的火控系統(tǒng)中。

表3 SylixOS下的解算時間對比Tab.3 Comparison of calculation times based on SylixOS ms