基于PSO優化的激光末制導炮彈諸元解算方法

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

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

0 引言

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

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

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

2)需根據當地氣象條件，查詢預先裝訂的射表，修正上述各影響因素。但這種方法解算精度較差，降低作戰效能。

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

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

1 彈丸的彈道模型

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

1.1 氣動力及作用力分析

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

表1 氣動布局

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

圖1 彈目之間角度關系

1)出炮口無控段：出炮口后末制導炮彈以低速旋轉的無控方式飛行，直至發動機工作；2)助推發動機工作段：發動機開始工作使炮彈的速度增加，從而延長飛行距離；3)助推發動機結束工作后無控段：助推發動機工作結束后，炮彈仍以低速旋轉的無控方式飛行；4)慣性制導段：陀螺解鎖后，慣性系統開始工作，為增加射程及末制導距離，此階段通過重力補償的方式使炮彈基本按直線飛行；5)末制導段：激光開始照射發出信號，當炮彈導引頭接收到目標反射的信號時，炮彈開始轉入末制導段，直至命中目標。

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

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

(1)

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

=+(-)

(2)

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

=

(3)

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

=

(4)

式中：為上下舵的舵偏轉角，沿右偏為正。

1.2 彈道方程建立

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

(5)

式中：為炮彈實時質量;為炮彈飛行速度；為發動機推力，由氣動數據進行插值得到;、分別為攻角、側滑角；為彈道傾角;為速度傾斜角;為彈道偏角。繞質心轉動動力學方程為：

(6)

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

(7)

式中：、、分別為飛行器質心的位置坐標。繞質心轉動運動學方程為：

(8)

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

(9)

式中：為質量秒流量。角度關系為：

(10)

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

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

圖2 末制導原理框圖

1.3 氣象模型

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

表2 氣象數據

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

(11)

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

(12)

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

(13)

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

2.1 模型方程解算

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

(14)

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

(15)

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

2.2 彈道模型校驗

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

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

步驟1：輸入初始數據。參數的初估值，狀態初值；待估計參數的上限和下限值；動力學系統試驗的觀測量和控制量輸入的實測值，即經預處理后的雷達數據。

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

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

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

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

3 基于粒子群的諸元解算

3.1 參數分析

分析表尺、程裝等參數變化對射程將產生影響。采用控制變量法，分別對遠區及近區兩種運動狀態進行分析。

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

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

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

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

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

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

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

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

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

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

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

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

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

·rand(,)·(,-,)

(16)

=,+,

(17)

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

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

(18)

圖6 解算流程

4 仿真與飛行試驗驗證

為驗證經PSO優化后模型的準確性和可行性，試驗分兩部分進行。可行性是指驗證彈道解算模型能否應用在實際中，是驗證中最為重要的一環。以2021年9月7日在某靶場一次飛行實驗為例，進行驗證。

4.1 仿真驗證

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

圖7 表尺及飛行時間的預測值

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

圖8 飛行彈道

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

4.2 實際打靶驗證

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

表3 射擊諸元

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

圖9 實際中靶效果圖

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

5 結論

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

該模型對實際工程應用提供了參考。