禁飛區(qū)影響下的空天飛機可達區(qū)域計算方法

章吉力，周大鵬，楊大鵬，劉然，劉凱,*

1. 大連理工大學 航空航天學院，大連 116024

2. 航空工業(yè)沈陽飛機設計研究所，沈陽 110035

空天飛機是一種能夠快速往返于稠密大氣層、臨近空間和軌道空間的可重復使用飛行器[1-2]，既可以像普通飛機一樣水平起降，又具備臨近空間持續(xù)機動飛行能力，具有重要的政治價值、經濟價值和軍事價值。

可達區(qū)域求解[1-2]是對空天飛機的覆蓋范圍和飛行能力進行評估的重要方法，基于可達區(qū)域求解的結果進行任務的繼續(xù)與中止、目標的決策與切換能夠大幅提升飛行器自主性。可達區(qū)域的求解方法已受到國內外眾多學者的關注。

可達區(qū)域的求解方法有很多種，主要有常值傾側角法[3-4]、剖面平移設計法[5]、優(yōu)化法[6-14]、剖面參數規(guī)劃法[15-17]等。由于常值傾側角法解算出的區(qū)域存在一定誤差，與實際區(qū)域相比偏小，但同時也具有簡單易行、穩(wěn)定性高、解算速度快的特性，因而成為工程中常用的方法。剖面平移設計法首先將再入過程約束轉化到H-V剖面，形成再入走廊，之后沿著走廊下邊界進行再入剖面規(guī)劃，形成航程最短的極限剖面，再將該剖面向上平移，形成航程適中的剖面，最后以零傾側角飛行，形成航程最長的剖面[5]，將這一系列剖面的終點連接起來，即得到可達區(qū)域。

優(yōu)化法是現有最為精確的解算方法，但解算速度相對較慢。文獻[7-8]基于目標函數變更求解構成可達區(qū)域邊界的軌跡。優(yōu)化指標涉及最大縱程、最小縱程，及若干中間量縱程填補邊界曲線，最后把這些彈道的落點連起來，即得到可達區(qū)域。文獻[13]還考慮到優(yōu)化得到的控制指令的平滑性，認為直接使用偽譜法優(yōu)化傾側角和攻角指令得到的剖面可行性差，使用控制變量的二階導數作為優(yōu)化變量能夠得到更平滑的控制指令。剖面參數規(guī)劃法介于優(yōu)化法和剖面平移設計法之間，首先確定一個基準剖面，剖面上有許多待求參數，通過設置與優(yōu)化法類似的特殊優(yōu)化指標，來確定這些待求參數，文獻[15]采用粒子群優(yōu)化的方式來求解這些待定參數，文獻[16]利用差分進化算法和傾側角插值求解待定參數，權衡兼顧了可達區(qū)域精度和求解效率。

還有學者通過對可達區(qū)域特性的分析，以解析和近似的方式來求解可達區(qū)域，文獻[18]考慮到可達區(qū)域同橢圓的近似性，以兩個最大橫程終點間線段為短軸，以最大縱程終點到短軸距離為半長軸構建半橢圓可達區(qū)域。文獻[19]通過解析同倫法，在求解4個子優(yōu)化問題的基礎上，通過構造合適的同倫參數延拓出縱向可達區(qū)。也有學者將可達區(qū)域求解應用于某些特殊場景。文獻[20]考慮了在具有路徑點約束情況下的可達區(qū)域估計方法。文獻[21]通過可達區(qū)域求解結果指導在軌星座設計方法，構建目標全覆蓋的空間戰(zhàn)略武器星座。

總的來說，現有的研究中對可達區(qū)域的求解方法鮮有考慮禁飛區(qū)的影響。主要存在以下難點：一方面，禁飛區(qū)與飛行器相對位置的不確定性引發(fā)禁飛區(qū)對可達區(qū)域影響的不確定性，另一方面，禁飛區(qū)的引入會對可達區(qū)域的形狀產生較大不規(guī)則改變，難以通過一個統(tǒng)一的方式或算法對影響后的可達區(qū)域進行描述。以上難點導致在現有方法中往往從制導層面完成禁飛區(qū)規(guī)避，在靠近禁飛區(qū)時添加額外的規(guī)避制導律，這局限了空天飛機的在線自主決策能力。

