摘要： 為快速、精確地求解高原環(huán)境下的航空炸彈彈道，在經典彈道模型的基礎上，面向高原環(huán)境提出一種基于改進Runge-Kutta的航空無控炸彈彈道解算（improved Runge-Kutta based uncontrolled bombs trajectory calculation， IRK-UBTC）方法。首先，針對炸彈下落過程中時刻變化的環(huán)境參數(shù)，提出了基于起飛點的空氣比重函數(shù)估算方法。然后，引入戰(zhàn)斗機所處環(huán)境的溫度、濕度、緯度、科爾奧利力等信息對經典炸彈質心運動微分方程組進行重構。最后，在三階Runge-Kutta算法的基礎上，提出了帶誤差控制機制的微分方程組實時解算方法。實驗結果表明，該方法在1.4 GHz主頻的嵌入式計算機中的解算周期小于5 ms，且最大射程誤差小于0.2%，能夠被有效應用于高原地區(qū)航彈彈道的實時解算。

關鍵詞： 航空無控炸彈; 彈道解算; 外彈道方程; 數(shù)值計算

中圖分類號： V 271.4+4

文獻標志碼： A

DOI：10.12305/j.issn.1001-506X.2024.06.25

Fast calculation method of aviation bomb trajectory adapted to plateau environment

ZHANG Baichuan， BI Wenhao^*， ZHANG An， LI Minghao

（School of Aeronautics， Northwestern Polytechnical University， Xi’an 710072， China）

Abstract： In order to solve the trajectory of aviation bombs in the plateau environment fast and accurately， based on the traditional trajectory model， an method based on improved Runge-Kutta based uncontrolled bombs trajectory calculation （IRK-UBTC） is proposed. Firstly， an estimation method of the air specific gravity function based on the take-off point is proposed considering the changing environmental parameters during the bomb’s falling process. Then， the temperature， humidity， latitude， Colorili force and other information of the fighter’s environment are introduced to the system of the traditional differential equations of the bomb’s center motion. Finally， based on the third-order Runge-Kutta algorithm， a real-time solution method for differential equations with error control mechanism is proposed. Experiment results show that the calculation period of the proposed method is less than 5 ms in an embedded computer with a main frequency of 1.4 GHz， and the maximum result error is less than 0.2%， which can be applied to real-time calculation of bomb trajectory in plateau areas effectively.

Keywords： aerial uncontrolled bomb; bomb trajectory solution; exterior ballistic equation; numerical calculation

0 引 言

在軍事領域中，航空無控炸彈因其體積小、重量輕、價格低廉等特點，成為了空對地作戰(zhàn)中的重要武器之一^［1-4］。為使航空炸彈適應現(xiàn)代飛機高空高速性能，加大其在戰(zhàn)場中的有效作用距離，需要針對多樣的作戰(zhàn)環(huán)境，從航空火力控制方面對彈道解算和瞄準進行研究。其中，彈道解算方法在機載火控系統(tǒng)中占據了重要的地位。戰(zhàn)斗機執(zhí)行空對地轟炸任務時，通常采用連續(xù)計算投放點（continuously computed release point， CCRP）^［5-6］或連續(xù)計算命中點（continuously computed impact point， CCIP）^［7］的瞄準方式。這類瞄準方式均以無風條件下的標準彈道為基準，綜合考慮風速、俯沖角等因素對投放點或命中點進行實時計算，為飛行員提供轟炸引導。其中，標準彈道的準確性決定了轟炸引導的精度，這就需要機載火控計算機能夠基于航空炸彈的氣動特性和戰(zhàn)斗機當前所處的大氣環(huán)境，不斷地對彈道進行實時解算。目前的機載彈道實時解算通常采用彈道表擬合法，即首先在地面高性能計算機中，根據建立的彈丸模型和彈道解算方法生成海量離散初始條件的彈道表;然后通過特定擬合算法形成連續(xù)空間下的彈道表達式;最后將彈道表裝載到機載火控計算機中，根據實際投放時的氣象條件添加彈道偏置，實時計算炸彈射程等信息。

在過去的幾十年里，眾多的國內外學者在彈丸模型和彈道解算方面發(fā)表了大量的研究成果。在彈丸模型方面，文獻［8］在彈丸平方阻力定律、西亞切阻力定律等基礎上，討論了速度函數(shù)與炸彈速度和空氣阻力加速度的關系，以及阻力系數(shù)和馬赫數(shù)之間的關系。文獻［9］通過建立彈體坐標系和地面坐標系的關系，在彈體坐標系下采用任意拉格朗日-歐拉控制方程對彈道進行求解。文獻［10］建立了基于剛體模型的彈丸質心運動模型，研究了噴流對彈道的影響規(guī)律。文獻［11］在橫風的干擾下，研究了固定舵二維修正彈的外彈道模型。文獻［12］將彈丸視為一個質點，并在攻角為零的條件下對彈丸的外彈道進行仿真，測得在彈丸飛行中的速度曲線以及加速度曲線。

