章吉力,周大鵬,楊大鵬,劉然,劉凱,*
1. 大連理工大學 航空航天學院,大連 116024
2. 航空工業沈陽飛機設計研究所,沈陽 110035
空天飛機是一種能夠快速往返于稠密大氣層、臨近空間和軌道空間的可重復使用飛行器[1-2],既可以像普通飛機一樣水平起降,又具備臨近空間持續機動飛行能力,具有重要的政治價值、經濟價值和軍事價值。
可達區域求解[1-2]是對空天飛機的覆蓋范圍和飛行能力進行評估的重要方法,基于可達區域求解的結果進行任務的繼續與中止、目標的決策與切換能夠大幅提升飛行器自主性。可達區域的求解方法已受到國內外眾多學者的關注。
可達區域的求解方法有很多種,主要有常值傾側角法[3-4]、剖面平移設計法[5]、優化法[6-14]、剖面參數規劃法[15-17]等。由于常值傾側角法解算出的區域存在一定誤差,與實際區域相比偏小,但同時也具有簡單易行、穩定性高、解算速度快的特性,因而成為工程中常用的方法。剖面平移設計法首先將再入過程約束轉化到H-V剖面,形成再入走廊,之后沿著走廊下邊界進行再入剖面規劃,形成航程最短的極限剖面,再將該剖面向上平移,形成航程適中的剖面,最后以零傾側角飛行,形成航程最長的剖面[5],將這一系列剖面的終點連接起來,即得到可達區域。
優化法是現有最為精確的解算方法,但解算速度相對較慢。文獻[7-8]基于目標函數變更求解構成可達區域邊界的軌跡。優化指標涉及最大縱程、最小縱程,及若干中間量縱程填補邊界曲線,最后把這些彈道的落點連起來,即得到可達區域。文獻[13]還考慮到優化得到的控制指令的平滑性,認為直接使用偽譜法優化傾側角和攻角指令得到的剖面可行性差,使用控制變量的二階導數作為優化變量能夠得到更平滑的控制指令。剖面參數規劃法介于優化法和剖面平移設計法之間,首先確定一個基準剖面,剖面上有許多待求參數,通過設置與優化法類似的特殊優化指標,來確定這些待求參數,文獻[15]采用粒子群優化的方式來求解這些待定參數,文獻[16]利用差分進化算法和傾側角插值求解待定參數,權衡兼顧了可達區域精度和求解效率。
還有學者通過對可達區域特性的分析,以解析和近似的方式來求解可達區域,文獻[18]考慮到可達區域同橢圓的近似性,以兩個最大橫程終點間線段為短軸,以最大縱程終點到短軸距離為半長軸構建半橢圓可達區域。文獻[19]通過解析同倫法,在求解4個子優化問題的基礎上,通過構造合適的同倫參數延拓出縱向可達區。也有學者將可達區域求解應用于某些特殊場景。文獻[20]考慮了在具有路徑點約束情況下的可達區域估計方法。文獻[21]通過可達區域求解結果指導在軌星座設計方法,構建目標全覆蓋的空間戰略武器星座。
總的來說,現有的研究中對可達區域的求解方法鮮有考慮禁飛區的影響。主要存在以下難點:一方面,禁飛區與飛行器相對位置的不確定性引發禁飛區對可達區域影響的不確定性,另一方面,禁飛區的引入會對可達區域的形狀產生較大不規則改變,難以通過一個統一的方式或算法對影響后的可達區域進行描述。以上難點導致在現有方法中往往從制導層面完成禁飛區規避,在靠近禁飛區時添加額外的規避制導律,這局限了空天飛機的在線自主決策能力。
本文提出一種禁飛區影響下的再入可達區域快速分類與求解方法,直接在評估點基于飛行狀態和禁飛區信息求解出禁飛區影響下的可達區域,進而通過射線法判斷備選目標點是否位于可達區域以內,對位于可達區域以內的可行目標點,應用分段預測—校正制導方法完成高精度導引。
空天飛機的再入過程如圖1所示,從再入點開始,到能量管理段(Terminal area energy management,TAEM)交接點為止。

圖1 空天飛機再入過程Fig.1 Reentry phase of aerospace plane
空天飛機再入動力學模型在彈道坐標系中建立,采用傾斜轉彎模式(Bank To Turn, BTT),整個再入過程均為無動力狀態,在考慮球形大地和地球自轉的基礎上,三維質點動力學模型為
(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)
高度和速度又可以用近似能量的方式來進行統一表征:
(8)
所以式(7)表示的終端狀態約束還可以采用能量形式表達為
ef=eTAEM
(9)
式中:ef為終端能量;eTAEM為相應的能量約束。


