徐 超，王治霖，王武剛，周永佳，白 嬋

(1 北京理工大學(xué)宇航學(xué)院，北京 100081;2 西北工業(yè)集團有限公司，西安 710043)

0 引言

為保證末制導(dǎo)炮彈命中率，需要提供較高精度的射擊諸元，即對射表精度有較高要求。諸元解算指炮彈發(fā)射前，根據(jù)實際發(fā)射條件和武器系統(tǒng)狀態(tài)計算出一組準確的諸元值，裝訂在武器系統(tǒng)上使其滿足要求，精確完成打擊任務(wù)。射擊諸元主要包括：裝藥號、射角(以下稱表尺)、發(fā)射方位角(以下稱方位)、機械陀螺工作時間(以下稱程裝)、激光照射起始時間(以下稱延遲)等。諸元解算時，需考慮導(dǎo)引頭視場范圍、過載能力、末制導(dǎo)段速度、射程等多約束條件；需考慮氣象條件(包括溫度、地面氣壓、橫風、縱風等)、馬格努斯效應(yīng)、藥溫等因素；需考慮延時、程裝、表尺、方位等多輸出問題。因此，可將諸元解算看成是一種具有多輸入多輸出的參數(shù)優(yōu)化問題，需要空間和時間雙維度耦合迭代的復(fù)雜尋優(yōu)計算。

目前，諸元解算存在以下問題：

1)彈道解算過程主要通過彈道積分的方法獲得導(dǎo)彈落點，涉及六自由度非線性微分方程的數(shù)值求解，導(dǎo)致彈道計算耗時較長，對作戰(zhàn)任務(wù)中的時敏目標精確打擊極為不利。

2)需根據(jù)當?shù)貧庀髼l件，查詢預(yù)先裝訂的射表，修正上述各影響因素。但這種方法解算精度較差，降低作戰(zhàn)效能。

針對上述問題，Ollerenshaw等基于彈丸線性理論，簡化了六自由度方程，推導(dǎo)了控制力作用下的彈丸轉(zhuǎn)向幅值和角度計算公式，結(jié)果與六自由度方程結(jié)果接近，但沒有考慮風的影響且是平射彈道。Hainz等提出了基于彈丸修正線性理論的快速彈道預(yù)測方法，相較于非線性六自由度、三自由度和質(zhì)點彈道數(shù)值積分進行彈道預(yù)測方法，具有計算機占用資源少及計算精度高等優(yōu)點,但為保證全彈道預(yù)測精度，彈丸線性理論在計算彈道時需實時更新迭代多次，使得解算復(fù)雜、數(shù)據(jù)量大，降低了彈道解算的實時性。丁天寶等創(chuàng)建了適應(yīng)寬海拔的彈道解算新方法。考慮非線性姿態(tài)運動，建立了阻力系數(shù)隨海拔高度變化的模型，研究結(jié)果表明新型彈道解算方法具有較高精度，適用于寬海拔作戰(zhàn)。趙東華等研究應(yīng)用二分法求根的思想求解外彈道方程組，從而決定了火炮射擊諸元的算法,但若想得到較高精度的諸元，需要積分的次數(shù)會隨之增加，計算量依舊較大。陳瑞軍等針對研究初期確立的基本諸元計算方法不能滿足有效攻擊區(qū)戰(zhàn)技指標及計算速度慢等問題，提出了將有效攻擊區(qū)中心為表載射程及簡化的彈道模型用于諸元計算的方法,但有效攻擊區(qū)中心需要計算諸多彈道，再計算捕獲區(qū)域的中心，若區(qū)域為不規(guī)則，精度可能會較差。

針對激光末制導(dǎo)炮彈諸元精確解算的需求，文中提出了一種基于PSO(particle swarm optimization)算法優(yōu)化的基本諸元解算方法。首先，分析作用力和氣動力，建立末制導(dǎo)炮彈的六自由度動力學(xué)和運動學(xué)方程組，構(gòu)建仿真模型。介紹了大氣模型，并分析高度、虛溫、氣壓、風速、風向數(shù)據(jù)對空氣密度、聲速、馬赫數(shù)等影響；其次,應(yīng)用龍格-庫塔法對六自由度動力學(xué)和運動學(xué)方程組進行解算。給出氣動參數(shù)辨識方法，進行彈道校驗；再次,分析了表尺、程裝等參數(shù)變化時，對射程和射高的影響，將其作為變量，應(yīng)用粒子群算法優(yōu)化得到最優(yōu)解，作為射表相應(yīng)的諸元解算結(jié)果。最后，將模型經(jīng)PSO優(yōu)化后進行仿真和實際打靶驗證，結(jié)果表明，無論是仿真驗證還是實際打靶試驗，這種方法精度更高，耗時更短，可滿足實戰(zhàn)化需求。

