基于變步長龍格庫塔法的高炮C-RAM解命中方法

劉新宇，舒立鵬，林智偉，吳 曄，朱柏飛，唐 旭

(西北機電工程研究所， 陜西 咸陽 712099)

1 引言

C-RAM(counter-rocket artillery mortar)是一種使用我方防空火力主動攔截來襲的火箭彈、榴彈和迫擊炮彈，以達到防御敵方地面曲射火力打擊的主動防護系統[1]。高炮C-RAM火控系統的任務是根據雷達對來襲RAM目標的跟蹤數據估計出其實時位置、速度并辨識其彈道系數，根據以上信息通過解命中計算得出目標的未來命中點并射擊，從而實現對RAM目標的攔截。但在火控處理時，由于RAM目標的運動模式的非線性性[2]，傳統的高炮解命中方法因假定模型以及計算方式的問題無法得到正確的解命中結果。由此需要設計一種基于RAM目標運動方程假定，通過積分方法計算未來點的高炮C-RAM解命中方法，以實現高炮C-RAM作戰的需要。

2 RAM目標運動特性

由于在實際進行C-RAM作戰時目標彈丸轉速、攻角、彈翼升力等狀態無法測得，這些狀態對于彈丸運動軌跡在短時間內的影響較小，在進行短時間的外彈道外推時相對于彈道系數估計誤差造成的未來點偏移其影響基本可以忽略不計，故本文采用彈丸質心運動模型作為外推RAM目標外彈道軌跡的公式，將除去空氣阻力與重力外的其他力對彈丸運動的影響看作是彈道系數誤差所造成的影響的一部分或是將其忽略。基于C-RAM作戰時一般在RAM目標的外彈道末端進行攔截，此時目標速度通常處于亞音速狀態，彈道系數值的變化幅度較小，故可以將彈道系數視為常值進行處理，并假設氣象條件是標準的，無風雨干擾。彈丸質心運動模型示意圖如圖1。

圖1 彈丸質心運動模型示意圖

彈丸質心運動方程組如式(1)所示:

(1)

式中:x、y、z、vx、vy、vz分別為彈丸的位置坐標與速度分量;c為彈道系數；g為重力加速度；cs為聲速；H(z)為空氣密度特征函數；G(v,cs)為空氣阻力函數[5-7]。

3 RAM類目標火控解命中方法

3.1 龍格庫塔默森法

(2)

其中：

(3)

誤差估計值計算公式如式(4)所示：

(4)

此誤差估計值的精度為3階，而此方法的積分精度為4階，對于大部分使用情況已經能夠滿足精度要求。

3.2 命中方程求解方法

(5)

(6)

以上命中點積分方法使得在進行彈道外推時預測命中點不會超過實際命中點過多，從而最大程度上的減少積分次數。

3) 命中點水平距離和高度求解:

(7)

(8)

4) 彈丸飛行時間求解:

(9)

圖2 結合龍格庫塔積分方法的迭代法收斂曲線

3.3 固定步長龍格庫塔默森法下積分誤差估計值及諸元解算誤差的統計經驗

總結出以下結論：

1) 解命中誤差與積分步長成正比；

2) 對于擁有同一初速度、速度方向、高度、彈道系數但高炮與彈丸軌跡水平投影斜距不同的彈道，解算誤差除以彈丸飛行時間的值與積分步長為線性關系，且斜率相對一致；

3) 不同斜距的外彈道軌跡下，x軸與y軸位置坐標估計誤差與積分步長的3次方的比值與積分步長為斜率相對一致的線性關系；

4) 不同斜距的外彈道軌跡下，z軸位置坐標估計誤差與積分步長的2次方的比值與積分步長為斜率相對一致的線性關系；

5) 同一條彈道在進行誤差估計時在在積分步長一致的情況下估計誤差的變化較小。

各項結論所總結的關系如圖3—圖6所示。

圖3 x軸解算誤差與彈飛時間的比值曲線

圖4 x軸估計誤差與步長3次方的比值曲線

圖5 z軸估計誤差除以步長2次方值的曲線

圖6 同一彈道在定步長下不同次積分的估計誤差變化

