返回器飛行管道的快速預測算法

鄒 文,張國斌,豐志偉,涂國勇,祿曉飛,張青斌*,楊 濤

(1. 國防科技大學 空天科學學院, 湖南 長沙 410073; 2. 湖南工商大學 計算機學院, 湖南 長沙 410205; 3. 中國人民解放軍63620部隊, 甘肅 酒泉 732750)

返回式航天器在完成預定任務后,其整體或其中一部分需要再入大氣層并在地面安全著陸。髙精度的落點預報結果和合理的搜救策略可極大地提高航天器返回器的回收效率和安全性[1-2]。返回式航天器的返回器在傘降著陸階段受氣象風的影響很大,如遇到大的橫風,返回器著陸點將產生很大偏移,造成目標搜救散布區擴大、搜救時間延長。我國“神舟”飛船返回器落點預報中,通過建立飛船回收過程動力學模型,結合“神舟”任務返回段實測數據分析,整個傘降著陸段飛行時間誤差有時高達30%,影響了氣象風漂移修正量和落點預報的精度。此外,對于嫦娥五號(Chang′e-5, CE-5)等夜間再入的返回器,光學探測手段受到了嚴重限制,增加了回收程序啟動后返回器及其分離物與搜救直升機碰撞的風險[3-4]。本文在飛船降落傘回收系統動力學的研究基礎上,將針對返回器在著陸場的搜救任務,探索采用高效的數值手段預測返回器在不確定因素作用下的可能運動空間,減小飛船的搜救范圍,縮短搜救時間,以提高地面返回搜救效率,保證第一時間找到返回器,確保搜救人員安全。

返回器的傘降運動軌跡受多種不確定因素共同影響,包括:初始運動參數、狀態參數以及環境參數等,其中風場的影響最為顯著。目前在工程領域中廣泛采用蒙特卡羅(Monte Carlo, MC)方法研究返回式航天器傘降階段的不確定量化評估[5-6]。張青斌等[7]利用MC方法對返回器和傘艙蓋的運動軌跡進行了分析,評估了返回器和傘艙蓋相撞的風險。王慧娟[8]利用MC方法研究了火星再入自適應開傘控制方法。美國國家航空航天局在對返回器進入、下降和著陸問題的研究中也大量采用了MC方法進行分析[9-11]。近十幾年,出現了從概率守恒原理的角度出發,應用概率密度演化方法研究動力學系統隨機特征的方法[12-13]。在概率密度演化方法中,系統的狀態變量被視為隨機變量,其真實運動對應于隨機狀態空間中的一條軌跡,隨機空間中的每一點均有一概率密度值與之對應。隨機劉維爾方程(stochastic Liouville equation, SLE)是概率密度演化方法中最為經典也是應用最為廣泛的理論方法[14]。Halder等[15]和Trisolini等[16]利用SLE方程研究了飛行器的超音速再入問題,通過與MC方法的對比計算,展示了SLE方程較高的計算精度和計算效率。Leonard等[17-18]和Klein等[19]利用SLE方程對不確定條件下的最優空投釋放點、主傘開傘點進行了分析和優化設計。盡管SLE方程已經成功應用于上述工程問題中,但是該方法仍存在兩點不足[20]:一是對于一些穩定系統,伴隨時間演化隨機狀態變量所在的空間會迅速收縮,這將導致概率密度值迅速增大,繼而使得積分困難;二是隨著狀態空間的畸變,對高維隨機變量的邊緣概率密度積分也會出現較大的誤差。為了解決這些問題,Leonard[20]利用Koopman算子和Frobenius-Perron算子的對偶特性以及Koopman算子的后拉機制[21]計算不確定量的期望,這種方法將系統末狀態與初始概率密度關聯,從而規避了對高維離散點的概率密度積分,在貨物空投的優化應用中取得了較好的效果[22]。利用Koopman算子后拉機制進行期望值計算不僅在收斂速度上優于MC方法,并且規避了SLE方法所引起的數值計算問題。在航天領域中,返回器及開傘過程中的分離物的軌跡預測是搜索任務中的重要計算需求。比如,獵戶座飛船測試小組在對飛船空投任務的落點區域規劃任務中,先后開發了名為Sasquatch的落點預測軟件[23-24],該軟件可分別計算出返回器、傘艙蓋和傘包等部件的期望落點及可能范圍,結合Debris Tool軟件還可給出搜救直升機是否處于安全空域。總的來說,目前工程上大多采用MC方法進行不確定的量化分析,但是針對返回落地搜救任務,需要研發一種高效、實時算法。