本文提出一種禁飛區(qū)影響下的再入可達區(qū)域快速分類與求解方法，直接在評估點基于飛行狀態(tài)和禁飛區(qū)信息求解出禁飛區(qū)影響下的可達區(qū)域，進而通過射線法判斷備選目標點是否位于可達區(qū)域以內，對位于可達區(qū)域以內的可行目標點，應用分段預測—校正制導方法完成高精度導引。

1 無禁飛區(qū)的可達區(qū)域求解

1.1 再入動力學模型

空天飛機的再入過程如圖1所示，從再入點開始，到能量管理段(Terminal area energy management，TAEM)交接點為止。

圖1 空天飛機再入過程Fig.1 Reentry phase of aerospace plane

空天飛機再入動力學模型在彈道坐標系中建立，采用傾斜轉彎模式(Bank To Turn, BTT),整個再入過程均為無動力狀態(tài)，在考慮球形大地和地球自轉的基礎上，三維質點動力學模型為

(1)

式中：V為速度；γ為彈道傾角；ψ為彈道偏角；r為地心距，表示飛行器到地心的距離；λ為飛行器在地表投影點的經度；φ為飛行器在地表投影點的緯度；g為重力加速度,其中g0=9.806 7 m/s2；m為飛行器的質量；ωe為地球自轉的角速度；σ為傾側角；D為阻力；L為升力，具體計算方式為

(2)

式中：ρ為大氣密度，可以視為高度的函數；A為參考面積；CL為升力系數；CD為阻力系數。再入過程的控制變量一般只有迎角α和傾側角σ，在設計再入彈道的過程中，迎角α的值由事先設定好的α-V剖面給出：

(3)

典型的再入過程約束由式(4)～式(6)給出：

(4)

(5)

(6)

再入終端滿足地心矩(高度)、速度約束：

(7)

高度和速度又可以用近似能量的方式來進行統(tǒng)一表征：

(8)

所以式(7)表示的終端狀態(tài)約束還可以采用能量形式表達為

ef=eTAEM

(9)

式中:ef為終端能量；eTAEM為相應的能量約束。

圖2 約束轉化得到的傾側角幅值邊界Fig.2 Boundary of bank angle translated from constraints

傾側角幅值滿足約束：

|σ|≤|σ|max

(10)

基于對再入過程的約束與定義，本文后續(xù)的可達區(qū)域計算的終點為TAEM交接點而非降落點。

1.2 不考慮禁飛區(qū)的可達區(qū)域求解

常用的可達區(qū)域求解方法主要包括常值傾側角法和優(yōu)化法。由于可達區(qū)域求解本質是基于某個特定狀態(tài)點對飛行器的飛行能力進行評估，因此，后續(xù)稱該初始狀態(tài)點為評估點。

當不考慮禁飛區(qū)約束時，常值傾側角法是求解可達區(qū)域最簡單、快速、有效的方法。常值傾側角法以任一評估點出發(fā)，采用恒定的傾側角值作為指導指令，積分生成軌跡直至滿足終端約束。通過在[-|σmax|,|σmax|]之間以適當的間隔選取一系列指令值，可以對應生成一系列軌跡，在經度—緯度剖面內將這一系列軌跡的終點連接起來，即構成飛行器在當前初始狀態(tài)下的可達區(qū)域。利用常值傾側角法求解可達區(qū)域的示意圖如圖3所示。

圖3 常值傾側角法求解可達區(qū)域Fig.3 Reachable domain obtained by constant bank angle method

常值傾側角法的算法本質只有積分過程，因而算法穩(wěn)定性好，計算速度較快，但受限于制導指令的單一性，求解出的可達區(qū)域并不完全精確，與實際可達區(qū)域相比偏小。

