基于偽譜法的再入可達域影響因素分析

李兆亭， 周 祥， 張洪波， 湯國建

(國防科技大學 空天科學學院，長沙 410073)

再入可達域是指可重復使用飛行器在滿足約束條件下，到達著陸或者交班狀態下所能允許的區域范圍.它是飛行器縱向機動能力與橫向機動能力的綜合反映，體現了飛行器的整體機動性能.可達域求解可為飛行器軌跡規劃與制導、應急著陸點選擇、制導律性能評估提供依據，因此具有重要意義.

再入可達域的求解方法按照原理可以分為兩類：參數化方法和非參數化方法.參數化方法是指將再入軌跡規劃問題的狀態量或控制量等進行參數化，采用某些優化方法對參數化后的問題進行優化求解；非參數化方法則不包含此類參數化、問題優化過程，常見的基于剖面的可達域求解方法就屬于此類.如曾夕娟等[1]在“r-V”空間內引入約束的對數表達式簡化再入走廊的表達，將再入軌跡分為3段設計，引入準平衡滑翔條件獲取解析制導律以跟蹤所規劃軌跡，通過平移“邊界跟蹤段”軌跡使其布滿再入走廊，從而得到最終的可達域.常松濤等[2]在基于E剖面的再入走廊下，引入了傾側角反向策略，改進了阻力加速度曲線插值方式和最大可行阻力加速度曲線的生成方法，使計算得到的可達區域更精確、邊界更均勻、彈道更平滑.He等[3]推導得到一個具有不確定性效應的再入走廊，其中攻角剖面可調節，縱向阻力剖面設計為上、下擬合安全邊界的插值結果，橫向升力-阻力剖面則利用準平衡滑翔條件得到，最后通過比例積分微分(PID)跟蹤器來跟蹤剖面，得到可達域.Zhang等[4]提出一種基于幾何預測軌跡的可達域計算方法，在D-E走廊下，對橫向運動狀態、轉彎半徑和角度進行分析求解，得到D-E剖面的最遠可達邊界，并提出幾何預測的軌跡；通過重復計算所有D-E剖面的可達邊界，生成可達域.上述基于剖面的可達域求解方法通常需要人工設計相關剖面，而剖面參數決定了其可達域范圍，因此得到的可達域通常只是實際可達域的一部分.此外，Liu等[5]提出一種使用多項式來描述可達域邊界的方法，使用反向傳播神經網絡快速確定空對地飛行器的可達域.Zhang等[6]通過快速解析的方法求解3個不同的最優控制問題，得到可達域的內、外邊界.章吉力等[7]對考慮禁飛區條件下的空天飛機再入可達區域問題進行研究，并基于極限繞飛軌跡提出一種不限禁飛區位置的可達區域求解方法.吳楠等[8]基于平衡滑翔假設和最優化飛行準則，通過數值積分獲得最大縱程和橫程彈道，基于橢圓近似法利用3個末端點構建可達域橢圓邊界.

可達域的參數化求解方法包含的種類較多，常見方法主要包括基于虛擬目標點[9-11]、基于偽譜法[12-14]和基于凸優化方法[15]等.如趙吉松等[16]提出一種基于稀疏差分法和網格細化技術的快速、高精度軌跡求解方法，并利用可達域的求解驗證了方法的有效性.藺君等[17]對攻角、傾側角進行參數化，利用帶約束的差分進化算法求解滿足再入過程約束和終端約束的再入軌跡.梁巨平等[18]選取傾側角為時間的分段常值函數，采用遺傳算法來求解一系列的軌跡優化問題.趙江等[19]設計了一種參數化的傾側角剖面，利用約束粒子群優化算法求解滿足再入過程約束和末端約束的最優滑翔軌跡.相比基于剖面的求解方法，可達域的參數化求解方法不需要設計相關剖面，其可達域范圍更大且更接近實際值，但受限于具體算法的性能，其全局最優性仍難以保證.孫勇等[20]基于擬能量概念提出了多段軌跡優化方法，將軌跡按照擬能量進行分段，對各段分別進行優化，減少了每次優化的變量個數，計算效率得到較大提高.

偽譜法具有良好的適用性和強大的求解能力，基于偽譜法求解可達域的方案較多.如王濤等[21]提出一種基于Gauss偽譜法的再入可達域計算方法，采用固定的攻角剖面，僅對傾側角進行單變量尋優.李柯等[22]基于Radau偽譜法求解得到飛行器的可達域，著重分析了升阻比、終端速度以及終端傾角對可達域的影響.Wang等[23]在準平衡滑行條件下，使用高斯偽譜法計算得到了返回航天器的可達域.解永鋒等[24]通過定義加權的橫程、縱程組合性能指標函數，將可達域求解問題轉化為組合性能指標最優的控制問題，采用勒讓德偽譜法快速計算可達域.