圍繞返回器的落點預報和搜救問題,在給出返回器傘降階段動力學模型的基礎上,結合Koopman算子后拉機制、Halton采樣[25](Halton sampling, HS)和正交設計[26](orthogonal design, OD)給出了返回器飛行管道的計算方法,通過計算分析說明本文提出方法的高效性。以飛行管道計算數據為支撐,研究了搜救直升機的安全飛行空域判定辦法,以嫦娥五號返回器的搜救任務為例,對搜救直升機飛行空域和飛行路線進行了分析和優化設計。

1 返回器降落傘回收階段的動力學模型

1.1 降落傘回收過程簡述

返回器一般選取彈道-升力方式再入,在進入大氣層后通過控制升力方向來調整彈道,從而在一定程度上減小落點散布范圍。當返回器沿著再入彈道到達指定開傘高度后,系統開始啟動降落傘減速程序。嫦娥五號返回器回收過程如圖1所示,當返回器到達開傘高度后,減速系統大體按以下步驟實施:①傘艙蓋以一定速度與返回器彈射分離,同時將引導傘從傘艙中拉出并打開,在引導傘的牽引下又將減速傘拉出并充滿;②減速傘與返回器分離并拉出主傘,主傘首先呈收口狀充滿;③主傘解除收口,完全充滿;④返回器著陸,切斷與主傘的連接。在返回器的降落傘減速階段,搜救直升機必須足夠接近返回器,但是又必須處于足夠安全的區域,避免與返回器相碰。

圖1 嫦娥五號返回器回收過程Fig.1 Recovery process of the reentry capsule of CE-5

1.2 物傘系統質心運動動力學

物傘系統的空間位置是返回搜救任務中測量跟蹤力量部署的重要依據。為了方便,將再入前的軌道動力學計算結果的末端狀態和降落傘減速階段的初始狀態相銜接,本文在地心坐標系下建立物傘系統的傘降運動方程。

1.2.1 坐標系定義

(1)

(2)

(a) 地心坐標系與地面觀測系(a) Geocentric coordinate system and ground observation coordinate system

(b) 地面觀測系與速度坐標系(b) Ground observation coordinate system and velocity coordinate system

1.2.2 動力學方程

為方便搜救指揮和地面觀測,本文選取返回器質心地心距r、地心緯度φ和經度θ描述物傘系統的位置,以速度大小V、速度偏角ψ和速度傾角γ描述返回器質心速度矢量。返回器傘降階段的動力學方程[27]可以表示為

(3)

在返回器物傘系統的運動過程中,降落傘所產生的阻力遠大于系統的升力,因此開傘后的動力學方程中的升力加速度項L可設為0。物傘系統所產生加速度可表示為

(4)

式中,ρ為大氣密度,Cd為阻力系數,S為參考面積,m為系統當前質量,ma為系統當前的附加質量。阻力面積CdS和附加質量ma可用以描述系統當前所處的物理狀態。因此,采用合理的模型量化物傘系統在各個階段的阻力系數及附加質量的變化規律即可對整個回收過程進行數值仿真。

1.2.3 降落傘充氣模型