圖2 約束轉化得到的傾側角幅值邊界Fig.2 Boundary of bank angle translated from constraints
傾側角幅值滿足約束:
|σ|≤|σ|max
(10)
基于對再入過程的約束與定義,本文后續的可達區域計算的終點為TAEM交接點而非降落點。
常用的可達區域求解方法主要包括常值傾側角法和優化法。由于可達區域求解本質是基于某個特定狀態點對飛行器的飛行能力進行評估,因此,后續稱該初始狀態點為評估點。
當不考慮禁飛區約束時,常值傾側角法是求解可達區域最簡單、快速、有效的方法。常值傾側角法以任一評估點出發,采用恒定的傾側角值作為指導指令,積分生成軌跡直至滿足終端約束。通過在[-|σmax|,|σmax|]之間以適當的間隔選取一系列指令值,可以對應生成一系列軌跡,在經度—緯度剖面內將這一系列軌跡的終點連接起來,即構成飛行器在當前初始狀態下的可達區域。利用常值傾側角法求解可達區域的示意圖如圖3所示。

圖3 常值傾側角法求解可達區域Fig.3 Reachable domain obtained by constant bank angle method
常值傾側角法的算法本質只有積分過程,因而算法穩定性好,計算速度較快,但受限于制導指令的單一性,求解出的可達區域并不完全精確,與實際可達區域相比偏小。
優化法通過多次更改目標函數生成特定軌跡來求解可達區域。首先,以任一評估點為優化的狀態初值,設置目標函數為航程最大和航程最小,求解出對應的軌跡并記錄最大、最小航程的值Lmax和Lmin;之后,設置優化指標為橫向航程最大,并將航程固定為Lmax和Lmin的中間量:Li=Lmin+wi(Lmax-Lmin),對wi在[0,1]之間進行一系列的取值,得到對應的一系列優化軌跡。即求解系列優化問題:
(11)


圖4 優化方法求解可達區域Fig.4 Reachable domain obtained by optimization method
優化方法得到的邊界雖然相對更加準確,但由于其底層依賴于復雜的非線性規劃求解器,實際應用中,一旦引入禁飛區約束,難以保證算法在求解各條軌跡時均能夠快速收斂;同時,運算時間長是一般優化算法的通病,這也束縛了它的在線應用前景。綜上所述,對于禁飛區影響下的在線可達區域求解,單純地使用上述2種方法均很難達到目的。
對于空天飛機這類具有較明顯升力面結構的飛行器,其再入的可達區域是一個缺角的近似橢圓,形狀如圖5所示。而禁飛區(圖中陰影部分N處)的引入,實質上是在原有的可達區域中引入了一片不可達區域,但這片不可達區域的形狀是不確定的。


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

由此,Rd1 d2f在禁飛區的影響下,初始可達區域Rabc中有兩片區域變得不可達,分別是子區域右側的RcQ1c1eQ2c2c和子區域左側的Rd1 d2f。
針對圖5中的諸元,存在一些關系需要補充說明。因此作出如下假設:
假設1示意圖中的可達區域是“充分必要”的,它能完整地反映飛行器飛行能力,也即:任何區域內的點均是可達的,任何區域外的點均是不可達的。區域內的所有點(λt,φt)滿足式(12):
? (λt,φt)∈Rabc
?σ(t)
(12)


基于以上兩點假設,圖5中的諸元滿足推論1和推論2。
推論1Ra1b1c1?Rabc,Ra2b2c2?Rabc
由禁飛區帶來的影響,只會導致原先可到達的區域,現在不可到達了,而不會使得原先不可到達的區域變得可到達。因此,從T1和T22個子評估點出發得到的子區域Ra1b1c1和Ra2b2c2均應位于原可達區域abc以內。
推論2Ra1b1c1和Ra2b2c2同Rabc相切,切點是Q1和Q2




圖6 簡化的“充分”可達區域Fig.6 Simplified "sufficient" reachable domain
1) 軌跡從禁飛區下(上)側經過,與禁飛區相切。
2) 軌跡終點是該側子區域的上(下)邊界點。
稱這樣的軌跡為極限繞飛軌跡。極限繞飛軌跡在切點前的部分與相切軌跡重合,在切合點后的部分是以最大傾側角幅值飛行產生的軌跡。
由極限繞飛軌跡與該側的相切子區域圍成的可達區域在圖中用色塊標記。由于極限繞飛軌跡實際并未達到飛行器的飛行能力極限,同時該軌跡下側的區域均在飛行能力覆蓋范圍內,因此得到的可達區域是“充分”的,即可達區域以內的點均是可以到達的。
對于子區域數量可能隨著禁飛區位置的變化發生改變的問題,在求解子區域之前,引入禁飛區分類策略,針對存在0個、1個、2個子區域的情況,實施分類求解。
禁飛區影響下的可達區域求解算法流程如圖7所示。