本文繼承發展了文獻[21]的求解思路，將攻角、傾側角同時作為控制量代入到優化問題中，采用偽譜法進行軌跡規劃，得到的可達域范圍更廣且更接近實際值.此外，描述了基于偽譜法的可達域快速求解方法，并對偽譜法進行了簡要闡述；然后基于上述方法對影響可達域的相關因素進行仿真研究.仿真結果表明，飛行器質量、氣動參考面積、大氣密度等會對短縱程軌跡產生明顯影響；同時，升阻比大小與可達域范圍成正相關.

1 再入飛行器及其可達域描述

介紹再入飛行器的動力學模型，并在換極坐標系下給出利用偽譜法求解可達域的思路.

1.1 動力學方程及約束

在求解可達域時，通常將飛行器的運動轉換到換極坐標系[25]中.在換極坐標系中，飛行器起點的經緯度均為0°，終點的經度和緯度分別描述了再入縱程和再入橫程.利用這些特性，可方便地標定偽譜法中的參數范圍，簡化彈道規劃算法.

在換極坐標系下，可以推導得到飛行器再入運動方程：

(1)

(2)

式中：ρ為大氣密度；M為飛行器質量；S為參考面積；CL和CD分別為升力系數和阻力系數，均與攻角α有關.同時，再入飛行過程中還需考慮多種過程約束，如熱流密度約束、過載約束、動壓約束等，具體為

(3)

在本文中，目標飛行器為CAV-H，其相關參數和氣動數據可參考文獻[26].

1.2 換極坐標系下的可達域求解

為了加快優化速度和提高優化質量，文獻[21]提出了一種基于Gauss偽譜法的再入可達域計算方法，采用固定的攻角剖面，將攻角作為狀態變量，僅對傾側角進行單變量尋優.本文延續該求解思路，但與之不同的是，攻角剖面并不是固定的，即攻角、傾側角將同時作為待優化變量.

再入可達域問題可表述為求解不同縱程條件下的最大橫程.在換極坐標系內，落點的橫程可用緯度表征:

J=min(cosφpf)

(4)

式中：φpf為換極坐標系中終點的緯度，如圖1所示.

圖1中，Lset為設定的縱程，λpf為換極坐標系中終點的經度，二者數值相等.不妨將Lset視為終端約束，通過改變Lset值得到不同縱程條件下的最大橫程軌跡.同時將換極坐標系下，φpf>0°的所有軌跡稱之為北半球軌跡，反之稱之為南半球軌跡.在進行偽譜法的參數設置時，考慮到地球自轉的影響，即微分方程(1)中的科氏加速度項和牽連加速度項均不為0，不妨將把北半球軌跡的緯度變化范圍設置為 [-2°, 90°]，將南半球軌跡的緯度變化范圍設置為 [-90°, 2°].

圖1 不同縱程下的最大橫程Fig.1 Maximum cross-ranges of different longitudinal ranges

求解上述不同縱程下的最大橫程問題，即可得到有效的可達域.但實際仿真中，同樣需要確定最大、最小縱程(Lmax和Lmin)，以便確定Lset的取值范圍.因此，最終的規劃流程如圖2所示.首先求解最大、最小縱程，標定Lmax和Lmin的具體數值；然后從縱程Lmax-ΔL開始，分別求解兩個最大橫程問題，依次減小縱程ΔL，直到達到Lmin；最后即可得到飛行器的再入可達域.

圖2 數值求解可達域的一般流程Fig.2 General procedure for numerically solving footprint

2 偽譜軌跡規劃方法

介紹偽譜法的原理，描述將原問題轉化為非線性規劃(NLP)問題的過程，給出偽譜法求解問題的一般流程.

2.1 偽譜離散化方法

偽譜法作為直接法中的典型代表，其原理與直接法一致：基于離散化一組節點上的控制變量或狀態變量，將最優控制問題轉換為NLP問題.不同的是，偽譜法使用全局多項式同時對狀態量和控制量進行參數化，并通過正交配點來近似微分代數方程.其一般流程如下所示.

首先引入新的時間變量τ，對原問題做如下轉換：

(5)

式中：t0為初始時間；tf為終端時間.轉換后得到的狀態量、控制量均采用拉格朗日多項式進行逼近.

(6)

式中：τ∈[-1, +1]；x為實際的狀態量；X為多項式擬合的狀態量；lj(τ),j=1, 2, …,Nk+1為拉格朗日多項式；τ1,τ2, …,τNk為配點，通常可以選擇Legendre-Gauss(LG)、Legendre-Gauss-Radau(LGR)、Legendre-Gauss-Lobatto(LGL)等正交配點[27].對于拉格朗日多項式，顯然有：