降落傘充氣過程是銜接各個工作階段的重要過程,也是物傘系統工作中最為復雜的過程。降落傘充氣展開涉及柔性織物與氣流的流固耦合作用、柔性體的幾何非線性運動以及柔性體的摩擦接觸,是一個具有強非線性的物理過程。依靠仿真計算模擬這一現象極為耗時,并且難以準確驗證仿真結果的正確性。在工程實踐中,利用大量試驗數據總結出的經驗公式已經在工程實踐中取得了大量的應用,也被證明是一種簡單可靠的方法[28]。本文采用如式(5)所示近似模型描述降落傘在充氣過程中的阻力面積的變化規律[28]。

(5)

式中:(CdS)1為初始充氣時刻tf1時的阻力面積;(CdS)2為充氣完成時刻tf時的阻力面積;n為充氣指數,一般可根據給定的降落傘型號查表得到。

不同于大氣中飛行的一般飛行器,降落傘的附加質量會對物傘系統動力學特性產生顯著的影響。工程中附加質量的簡化計算公式[29]為

(6)

式中,Ka為附加質量系數,ρ∞為傘質心處的大氣密度,D0為降落傘的名義直徑。

2 返回器飛行管道計算方法

在著陸場返回器搜救任務中,搜救直升機飛行引導系統需要和指揮控制系統進行高頻、及時的信息交互,為搜救直升機提供安全空域引導和前往區域坐標。返回器物傘系統的自身模型和環境參數具有諸多不確定性,要在指揮控制系統信息更新時間間隔內獲取高維隨機空間的概率分布信息是十分困難的。因此,本文采取的研究路線是采取高效的隨機空間采樣方法和期望彈道算法,在較少的采樣計算次數下計算期望彈道和包絡返回器及其分離物可能存在空間的飛行管道。盡管計算結果缺乏返回器的時空概率分布信息,但在搜救直升機安全空域引導和搜救方向引導方面是能夠滿足工程需求的。

2.1 基于Koopman算子理論計算期望彈道

設一動力學系統的狀態空間為Ω,Ω?d,d為狀態變量的維數。狀態變量在空間Ω中的概率密度分布函數為f,f∈L1且f≥0?？疾鞙y度空間(Ω,A,μ),其中A是σ代數[21],μ為概率度量,定義為

(7)

式中,A為σ代數子集,A∈A。

系統還存在觀測函數g:Ω→,g∈L∞。系統從初始時刻到T時刻的動力學映射表示為S:Ω→Ω,系統初始時刻狀態向量為x0,對應在空間Ω中的概率密度為f0,系統T時刻的狀態向量為xT,對應的空間Ω中的概率密度為fT。定義xT=S(x0),Koopman算子U:L∞→L∞定義為

g(xT)=Ug(x0)

(8)

Frobenius-Perron算子P:L1→L1定義為

fT=Pf0

(9)

其中,Koopman算子和Frobenius-Perron算子間的對偶關系[21]表示為

〈PfT,g〉=〈f0,Ug〉

(10)

式中,〈·〉表示內積運算,即

(11)

需要指出,式(10)成立要求觀測函數g在狀態空間中具有緊支性。

從初始時刻到T時刻,在系統動力學演化不加入任何新的隨機因素的情況下,有

(12)

式中,S-1(A)為σ代數子集A的逆映射集,若動力學映射S滿足可逆和可微條件,則進一步有

(13)

式中,|dS-1(x)/dx|是動力學逆映射S-1(x)雅可比矩陣的行列式。對于一些簡單的動力學系統,利用式(13)求取Frobenius-Perron算子作用下狀態空間概率密度方程的演化是可行的。對于為常微分方程組控制的復雜動力學系統,利用式(13)求取概率密度演化方程則是非常困難的,此時可利用Koopman算子和Frobenius-Perron算子間的對偶關系式(10)推導出經典的SLE方程,有關SLE方程求解系統概率密度函數演化可以參考文獻[13,15-16]。

系統觀測量g在時刻T的期望值可表示為