在彈道解算方面，通?；诮⒑玫膹椡柽\動模型，采用Runge-Kutta方法^［13-16］、Adams預測-校正方法^［17-20］等生成彈道表，并使用特定的算法進行擬合，用于機載火控計算機對炸彈落點的實時解算。文獻［21］利用三次樣條曲線方程對炸彈的特征進行擬合，用于優(yōu)化射表。經過實彈檢驗，該方法具有較高的精度和較快的收斂性。文獻［22］通過對彈丸散布的研究，提出彈藥的一維射程修正技術，通過對彈丸加裝阻尼環(huán)的方法，人為地增加彈丸的阻尼系數(shù)，使彈丸在不同環(huán)境下的射程能夠匹配射表。

文獻［23］利用改進擴展卡爾曼濾波（extended Kalman filter， EKF）算法，對炸彈的射程誤差進行修正，提高了對彈丸飛行軌跡預測的準確率。文獻［24］利用圖形處理器單元（graphics processing unit， GPU）的并行計算特性，利用統(tǒng)一計算設備架構（compute unified device architecture， CUDA），基于四階Runge-Kutta算法提出一種大規(guī)模彈道快速計算的方法。經仿真實驗，該方法能夠同時處理一百余條彈道數(shù)據，但其求解精度嚴格依賴于求解步長^［25］。彈丸質心運動微分方程的在線解算方法也在工程中被廣泛研究，理論上越小的積分步長就能帶來越高的求解精度，然而與此同時也大幅增加了解算時間^［26］。

可以看出，目前的研究均集中于在彈丸模型中對環(huán)境參數(shù)進行簡化，并利用如Runge-Kutta等彈道解算方法計算獲得彈道表。然而在面對高原地區(qū)的彈道解算任務時，現(xiàn)有的方法多以平原地區(qū)的彈道表和實驗結果為基礎，通過特定規(guī)律轉換得出^［27］解算結果。這種方法存在以下缺陷：首先，高原環(huán)境中大氣條件更加復雜，目前的環(huán)境參數(shù)均是基于經驗公式迭代，在不同季節(jié)中可能出現(xiàn)計算出落點誤差較大的情況^［28］;并且，由于缺乏對具體氣象條件的分析，以及雷諾數(shù)變化對彈丸阻力的影響，利用傳統(tǒng)彈道表轉換方法難以在高原地區(qū)解算出較為準確的炸彈落點^［29］。此外，高原環(huán)境下的彈道模型需要裝填的參數(shù)較平原彈道更多，所生成的彈道表會出現(xiàn)維度災難的問題。

因此，面對現(xiàn)代化戰(zhàn)場中高原投彈環(huán)境的需求，亟需構建更為精確的彈道模型。本文通過引入不同環(huán)境下的氣壓、溫度、重力、科爾奧利力等因素，提高高原環(huán)境下彈道解算結果的精確度。本文針對高原環(huán)境中的氣象特性以及具體應用場景的解算周期需求，提出一種基于改進Runge-Kutta的航空無控炸彈彈道解算（improved Runge-Kutta based uncontrolled bombs trajectory calculation， IRK-UBTC）方法，以快速求解出高精度的炸彈彈道。本文的主要創(chuàng)新點包括：

（1） 提出基于起飛點的空氣比重函數(shù)估算方法，以實時獲取下落過程中所處環(huán)境的大氣參數(shù);

（2） 在三階Runge-Kutta算法的基礎上，基于步長誤差控制基準公式提出IRK-UBTC方法，能夠在機載火控計算機上實時對彈道模型進行解算，有效地規(guī)避了彈道表擬合法所可能產生的維度災難。

本文結構如下：第1節(jié)介紹了經典炸彈質心運動微分方程;第2節(jié)提出基于改進Runge-Kutta的IRK-UBTC方法;第3節(jié)通過實驗驗證了本文所提模型和方法的有效性及優(yōu)越性;第4節(jié)為總結，并提出下一步的研究方向。

1 炸彈彈道解算數(shù)學模型

由于航空無控炸彈在彈體設計的過程中已解決飛行穩(wěn)定性問題，即章動角很小，本文的研究近似認為作用在彈丸上的空氣阻力的方向與彈軸方向一致且與速度方向相反，并基于此條件對彈丸的質心運動進行分析。

在地理坐標系O_ex_ey_ez_e下對彈丸運動模型進行描述，航空炸彈的質心運動微分方程組如下所示：