圖7 禁飛區影響下的可達區域求解算法流程Fig.7 Flow chart of algorithm for obtaining the influenced reachable domain
在2.1節中,給出了一般情形下禁飛區對可達區域的影響特性,在實際情況中,受禁飛區與飛行軌跡相對位置的影響,對于2.1節圖5給出的諸元,可能并不總是全部存在,比如有時飛行器可能只能從禁飛區的某一側通過,則對應的子區域也只有一個。為此,需要對這些情況進行分類討論。圖8給出了不同類的禁飛區位置示意圖。

圖8 不同類的禁飛區位置示意圖Fig.8 Diagram of different classes of no-fly zones
需要指出,本節和后續2.3節實現的是實際可能面對的多種復雜情形下圖5中子區域Ra1b1c1和Ra2b2c2的求解。
對禁飛區進行分類的標準遵循禁飛區與3條特殊軌跡的相對位置。Lmax為最大航程軌跡,+Lmin和-Lmin分別為下側和上側的最小航程軌跡。
對于經度—緯度剖面內的任一條軌跡,可用矩陣[λi,φi]表示:
(13)
定義軌跡距離S為軌跡到禁飛區圓心的最短距離,由式(14)給出:
(14)
式中:λn為禁飛區圓心的經度;φn為禁飛區圓心的緯度。軌跡同禁飛區的關系標識由式(15)給出:
(15)
式中:1為相交;0為相離;Rn為禁飛區的半徑。由位置關系(相交或相離) 可將禁飛區分為4個大類,8個小類。如表1所示。

表1 基于同Lmax,+Lmin&-Lmin相對位置關系的禁飛區分類
對于禁飛區同3條特殊軌跡均相離的情況,仍然存在第0類,1-B類和1-C類3種情形,可以按照式(16)原則進一步劃分:
(16)
式中:N表示禁飛區;[λmax,φmax]為Lmax對應的軌跡點經緯度;φmax([λmax,φmax]|λn)為Lmax在λ=λn處的插值緯度。
對第一大類,禁飛區影響下的可達區域存在2個子區域,即對應2.1節圖5中的Ra1b1c1和Ra2b2c2;第二大類和第三大類,禁飛區影響下的可達區域僅存在一個子區域,對于第零大類和第四大類,禁飛區影響下的可達區域不存在子區域。更多細節在表2中展示。

表2 不同分類的禁飛區規避策略
綜合考量效率和穩定性,本節基于常值傾側角對應的一系列軌跡,采用二分法搜索相切軌跡,再從相切點(子評估點)T1,T2出發計算子區域。
由常值傾側角法計算軌跡落點的方法如式(17)所示:
(17)
式中:σc的取值為
(18)
依據待搜索軌跡所在的大致區域,算法由兩對共4個二分法搜索環節構成,4個環節分別稱Up1,Up2,Down1和Down2。以0°傾側角對應軌跡Lmax為邊界,Up區為上半區,Down區為下半區;1代表禁飛區左側,2代表禁飛區右側。圖9給出了搜索區域以及搜索的目標軌跡。

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

(20)
相切軌跡求解流程見①~③
① 在單個搜索區內找到一條與禁飛區相交的軌跡[λstart,φstart],軌跡對應的傾側角值為σstart,作為初始搜索區間的邊界。其中,Up1和Down1初始搜索區間在[λstart,φstart]左側,Up2和Down2初始搜索區間在[λstart,φstart]右側。
② 定義軌跡距離S為目標函數,對任意的傾側角值σu,有其對應軌跡[λu,φu]及最短距離值S(σu)
(21)
(22)
4個環節的主要差異在初始搜索區間,如表3所示。

表3 不同環節搜索區域和初始區間
依據分類結果,需要分別采用不同的環節,在不同區間搜索相切軌跡,具體如表4所示。

表4 各分類應用環節

需要指出的是,對于部分升阻比偏高的飛行器,其地面軌跡在末段可能出現繞回現象,即在低速、大傾側時,彈道偏角在短時間內就會發生較大改變。本文給出的方法適用于不發生繞回現象的情況。
在求解得到可達區域以后,若想要辨明某個目標點是否位于可達區域內,還需要設計在線判定算法。
對于不考慮禁飛區的情形,可達區域的邊界相對簡單,即第2章中所述區域Rabc。在仿真與實際應用中,根據2.3節提到的可達區域求解方法,Rabc并非是由解析方程確定的橢圓,而是通過若干個邊界點的連接來確定的凹多邊形。對于平面內任意多邊形區域,均可以采用射線法判斷來點是否位于區域內,如圖10所示。