E[g]=〈fT,g〉=〈Pf0,g〉

(14)

利用Koopman算子和Frobenius-Perron算子間的對偶關系,觀測量g在時刻T的期望值還可以表示為

E[g]=〈f0,Ug〉

(15)

從式(14)右端來看,觀測量g在T時刻的期望值為Frobenius-Perron算子作用下的概率密度函數與觀測量的內積;從式(15)來看,觀測量g在T時刻的期望值也等于系統T時刻的觀測量與初始概率密度的內積,注意到在式(15)中,Koopman算子作用的觀測量與初始狀態的概率密度函數相關聯,即將當前觀測量拉回至與初始概率密度函數相關聯,這一機制被稱為Koopman算子的后拉機制。一般而言,觀測量是離散形式的,因此,對觀測量的期望計算涉及高維空間的離散數據的積分問題。對該類問題的處理可以采用Binning方法[20]、Lobachevsky樣條曲線[30]、Radial Basis Function方法[31]等進行計算?？紤]到Halton采樣的空間均勻性,為簡化計算,本文中利用Binning方法計算觀測量的期望值,即

(16)

2.2 返回器飛行管道計算流程

返回器再入過程中包含多種不確定因素,具體可分為環境不確定因素和物-傘系統自身不確定因素。環境不確定因素包括風速、風向和大氣密度,物-傘系統自身不確定因素包括初始位置(地心距、經度,地心緯度)、初始速度(速度大小、速度傾角、速度偏角)、阻力系數和返回器質量。本文中認為上述11個不確定因素服從正態分布且相互獨立,聯合概率密度分布函數f為上述11個不確定因素的函數。

在上述11個隨機變量組成的高維隨機空間中進行全面采樣的計算成本是十分昂貴的,即使每個隨機因素選取3個水平,采樣數也會超過17萬次。另一方面,根據搜救指揮總體任務需求,飛行管道的單次更新一般需要在3 s內完成,以便根據實時測量數據更新計算結果。本文中利用C++語言編寫動力學計算程序,利用變步長龍格庫塔法進行動力學積分[32],在一臺CPU時鐘頻率為2.5 GHz的電腦上,一條彈道的平均計算時間約為2.9 ms,在保證任務需求前提下,采樣計算的彈道總數最大不能超過1 020。由于計算時間的限制,必須對隨機空間的采樣策略進行合理的設計,以保證通過較少次數的采樣計算得到精度較高的期望彈道和飛行管道。以返回器為例,對飛行管道的計算可分為如下3個步驟。

重復以上步驟,可以計算出傘艙蓋和減速傘等返回器分離物的期望彈道和飛行管道信息。在實際應用中,采樣點個數N應根據硬件能力和系統更新頻率要求進行限制,在本文中采樣點個數設為500。

以返回器為例,判定搜救直升機是否位于遠離返回器飛行管道的一般流程為:

1)輸入當前時間和位置,根據當前高度分別計算該高度下返回器的安全時間以及飛行管道半徑(安全半徑)。

2)判定直升機是否在安全時間之后位于該高度,若是,直升機是安全的,否則轉步驟3。步驟2稱為時間安全判定。

3)計算得到直升機與返回器期望位置間的距離,判定該距離是否大于當前高度的管道半徑,若是,則直升機是安全的,否則,當前位置存在風險,應沿返回器期望位置所在方向的相反向飛離。步驟3稱為空間安全判定。

上述搜救直升機是否處于安全空域的判定流程如圖3所示。重復上述過程可判定直升機是否位于傘艙蓋和減速傘的安全空域。

圖3 安全空域的判定流程Fig.3 Determining process of the safe airspace

3 返回器搜救任務不確定性分析