優(yōu)化法通過多次更改目標函數生成特定軌跡來求解可達區(qū)域。首先，以任一評估點為優(yōu)化的狀態(tài)初值,設置目標函數為航程最大和航程最小，求解出對應的軌跡并記錄最大、最小航程的值Lmax和Lmin；之后，設置優(yōu)化指標為橫向航程最大，并將航程固定為Lmax和Lmin的中間量：Li=Lmin+wi(Lmax-Lmin)，對wi在[0,1]之間進行一系列的取值，得到對應的一系列優(yōu)化軌跡。即求解系列優(yōu)化問題：

(11)

圖4 優(yōu)化方法求解可達區(qū)域Fig.4 Reachable domain obtained by optimization method

優(yōu)化方法得到的邊界雖然相對更加準確，但由于其底層依賴于復雜的非線性規(guī)劃求解器，實際應用中，一旦引入禁飛區(qū)約束，難以保證算法在求解各條軌跡時均能夠快速收斂；同時，運算時間長是一般優(yōu)化算法的通病，這也束縛了它的在線應用前景。綜上所述，對于禁飛區(qū)影響下的在線可達區(qū)域求解，單純地使用上述2種方法均很難達到目的。

2 禁飛區(qū)影響下的可達區(qū)域求解

2.1 禁飛區(qū)對可達區(qū)域的影響特性分析

對于空天飛機這類具有較明顯升力面結構的飛行器，其再入的可達區(qū)域是一個缺角的近似橢圓，形狀如圖5所示。而禁飛區(qū)(圖中陰影部分N處)的引入，實質上是在原有的可達區(qū)域中引入了一片不可達區(qū)域，但這片不可達區(qū)域的形狀是不確定的。

圖5 禁飛區(qū)對可達區(qū)域的影響Fig.5 Influence of no-fly zone on reachable domain

由此，Rd1 d2f在禁飛區(qū)的影響下，初始可達區(qū)域Rabc中有兩片區(qū)域變得不可達，分別是子區(qū)域右側的RcQ1c1eQ2c2c和子區(qū)域左側的Rd1 d2f。

針對圖5中的諸元，存在一些關系需要補充說明。因此作出如下假設：

假設1示意圖中的可達區(qū)域是“充分必要”的，它能完整地反映飛行器飛行能力，也即：任何區(qū)域內的點均是可達的，任何區(qū)域外的點均是不可達的。區(qū)域內的所有點(λt，φt)滿足式(12):

? (λt，φt)∈Rabc

?σ(t)

(12)

基于以上兩點假設，圖5中的諸元滿足推論1和推論2。

推論1Ra1b1c1?Rabc,Ra2b2c2?Rabc

由禁飛區(qū)帶來的影響，只會導致原先可到達的區(qū)域，現在不可到達了，而不會使得原先不可到達的區(qū)域變得可到達。因此，從T1和T22個子評估點出發(fā)得到的子區(qū)域Ra1b1c1和Ra2b2c2均應位于原可達區(qū)域abc以內。

推論2Ra1b1c1和Ra2b2c2同Rabc相切，切點是Q1和Q2

圖6 簡化的“充分”可達區(qū)域Fig.6 Simplified "sufficient" reachable domain

1) 軌跡從禁飛區(qū)下(上)側經過，與禁飛區(qū)相切。

2) 軌跡終點是該側子區(qū)域的上(下)邊界點。

稱這樣的軌跡為極限繞飛軌跡。極限繞飛軌跡在切點前的部分與相切軌跡重合，在切合點后的部分是以最大傾側角幅值飛行產生的軌跡。

由極限繞飛軌跡與該側的相切子區(qū)域圍成的可達區(qū)域在圖中用色塊標記。由于極限繞飛軌跡實際并未達到飛行器的飛行能力極限，同時該軌跡下側的區(qū)域均在飛行能力覆蓋范圍內，因此得到的可達區(qū)域是“充分”的，即可達區(qū)域以內的點均是可以到達的。

對于子區(qū)域數量可能隨著禁飛區(qū)位置的變化發(fā)生改變的問題，在求解子區(qū)域之前，引入禁飛區(qū)分類策略，針對存在0個、1個、2個子區(qū)域的情況，實施分類求解。

禁飛區(qū)影響下的可達區(qū)域求解算法流程如圖7所示。

圖7 禁飛區(qū)影響下的可達區(qū)域求解算法流程Fig.7 Flow chart of algorithm for obtaining the influenced reachable domain

