劉 君，陳克俊，謝 愈，湯國(guó)建

(國(guó)防科學(xué)技術(shù)大學(xué)航天與材料工程學(xué)院，長(zhǎng)沙 410073)

0 引言

近年來(lái)，隨著反導(dǎo)技術(shù)的提高，傳統(tǒng)彈道導(dǎo)彈的生存能力受到嚴(yán)重挑戰(zhàn)，由于助推滑翔飛行器飛行軌跡難以預(yù)測(cè)，容易突防，因此助推滑翔飛行器日益成為各國(guó)研究的熱點(diǎn)。

發(fā)射諸元是指導(dǎo)彈從發(fā)射點(diǎn)到目標(biāo)點(diǎn)的一條飛行軌跡的特征參數(shù)。助推滑翔飛行器采用非彈道式飛行軌跡，同彈道導(dǎo)彈相比，飛行軌跡復(fù)雜、發(fā)射諸元數(shù)目眾多，且發(fā)射諸元方程組為不定方程組，而傳統(tǒng)的諸元計(jì)算方法如模擬打靶法、區(qū)間插值法、參數(shù)擬合法等[3～6]無(wú)法解決不定方程組問(wèn)題。因此，助推滑翔飛行器發(fā)射諸元的快速計(jì)算是一個(gè)全新的課題。

文中以助推滑翔飛行器發(fā)射諸元計(jì)算方法研究為背景，建立了助推滑翔飛行器動(dòng)力學(xué)模型，討論了助推滑翔飛行器發(fā)射諸元的選取及計(jì)算難點(diǎn)，通過(guò)分析助推滑翔飛行器發(fā)射諸元對(duì)飛行器飛行軌跡的影響，提出了助推滑翔飛行器發(fā)射諸元計(jì)算策略，采用牛頓迭代模擬打靶法解決助推滑翔飛行器的諸元快速計(jì)算問(wèn)題。

1 助推滑翔飛行器動(dòng)力學(xué)模型

助推滑翔飛行器由助推火箭和滑翔飛行器兩部分構(gòu)成。助推火箭首先將滑翔飛行器運(yùn)送到高空，然后與之分離，緊接著滑翔飛行器再入大氣層后，由氣動(dòng)力實(shí)現(xiàn)變軌進(jìn)入無(wú)動(dòng)力滑翔飛行階段。助推滑翔飛行器飛行軌跡可劃分為助推段、變軌段和滑翔段，

其中變軌段和滑翔段也可統(tǒng)稱為滑翔段。在發(fā)射坐標(biāo)系o－xyz中，助推滑翔飛行器的質(zhì)心運(yùn)動(dòng)數(shù)學(xué)模型[2]可表示為:

式中:V為飛行器的飛行速度;ρ為飛行器的位置坐標(biāo);m為飛行器的質(zhì)量;P為發(fā)動(dòng)機(jī)推力;Fc為控制力;Fa為作用在飛行器上的空氣動(dòng)力;g為地球引力加速度;ae為牽連加速度;ac為哥氏加速度;ak為附加哥氏加速度。各加速度描述為:

牽連加速度:ae=ωe× ωe×( r)

哥氏加速度:ac=2ωe×V

式中:˙m為發(fā)動(dòng)機(jī)秒耗量;ωT為飛行器相對(duì)于慣性坐標(biāo)系的轉(zhuǎn)動(dòng)角速度矢量;ρe為質(zhì)心到噴管出口中心點(diǎn)距離。

助推段結(jié)束后，滑翔飛行器的飛行軌跡及姿態(tài)調(diào)整均通過(guò)氣動(dòng)力作用實(shí)現(xiàn)，在此階段飛行器質(zhì)心運(yùn)動(dòng)數(shù)學(xué)模型式(1)中發(fā)動(dòng)機(jī)推力P和附加哥氏加速度ak均為零。

2 助推滑翔飛行器發(fā)射諸元