本節以CE-5返回器在不確定條件下的安全空域計算問題和搜救引導問題為例進行分析。安全空域計算結果確定了每一時刻的最低安全高度和返回器及其分離物的管道,為空中和地面搜救人員提供安全指引,同時給出返回器、傘艙蓋以及減速傘的期望落點,引導搜救力量對返回器進行回收。CE-5返回器的主要標稱參數如表1所示。利用探空氣球測量的著陸場風場數據作為仿真風場輸入計算返回器的飛行彈道,不同高度上的風速大小如圖4所示。系統的不確定性參數均視為以名義值為期望值的正態分布,參考已有相關研究[9]和CE-5返回器的自身特點,其3倍標準差的范圍如表2所示。

表1 物傘系統參數名義值

圖4 名義風場Fig.4 Nominal wind field

表2 回收系統不確定參數及取值范圍

不同采樣方法得到的在地面處的飛行管道半徑如圖5所示,其中HS方法和MC方法的采樣次數均為10萬次,OD方法采用標準正交表L12(211)進行采樣,從圖5中結果可見,HS方法得到的管道半徑要明顯大于MC方法得到的計算結果,正交設計采樣得到的最大半徑又要大于HS方法計算得到的管道半徑,極差分析得到的正交設計最優方案(optimized plan of orthogonal design, OPOD)的管道半徑又略有增大。通過對正交設計采樣計算結果進行極差分析得到使飛行管道半徑最大的隨機參數取值方案,方案如表3所示,其中,ψw表示風向偏航角,Vw表示風速大小。在飛行管道半徑的試驗方案中,返回器的質量和速度都取最小值,其他因素均取最大值,這與返回器傘降過程的物理特性是一致的:返回器質量最小,初速度最小,風速和空氣密度最大,能夠使得返回器達到最大的風飄距離。

圖5 不同方法計算得到的飛行管道半徑Fig.5 Fight pipeline radius calculated by different methods

表3 正交試驗分析得到的偏離期望落點最大的方案Tab.3 Scheme with the largest deviation from the expected landing point obtained by orthogonal test analysis

期望彈道是返回器飛行管道計算的必要結果,由于更新時間的需要,期望彈道的計算量必須限制在一定的數量下。對比HS方法和MC方法計算期望落點東向坐標的收斂過程,如圖6所示,從圖中可見,HS方法的收斂速度要遠優于MC方法,在采樣次數小于1 000的限制下,MC方法計算的期望落點坐標仍處于不穩定的波動狀態,而HS方法已經收斂至落點坐標25 m以內的狀態。顯然,HS方法能夠更快地收斂至期望落點。

圖6 落點坐標期望值隨采樣次數增加的收斂曲線Fig.6 Convergence curve of the expected landing point coordinate with increasing of sampling times

在給定輸入下計算得到的返回器、傘艙蓋和減速傘的飛行管道分別如圖7(a)～(c)所示。圖中直角坐標系為原點在開傘點地面投影點的天東北坐標系,名義彈道指按名義參數計算得到的彈道。

(a) 返回器飛行管道(a) Flight pipeline of the reentry capsule

(b) 傘艙蓋飛行管道(b) Flight pipeline of the parachute cap

(c) 減速傘飛行管道(c) Flight pipeline of the drogue parachute

從飛行管道的計算結果可見,由于當地風場風向偏東,返回器、傘艙蓋和減速傘均明顯向東偏移。因為質量偏小,傘艙蓋和減速傘受風場影響相比于返回器更明顯,落地時刻傘艙蓋和減速傘主傘包組合體的可能存在半徑相比于返回器分別增大了25%和66%,計算結果反映了實際搜救任務中傘艙蓋和減速傘難以找到的實際問題。通過對數據的處理,不同高度h下的返回器飛行管道半徑R和安全時間Ts如表4所示,表中的安全時間從降落傘打開指令發出開始計時,飛行管道半徑和安全時間的變化趨勢如圖8所示,伴隨高度的增加,飛行管道半徑和安全時間呈減小趨勢,這與系統基本特性是相符的。