(7)

對拉格朗日多項式進行微分，可以得到微分矩陣算子：

(8)

DX-f(X)=0

(9)

偽譜法中對于動力學微分方程的處理如式(9)所示，其作為等式約束與其他約束一起組成最終的NLP問題.可以看到，由于引入了微分方程的右側項f(X)，此約束一般而言是非線性約束，所以最終形成的離散問題，也就是NLP問題.關于再入過程中的其他約束，如式(3)所示，可以通過狀態量、控制量的擬合式(6)進行顯示表達.

最終，原再入軌跡規劃問題就被轉化為NLP問題.實際上，考慮到軌跡規劃問題的復雜性，求解時不僅需要增加多項式的階數，通常也需要增加多項式的個數，即為常用的ph網格細化方法，將歸一化后的總時間[-1, 1]劃分為多個區間，每一個區間對應一個拉格朗日多項式.

2.2 NLP問題及其求解流程

利用上述偽譜離散化方法，可以將原問題轉化為NLP問題.通常可以將NLP問題描述為

minΦ(Z)

(10)

s.t.Fmin≤F(Z)≤Fmax

Zmin≤Z≤Zmax

式中：Φ為性能指標函數；Z為自變量，通常包括狀態量、控制量、時間變量、參數變量等；Zmax和Zmin分別為上述各變量的上、下限；F為等式和不等式約束，通常包含關于動力學微分方程的等式約束、動壓約束、過載約束、熱流約束等；Fmax和Fmin分別對應約束的上、下限，特別地，等式約束對應的上、下限均為0.

得到如式(10)所示的NLP問題后，可以采用通用的非線性規劃求解器對問題進行求解，如稀疏非線性優化器(SNOPT)[28]或內點優化器(IPOPT)[29]等.為使NLP問題能夠被求解，需要即時計算一些關鍵信息，如梯度、雅可比矩陣、海森矩陣等.此類信息的解算速度嚴重影響了最終的計算時間，為此，文獻[30]推導出了NLP問題的梯度、雅可比矩陣和拉格朗日海森矩陣的明確表達式，確定了NLP問題的稀疏結構，可以顯著提高NLP問題的計算效率和可靠性.

在NLP問題被求解后，通常來說此時的解很難滿足期望的精度要求.因此，需要對網格進行細化，在更細密的網格下離散原問題，并求解得到更大規模的NLP問題.如此，上述網格被不斷細化，直到滿足迭代截止條件.

綜上所述，偽譜法求解優化問題的流程如圖3所示.圖中，原問題首先在初始網格條件下進行離散化得到NLP問題，然后調用相關求解器求解此類問題，最后不斷加密網格并求解離散化得到的NLP問題，直到滿足截止條件.需要注意的是，由于稀疏網格條件下NLP問題與原問題不等價，NLP問題的收斂容限可以適當增加以減少計算時間、提高求解效率.

圖3 偽譜法求解優化問題的一般流程Fig.3 General flow of pseudospectral method for solving optimization problems

3 仿真分析

將上述偽譜法應用于圖3所示的可達域求解方案中，即可快速求解再入飛行器的可達域.

3.1 標準情況下的可達域求解

飛行器在初始時刻的經度、緯度被設定為(0°，0°)，仿真采用的大氣環境為指數型大氣模型，其他參數如表1所示.仿真計算機硬件條件為AMD Ryzen 7 5800H@3.2 GHz，偽譜法求解器SWIFT由本小組獨立開發.

表1中給出了偽譜法中一些參數的設置范圍.其中，h為飛行器的滑翔高度，φS為求解換極坐標系下南半球的飛行軌跡時的緯度，φN為求解換極坐標系下北半球的飛行軌跡時的緯度.為保證問題的全局最優性，部分參數范圍取得較寬泛，如速度方位角、速度傾角等.

表1 偽譜法求解可達域時的部分參數設置Tab.1 Parameters for solving footprint by pseudospectral method

此外，一些其他參數設置如下：初始網格設置為2×2，即將總時間分為2段，每一段的軌跡采用2階多項式進行逼近；最大迭代次數設置為5；原問題期望的收斂容限設置為10-4；最大飛行時間設置為 4 500 s；過載不超過3 g，動壓不超過100 kPa，熱流不超過150 kW/m2.最終得到的可達域仿真結果如圖4所示.