助推滑翔飛行器飛行軌跡由助推段、變軌段和滑翔段的飛行特征參數(shù)共同確定，這些飛行特征參數(shù)均是助推滑翔飛行器的發(fā)射諸元。

圖1 助推滑翔飛行器助推段飛行程序角

2.1 助推段發(fā)射諸元

助推滑翔飛行器助推段飛行軌跡與傳統(tǒng)彈道導(dǎo)彈主動(dòng)段相似，可采用相同的飛行程序角模型[2](見(jiàn)圖 1):飛行器垂直起飛后不久，飛行器在稠密大氣層內(nèi)實(shí)行有攻角轉(zhuǎn)彎，在飛行速度接近跨音速段之前，使攻角收縮為零;飛行速度達(dá)到音速后在重力作用下實(shí)現(xiàn)轉(zhuǎn)彎;當(dāng)飛行器飛出稠密大氣層后，在飛行軌跡設(shè)計(jì)時(shí)，通過(guò)程序角速度˙φ確保飛行器達(dá)到預(yù)定的再入狀態(tài);為提高飛行器的穩(wěn)定性，在發(fā)動(dòng)機(jī)啟動(dòng)、關(guān)機(jī)以及級(jí)間分離前后，要求飛行器保持常值程序角飛行，并盡量減小飛行攻角。據(jù)此，飛行程序角模型可取為:其中:t1、t2分別為攻角轉(zhuǎn)彎起始和結(jié)束時(shí)間;t3、t4分別為程序角速度轉(zhuǎn)彎起始和結(jié)束時(shí)間;t5為發(fā)動(dòng)機(jī)關(guān)機(jī)時(shí)間;αm為最大飛行攻角;Km為攻角參數(shù)。對(duì)于確定的某型飛行器而言，t1、t2、t3、t4、t5均可以取為常數(shù)。同時(shí)，由于地球形狀不規(guī)則以及自轉(zhuǎn)影響，飛行器的飛行軌跡還與發(fā)射點(diǎn)經(jīng)緯度(λ0，φ0)、發(fā)射方位角A0有關(guān)。因此，助推滑翔飛行器助推段發(fā)射諸元包括發(fā)射點(diǎn)經(jīng)緯度(λ0，φ0)、發(fā)射方位角 A0、最大飛行攻角αm、程序角速度˙φ以及攻角參數(shù)Km。

2.2 滑翔段發(fā)射諸元

助推段結(jié)束后滑翔飛行器的飛行軌跡與彈道導(dǎo)彈被動(dòng)段完全不同，其飛行軌跡可細(xì)分為:變軌段與滑翔段，其每一階段的飛行軌跡均有其獨(dú)特的飛行特征參數(shù)。

滑翔飛行器變軌段屬于無(wú)動(dòng)力飛行，其目的在于在末端獲得滿足滑翔的初始條件。飛行器的初始滑翔位置并不固定，初始滑翔條件一般是要求高度H、速度傾角滿足某特定條件，變軌段飛行器需要充分利用氣動(dòng)力作用進(jìn)行飛行軌跡調(diào)整，因此可將飛行程序攻角αc作為發(fā)射諸元。達(dá)到初始滑翔條件后飛行器進(jìn)入無(wú)動(dòng)力滑翔階段，此階段諸元與制導(dǎo)方法有關(guān)，制導(dǎo)標(biāo)準(zhǔn)參數(shù)可視作發(fā)射諸元，這里為了簡(jiǎn)便的說(shuō)明發(fā)射諸元計(jì)算的原理，不再增加諸元，制導(dǎo)標(biāo)準(zhǔn)參數(shù)視為固定值。當(dāng)飛行器滿足終端條件，計(jì)算得到相對(duì)落點(diǎn)偏差(ΔL，ΔZ)。變軌段初始滑翔條件高度H和速度傾角約束實(shí)際是一個(gè)互相制約的條件，這里選取滑翔起點(diǎn)高度H作為約束條件。