式中：C代表航空炸彈的彈道系數(shù);H_τ代表空氣比重函數(shù);G代表炸彈飛行過程中的空氣阻力函數(shù);g代表重力加速度;x，y，z為炸彈質心在坐標系O_ex_ey_ez_e中的空間坐標;v_x，v_y，v_z為炸彈質心速度在坐標系O_ex_ey_ez_e 3個軸上的分量。

在航空炸彈彈道解算模型中選取了地理坐標系作為基準，而在投放瞄準的過程中則需要根據下式轉換為戰(zhàn)斗機坐標系下炸彈的射程：

x_f=x·cos θ+y·sin θ

y_f=－x·sin θ+y·cos θ（2）

式中：x_f代表戰(zhàn)斗機坐標系下的炸彈射程;y_f代表戰(zhàn)斗機坐標系下的炸彈橫向偏差;θ為戰(zhàn)斗機航向角。

2 基于改進Runge-Kutta的彈道解算

為克服傳統(tǒng)彈道表擬合法難以應對動態(tài)環(huán)境參數(shù)的問題，同時規(guī)避過多環(huán)境參數(shù)導致的維度災難，本文提出IRK-UBTC方法。首先分析了動態(tài)環(huán)境參數(shù)下的炸彈運動方程，提出基于起飛點的空氣比重函數(shù)估算方法;然后根據炸彈下落過程中的受力分析，對質心運動微分方程進行了重構;最后提出基于改進Runge-Kutta的微分方程組實時解算方法。所提出的方法通過對迭代過程中單步誤差進行嚴格的控制，兼顧了彈道求解的速度和精度。

2.1 基于起飛點的空氣比重函數(shù)估算方法

在航空炸彈的下落過程中，其所處環(huán)境的參數(shù)（大氣密度、濕度、溫度、重力加速度等）在實時改變，同時也動態(tài)地影響著炸彈下落的軌跡。然而，彈藥上通常未裝備環(huán)境信息傳感器，無法實時反饋當前所處的大氣環(huán)境狀態(tài)，從而難以應對動態(tài)環(huán)境對彈道產生影響。

針對上述問題，本節(jié)通過對大氣靜力方程進行分析，結合外彈道學理論及各項經驗參數(shù)，基于戰(zhàn)斗機起飛點的大氣狀態(tài)提出空氣比重函數(shù)H_τ（z）的估算方法，如下所示：

H_τ（z）=h（z）·T_v，ONh_ON·T_v（z）（3）

式中：h（z）和T_v（z）分別代表炸彈所處環(huán)境的動態(tài)壓強函數(shù)（單位：mmHg）和動態(tài)修正溫度函數(shù)（單位：K）;h_ON和T_v，ON分別代表戰(zhàn)斗機起飛時的壓強和修正溫度。

由于對流層中的大氣壓強隨海拔高度的增加而以指數(shù)形式降低，本文基于此規(guī)律提出動態(tài)壓強函數(shù)的估算方法，如下所示：

h（z）=expln（h₀/h_ON）H₀－H_ON（H₀－H_ON－z）+ln（h_ON）·

（1+δ_h·（H₀－z））（4）

式中：H₀和H_ON分別代表戰(zhàn)斗機投放炸彈和戰(zhàn)斗機起飛時的海拔高度;δ_h代表壓強補償系數(shù);z代表炸彈相對戰(zhàn)斗機的下落距離。其中，壓強補償系數(shù)δ_h的計算方法如下所示：

δ_h=h_ONh_std×（1－2.190 5×10^－5×H_ON）^5.4－1（5）

式中：h_std代表航空兵標準氣象條件的大氣壓強，通常取750 mmHg。

動態(tài)修正溫度函數(shù)T_v（z）同樣隨環(huán)境參數(shù)而變化，主要影響因素為大氣溫度、濕度和壓強，其估算方法如下所示：

T_v（z）=T（z）1－0.375·e（z）/h（z）（6）

式中：T（z）和e（z）分別代表大氣溫度和水蒸氣分壓。在對流層中，大氣溫度與海拔高度呈線性關系，則基于起飛點的大氣溫度估算方法如下所示：

式中：T₀和T_ON分別代表戰(zhàn)斗機投放炸彈和戰(zhàn)斗機起飛時的大氣溫度;δ_T代表溫度補償系數(shù)，計算方法如下所示：

δ_T=T_ON/（T_std－6.328×10^－3×H_ON－1）（8）

式中：T_std代表航空兵標準氣象條件的大氣溫度，通常取288.15 K。

式（6）中的水蒸氣分壓為濕空氣中水蒸氣所形成的壓力，在濕度為RH的空氣中，水蒸氣分壓e的計算方法如下所示：

e（z）=RH（z）×4.579×10^{7.45×（T（z）－273.15）T（z）－38.15}（9）

與大氣溫度相似，在炸彈下落過程中，空氣濕度也可以近似認為與所處海拔高度呈線性關系。因此，基于起飛點的大氣濕度估算方法如下所示：

