再入飛行器禁飛區(qū)規(guī)避近似解析制導(dǎo)方法

趙亮博，朱廣生，莊 凌，宋加洪，盧寶剛

（1.北京航天長(zhǎng)征飛行器研究所，北京，100076；2.中國(guó)運(yùn)載火箭技術(shù)研究院，北京，100076）

0 引 言

再入飛行器的機(jī)動(dòng)制導(dǎo)問(wèn)題，一直是再入飛行器制導(dǎo)控制領(lǐng)域的研究熱點(diǎn)[1～3]。飛行器再入制導(dǎo)過(guò)程具有非線性、時(shí)變性和不確定性等特點(diǎn)，其飛行路徑除了需要滿足過(guò)載、動(dòng)壓、熱流率約束外，還需要滿足由自然、軍事等因素造成的禁飛區(qū)約束。

禁飛區(qū)約束是一種復(fù)雜的路徑約束，現(xiàn)有禁飛區(qū)規(guī)避方法主要集中在離線最優(yōu)軌跡解算中，如經(jīng)典的數(shù)值求解方法偽譜法[4]，通過(guò)不同配點(diǎn)和多項(xiàng)式插值方式實(shí)現(xiàn)全局軌跡的優(yōu)化求解。但是離線求解方法解算效率低，無(wú)法滿足實(shí)時(shí)制導(dǎo)需求。隨著探測(cè)技術(shù)發(fā)展，對(duì)禁飛區(qū)域的實(shí)時(shí)或近似實(shí)時(shí)探測(cè)成為可能，因此需要研究實(shí)現(xiàn)禁飛區(qū)在線規(guī)避的制導(dǎo)算法。文獻(xiàn)[5]提出了一種考慮禁飛區(qū)的預(yù)測(cè)校正制導(dǎo)方法，利用航向角誤差走廊的形式控制飛行器的側(cè)向運(yùn)動(dòng)，當(dāng)禁飛區(qū)等約束急劇變化時(shí)預(yù)測(cè)校正的制導(dǎo)方法易產(chǎn)生性能退化。文獻(xiàn)[6]利用一種基于三觸角預(yù)測(cè)方法的側(cè)滑角瞬變策略，實(shí)現(xiàn)了禁飛區(qū)的規(guī)避，但機(jī)動(dòng)軌跡的生成依賴于飛行器動(dòng)力學(xué)方程積分而來(lái)的3個(gè)“觸角”，計(jì)算代價(jià)大。

基于Dubins曲線的路徑規(guī)劃方法在機(jī)器人、無(wú)人機(jī)的軌跡規(guī)劃中得到了廣泛應(yīng)用[7～9]。文獻(xiàn)[10]第1次提出Dubins路徑的概念，在轉(zhuǎn)彎半徑對(duì)車輛運(yùn)動(dòng)影響的分析上，設(shè)計(jì)了滿足約束條件的最短路徑并進(jìn)行了理論證明，即滿足曲率限制的最短曲線由圓弧與直線構(gòu)成。文獻(xiàn)[11]設(shè)計(jì)了基于 Dubins曲線的高超聲速飛行器離線路徑規(guī)劃方法，但未結(jié)合具體的飛行器轉(zhuǎn)彎能力和終端約束進(jìn)行分析。

本文在對(duì)再入飛行器轉(zhuǎn)彎能力分析的基礎(chǔ)上，得出飛行器轉(zhuǎn)彎半徑和傾側(cè)角的近似解析關(guān)系，結(jié)合Dubins曲線的路徑規(guī)劃方法生成禁飛區(qū)規(guī)避的制導(dǎo)策略。同時(shí)為修正規(guī)避引起的終端誤差，建立基于能量運(yùn)動(dòng)模型來(lái)解析預(yù)測(cè)航程，通過(guò)割線法對(duì)傾側(cè)角指令進(jìn)行迭代修正，實(shí)現(xiàn)航程校正控制。通過(guò)速度傾角的反饋控制設(shè)計(jì)，滿足終端高度約束。最后基于CAV-H的數(shù)學(xué)模型仿真驗(yàn)證了該算法在禁飛區(qū)繞飛規(guī)避的有效性，滿足終端約束條件且計(jì)算效率高、實(shí)時(shí)性好。