綜上所述，助推滑翔飛行器發(fā)射諸元除了具有與彈道導(dǎo)彈相同的助推段諸元:發(fā)射方位角A0、最大程序攻角αm、程序角速度˙φ、攻角參數(shù)Km，還增加了滑翔段諸元程序攻角αc等，發(fā)射諸元函數(shù)方程如式(3)所示:

由此可見(jiàn):1)助推滑翔飛行器發(fā)射諸元眾多;2)發(fā)射諸元函數(shù)方程(3)只有3個(gè)終端約束條件卻有5個(gè)未知量，是一個(gè)典型的非線性不定方程組，具有無(wú)窮多解。這些都為助推滑翔飛行器發(fā)射諸元的計(jì)算制造了困難。

3 諸元計(jì)算策略

按照常規(guī)的方法，顯然非線性不定方程組(3)很難求解。為了尋求助推滑翔飛行器發(fā)射諸元的計(jì)算方法，首先分析發(fā)射諸元對(duì)飛行器飛行軌跡的影響。為了研究方便，文中采用影響因子E來(lái)分析發(fā)射諸元對(duì)彈道的影響，影響因子E定義如下:

其中:?X/X為自變量的變化比值，?Y/Y為?X引起的因變量變化比值。影響因子E的值越大表明自變量對(duì)因變量的影響越大。

圖2 飛行軌跡示意圖

圖2 為助推滑翔飛行器飛行軌跡示意圖，N為北極，O為地心，飛行器從發(fā)射點(diǎn)L起飛，在C點(diǎn)降落。弧線LC對(duì)應(yīng)的地心角β稱之為射程角，弧線LC與弧線LN構(gòu)成的球面角Ad稱

之為落點(diǎn)方位角。助推滑翔飛行器的飛行軌跡可以由射程角β、落點(diǎn)方位角Ad、滑翔起點(diǎn)高度H表示。

1)發(fā)射方位角對(duì)彈道的影響

由表1分析可知，發(fā)射方位角對(duì)射程角β、落點(diǎn)方位角Ad、滑翔起點(diǎn)高度H均有較大的影響。

表1 發(fā)射方位角影響因子

2)飛行程序?qū)椀赖挠绊?/p>

由表2的影響因子分析可知，主動(dòng)段飛行程序參數(shù)最大飛行攻角αm、程序角速度˙φ、攻角參數(shù)Km對(duì)射程角β和滑翔起點(diǎn)高度H有較大影響，而對(duì)落點(diǎn)方位角Ad的影響較小，變軌段飛行程序參數(shù)程序攻角αc均對(duì)滑翔起點(diǎn)高度H有較大影響，而對(duì)其它兩個(gè)彈道參數(shù)影響較小。

表2 飛行程序影響因子

綜上所述，飛行器射程的遠(yuǎn)近主要是由發(fā)射角A0、最大飛行攻角αm、程序角速度˙φ所主導(dǎo)，而飛行器落點(diǎn)的方位角Ad則主要是由發(fā)射角A0決定的。因此，根據(jù)發(fā)射諸元對(duì)飛行軌跡的影響大小可將諸元分成兩組:一組為發(fā)射方位角A0、最大飛行攻角αm、程序角速度˙φ;另一組為發(fā)射方位角A0、攻角參數(shù)Km以及程序攻角αc。這樣發(fā)射諸元函數(shù)方程組可分解為以下兩個(gè)恰定方程組:

所以助推滑翔飛行器發(fā)射諸元的可分為如下兩步:第一步，固定參數(shù)Km與程序攻角αc，放寬落點(diǎn)與滑翔起點(diǎn)高度精度要求，利用牛頓迭代算法求解方程組(5)解出發(fā)射方位角A0粗略值以及最大飛行攻角αm、三級(jí)程序角速度˙φ3。第二步，在前步計(jì)算的基礎(chǔ)上固定最大飛行攻角αm、程序角速度˙φ，利用牛頓迭代算法求解方程組(6)解出最終的發(fā)射方位角A0精確值、攻角參數(shù)Km和程序攻角αc。

4 算法設(shè)計(jì)