2.2 禁飛區(qū)對可達區(qū)域的影響特性分類

在2.1節(jié)中，給出了一般情形下禁飛區(qū)對可達區(qū)域的影響特性，在實際情況中，受禁飛區(qū)與飛行軌跡相對位置的影響，對于2.1節(jié)圖5給出的諸元，可能并不總是全部存在，比如有時飛行器可能只能從禁飛區(qū)的某一側通過，則對應的子區(qū)域也只有一個。為此，需要對這些情況進行分類討論。圖8給出了不同類的禁飛區(qū)位置示意圖。

圖8 不同類的禁飛區(qū)位置示意圖Fig.8 Diagram of different classes of no-fly zones

需要指出，本節(jié)和后續(xù)2.3節(jié)實現的是實際可能面對的多種復雜情形下圖5中子區(qū)域Ra1b1c1和Ra2b2c2的求解。

對禁飛區(qū)進行分類的標準遵循禁飛區(qū)與3條特殊軌跡的相對位置。Lmax為最大航程軌跡，+Lmin和-Lmin分別為下側和上側的最小航程軌跡。

對于經度—緯度剖面內的任一條軌跡，可用矩陣[λi,φi]表示：

(13)

定義軌跡距離S為軌跡到禁飛區(qū)圓心的最短距離，由式(14)給出：

(14)

式中:λn為禁飛區(qū)圓心的經度；φn為禁飛區(qū)圓心的緯度。軌跡同禁飛區(qū)的關系標識由式(15)給出：

(15)

式中:1為相交；0為相離；Rn為禁飛區(qū)的半徑。由位置關系(相交或相離) 可將禁飛區(qū)分為4個大類，8個小類。如表1所示。

表1 基于同Lmax,+Lmin&-Lmin相對位置關系的禁飛區(qū)分類

對于禁飛區(qū)同3條特殊軌跡均相離的情況，仍然存在第0類，1-B類和1-C類3種情形，可以按照式(16)原則進一步劃分：

(16)

式中:N表示禁飛區(qū)；[λmax，φmax]為Lmax對應的軌跡點經緯度；φmax([λmax,φmax]|λn)為Lmax在λ=λn處的插值緯度。

對第一大類，禁飛區(qū)影響下的可達區(qū)域存在2個子區(qū)域，即對應2.1節(jié)圖5中的Ra1b1c1和Ra2b2c2；第二大類和第三大類，禁飛區(qū)影響下的可達區(qū)域僅存在一個子區(qū)域，對于第零大類和第四大類，禁飛區(qū)影響下的可達區(qū)域不存在子區(qū)域。更多細節(jié)在表2中展示。

表2 不同分類的禁飛區(qū)規(guī)避策略

2.3 相切軌跡搜索與子區(qū)域生成

綜合考量效率和穩(wěn)定性，本節(jié)基于常值傾側角對應的一系列軌跡，采用二分法搜索相切軌跡，再從相切點(子評估點)T1，T2出發(fā)計算子區(qū)域。

由常值傾側角法計算軌跡落點的方法如式(17)所示：

(17)

式中:σc的取值為

(18)

依據待搜索軌跡所在的大致區(qū)域，算法由兩對共4個二分法搜索環(huán)節(jié)構成，4個環(huán)節(jié)分別稱Up1，Up2，Down1和Down2。以0°傾側角對應軌跡Lmax為邊界，Up區(qū)為上半區(qū)，Down區(qū)為下半區(qū)；1代表禁飛區(qū)左側，2代表禁飛區(qū)右側。圖9給出了搜索區(qū)域以及搜索的目標軌跡。

圖9 搜索區(qū)域示意圖Fig.9 Diagram of search area

(19)

(20)

相切軌跡求解流程見①～③

① 在單個搜索區(qū)內找到一條與禁飛區(qū)相交的軌跡[λstart，φstart]，軌跡對應的傾側角值為σstart，作為初始搜索區(qū)間的邊界。其中，Up1和Down1初始搜索區(qū)間在[λstart，φstart]左側，Up2和Down2初始搜索區(qū)間在[λstart，φstart]右側。