RH（z）=RH₀－RH_ONH₀－H_ON（H₀－H_ON－z）+RH_ON（10）

在式（4）～式（10）的基礎上，可在炸彈下落過程中對大氣溫度、壓強和濕度進行實時估算，并獲得相比航空兵標準氣象條件更準確的空氣比重函數(shù)H_τ，為航空無控炸彈彈道的精確解算提供理論基礎。

2.2 質心運動微分方程重構

由于地球自轉的影響，炸彈在下落時除了受到空氣阻力加速度J和地球重力加速度g的作用，還受到由科爾奧利力產生的科氏加速度j_e的作用^［30］，如圖1所示。

炸彈下落過程中的阻力函數(shù)G是計算阻力加速度J的重要因素之一^［31］，其計算方法如下所示：

G（z）=4.737×10^－4v_zC_x0v_z6.410 5×g（φ）·T_v×H_τ（z）（11）

式中：v_z代表炸彈下落高度為z時的速度;C_x0為阻力定律函數(shù)，通常通過查詢阻力定律表獲得;φ為投彈點所在地理位置的緯度。

由于地球不是一個標準的球體，在相同的海拔高度下，緯度越大，重力加速度就越大。因此，基于實驗人員測量所得的海平面緯度與重力加速度之間的關系，通過三次擬合得到重力加速度與緯度的映射方程，如下所示：

g（φ）=－1.49×10^－7×φ³+

2.02×10^－5×φ²－3.55×10^－5×φ+9.78（12）

在無控炸彈下落過程中，重力加速度g在實時變化，而投放點到命中點的緯度變化可以忽略不計，因此重力加速度與緯度的映射方程可簡化為如下所示：

式中：φ₀代表戰(zhàn)斗機投彈時所處的緯度; R₀代表地球半徑，通常取6 371 km。

為進一步提高對炸彈彈道的解算精度，本文在彈道模型中綜合考慮地球自轉產生的科氏加速度對炸彈的作用。定義參考系為地球表面，動點為炸彈質心，則在炸彈下落過程中的科氏加速度j_e如下所示：

j_e=2|Ω||v|sin（Ω，v）（14）

式中：sin（Ω，v）為地球自轉角速度方向（沿地軸指向正北）與炸彈速度方向夾角的正弦值;Ω代表地球自轉角速度，方向自西向東，取值為7.292×10^-5;j_e方向由右手定則確定。

綜上，在定義阻力函數(shù)、動態(tài)重力加速度和科氏加速度后，重構的質心運動微分方程如下所示：

2.3 基于改進Runge-Kutta的實時彈道解算方法

在求解第2.2節(jié)中所重構的質心運動微分方程時，需要考慮起飛點環(huán)境參數(shù)數(shù)據、投彈點環(huán)境參數(shù)數(shù)據、目標高度、相對高度等共12項參數(shù)。若在地面高性能計算機中以12個維度解算彈道表并保存，會出現(xiàn)彈道表過大的情況，形成維度災難。因此，針對第2.2節(jié)中所重構的彈丸質心運動微分方程組，需要研究機載火控計算機中的實時解算方法。

在彈道學中對微分方程組的求解通常采用高階Runge-Kutta算法。對于式（15）中速度-時間微分方程dv/dt=f（t，v）， v∈{v_x，v_y，v_z}，經典的三階Runge-Kutta算法可表示為

式中：Δ代表迭代步長。

在迭代過程中，為取得更小的截斷誤差，需要將步長盡可能地縮小。然而隨著步長的縮小，在一定求解范圍內所需要完成的步數(shù)便會隨之增加，導致計算量大幅增加。對于航空無控炸彈的彈道解算場景，通常需要取算法的迭代步長小于5 ms。此時，常規(guī)的三階Runge-Kutta算法已無法滿足機載火控計算機對實時性的需求。

因此，本節(jié)結合彈丸下落時的運動特性，提出一種改進的Runge-Kutta微分方程實時解算方法IRK-UBTC。基于給定的總體誤差條件ε，動態(tài)調整微分方程的解算步長，從而以最低的時間復雜度為代價，實時求解炸彈彈道。

炸彈的自由下落過程可總結為3個階段，如圖2所示。

（1） 減速階段：炸彈投放初期，由于炸彈水平方向初速度較大，從而空氣阻力加速度較大，因此合速度在投放初期呈現(xiàn)降低的趨勢。

（2） 加速階段：經過一段時間后，在重力加速度作用下，重力主導的Z_e方向速度分量顯著增大，合速度呈增加的趨勢。

（3） 平緩及下降階段：炸彈在加速下落過程中，空氣阻力逐漸增大，與重力相抵，Z_e軸方向加速度逐漸變小，炸彈的飛行速度趨于平緩后逐步下降。