對不同彈道進行仿真后確認了以上統計經驗具有廣泛性，故可以將實際誤差、誤差估計值、積分步長之間的關系總結為如式(10)所示的公式。

(10)

式(10)中：Cx、Cy、Cz為各軸的誤差換算比值常數；Ex、Ey、Ez為各軸位置估計誤差；ht為積分步長；Ex_ture、Ey_ture、Ez_ture為實際位置解算誤差；tf為彈丸飛行時間。

在對多種彈丸的彈道進行大量仿真驗證后發現，在不同彈道系數及彈丸速度下，誤差換算比值常數Cx、Cy、Cz是不同的。對一飛行軌跡起點固定的彈丸進行定步長解命中仿真，設定彈道系數為0.5，改變每次計算時x軸的初始速度，x軸誤差換算比值常數隨積分步長與x軸速度變化曲面如圖7。積分步長為0.05 s時，x軸誤差換算比值常數在不同彈道系數下隨x軸速度變化的曲面如圖8。

圖7 不同速度、積分步長下的誤差換算比值常數變化曲面

圖8 不同速度、彈道系數下的誤差換算比值常數變化曲面Fig.8 Variation curve of error conversion ratio constant under different speeds and ballistic coefficient

從圖7、圖8的結果中可以看出該比值常數在不同速度以及彈道系數下是不同的，變化規律也呈非線性，但不隨積分步長變化。在使用式(10)估計實際解算誤差時，比值常數的取值將決定估計出的實際解算誤差準確與否。

對于比值常數的取值問題，在實際應用時的解決方式有3種：

1) 通過大量的實驗來制作對應不同彈道系數、速度的比值常數表，在計算時根據目標的彈道系數以及速度從表中取值；

2) 將比值變化曲線近似線性化，在計算時根據近似函數得出比值用于計算；

3) 取一固定值，保證大部分高概率情況的解算效果，犧牲小部分低概率情況的解算精度和速度。

(11)

其中：

(12)

其誤差估計值計算公式如式(13)所示：

(13)

夏普法中K5正好是下一次積分計算時的K1，因此只需要在第一步多算一次K1，之后每步積分只需要計算4次右端函數。其誤差換算公式如式(14)所示，形式與默森法基本相同，但積分步長次方值不同。

(14)

根據式(14)計算出的夏普法誤差換算比值常數變化曲面如圖9、圖10。可以看出相比于默森法的比值常數變化曲線，夏普法的比值常數在速度和彈道系數較大時，同一速度、彈道系數下其不同積分步長對應比值常數變化較大。故在速度和彈道系數較大的情況下，若要進行最優步長控制，夏普法的步長控制難度要大于默森法，故相比于默森法，夏普法并不適用于C-RAM的解命中計算。

圖9 彈道系數=0.2，不同速度、積分步長下夏普法的誤差換算比值常數變化曲面

圖10 彈道系數=0.5，不同速度、積分步長下夏普法的誤差換算比值常數變化曲面

3.4 基于最優步長控制策略的龍格庫塔默森法

具體的最優步長控制策略為：

① 首先以一較為合適的固定步長hm進行第一次解命中，求得彈丸飛行時間tf。

② 給定誤差限制ε0，以固定步長hm作為第一次積分步長，以上一次解命中的彈飛時間tf作為本次解算參考彈飛時間，求若以本次積分步長進行解命中時各軸諸元解算的誤差值εx、εy、εz，公式為式(15)所示。

(15)

若εx≤ε0且εy≤ε0且εz≤ε0，則本次積分成功，通過式(16)確定下一步步長hm+1。

(16)

若εx>ε0或εy>ε0或εz>ε0，則本次積分誤差過大，同樣通過式(16)求出一個積分步長hm+1。為了盡量不重復進行積分運算，減小計算時間，在一次諸元解算中除第一次積分外均無需重新積分。在同一次解算中同一積分步長造成的誤差相差不大，除第一次積分外若出現計算誤差值大于誤差限制值的情況，只需要減小下一次積分步長保證下一次積分的精度即可。若對解算精度的要求較高，則只需要將誤差限制值設置的比理想誤差值稍小就能控制每步積分的誤差在范圍內。