模擬打靶法基本原理如下:由發(fā)射點(diǎn)和目標(biāo)點(diǎn)位置確定初始積分條件，對(duì)于給定發(fā)射方位角A、最大飛行攻角αm積分彈道求得落點(diǎn)的位置坐標(biāo)，與目標(biāo)點(diǎn)位置坐標(biāo)比較確定落點(diǎn)計(jì)算誤差，以此誤差修正A、αm，再次積分彈道，求得新落點(diǎn)，經(jīng)過(guò)如此若干次反復(fù)計(jì)算，當(dāng)落點(diǎn)計(jì)算誤差達(dá)到要求時(shí)，最終可求得給定發(fā)射點(diǎn)和目標(biāo)點(diǎn)的一條彈道。

文中采用牛頓迭代模擬打靶算法解算發(fā)射諸元[1]。為表述方便，可用向量Y表示助推滑翔飛行器飛行軌跡參數(shù)，向量X表示發(fā)射諸元，前文描述的諸元函數(shù)關(guān)系方程式(5)、(6)可改寫為如下通用形式:

將式(7)在初始狀態(tài)參數(shù)X0處線性展開(kāi)得:

可進(jìn)一步整理式(8)得到:

圖3 計(jì)算流程圖

5 算例仿真分析及結(jié)論

給定表3中4組初始諸元計(jì)算條件，采用文中設(shè)計(jì)的助推滑翔飛行器諸元計(jì)算方法計(jì)算得到的發(fā)射諸元如表4所示，將這些諸元代入飛行仿真程序，得到飛行器落點(diǎn)及滑翔起點(diǎn)高度偏差如表5所示，圖4為表3中4組飛行軌跡圖。

圖4 助推滑翔飛行器飛行軌跡

表3 初始條件

表4 諸元計(jì)算結(jié)果

表5 仿真結(jié)果

文中基于發(fā)射諸元對(duì)助推滑翔飛行器飛行軌跡的影響分析，提出了助推滑翔飛行器諸元計(jì)算策略，并采用牛頓迭代模擬打靶法設(shè)計(jì)了諸元計(jì)算算法。實(shí)際算例中飛行器落點(diǎn)及滑翔起點(diǎn)高度偏差很小，計(jì)算迭代次數(shù)較少，能很好的滿足助推滑翔飛行器快速發(fā)射的要求。因此，文中設(shè)計(jì)的諸元計(jì)算方法是可行的。由于文中待解函數(shù)沒(méi)有解析表達(dá)式，牛頓迭代諸元計(jì)算算法的收斂性有待數(shù)學(xué)上進(jìn)一步論證，但現(xiàn)實(shí)計(jì)算證明是可以收斂的，收斂速度和收斂性與迭代初值的選取有關(guān)，其具體情況有待進(jìn)一步研究。

[1]陳克俊.載入飛船上升段軌道的 Newton迭代設(shè)計(jì)法[J].國(guó)防科技大學(xué)學(xué)報(bào)，1992，14(2):66 －71.

[2]賈沛然，陳克俊，何力.遠(yuǎn)程火箭彈道學(xué)[M].長(zhǎng)沙:國(guó)防科技大學(xué)出版社，1993:217－220.

[3]郝佳新，丁振國(guó)，陳克俊.彈道導(dǎo)彈基本諸元的快速裝訂算法研究[J].計(jì)算機(jī)仿真，2008，25(4):25 －29.

[4]雍恩米，唐國(guó)金.基于攝動(dòng)制導(dǎo)的彈道導(dǎo)彈發(fā)射諸元的仿真算法[J].系統(tǒng)仿真學(xué)報(bào)，2005，17(5):1048 －1051.

[5]王海麗，陳磊，胡小平.彈道基本諸元的快速裝訂算法研究[J].國(guó)防科技大學(xué)學(xué)報(bào)，1999，21(2):5 －8.

[6]韓東，谷宏強(qiáng)，李洪儒.彈道導(dǎo)彈射程的擬合算法[J].軍械工程學(xué)院學(xué)報(bào)，2002，14(3):26 －28.