表4 返回器在不同高度下的飛行管道半徑和安全時間Tab.4 Flight pipeline radius and safe time of the reentry capsule under different heights

(a) 飛行管道半徑(a) Flight pipeline radius

(b) 安全時間(b) Safe time

以搜救直升機在距離地面2 km高度工作為例,對返回器的搜救任務進行說明。一方面,根據計算結果,當搜救直升機在名義開傘指令下達后的489 s、551 s和987 s以后出現在距離地面2 km高度上,將幾乎不可能與返回器、傘艙蓋和減速傘發生空中碰撞。另一方面,在考慮隨機因素擾動的情況下,返回器最晚會在名義開傘指令下達后的667 s落地。地面2 km高度上的飛行管道截面如圖9所示,返回器期望落點包含于3個飛行管道之內,因此,要保證在返回器落地時,搜救直升機出現在落地點上空,搜救直升機必須承擔一定與傘艙蓋或減速傘相碰撞的風險。在保證搜救直升機不與返回器相撞的前提下,應盡可能合理地規劃直升機飛行路線以降低搜救直升機與返回器分離物相撞的風險。

圖9 地面2 km高度上的飛行管道截面Fig.9 Flight pipeline cross section at 2 km above ground

搜救直升機的搜索流程按如下步驟實施:①在開傘指令下達后的489 s之前,可由安全判定流程使搜救直升機位于2 km高度的飛行管道之外;②489 s以后,從返回器、傘艙蓋和減速傘飛行管道最外邊界選取一個最優出發點開始進行搜救。

為定量刻畫搜救直升機與返回器分離物相撞的風險,定義碰撞風險指標為

(17)

其中,s為飛行路線,y和z分別為飛行路線點在天東北坐標系下的東向和北向坐標,y2和z2為傘艙蓋在2 km高度的期望坐標,y3和z3為減速傘在2 km高度的期望坐標。設搜救直升機的飛行速度為260 km/h,從出發點至期望落點,飛行路線為直線。通過計算得到從相對期望落點不同方向管道邊界點進行搜救的風險指標數值和達到期望落點所需時間的曲線,如圖10所示。最優進入點及搜索方向標注于圖9中,飛行方向為北偏東1.6°,搜救出發點距離期望落點4.2 km, 所需飛行時間為59 s。當搜救直升機達到期望落點時,名義狀態下,返回器仍需24 s才能落地。相比于最優進入點,若搜救直升機從位于期望落點正東方向的管道邊界點進行搜索,則對應的風險指標將增大116%,搜索所需時間將增大2.4倍,可見,搜索路線和搜索出發點的優化設計對于搜救任務的安全性和快速性均有較大的提升。

圖10 期望落點不同方向搜索出發點對應的 風險指標和期望搜索時間Fig.10 Risk index and expected search time corresponding to the departure point in different directions of the expected landing point

4 結論

在航天器返回搜救任務中,對著陸器的飛行管道進行高效、準確的預報分析具有非常重要的工程應用價值。本文為了將仿真算法融合到搜救指揮系統中,重點研究了返回器減速過程的不確定量化的分析算法,主要結論如下:

1) 基于Koopman算子理論的快速預測算法可以應用到航天器返回過程的不確定量化分析之中。Koopman算子將初始概率密度值與當前狀態關聯,避免了概率積分和高維離散點數值積分所帶來的煩瑣步驟和數值誤差。利用Halton采樣方法在隨機空間采點的均勻性和正交設計方法,可進一步縮減計算量,以較少的采樣點得到系統響應邊界。

2) 本文的仿真算法已應用到CE-5返回器的搜救任務之中。設計了搜救直升機的安全空域判定流程,用于搜救直升機的實時安全飛行空域引導和搜救目標指引。通過數據可視化處理,增強了搜救態勢的直觀性,較好地解決了搜救安全性和搜救及時性的矛盾。

本文的研究方法也可推廣應用到降落傘空投試驗等相關領域。