② 定義軌跡距離S為目標函數，對任意的傾側角值σu，有其對應軌跡[λu，φu]及最短距離值S(σu)

(21)

(22)

4個環(huán)節(jié)的主要差異在初始搜索區(qū)間，如表3所示。

表3 不同環(huán)節(jié)搜索區(qū)域和初始區(qū)間

依據分類結果，需要分別采用不同的環(huán)節(jié)，在不同區(qū)間搜索相切軌跡，具體如表4所示。

表4 各分類應用環(huán)節(jié)

需要指出的是，對于部分升阻比偏高的飛行器，其地面軌跡在末段可能出現繞回現象，即在低速、大傾側時，彈道偏角在短時間內就會發(fā)生較大改變。本文給出的方法適用于不發(fā)生繞回現象的情況。

3 可達性判定與制導

3.1 可達性判定方法

在求解得到可達區(qū)域以后，若想要辨明某個目標點是否位于可達區(qū)域內，還需要設計在線判定算法。

對于不考慮禁飛區(qū)的情形，可達區(qū)域的邊界相對簡單，即第2章中所述區(qū)域Rabc。在仿真與實際應用中，根據2.3節(jié)提到的可達區(qū)域求解方法，Rabc并非是由解析方程確定的橢圓，而是通過若干個邊界點的連接來確定的凹多邊形。對于平面內任意多邊形區(qū)域，均可以采用射線法判斷來點是否位于區(qū)域內，如圖10所示。

圖10 射線法示意圖Fig.10 Diagram of ray method

所謂射線法，是從被測點向任意方向作一條射線，判斷射線與多邊形邊界的交點數量。如果交點的數量為奇數，則被測點在多邊形內；如果交點的數量為偶數，則被測點在多邊形以外。應用該方法即可判定是否有目標點位于可達區(qū)域內。

3.2 考慮禁飛區(qū)規(guī)避的分段制導方法

預測校正方法是近年來熱門的在線再入制導方法，許多學者進行了相關研究[24-25]。本文采用一種分段預測—校正制導方法，具體邏輯詳見文獻[26]。

該方法與傳統(tǒng)預測—校正制導的主要區(qū)別有兩點：一是在縱向制導環(huán)引入分段目標函數，如下所示：

(23)

式中：Sp為當前點到預測落點的距離；S0為當前點到目標點的距離；Sfp為目標點到預測落點的距離。各項參數的具體計算方法也可參見文獻[26]。

二是引入了禁飛區(qū)規(guī)避邏輯，當飛行器抵近禁飛區(qū)時，采取傾側角反轉或增大傾側角幅值的策略來完成規(guī)避。

4 仿真校驗

4.1 禁飛區(qū)影響下的可達區(qū)域分類仿真

綜合考量效率和穩(wěn)定性，本文采用常值傾側角法在評估點R和子評估點T1，T2計算可達區(qū)域和子區(qū)域，需要注意，由于常值傾側角法求解可達區(qū)域與實際的可達區(qū)域相比相對偏小，仿真中實際計算得到的區(qū)域無法嚴格滿足推論1和推論2，最終所得影響下可達區(qū)域與實際影響下可達區(qū)域相比也相對偏小，即得到禁飛區(qū)影響下的可達區(qū)域是“充分”的結果。

以空天飛機的再入階段為例進行仿真分析，仿真對象為某具有升力結構的空天飛機，再入點的初始狀態(tài)由表5給出。

表5 縱向初始條件

首先，沿最大航程軌跡Lmax布置大量禁飛區(qū)，禁飛區(qū)的半徑均為200 km。布置結果如圖11所示，圖中，實線邊界為不考慮禁飛區(qū)的情形下得到的初始可達區(qū)域(即Rabc)，虛線為最大航程軌跡Lmax，即0°傾側角對應軌跡。

采用本文所述算法對圖11所有禁飛區(qū)進行分類，分類結果如圖12所示。

圖11 禁飛區(qū)分布示意圖Fig.11 Distribution diagram of no-fly zones

圖12 禁飛區(qū)分類結果Fig.12 Classification results of no-fly zones