1 飛行器運(yùn)動(dòng)建模

1.1 再入運(yùn)動(dòng)方程

本文研究的飛行器采用無(wú)動(dòng)力滑翔的方式再入，采取傾斜轉(zhuǎn)彎（Bank To Turn，BTT）的機(jī)動(dòng)策略。為研究飛行器在大氣層內(nèi)的運(yùn)動(dòng)規(guī)律、分析飛行器轉(zhuǎn)彎能力、設(shè)計(jì)制導(dǎo)策略，在不考慮地球自轉(zhuǎn)影響的情況下，建立飛行器再入運(yùn)動(dòng)方程組為

式中r,λ,φ,V,θ,ψ分別為地心距、經(jīng)度、緯度、速度、速度傾角和航向角；σ為傾側(cè)角；g為重力加速度；L，D分別為氣動(dòng)升力和阻力，且：

式中ρ為大氣密度；Sref為參考面積；CL，CD分別為升力系數(shù)和阻力系數(shù)。

由于飛行結(jié)束時(shí)間不確定，給航程解析預(yù)測(cè)帶來(lái)了困難。為此，定義如下形式的無(wú)量綱變量e能量：

式中為無(wú)量綱的地心距，無(wú)量綱化參數(shù)為地球半徑R0；為無(wú)量綱速度，無(wú)量綱化參數(shù)為，g0為海平面重力加速度。為解析預(yù)測(cè)航程需要，建立如下歸一化的基于能量的運(yùn)動(dòng)方程：

式中為無(wú)量綱阻力加速度，無(wú)量綱系數(shù)為g0。

1.2 再入約束建模

本文考慮的終端約束有高度約束、速度約束及終端位置約束，以能量為自變量時(shí)表達(dá)式如下：

式中ef,rf,Vf,λf,φf(shuō)分別為以能量為自變量描述的終端能量、地心距、速度、經(jīng)度和緯度。

由于地理、軍事等因素會(huì)造成飛行器不能在某些區(qū)域飛行，即禁飛區(qū)，通常用無(wú)限高圓柱表示。用經(jīng)緯度及禁飛區(qū)半徑rN描述禁飛區(qū)如下：

式中λz,φz分別為禁飛區(qū)中心的經(jīng)度和緯度。

函數(shù)s((λz,φz), (λ,φ) )表示兩點(diǎn)間的大地線長(zhǎng)度：

2 制導(dǎo)策略

本文將制導(dǎo)策略分為禁飛區(qū)規(guī)避制導(dǎo)和完成規(guī)避后的校正策略兩部分設(shè)計(jì)，具體算法流程如圖1所示。

圖1 制導(dǎo)算法流程 Fig.1 Guidance Algorithm Flow

2.1 轉(zhuǎn)彎能力分析

為設(shè)計(jì)禁飛區(qū)規(guī)避策略，需要對(duì)飛行器的轉(zhuǎn)彎能力進(jìn)行分析。對(duì)于采用BTT機(jī)動(dòng)方式的飛行器，通過(guò)傾斜升力面的方式使升力產(chǎn)生側(cè)向分量，提供轉(zhuǎn)彎的需用向心力。

假設(shè)轉(zhuǎn)彎過(guò)程為無(wú)側(cè)滑的平衡飛行狀態(tài)，即轉(zhuǎn)彎過(guò)程中速度傾角近似為零：

轉(zhuǎn)彎過(guò)程中的需用向心力由升力的分量Lsinσ提供，定義轉(zhuǎn)彎半徑為rz，可得：

聯(lián)立式（8）、式（9）可得：

在水平轉(zhuǎn)動(dòng)平面內(nèi)忽略離心力項(xiàng)，可得機(jī)動(dòng)轉(zhuǎn)彎半徑的近似解析解：