基于上述分析，在對炸彈質心運動微分方程組的求解過程的初期，速度變化程度較為劇烈，需要盡可能減小解算步長;而在炸彈下落末端，適當增加求解步長不會造成過多的精度損失。

因此，定義微分方程組動態(tài)迭代步長Δ，如下所示：

Δ=［a+b（H₀－H_ON）］·exp{z/c}（17）

式中：a和b為步長基準控制參數(shù);c為步長增量控制參數(shù)。

同時，為保證在所有投彈條件下均滿足給定的總體誤差范圍，本文進一步給出了基于截斷誤差的步長控制策略。首先，從時間點t_i出發(fā)，將以步長Δ求取的速度近似值記為v^Δ_i+1;然后，將以步長Δ/2求取的速度近似值記為v^Δ/2_i+1。由于三階Runge-Kutta公式對于步長的截斷誤差存在O（Δ⁴），故有：

式中：v［t_i+1］為t+（i+1）Δ時刻的速度真實值;k為某特定常數(shù)。由式（18）和式（19）以及總體誤差條件ε可得：

v［t_i+1］－v^Δ/2_i+1v［t_i+1］≈17v^Δ/2_i+1－v^Δ_i+1v［t_i+1］≤ε（20）

|v^Δ/2_i+1－v^Δ_i+1|max{|v^Δ_i+1|，|2v^Δ/2_i+1|}≤7ε（21）

則式（21）為步長誤差控制基準。進一步地，針對式（15）中的后3項對炸彈質心的迭代公式，由于方程中變量與步長無關，后3項的求解可利用Eula法進行簡化，即

p_i+1=p_i+Δ·f（t_i，p_i）（22）

式中：p∈{x，y，z}為炸彈在3個分方向的位移。

綜上，定義步長誤差控制周期η，本文所提出的IRK-UBTC方法流程圖如圖3所示，共分為11個步驟。

步驟 1 在戰(zhàn)斗機起飛前裝載起飛點的各項環(huán)境參數(shù)，包括海拔高度H_ON、大氣壓強h_ON、大氣溫度T_ON、修正大氣溫度τ_ON，以及大氣濕度RH_ON。

步驟 2 輸入炸彈投放點的各項環(huán)境參數(shù)及目標相對信息，包括海拔高度H₀、大氣壓強h₀、大氣溫度T₀、修正大氣溫度τ₀、大氣濕度RH₀、所處緯度φ、目標海拔高度H_T以及彈道系數(shù)C，并定義迭代次數(shù)i=0。

步驟 3 根據式（17）計算迭代步長Δ。

步驟 4 根據式（4）～式（10）估算炸彈當前所處環(huán)境的各項環(huán)境參數(shù)。

步驟 5 針對式（15）炸彈質心運動微分方程組的前3項，利用改進Runge-Kutta微分方程實時解算方法，以步長Δ求解出v^Δ_x，i+1、v^Δ_y，i+1與v^Δ_z，i+1。

步驟 6 判斷迭代次數(shù)i是否為η的整數(shù)倍，若是則繼續(xù)執(zhí)行步驟7，否則跳轉至步驟10。

步驟 7 利用改進Runge-Kutta微分方程實時解算方法，以步長Δ/2求解出v^Δ/2_x，i+1、v^Δ/2_y，i+1與v^Δ/2_z，i+1。

步驟 8 對于所有（v^Δ/2_i+1，v^Δ_i+1）∈{（v^Δ/2_x，i+1，v^Δ_x，i+1）， （v^Δ/2_y，i+1，v^Δ_y，i+1），（v^Δ/2_z，i+1，v^Δ_z，i+1）}，若滿足步長誤差控制基準式（21），則跳轉至步驟10，否則繼續(xù)執(zhí)行步驟9。

步驟 9 將步長折半后迭代，即Δ′=Δ/2。

步驟 10 根據式（22）所述Eula法求取下一時刻的炸彈位移（x_i+1，y_i+1，z_i+1），并將迭代次數(shù)i增加1。

步驟 11 判斷炸彈下落距離z是否達到目標點位置，若是則結束，否則跳轉至步驟3繼續(xù)迭代。

3 數(shù)值分析

本節(jié)在三維環(huán)境中對上文所提出的航空炸彈彈道解算方法進行數(shù)值分析。炸彈及仿真環(huán)境默認參數(shù)如表1所示。本節(jié)中所有仿真代碼均基于C99標準實現(xiàn)，分別在Windows環(huán)境下（Intel Core^TM2 Duo CPU E7 500 @ 2.93 GHz，2 GB RAM）和Linux環(huán)境下（RaspberryPi-3B+，BCM288 37B0 SoC ARM 64-bit CortexTM-A53 CPU 1.4 GHz，1 GB RAM）進行測試。