1 彈丸的彈道模型

末制導(dǎo)炮彈由舵機、導(dǎo)引頭、發(fā)動機、慣性陀螺、尾翼等分部件組成，包含氣動、動力、控制、制導(dǎo)等系統(tǒng)。由于各系統(tǒng)部件參數(shù)、氣象測量、發(fā)射平臺均存在誤差，若要在各種誤差下仍能精確命中目標，需建立更加精確的彈丸模型。

1.1 氣動力及作用力分析

末制導(dǎo)炮彈在飛行過程中，不同階段其氣動布局均有變化，主要氣動布局如表1所示。

表1 氣動布局

根據(jù)不同飛行模式代入對應(yīng)參數(shù)，使模型更加準確。為得到準確的飛行數(shù)學(xué)模型，根據(jù)表1中不同的氣動布局進行理論值計算及風洞試驗，得出各狀態(tài)下的阻力系數(shù)、升力系數(shù)、傳遞比、滾轉(zhuǎn)力矩系數(shù)等氣動參數(shù)。對于末制導(dǎo)炮彈，一般情況下有兩種運動狀態(tài)：一種是用于近距離射擊，稱為近區(qū)彈道；另一種是用于遠距離射擊，稱為遠區(qū)彈道。以遠區(qū)為例，末制導(dǎo)炮彈的彈道組成如圖1所示，分為以下5段：

圖1 彈目之間角度關(guān)系

1)出炮口無控段：出炮口后末制導(dǎo)炮彈以低速旋轉(zhuǎn)的無控方式飛行，直至發(fā)動機工作；2)助推發(fā)動機工作段：發(fā)動機開始工作使炮彈的速度增加，從而延長飛行距離；3)助推發(fā)動機結(jié)束工作后無控段：助推發(fā)動機工作結(jié)束后，炮彈仍以低速旋轉(zhuǎn)的無控方式飛行；4)慣性制導(dǎo)段：陀螺解鎖后，慣性系統(tǒng)開始工作，為增加射程及末制導(dǎo)距離，此階段通過重力補償?shù)姆绞绞古趶椈景粗本€飛行；5)末制導(dǎo)段：激光開始照射發(fā)出信號，當炮彈導(dǎo)引頭接收到目標反射的信號時，炮彈開始轉(zhuǎn)入末制導(dǎo)段，直至命中目標。

對于近區(qū)彈道，由于助推發(fā)動機不工作，同時也不需要滑翔增大射程，因此近區(qū)彈道僅由自由飛行段和末制導(dǎo)段組成。

末制導(dǎo)炮彈在飛行過程中所受的外力主要是空氣動力，可分解為阻力、升力和側(cè)向力，其表達式為：

(1)

式中：，，分別代表阻力、升力、側(cè)向力;、、為升力系數(shù)、阻力系數(shù)、側(cè)向力系數(shù);為動壓頭；為參考面積。由計算和風洞試驗得出，發(fā)動機推力可表示為：

=+(-)

(2)

式中：為每秒燃料的消耗量;為燃氣在噴口噴出的速度;為發(fā)動機噴口截面積;為噴口截面處燃氣流壓強;為噴口周圍的大氣靜壓強。此外對于舵機控制力，以有兩對舵為例，即上下舵(2舵、4舵)和左右舵(1舵、3舵)。對于左右舵：

=

(3)

式中：為舵片效率;為左右舵的舵偏轉(zhuǎn)角，沿上偏為正。對于上下舵：

=

(4)

式中：為上下舵的舵偏轉(zhuǎn)角，沿右偏為正。

1.2 彈道方程建立

彈丸的剛體六自由度彈道方程能夠準確表示彈丸的飛行運動狀態(tài)，建立運動學(xué)模型和動力學(xué)模型。質(zhì)心運動動力學(xué)方程為：

(5)

式中：為炮彈實時質(zhì)量;為炮彈飛行速度；為發(fā)動機推力，由氣動數(shù)據(jù)進行插值得到;、分別為攻角、側(cè)滑角；為彈道傾角;為速度傾斜角;為彈道偏角。繞質(zhì)心轉(zhuǎn)動動力學(xué)方程為：