圖4(a)～(c)中， 紅色、綠色、藍色曲線分別代表最大縱程軌跡、北半球最小縱程軌跡、南半球最小縱程軌跡，其中最大縱程約為136.4°，南北半球最小縱程約為3.62°；32條黑色曲線代表不同縱程條件下的最值橫程軌跡，其中所有軌跡的最值橫程約為±54.52°.將所有軌跡的終點連接在一起，即可形成最終的飛行器可達域.可見，標準情況下此類飛行器的可達域接近于橢圓形.同時給出文獻[21]中帶有攻角剖面的可達域計算結果，如圖4(d)中黃色實線所示，可見受攻角剖面的限制，這種情況下得到的可達域相對較小.此時，最大縱程約為112.6°，相比136.4°，損失約17.45%.最值橫程約為±41.53°，相比±54.52°，損失23.83%.

整體而言，對于構成可達域的圖4(a)～(c)中的35條軌跡，其速度、高度、速度方位角等隨時間變化較為規律.對于速度而言，縱程越大，速度減小得越慢.同時從4(b)也可以看到，速度減小變慢的原因在于飛行器的滑翔高度較高；從4(c)可知，整個可達域條件下的速度方位角在-90°～270°之間變化.

圖4 各條軌跡的狀態量變化Fig.4 State change for each trajectory

給出上述35條規劃軌跡的控制量變化，如圖5所示.其中，攻角的變化總是以一個較大的值開始，逐步過渡到最大升阻比攻角附近，最后階段受表1中終端速度傾角范圍設置的影響，攻角快速拉升到30°附近，以抬高速度傾角至設定范圍.傾側角的變化則明顯分為南、北半球兩部分.當規劃北半球軌跡時，傾側角多數情況下小于0°；南半球軌跡則相反.同時可以看到，縱程越大的軌跡，起點處的傾側角絕對值越小，相對時間的平均斜率的絕對值也越小；而縱程越小的軌跡，起點處的傾側角絕對值越大，相對時間的平均斜率的絕對值也越大.

圖5 各條軌跡的控制量變化Fig.5 Control changes for each trajectory

3.2 參數拉偏條件下的可達域

研究某些參數對可達域的影響，將某些參數進行拉偏，在其他參數不變的條件下求解可達域.一些拉偏情形如表2所示；仿真得到的可達域隨拉偏參數的變化如圖6所示.

表2 5種參數拉偏情形Tab.2 Five cases of parameter deviation

圖6(a)中，在質量拉偏±5%之間，飛行器的可達域基本沒有變化；而當質量拉偏+10%時，短縱程軌跡的某些約束達到上限，其最小縱程有了明顯變化，但最大縱程基本沒有變化.同樣的現象出現在圖6(b)和6(c)中，氣動參考面積和大氣密度在±5%之間波動時，可達域不發生變化；當拉偏在 -8% 時，其最小縱程有了明顯的損失，而最大縱程不變.

圖6 可達域隨參數拉偏情況的變化Fig.6 Variation of footprint with parameter deviation

結合式(1)～(3)可知，飛行器的質量、氣動參考面積、大氣密度僅出現在式(2)和式(3)中，并不顯含在動力學微分方程(1)中.并且質量、面積、大氣密度的變化僅會導致升力、阻力的同步變化，因此對于飛行器某一時刻的飛行狀態而言，通過適當改變攻角、傾側角，就可以使得式(1)中的左側微分項完全一致，則在一定范圍內，可達域是不變的.而超出此范圍后，受到熱流、過載、動壓等約束的影響，部分軌跡，特別是短縱程軌跡，可能存在達到約束上界的情況，因此必然要通過增加縱程以降低過程約束.

此外，圖6(d)為升力和阻力系數拉偏條件下的可達域計算結果.其中，實線代表升力系數拉偏的可達域，虛線代表阻力系數拉偏的可達域.可見，升力系數的增加和阻力系數的減小均會導致可達域范圍增加，反之亦然.同時，升力系數拉偏+10%和阻力系數拉偏-10% 的結果一致，升力系數拉偏-10% 和阻力系數拉偏 +10% 的結果也一致，說明升力、阻力系數的變化對可達域的影響是等價的.更準確地，升阻比的大小影響了可達域的大小：升阻比越大，可達域越大；升阻比越小，可達域越小.

4 結語

本文描述了基于偽譜法的可達域快速求解方法，并基于此方法對影響可達域的相關因素進行了研究與仿真.仿真結果表明，在給定參數條件下，飛行器質量、氣動參考面積、大氣密度等在±5%范圍內不會導致可達域的改變；超出一定范圍后，如飛行器質量增加10%、氣動參考面積減少8%、大氣密度減少8%等，均會對短縱程軌跡產生明顯影響，且體現在可達域的左半部分，其右半部分不受影響；升阻比對可達域的影響較大，其大小與可達域范圍成正相關.