圖10 射線法示意圖Fig.10 Diagram of ray method
所謂射線法,是從被測點向任意方向作一條射線,判斷射線與多邊形邊界的交點數量。如果交點的數量為奇數,則被測點在多邊形內;如果交點的數量為偶數,則被測點在多邊形以外。應用該方法即可判定是否有目標點位于可達區域內。
預測校正方法是近年來熱門的在線再入制導方法,許多學者進行了相關研究[24-25]。本文采用一種分段預測—校正制導方法,具體邏輯詳見文獻[26]。
該方法與傳統預測—校正制導的主要區別有兩點:一是在縱向制導環引入分段目標函數,如下所示:
(23)
式中:Sp為當前點到預測落點的距離;S0為當前點到目標點的距離;Sfp為目標點到預測落點的距離。各項參數的具體計算方法也可參見文獻[26]。
二是引入了禁飛區規避邏輯,當飛行器抵近禁飛區時,采取傾側角反轉或增大傾側角幅值的策略來完成規避。
綜合考量效率和穩定性,本文采用常值傾側角法在評估點R和子評估點T1,T2計算可達區域和子區域,需要注意,由于常值傾側角法求解可達區域與實際的可達區域相比相對偏小,仿真中實際計算得到的區域無法嚴格滿足推論1和推論2,最終所得影響下可達區域與實際影響下可達區域相比也相對偏小,即得到禁飛區影響下的可達區域是“充分”的結果。
以空天飛機的再入階段為例進行仿真分析,仿真對象為某具有升力結構的空天飛機,再入點的初始狀態由表5給出。

表5 縱向初始條件
首先,沿最大航程軌跡Lmax布置大量禁飛區,禁飛區的半徑均為200 km。布置結果如圖11所示,圖中,實線邊界為不考慮禁飛區的情形下得到的初始可達區域(即Rabc),虛線為最大航程軌跡Lmax,即0°傾側角對應軌跡。
采用本文所述算法對圖11所有禁飛區進行分類,分類結果如圖12所示。

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

圖12 禁飛區分類結果Fig.12 Classification results of no-fly zones
仿真結果表明,算法能夠穩定高效地完成對經度/緯度剖面內的禁飛區的分類,用以支撐后續可達區域求解過程。
對于不同分類的禁飛區,需要采用不同的搜索環節搜索相切軌跡和進行子區域求解。在圖12中任意選擇每一類的禁飛區進行可達區域求解,選取的禁飛區圓心位置在表6中給出。最終得到禁飛區影響下的可達區域如圖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 各類禁飛區圓心位置
圖中,大橢圓邊界為初始可達區域(即圖5中的Rabc),小橢圓邊界為子區域(即圖5中的Ra1b1c1和Ra2b2c2),圓形為禁飛區,禁飛區影響下的可達區域在圖中用色塊填充。
仿真結果表明,算法對不同位置的禁飛區實現了有效分類,并基于分類結果求解得到了相切軌跡、極限繞飛軌跡和子區域,最終獲得禁飛區影響下的可達區域,對于經度—緯度剖面內任意位置、任意分類的圓形禁飛區,在不超過飛行器物理能力的前提下,算法均能夠穩定地實現可達區域求解。
結合本文仿真對象實際情況,可以給出本文所述方法的一個初步適用范圍,對升阻比在2以下的飛行器,取末端能量條件HTAEM=21.9 km,VTAEM=764 m/s時,能夠應用上述方法實現飛行能力評估。
以圖13(a)情形1-A的禁飛區為基礎進行本節仿真,初始條件與表5中完全一致。在圖13(a)所示的區域中選擇了兩個位于可達區域內的目標點:
(λf1,φf1)=(90°,45°)
(λf2,φf2)=(80°,35°)
仿真結果如圖16和圖17所示。從圖中可以發現,對于選取的2個目標點,算法可以有效實現將飛行器向目標點導引,并且飛行軌跡不經過禁飛區。

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

圖17 面向目標點2導引的飛行軌跡Fig.17 Trajectory of guidance to Target 2
表7數據顯示,相比于單段目標函數預測—校正制導方法[23],引入分段目標函數后的制導精度有較明顯提升。

表7 落點誤差對比
1) 分析了一般情形下禁飛區對可達區域的影響特性;結合實際飛行中可能存在的禁飛區與飛行器的相對位置情況,提出借助最大航程軌跡Lmax、上側最小航程軌跡-Lmin和下側最小航程軌跡+Lmin,對經度/緯度剖面內禁飛區存在位置展開分類。
2) 依據分類結果,分別求解相切軌跡、極限繞飛軌跡和子區域;通過子區域和極限繞飛軌跡確定可達區域與不可達區域的邊界線,得到禁飛區影響下的可達區域。
3) 介紹了判定目標點是否位于可達區域以內的射線法和可以應用于可達區域以內目標點的分段預測校正制導方法。
4) 通過數值仿真驗證了提出算法的有效性和穩定性。對于經度/緯度剖面內散布的禁飛區,算法均能實現可達區域的求解。對確定位于可達區域以內的目標點,通過應用分段預測校正方法,能夠實現導引并滿足終端精度要求。