(6)

式中：、、分別為飛行器在彈體坐標各軸的轉(zhuǎn)動慣量;、、分別為飛行器受外力對質(zhì)心產(chǎn)生的力矩在單體各軸上的投影;、、分別為飛行器繞質(zhì)心旋轉(zhuǎn)的俯仰、偏航、滾轉(zhuǎn)角速度。質(zhì)心運動方程為：

(7)

式中：、、分別為飛行器質(zhì)心的位置坐標。繞質(zhì)心轉(zhuǎn)動運動學(xué)方程為：

(8)

式中：?、分別為飛行器的俯仰角、偏航角。發(fā)動機工作時，炮彈質(zhì)量變化方程為：

(9)

式中：為質(zhì)量秒流量。角度關(guān)系為：

(10)

對于慣導(dǎo)段，其彈道為重力補償彈道，即重力沿垂直于速度方向的法向分量與攻角產(chǎn)生的補償升力之間存在平衡關(guān)系。這種平衡關(guān)系可保證飛行時其縱軸對地平面的傾斜角基本不變，增大射程。

在末制導(dǎo)段上，延時開始時刻，導(dǎo)引頭開始工作，測量出彈目線旋轉(zhuǎn)角速度，再根據(jù)過重力補償?shù)谋壤龑?dǎo)引制導(dǎo)律給出對舵機的控制信號，修正縱向和橫向偏差，直至命中目標。末制導(dǎo)原理框圖如圖2所示。

圖2 末制導(dǎo)原理框圖

1.3 氣象模型

對于氣象數(shù)據(jù)，由氣象雷達車測量地面和高空的氣象條件進行采集，包括氣溫、氣壓、風速、風向等。以2021年9月7日在某靶場一次飛行實驗為例，以200 m分層的氣象數(shù)據(jù)如表2所示。

表2 氣象數(shù)據(jù)

對表2進行插值，即可得到不同高度下的虛溫、氣壓、風速、風向數(shù)據(jù)。對于虛溫，需將其數(shù)據(jù)轉(zhuǎn)換為開氏溫度。對于不同表尺和程裝數(shù)據(jù)，彈道高度都不相同。因此在應(yīng)用此表格時，高度應(yīng)覆蓋彈道高。下面分析這些數(shù)據(jù)變化對空氣密度、聲速等的影響。首先分析對空氣密度的影響：

(11)

式中：為氣體常數(shù)，其值為287.14；為當前高度下的氣壓；為當前高度下的溫度；為當前高度下的空氣密度。其次，分析對聲速的影響：

(12)

式中：為當前高度下的聲速。最后，對于當前高度下的風速和風向，主要影響彈丸的速度。聲速和來流速度對當前高度的馬赫數(shù)影響，可表示為：

(13)

2 模型方程解算及彈道模型校驗

2.1 模型方程解算

對于上述六自由度微分方程，常用的微分方程的數(shù)值解法有歐拉法、龍格-庫塔法等。歐拉法的特點是簡單易行，但精度低。在同樣計算步長的條件下，龍格-庫塔法的計算精度要比歐拉法高，因此采用龍格-庫塔法計算六自由度微分方程。設(shè)有一階微分方程：

(14)

若已知時刻的參數(shù)值，則可用龍格-庫塔法求+1=+Δ時刻的+1的近似值，龍格-庫塔公式為：

(15)

四階龍格-庫塔法每積分一個步長，需要計算4次右端函數(shù)值，并將其線性組合求出被積函數(shù)的增量Δ。四階龍格-庫塔法除了精度較高外，還易于編制計算程序，改變步長方便，是一種自啟動的單步數(shù)值積分方法。因此，對于六自由度微分方程采用此方法進行解算。

2.2 彈道模型校驗

飛行試驗可獲取豐富的測量信息，包括地面雷達獲取的速度、位置信息以及彈上黑匣子慣性測量數(shù)據(jù)，利用這些數(shù)據(jù)可識彈道模型重要的氣動參數(shù)，實現(xiàn)彈道模型的高置信度校驗。

氣動參數(shù)是彈道模型構(gòu)建的基礎(chǔ)參數(shù)，常用的參數(shù)辨識方法有最小二乘法、最大似然、最小方差、最小風險等。文中應(yīng)用牛頓-拉夫遜算法求解氣動參數(shù)辨識問題，氣動參數(shù)辨識的具體過程如下：