表1 航空炸彈彈道仿真參數(shù)

Table 1 Air bomb trajectory simulation parameter參數(shù)取值炸彈彈道系數(shù)C0.462 859戰(zhàn)斗機起飛時海拔高度H_ON/m0戰(zhàn)斗機起飛時大氣壓強h_ON/（mmHg）750戰(zhàn)斗機起飛時大氣溫度T_ON/K288.15戰(zhàn)斗機起飛時大氣濕度RH_ON50%戰(zhàn)斗機投彈時海拔高度H₀/m2 000～10 000戰(zhàn)斗機投彈時大氣壓強h₀/（mmHg）288.15-6.328×10^-3×H₀戰(zhàn)斗機投彈時大氣溫度T₀/K760×

（1-2.190 5×10^-5×H₀）^5.4戰(zhàn)斗機投彈時大氣濕度RH₀50%戰(zhàn)斗機投彈時所處緯度φ30目標海拔高度H_t/m1 000戰(zhàn)斗機投彈時炸彈速度v₀/（km/h）800戰(zhàn)斗機投彈時航向角θ0步長基準控制參數(shù)及步長增量

控制參數(shù)（a， b， c）（0.02， 270， 15 000）

3.1 算法誤差與分析

本節(jié)在Linux環(huán)境下裝載相關參數(shù)，并在戰(zhàn)斗機投彈點海拔2 000～10 000 m處隨機取1 000點進行彈道仿真計算。記錄仿真計算結果和運行時間，并與標準彈道進行對比，結果如表2、圖4和圖5所示。

由實驗結果可知，本文所提出的IRK-UBTC方法射程誤差最大不超過1.162 0%，時間誤差最大不超過51 ms，且在RaspberryPi平臺上的解算周期小于5 ms。這是由于在解算過程中，步長誤差控制基準式（21）對期望誤差進行了嚴格的控制。相較于傳統(tǒng)的定步長和變步長三階Runge-Kutta算法，本文所提出的IRK-UBTC方法的射程誤差和下落時間誤差均相對較低。盡管變步長三階Runge-Kutta算法的運行速度較快，但精度仍低于IRK-UBTC方法。這是由于在誤差和解算時長允許范圍內，IRK-UBTC方法盡可能地增加解算步長，以降低解算周期。本文所提出的方法具有精度高、周期短的特點，能夠被部署在機載火控計算機中，對彈道進行實時解算。

3.2 IRK-UBTC方法精度分析

為對比IRK-UBTC方法與傳統(tǒng)Runge-Kutta算法^［31］的精度，本節(jié)對炸彈質心微分方程組進行解算。同時，取四階Runge-Kutta算法為基準，計算炸彈射程殘差曲線，如圖6所示。

由圖6可知，IRK-UBTC方法與標準四階Runge-Kutta算法相比，最大飛行距離殘差相差不超過0.1 m。計算其判定系數(shù)R²=（1-6×10^-9）≈1，證明利用IRK-UBTC方法解算獲得的彈道數(shù)據近似等于四階Runge-Kutta算法。這是由于在構建IRK-UBTC的變步長策略時，對炸彈運動初期變化劇烈的階段采用了較大的步長，而對炸彈的平緩飛行階段采用了較小的步長。因此，炸彈下落前段的解算精度不低于基準算法，同時省去了下落后段解算過程中的無效迭代。

3.3 高原彈道實驗分析

本節(jié)用某型炸彈在高原地區(qū)進行投放測試，使用本文算法及傳統(tǒng)算法對彈道進行解算，并對比射程誤差。對比結果如表3所示。

由表3可知，傳統(tǒng)算法與本文所提算法在高原地形環(huán)境下都可解算出較為準確的彈道。并且，在高原地形環(huán)境下，本文所提出的IRK-UBTC的解算精度相較于傳統(tǒng)算法精度更高。這是由高原地區(qū)空氣較為稀薄、空氣阻力較小、重力加速度較小等因素綜合導致的。本文提出的算法在對炸彈彈道進行迭代求解時，綜合考慮了下落過程中的高度、大氣壓強、溫度、濕度以及地球自轉所帶來的影響，對空氣阻力、重力以及科爾奧利力進行實時的更新，從而通過計算獲得相對更準確的炸彈彈道。同時，通過更改彈道系數(shù)及戰(zhàn)斗機上裝訂的投彈環(huán)境參數(shù)，即可適配不同規(guī)格的航空無控炸彈。與傳統(tǒng)算法相比，本文提出的IRK-UBTC算法更符合炸彈的實際下落場景，在工程上具有較高的可行性。

4 結 論

