返回器飞行管道的快速预测算法

邹 文,张国斌,丰志伟,涂国勇,禄晓飞,张青斌*,杨 涛

(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。考察测度空间(Ω,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返回器的搜救任务之中。设计了搜救直升机的安全空域判定流程,用于搜救直升机的实时安全飞行空域引导和搜救目标指引。通过数据可视化处理,增强了搜救态势的直观性,较好地解决了搜救安全性和搜救及时性的矛盾。

本文的研究方法也可推广应用到降落伞空投试验等相关领域。