因此，一定轉(zhuǎn)彎半徑rz下傾側(cè)角σ的近似解析解：

為驗(yàn)證該公式，針對(duì)CAV-H飛行器在高度50 km、速度 5000 m/s和高度 55 km、速度 5500 m/s兩組不同的工況進(jìn)行仿真驗(yàn)證。

圖2 工況一：高度 50 km，速度 5000 m/s Fig.2 Condition 1: height 50 km, velocity 5000 m/s

圖3 工況二：高度 55 km，速度 5500 m/s Fig.3 Condition 2: height 55 km, velocity 5500 m/s

該公式建立在平衡滑翔飛行的基礎(chǔ)上，可以用作該狀態(tài)下傾側(cè)角和轉(zhuǎn)彎半徑的解析計(jì)算。從圖2、圖3中可以看出，當(dāng)傾側(cè)角大于25°時(shí)解析解和數(shù)值解誤差較小，式（12）可為機(jī)動(dòng)策略的選取提供依據(jù)。

此外，在實(shí)際轉(zhuǎn)彎過(guò)程中由于升力在鉛垂方向上的分量減小，如果不對(duì)攻角加以調(diào)節(jié)會(huì)引起飛行高度的下降。為方便設(shè)計(jì)參考，本文給出了無(wú)攻角聯(lián)合調(diào)節(jié)時(shí)機(jī)動(dòng)轉(zhuǎn)彎引起的高度損失近似估算公式。

假設(shè)轉(zhuǎn)彎前后飛行器的高度分別為h1,h2，大氣密度分別為ρ1,ρ2，根據(jù)大氣密度的近似估算公式可得：

轉(zhuǎn)彎后為維持平衡飛行，推導(dǎo)出大氣密度變化規(guī)律：

因此，高度變化量Δh和傾側(cè)角間的近似解析關(guān)系為

該關(guān)系式可用做轉(zhuǎn)彎引起高度損失的近似估算，對(duì)機(jī)動(dòng)策略的選取具有參考意義。

2.2 基于Dubins曲線的禁飛區(qū)規(guī)避策略

飛行器由起始狀態(tài)調(diào)整至目標(biāo)狀態(tài)所經(jīng)過(guò)的路徑是經(jīng)計(jì)算選擇的一個(gè)起始轉(zhuǎn)向圓rz1與一個(gè)目標(biāo)轉(zhuǎn)向圓rz2之間通過(guò)合適的公切線相連得到的有向曲線，此有向曲線即是一條Dubins路徑[10]，如圖4所示。

圖4 Dubins路徑 Fig.4 Dubins Path

依據(jù) Dubins曲線的組成原理確定本文的繞飛策略，其中需要求解的變量包括 2次轉(zhuǎn)彎的轉(zhuǎn)彎半徑rz1,rz2及轉(zhuǎn)彎角度θ1,θ2。同時(shí)，針對(duì)轉(zhuǎn)彎過(guò)程中可能出現(xiàn)的高度損失，采用攻角和傾側(cè)角聯(lián)合調(diào)整的方法設(shè)計(jì)攻角大小。具體的規(guī)避流程如圖5所示。

圖5 基于Dubins曲線的規(guī)避流程 Fig.5 Avoiding Process Based on Dubins Curve

假設(shè)對(duì)禁飛區(qū)的探測(cè)距離為rdet，即在距離rdet時(shí)發(fā)現(xiàn)禁飛區(qū)。對(duì)于半徑為rz大小的禁飛區(qū)，Dubins曲線起始圓的最小轉(zhuǎn)彎半徑rz1計(jì)算公式為

由式（16）可得，飛行器探測(cè)到禁飛區(qū)的時(shí)機(jī)越早、探測(cè)距離越遠(yuǎn)，需用機(jī)動(dòng)轉(zhuǎn)彎半徑越大。聯(lián)合式（12）與式（16）即可求出規(guī)避禁飛區(qū)的需用傾側(cè)角σrz1。同時(shí)，根據(jù)相切點(diǎn)位置通過(guò)解析幾何的方式可求出轉(zhuǎn)彎角度1θ的大小。