本文在經典平原航空無控炸彈彈道模型的基礎上，通過引入溫度、濕度、戰(zhàn)斗機緯度、科爾奧利力等多項環(huán)境參數(shù)，面向高原環(huán)境建立了彈道質心運動微分方程組。然后，為解決模型參數(shù)過多導致的生成彈道表時的維度爆炸問題，本文在三階Runge-Kutta算法的基礎上提出了帶誤差控制機制的微分方程組實時解算方法。最后，通過在高原地區(qū)的炸彈投放實驗，驗證了所提出的IRK-UBTC算法在高原環(huán)境下具有彈道數(shù)據解算快、精度高的特點，在工程上具有可觀的應用前景。

參考文獻

［1］SENAY N. The strategic level optimization of air to ground missiles for turkish air force decision support system［R］. Ohio： Wright-Patterson Air Force Base， 2012.

［2］ABOUSEADA H A H A. Modeling， simulating and control of free falling bomb using PID［C］∥Proc.of the International Undergraduate Research Conference， 2022.

［3］KOZLOVSKA M， SYBKOVA H， OTAHAL P P. Radiation monitoring after experimental dirty bomb explosion［J］. Radiation Protection Dosimetry， 2023， 199（8/9）： 1012-1020.

［4］CARPENTER C， MONTGOMERY A H. The stopping power of norms： saturation bombing， civilian immunity， and US attitudes toward the laws of war［J］. International Security， 2020， 45（2）： 140-169.

［5］ZHANG A， XU H， BI W H， et al. Adaptive mutant particle swarm optimization based precise cargo airdrop of unmanned ae-rial vehicles［J］. Applied Soft Computing， 2022， 130： 109657.

［6］LEONARD A， ROGERS J， GERLACH A. Probabilistic release point optimization for airdrop with variable transition altitude［J］. Journal of Guidance， Control， and Dynamics， 2020， 43（8）： 1487-1497.

［7］ATANASOV M. Mathematical model for operation of aviation systems for delivery of special means to air and earth objects［J］. Aerospace Research Bulgaria， 2022， 1（34）： 138-148.

［8］王勇亮， 趙成仁， 盧穎. 炸彈空氣阻力加速度的仿真與實現(xiàn)［J］. 彈箭與制導學報， 2006（S1）： 251-252， 256.

WANG Y L， ZHAO C R， LU Y. Simulation and realization of bomb air drag force acceleration［J］. Journal of Projectiles， Rockets， Missiles and Guidance， 2006， 26（S1）： 251-252， 256.

［9］SAHU J. Time-accurate numerical prediction of free-flight aerodynamics of a finned projectile［J］. Journal of Spacecraft and Rockets， 2008， 45（5）： 946-954.

［10］SAHU J. Unsteady aerodynamic simulations of a finned projectile at a supersonic speed with jet interaction［C］∥Proc.of the 52nd AIA Aerospace Sciences Meeting， 2014： 3024

［11］郝永平， 陳闖， 張嘉易， 等. 固定舵二維修正彈外彈道仿真與動態(tài)模擬［J］. 兵工學報， 2018， 39（4）： 67-76.

HAO Y P， CHEN C， ZHANG J Y， et al. Trajectory and dynamic simulations of two-dimensional trajectory correction projectile with fixed canards［J］. Acta Armamentarii， 2018， 39（4）： 67-76.

［12］王莎莎， 張東東， 李新娥， 等. 基于FLUENT和ADAMS的外彈道加速度仿真［J］. 探測與控制學報， 2019， 41（5）： 47-52.

WANG S S， ZHANG D D， LI X E， et al. Exterior trajectory acceleration simulation based on FLUENT and ADAMS［J］. Journal of Detection amp; Control， 2019， 41（5）： 47-52.

［13］AHMADIANFAR I， HEIDARI A A， GANDOMI A H， et al. RUN beyond the metaphor： an efficient optimization algorithm based on Runge Kutta method［J］. Expert Systems with Applications， 2021， 181： 115079.

［14］LI J W， LI X， JU L L， et al. Stabilized integrating factor Runge-Kutta method and unconditional preservation of maximum bound principle［J］. SIAM Journal on Scientific Computing， 2021， 43（3）： 1780-1802.

［15］CHEN H， AHMADIANFAR I， LIANG G， et al. A successful candidate strategy with Runge-Kutta optimization for multi-hydropower reservoir optimization［J］. Expert Systems with Applications， 2022， 209： 118383.

［16］QIAO Z L， LI L， ZHAO X C， et al. An enhanced Runge-Kutta boosted machine learning framework for medical diagnosis［J］. Computers in Biology and Medicine， 2023， 160： 106949.

［17］DIETHELM K， FORD N J， FREED A D. A predictor-corrector approach for the numerical solution of fractional differential equations［J］. Nonlinear Dynamics， 2002， 29（1/4）： 3-22.