仿真結果表明，算法能夠穩(wěn)定高效地完成對經度/緯度剖面內的禁飛區(qū)的分類，用以支撐后續(xù)可達區(qū)域求解過程。

4.2 禁飛區(qū)影響下的可達區(qū)域求解仿真

對于不同分類的禁飛區(qū)，需要采用不同的搜索環(huán)節(jié)搜索相切軌跡和進行子區(qū)域求解。在圖12中任意選擇每一類的禁飛區(qū)進行可達區(qū)域求解，選取的禁飛區(qū)圓心位置在表6中給出。最終得到禁飛區(qū)影響下的可達區(qū)域如圖13～圖15所示。

圖13 情形1-A,1-B,1-C仿真實例Fig.13 Cases 1-A,1-B and 1-C simulation example

圖14 情形2-A,2-B仿真實例Fig.14 Cases 2-A and 2-B simulation example

圖15 情形3-A,3-B仿真實例Fig.15 Cases 3-A and 3-B simulation example

表6 各類禁飛區(qū)圓心位置

圖中，大橢圓邊界為初始可達區(qū)域(即圖5中的Rabc)，小橢圓邊界為子區(qū)域(即圖5中的Ra1b1c1和Ra2b2c2)，圓形為禁飛區(qū)，禁飛區(qū)影響下的可達區(qū)域在圖中用色塊填充。

仿真結果表明，算法對不同位置的禁飛區(qū)實現了有效分類，并基于分類結果求解得到了相切軌跡、極限繞飛軌跡和子區(qū)域，最終獲得禁飛區(qū)影響下的可達區(qū)域，對于經度—緯度剖面內任意位置、任意分類的圓形禁飛區(qū)，在不超過飛行器物理能力的前提下，算法均能夠穩(wěn)定地實現可達區(qū)域求解。

結合本文仿真對象實際情況，可以給出本文所述方法的一個初步適用范圍，對升阻比在2以下的飛行器，取末端能量條件HTAEM=21.9 km，VTAEM=764 m/s時，能夠應用上述方法實現飛行能力評估。

5.3 面向可達區(qū)域內目標點的規(guī)避制導仿真

以圖13(a)情形1-A的禁飛區(qū)為基礎進行本節(jié)仿真，初始條件與表5中完全一致。在圖13(a)所示的區(qū)域中選擇了兩個位于可達區(qū)域內的目標點：

(λf1，φf1)=(90°，45°)

(λf2，φf2)=(80°，35°)

仿真結果如圖16和圖17所示。從圖中可以發(fā)現，對于選取的2個目標點，算法可以有效實現將飛行器向目標點導引，并且飛行軌跡不經過禁飛區(qū)。

圖16 面向目標點1導引的飛行軌跡Fig.16 Trajectory of guidance to Target 1

圖17 面向目標點2導引的飛行軌跡Fig.17 Trajectory of guidance to Target 2

表7數據顯示，相比于單段目標函數預測—校正制導方法[23]，引入分段目標函數后的制導精度有較明顯提升。

表7 落點誤差對比

5 結 論

1) 分析了一般情形下禁飛區(qū)對可達區(qū)域的影響特性；結合實際飛行中可能存在的禁飛區(qū)與飛行器的相對位置情況，提出借助最大航程軌跡Lmax、上側最小航程軌跡-Lmin和下側最小航程軌跡+Lmin，對經度/緯度剖面內禁飛區(qū)存在位置展開分類。

2) 依據分類結果，分別求解相切軌跡、極限繞飛軌跡和子區(qū)域；通過子區(qū)域和極限繞飛軌跡確定可達區(qū)域與不可達區(qū)域的邊界線，得到禁飛區(qū)影響下的可達區(qū)域。

3) 介紹了判定目標點是否位于可達區(qū)域以內的射線法和可以應用于可達區(qū)域以內目標點的分段預測校正制導方法。

4) 通過數值仿真驗證了提出算法的有效性和穩(wěn)定性。對于經度/緯度剖面內散布的禁飛區(qū)，算法均能實現可達區(qū)域的求解。對確定位于可達區(qū)域以內的目標點，通過應用分段預測校正方法，能夠實現導引并滿足終端精度要求。