通過(guò)目標(biāo)圓的半徑rz2選取對(duì)應(yīng)傾側(cè)角σrz2大小時(shí)需進(jìn)行判斷，假設(shè)當(dāng)前飛行狀態(tài)下轉(zhuǎn)彎半徑rz對(duì)應(yīng)的傾側(cè)角為σrz，通過(guò)式（23）得出的傾側(cè)角大小為σcmd，則：

然后，通過(guò)判斷航向角誤差是否為0，即飛行器是否瞄準(zhǔn)目標(biāo)點(diǎn)，作為轉(zhuǎn)彎結(jié)束的標(biāo)志并進(jìn)入校正環(huán)節(jié)，不再計(jì)算目標(biāo)圓旋轉(zhuǎn)角度2θ。

由式（15）可得，轉(zhuǎn)彎時(shí)如果不對(duì)總升力進(jìn)行調(diào)節(jié)，會(huì)引起飛行器的高度損失。本文通過(guò)聯(lián)合調(diào)整攻角的方式來(lái)調(diào)整升力，補(bǔ)償高度損失，近似平衡滑翔飛行時(shí)存在如下關(guān)系式：

可求得該飛行條件下的需用升力，通過(guò)氣動(dòng)反插值的方式求出需用攻角。

需要指出的是，當(dāng)傾側(cè)過(guò)大時(shí)為了補(bǔ)償高度損失所需用的攻角較大，大攻角飛行時(shí)升阻比低會(huì)產(chǎn)生較大的速度損失。因此需設(shè)定合理的攻角聯(lián)合調(diào)節(jié)范圍［αlow,αmax］，針對(duì)CAV-H飛行器，本文選取的攻角調(diào)節(jié)區(qū)間為[8°，15°]。

2.3 航程解析預(yù)測(cè)校正控制

完成禁飛區(qū)規(guī)避后，飛行器高度和速度會(huì)產(chǎn)生相應(yīng)的損失。為滿足終端約束，需對(duì)航程及高度進(jìn)行校正。本文通過(guò)傾側(cè)角大小的調(diào)節(jié)，實(shí)現(xiàn)航程的校正。根據(jù)文獻(xiàn)[12]的推導(dǎo)，基于能量的運(yùn)動(dòng)模型進(jìn)行航程預(yù)測(cè)：

式中(ef)為從能量e到能量ef的歸一化預(yù)測(cè)航程，歸一化變量為地球半徑R0。此時(shí)，預(yù)測(cè)航程和待飛航程的差值：

式中s為當(dāng)前真實(shí)的待飛航程，可通過(guò)式（7）并結(jié)合當(dāng)前點(diǎn)和終點(diǎn)的經(jīng)緯度坐標(biāo)求出。

通過(guò)調(diào)節(jié)傾側(cè)角大小的方式實(shí)現(xiàn)對(duì)航程的在線校正，采用割線法不斷迭代傾側(cè)角大小：

2.4 終端高度控制

在再入飛行的大部分時(shí)間內(nèi)，飛行攻角處于最大升阻比大小附近，通常采用攻角剖面的形式設(shè)計(jì)攻角大小，同時(shí)實(shí)現(xiàn)滿足過(guò)程約束：

式中αm為較大攻角，使飛行器在初始階段快速達(dá)到平衡滑翔的飛行狀態(tài)；αL/D為最大升阻比攻角；VL/D和Vm為設(shè)計(jì)變量，實(shí)現(xiàn)對(duì)不同飛行階段攻角大小的控制。

在飛行末段，為實(shí)現(xiàn)滿足終端高度約束，不采用攻角剖面的形式。本文通過(guò)設(shè)計(jì)一種速度傾角反饋控制的策略調(diào)整攻角大小，期望速度傾角dθ設(shè)計(jì)為