［18］WU G C， WEI J L， LUO C， et al. Parameter estimation of fractional uncertain differential equations via Adams method［J］. Nonlinear Analysis： Modelling and Control， 2022， 27（3）： 413-427.

［19］ZABIDI N A， MAJID Z A， KILICMAN A， et al. Numerical solution of fractional differential equations with Caputo derivative by using numerical fractional predict-correct technique［J］. Advances in Continuous and Discrete Models， 2022， 2022（1）： 1-23.

［20］ODIBAT Z， ERTURK V S， KUMAR P， et al. Dynamics of generalized Caputo type delay fractional differential equations using a modified predictor-corrector scheme［J］. Physica Scripta， 2021， 96（12）： 125213.

［21］QI Z， LYSTER D. Multi-spline technique for the extraction of drag coeffidents from radar data［J］. Journal of Beijing Institute of Technology， 1994， 3（1）： 33-42.

［22］李杰， 馬寶華. 迫擊炮彈一維射程修正引信技術研究［J］. 兵工學報， 2001， 22（4）： 553-555.

LI J， MA B H. Research on one-dimensional range correction fuze technology of mortar shells［J］. Acta Armamentarii， 2001， 22（4）： 553-555.

［23］雷曉云， 張志安， 杜忠華. 基于改進無跡卡爾曼濾波的彈道射程修正算法研究［J］. 兵工學報， 2018， 39（9）： 1701-1710.

LEI X Y， ZHANG Z A， DU Z H. Ballistic range correction algorithm based on an improved unscented Kalman filter［J］. Acta Armamentarii， 2018， 39（9）： 1701-1710.

［24］左軍濤， 朱恩成， 黃四牛， 等. 基于GPU的彈道快速計算方法［J］. 火力與指揮控制， 2012， 37（9）： 193-197.

ZUO J T， ZHU E C， HUANG S N， et al. Fast calculation method base on GPU for trajectory［J］. Fire Control amp; Command Control， 2012， 37（9）： 193-197.

［25］項帆， 王雨時， 聞泉， 等. 無控彈丸剛體外彈道學應用綜述［J］. 探測與控制學報， 2021， 43（4）： 14-26.

XIANG F， WANG Y S， WEN Q， et al. Summary of uncontrolled projectile rigid exterior ballistics application［J］. Journal of Detection amp; Control， 2021， 43（4）： 14-26.

［26］楊青， 李若， 蔡振寧. 六自由度外彈道方程組的快速數(shù)值方法［J］. 高等學校計算數(shù)學學報， 2014， 36（3）： 253-270.

YANG Q， LI R， CAI Z N. A fast numercial method for the six degrees of freedom model in the exrernal ballistics［J］. Journal of Computational Mathematics of Colleges and Universities， 2014， 36（3）： 253-270.

［27］丁天寶， 何朝， 王良明， 等. 高速旋轉炮彈寬海拔彈道解算方法［J］. 兵工學報， 2021， 42（1）： 209-213.

DING T B， HE Z， WANG L M， et al. Calculation method of firing trajectory of high spinning projectile adapted to wide altitude［J］. Acta Armamentarii， 2021， 42（1）： 209-213.

［28］趙靜， 杜忠華， 趙永平， 等. 基于三次B樣條曲線的一維彈道修正彈空氣阻力系數(shù)擬合［J］. 火力與指揮控制， 2015， 40（4）： 123-126.

ZHAO J， DU Z H， ZHAO Y P， et al. Simulation for air resistance coefficient of one-dimension trajectory projectile based on cubic B-spline curve［J］. Fire Control amp; Command Control， 2015， 40（4）： 123-126.

［29］FENRICK W J. Targeting and proportionality during the NATO bombing campaign against Yugoslavia［J］. European Journal of International Law， 2001， 12（3）： 489-502.

［30］MAHMOUDIAN M， FILHO J， MELICIO R， et al. Three-dimensional performance evaluation of hemispherical coriolis vibratory gyroscopes［J］. Micromachines， 2023， 14（2）： 254.

［31］浦發(fā). 航空炸彈的標準落下時間和阻力定律及彈道系數(shù)關系問題的探討［J］. 兵工學報， 1985（3）： 1-7.

PU F. Discussion on the relationship between the standard falling time and the drag law and the ballistic coefficient of air bombs［J］. Acta Armamentarii， 1958（3）： 1-7.

作者簡介

張百川（1996—），男，博士研究生，主要研究方向為彈道解算、態(tài)勢預測、任務分配。

畢文豪（1986—），男，副研究員，博士，主要研究方向為復雜系統(tǒng)建模、仿真與效能評估。

張 安（1962—），男，教授，博士，主要研究方向為火力控制、復雜系統(tǒng)建模。

李銘浩（1997—），男，博士研究生，主要研究方向為火力控制、多無人機協(xié)同控制。