步驟1：輸入初始數(shù)據(jù)。參數(shù)的初估值，狀態(tài)初值；待估計參數(shù)的上限和下限值；動力學(xué)系統(tǒng)試驗的觀測量和控制量輸入的實測值，即經(jīng)預(yù)處理后的雷達數(shù)據(jù)。

步驟3：準則函數(shù)的計算。利用模型輸出值和實測數(shù)據(jù)，求得準則函數(shù)。

步驟4：噪聲協(xié)方差矩陣計算。利用模型輸出值和實測數(shù)據(jù)求得該矩陣。

步驟5：待識別參數(shù)增量Δ的計算。通過牛頓-拉夫遜算法的迭代公式求得本次迭代增量。并通過靈敏度公式計算出靈敏度。

在上述過程中，待辨識參數(shù)的迭代初值通常取過去的試驗值或理論計算值。的上限和下限按待辨識參數(shù)的物理意義，結(jié)合實際情況給定，收斂指標則根據(jù)精度要求進行取值。

3 基于粒子群的諸元解算

3.1 參數(shù)分析

分析表尺、程裝等參數(shù)變化對射程將產(chǎn)生影響。采用控制變量法，分別對遠區(qū)及近區(qū)兩種運動狀態(tài)進行分析。

對于近區(qū)彈道，由于助推發(fā)動機不工作，同時也不需要滑翔增大射程，因此只考慮表尺變化時對射程的影響。選取初速為326 m/s的炮彈，無程裝，對表尺為300 mil、350 mil、400 mil進行仿真,仿真結(jié)果如圖3所示。

圖3 表尺變化對初速為326 m/s無程裝炮彈射程的影響

由圖3可知，對于近區(qū)彈道，當表尺增大時，射程及射高都同時增大。對于遠區(qū)彈道，分析表尺變化時對射程的影響。選取初速為760 m/s的炮彈，固定程裝為110分化，對表尺為480 mil、500 mil、520 mil進行仿真，仿真結(jié)果如圖4所示。

圖4 表尺變化對初速760 m/s固定程裝110分化炮彈射程的影響

由圖4可知，與近區(qū)彈道相同，當表尺增大時，遠區(qū)射程及射高都同時增大。分析程裝變化對射程的影響，選取初速為760 m/s的炮彈，固定表尺為500 mil，對表尺為105分化、110分化、115分化進行仿真，仿真結(jié)果如圖5所示。

圖5 程裝變化對初速760 m/s固定表尺500 mil炮彈射程的影響

由圖5可知，當程裝時間增大時，射程增大，但對射高并無明顯影響。由上述分析可知，無論近區(qū)還是遠區(qū)彈道，不同參數(shù)對射程均有不同影響。同時，對射高有一定影響，再次證明了氣象數(shù)據(jù)中高度需覆蓋射高，否則會對彈道產(chǎn)生較大影響。

3.2 基于粒子群的諸元解算方法

諸元解算要考慮的因素較多，包括氣象條件、藥溫、馬格努斯效應(yīng)等；同時，所考慮的約束條件較多，包括期望射程、程裝與慣性陀螺工作的角度范圍、導(dǎo)引頭有限視場角、彈體機動過載能力、末制導(dǎo)段速度等均應(yīng)滿足彈丸的使用要求；諸元解算的輸出量較多，包括表尺、方位、程裝、延時等。

對于多輸入多輸出的參數(shù)優(yōu)化問題，粒子群算法具有參數(shù)少、執(zhí)行效率較高等優(yōu)點，可以在有限計算量下得到全局最優(yōu)解。粒子群算法所有的粒子都有一個由被優(yōu)化的函數(shù)決定的適應(yīng)度值，每個粒子速度決定他們的飛行的方向和距離，粒子群們就追隨當前的最優(yōu)粒子在解空間中搜索。

粒子群初始化為一群隨機粒子，然后通過迭代找到最優(yōu)解。在每一次迭代中，粒子通過跟蹤兩個“極值”來更新。第一個就是粒子本身所找到的最優(yōu)解，這個解叫做個體極值；另一個極值是整個種群目前找到的最優(yōu)解，這個極值是全局極值。另外也可以不用整個種群而只是用其中一部分作為粒子的鄰居，那么在所有鄰居中的極值就是局部極值。