設(shè)計(jì)線性反饋跟蹤控制律：

聯(lián)合式（1）可求得需用升力Lc：

同樣，通過(guò)氣動(dòng)反插值的方式求出需用攻角：

對(duì)于航向的控制采用傳統(tǒng)航向角誤差走廊的形式實(shí)現(xiàn)，在此不做贅述。

3 仿真結(jié)果及分析

本文采用CAV-H飛行器的數(shù)學(xué)模型，進(jìn)行仿真驗(yàn)證[14]，仿真初始及終端條件見(jiàn)表1。

表1 初始條件及終端條件 Tab.1 Initial Conditions and Terminal Conditions

假設(shè)飛行器探測(cè)距離為2000 km，探測(cè)到禁飛區(qū)中心經(jīng)緯度為（30°，40°），禁飛區(qū)半徑500 km。通過(guò)Dubins曲線的路徑規(guī)劃策略生成的起始圓起點(diǎn)經(jīng)緯度（20.92°，23.71°），旋轉(zhuǎn)角度θ1=28.07°，旋轉(zhuǎn)半徑rz1=3750 km，對(duì)應(yīng)的近似解析傾側(cè)角大小σrz1=46.99°；校正制導(dǎo)起點(diǎn)經(jīng)緯度（50.34°，48.80°），具體仿真結(jié)果如圖6～9所示。

圖6 CAV-H規(guī)避禁飛區(qū)三維軌跡 Fig.6 Three-dimensional Trajectory Formed by CAV-H Avoiding the No-fly Zone

圖7 CAV-H規(guī)避禁飛區(qū)星下點(diǎn)軌跡 Fig.7 Subsatellite Track Formed by CAV-H Avoiding the No-fly Zone

由圖6、圖7可得，飛行器能夠有效實(shí)現(xiàn)禁飛區(qū)的規(guī)避，其終端脫靶量611.80 m，精度較高。

圖8為飛行器的高度、速度曲線，終端高度誤差15.42 m、速度誤差14.03 m/s，均在合理偏差范圍內(nèi)，且高度及速變化平穩(wěn)。

圖8 高度和速度曲線 Fig.8 Height and Velocity Curve

圖9為飛行器的攻角和傾側(cè)角曲線。

圖9 攻角和傾側(cè)角曲線 Fig.9 Angle of Attack and Bank Curve

飛行器分別在519.8 s和809.8 s進(jìn)行了機(jī)動(dòng)轉(zhuǎn)彎，實(shí)現(xiàn)了禁飛區(qū)的繞飛規(guī)避和目標(biāo)的重新瞄準(zhǔn)。在接近目標(biāo)時(shí)，攻角逐漸增大校正終端高度，同時(shí)連續(xù)調(diào)整傾側(cè)符號(hào)修正航向誤差。

4 結(jié) 論

本文提出了一種基于近似解析解的制導(dǎo)方法，能夠有效實(shí)現(xiàn)禁飛區(qū)的繞飛規(guī)避，同時(shí)滿足終端航程及高度約束，該方法具有以下特點(diǎn)：

a）對(duì)飛行器的轉(zhuǎn)彎能力進(jìn)行了分析，得到了機(jī)動(dòng)轉(zhuǎn)彎半徑和傾側(cè)角的近似解析關(guān)系；

b）結(jié)合 Dubins曲線的路徑規(guī)劃方法設(shè)計(jì)規(guī)避策略，結(jié)構(gòu)簡(jiǎn)單、效率高；

c）通過(guò)對(duì)航程的解析預(yù)測(cè)和速度傾角的反饋設(shè)計(jì)，實(shí)現(xiàn)了終端航程及高度的控制，計(jì)算效率高、誤差小。

通過(guò)數(shù)學(xué)仿真驗(yàn)證了算法的有效性，得到的近似解析關(guān)系對(duì)機(jī)動(dòng)制導(dǎo)算法的設(shè)計(jì)也具有參考價(jià)值，工程實(shí)用性強(qiáng)。