③設置最小步長hmin，一旦積分步長小于此步長值則不再適用步長控制，此舉是為了避免計算機在步長過小時由于舍入誤差造成數值偏差比過大，從而計算出一個過大的誤差值，引起步長振蕩導致解算失敗。

4 實例仿真

為驗證基于統計特征的最優步長控制策略是否能夠有效控制積分步長使解命中誤差保持在可接受范圍，對某2個RAM目標進行解命中計算，假設目標一被探測到時的彈道起點為x=500 m、y=1 000 m、z=1 000 m,高炮位于坐標原點，彈丸速度向量為x軸方向-250 m/s、y軸方向-250 m/s、z軸方向-90 m/s，彈道系數為0.2。

目標二被探測到時的彈道起點為x=500 m、y=1 000 m、z=1 500 m,高炮位于坐標原點，彈丸速度向量為x軸方向-200 m/s、y軸方向-150 m/s、z軸方向-150 m/s，彈道系數為0.3。

假設解命中過程中所得到的目標運動狀態以及彈道系數均為理想值，每次諸元解算間隔為0.01 s，共解算600次。基于變步長默森法的解命中法為方法一，基于變步長夏普法的解命中法為方法2，文獻[2]中的定步長解命中法為方法3，變步長法設置誤差限制為各軸位置誤差小于0.1 m，方法1比值常數取值為Cx=Cy=14 000、Cz=110；方法2比值常數取值為Cx=Cy=2 200、Cz=20；定步長法設置積分步長為0.05 s。

分別使用3種解命中方法對目標一及目標二進行仿真。諸元解算后進行逆解得到的解算誤差如圖11—圖13所示。

圖11 x軸位置解算誤差圖

圖12 y軸位置解算誤差圖

圖13 z軸位置解算誤差圖

從仿真結果中可以看出，相對于定步長的C-RAM解命中的方法，變步長法能夠使各軸方向上的解算誤差基本保持在設定的誤差范圍內。但由于比值常數使用了取固定值的方法，導致對誤差的控制精度較差，解算誤差超出設置范圍的情況，若是采取使用比值常數表或者擬合比值變化曲線的方法，將能進一步的提升誤差控制的精度。

在變步長法中，相比默森法的誤差控制效果，夏普法在對目標一解命中時效果與默森法接近，對目標二解命中時誤差小于默森法。但誤差較小并不意味著控制效果好，夏普法與默森法的積分精度同為四階且誤差估計精度也同為三階，在相同的精度條件下，更小的解算誤差代表消耗了更多的解算時間，證明了在積分步長控制效果上夏普法要弱于默森法。且仿真中發現對夏普法應用本文中的最優步長控制時不能將上一次計算出的K5代入下一次計算中，否則將會因為步長的變化導致誤差估計值的計算出現錯誤，進而無法有效控制積分步長，故在每步計算量上夏普法與默森法也基本相同，使其不具備計算量的優勢。

變步長默森法積分步長的控制效果如圖14、圖15所示。

圖15 變步長法解算目標二時積分步長控制云圖

由圖14、圖15所展示的積分步長控制情況中可以看出，本文中的變步長C-RAM解命中方法能夠實現對積分步長的最優控制。其能夠實時的調節積分步長的大小，在積分誤差較大的情況下減小步長保證精度，積分誤差較小時增加步長減少積分次數。實現在保證解命中精度的情況下盡可能增加積分步長的長度，減少計算時間。

5 結論

1) 在基于彈丸質心運動模型模擬理想RAM目標外彈道軌跡的情況下，基于變步長龍格庫塔法的高炮C-RAM解命中方法可以完成高炮C-RAM解命中任務。

2) 基于C-RAM解命中誤差與積分估計誤差規律應用的最優步長控制策略能夠完成對積分步長的最優控制，相比使用定步長積分法進行C-RAM解命中計算，該方法可控制目標在不同距離、速度、彈道系數情況下的諸元解算誤差。

3) 在同樣應用最優步長控制方法同樣精度的變步長積分法中，默森法相比于夏普法步長控制效果更好，更適合用于C-RAM解命中計算。

若能夠在工程應用方面解決誤差換算比值常數表的編制問題，則該方法便可以工程化應用于高炮防空系統中，可提高高炮的C-RAM作戰能力。