粒子位置可以用一個維的向量來表示，第個粒子位置表示為：=(,1,…,,,…,,)，第個粒子的局部極值位置為=(,1,…,,,…,,)，全局最優(yōu)解的位置=(,1,…,,,…,,)第個粒子的速度為=(,1,…,,,…,,)。在找到粒子的位置和速度，以及局部極值和全局最優(yōu)解的位置之后，粒子根據(jù)下面兩個公式來更新自己的位置和速度：

=·,+·rand(,)·(,-,)+

·rand(,)·(,-,)

(16)

=,+,

(17)

式中:,是粒子的第維速度；是粒子速度的慣性系數(shù)；、是兩個非負常數(shù)，稱為學(xué)習因子，分別調(diào)節(jié)向個體最優(yōu)解和全局最優(yōu)解方向飛行的步長；rand(,)是隨機函數(shù)，產(chǎn)生在區(qū)間[0，1]的隨機數(shù)；,是粒子的第維位置；,是粒子的局部極值在第維的位置；,是粒子的全局最優(yōu)解在第維的位置。

由于表尺、程裝等參數(shù)變化對射程有較大影響，因此將具有決定性作用的表尺、程裝等作為變量，為其初始配置位置與速度，生成一群隨機粒子。根據(jù)各項約束條件及期望射程，設(shè)計相應(yīng)的復(fù)合適應(yīng)度函數(shù)。設(shè)定學(xué)習因子，通過速度和位置的更新得到新的粒子群，利用經(jīng)校驗的彈道模型，計算各粒子的適應(yīng)度值并選擇當前的個體極值，結(jié)合全局極值，為下一次的粒子更新提供依據(jù)。經(jīng)過多次迭代，最優(yōu)粒子即為最終的最優(yōu)解，作為諸元取值。適應(yīng)度函數(shù)可表示為：

(18)

圖6 解算流程

4 仿真與飛行試驗驗證

為驗證經(jīng)PSO優(yōu)化后模型的準確性和可行性，試驗分兩部分進行?？尚行允侵蛤炞C彈道解算模型能否應(yīng)用在實際中，是驗證中最為重要的一環(huán)。以2021年9月7日在某靶場一次飛行實驗為例，進行驗證。

4.1 仿真驗證

此次飛行試驗的期望飛行距離為5 020 m、5 508 m、5 975 m，為近區(qū)彈道，根據(jù)各裝藥號的射程覆蓋，3次試驗均選擇初始速度為326 m/s的炮彈。采用氣象雷達車測量了地面和高空的氣象條件，包括氣溫、氣壓、風速、風向等。根據(jù)雷達測量，炮位基準位置為1 640 m。以期望射程5 000 m為例，將氣象及期望射程等參數(shù)輸入到模型中，通過粒子群算法優(yōu)化后，輸出的表尺及飛行總時間(可求出程裝、延時)值如圖7所示。

圖7 表尺及飛行時間的預(yù)測值

如圖7可知，表尺為322.933 mil，飛行時間為20.049 9 s。通常情況下射角為整數(shù)值，飛行時間取兩位小數(shù)，因此射角為323 mil，飛行時間為20.05 s。同樣的方法計算得出5 500 m時，射角為379 mil，飛行時間為23.12 s；在6 000 m時，發(fā)射角為449 mil，飛行時間為26.75 s。無控飛行彈道如圖8所示。

圖8 飛行彈道

從圖8可以看出，實際射程與期望射程的誤差滿足諸元解算的要求。

4.2 實際打靶驗證

在某靶場開展激光末制導(dǎo)炮彈飛行試驗，對諸元優(yōu)化方法進行進一步驗證。在3個射程下，通過文中的優(yōu)化方法，得到相應(yīng)的諸元參數(shù)，如表3所示。

表3 射擊諸元

根據(jù)上述諸元進行射擊，均命中目標，圖9為實際中靶效果圖。

圖9 實際中靶效果圖

從圖9可看出，基于粒子群優(yōu)化的射擊諸元均可命中目標，脫靶量較小。通過實際飛行實驗，進一步驗證了這種方法的有效性。

5 結(jié)論

通過龍格-庫塔法對彈道進行解算，給出氣動參數(shù)辨識方法分析了表尺、程裝等參數(shù)變化對射程將產(chǎn)生影響；應(yīng)用粒子群算法優(yōu)化得到最優(yōu)解，作為射表相應(yīng)的諸元取值。通過仿真和實際打靶驗證說明了彈丸模型經(jīng)粒子群算法優(yōu)化后準確性高、可行性強。

該模型對實際工程應(yīng)用提供了參考。