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

邹文; 张国斌; 丰志伟; 涂国勇; 禄晓飞; 张青斌; 杨涛; ZOU Wen; ZHANG Guobin; FENG Zhiwei; TU Guoyong; LU Xiaofei; ZHANG Qingbin; YANG Tao

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

doi: 10.11887/j.cn.202401003

邹文^1,2 ，张国斌¹ ，丰志伟¹ ，涂国勇³ ，禄晓飞³ ，张青斌¹ ，杨涛¹

1. 国防科技大学空天科学学院, 湖南长沙 410073

2. 湖南工商大学计算机学院, 湖南长沙 410205

3. 中国人民解放军63620部队，甘肃酒泉 732750

基金项目: 国家自然科学基金资助项目（11902331）

详细信息

作者简介

邹文（1981—），女，湖南长沙人，博士研究生，E-mail: 8320763@qq.com

通讯作者

张青斌（1975—），男，山西大同人，研究员，博士，博士生导师，E-mail: qingbinzhang@sina.com

中图分类号: V555

文献标识码: A

文章编号: 1001-2486(2024)01-022-10

Fast prediction algorithm of flight pipeline of reentry capsule

ZOU Wen^1,2 ， ZHANG Guobin¹ ， FENG Zhiwei¹ ， TU Guoyong³ ， LU Xiaofei³ ， ZHANG Qingbin¹ ， YANG Tao¹

1. College of Aerospace Science and Engineering, National University of Defense Technology, Changsha 410073 , China

2. School of Computer Science, Hunan University of Technology and Business, Changsha 410205 , China

3. The PLA Unit 63620, Jiuquan 732750 , China

摘要

针对返回器回收任务中对安全空域和期望落点的计算需求，提出了基于Koopman算子的飞行管道快速预测算法，给出了搜救直升机安全飞行空域的判定流程。建立了物伞动力学模型，利用Halton采样方法从随机空间中均匀采点，计算得到多条可能弹道；采用Koopman算子的后拉机制，将初始概率密度值与当前状态关联，得到不确定条件下返回器及其分离部件的飞行管道和期望弹道。仿真结果表明，基于Koopman算子的飞行管道快速预测算法在收敛速度和精度上都要显著优于Monte Carlo方法；利用飞行管道计算结果对搜救直升机飞行路线进行规划后，碰撞风险最大降低54%且搜索时间减少70%。飞行管道预测算法已成功应用到嫦娥五号的回收任务中。

关键词

降落伞回收技术 / 不确定性分析 / Koopman算子 / 返回器 / 探月任务

Abstract

To meet the computational requirements of safe airspace and the expected landing point in the recovery mission of the reentry capsule, a fast prediction algorithm of the flight pipeline based on the Koopman operator approach was proposed for the reentry capsule and the determination process for safe flight airspace of search and rescue helicopters was provided. The body-parachute dynamic model was constructed. A group of discrete state points was uniformly selected from the random state space by using the Halton sampling method, and the multiple possible trajectory was calculated. Based on the back pulling mechanism of Koopman operator, the initial probability density value was associated with the current state to obtain the flight pipeline and desired trajectory of the reentry capsule and its separation parts under uncertain conditions.The simulation results show that the fast prediction algorithm of the flight pipeline based on the Koopman operator approach is significantly better than the Monte Carlo method in terms of convergence speed and accuracy. After using the flight pipeline calculation results to plan the flight route of the rescue helicopter, the collision risk is reduced by 54% at most and the corresponding search time is reduced by 70%. The proposed algorithm has been successfully applied to the Chang′e-5 recovery mission.

Keywords

parachute recovery technology / uncertainty analysis / Koopman operator / reentry capsule / lunar missions

1 返回器降落伞回收阶段的动力学模型 1.1 降落伞回收过程简述 1.2 物伞系统质心运动动力学 1.2.1 坐标系定义 1.2.2 动力学方程 1.2.3 降落伞充气模型 2 返回器飞行管道计算方法 2.1 基于Koopman算子理论计算期望弹道 2.2 返回器飞行管道计算流程 3 返回器搜救任务不确定性分析 4 结论

返回式航天器在完成预定任务后，其整体或其中一部分需要再入大气层并在地面安全着陆。髙精度的落点预报结果和合理的搜救策略可极大地提高航天器返回器的回收效率和安全性^[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.2 物伞系统质心运动动力学

物伞系统的空间位置是返回搜救任务中测量跟踪力量部署的重要依据。为了方便，将再入前的轨道动力学计算结果的末端状态和降落伞减速阶段的初始状态相衔接，本文在地心坐标系下建立物伞系统的伞降运动方程。

1.2.1 坐标系定义

返回器再入过程主要涉及的坐标系有：地心坐标系O_E-X_EY_EZ_E、地面观测系O-xyz以及速度坐标系o₁-x_vy_vz_v。地心坐标系O_E位于地球质心，O_EX_E指向某时刻起的零度子午线，O_EZ_E垂直于赤道平面指向北极，O_EY_E与O_EX_E和O_EZ_E构成右手直角坐标系。地面观测系原点O固定于地球上的某一点，Ox轴延长线交于地心O_E，方向背离地心，Oz轴指向北方，Oy轴与Ox轴和Oz轴构成右手直角坐标系。速度坐标系原点o₁位于系统质心，o₁y_v沿质心速度方向，o₁x_v在铅垂平面Oyz内，o₁z_v与o₁x_v和o₁y_v构成右手直角坐标系。地心坐标系与地面观测系的变换如图2（a）所示，从地心坐标系到地面观测系的旋转顺序为：首先绕Z_E轴旋转经度θ角；然后，绕新坐标系的Y^（1）_E轴的反向旋转地心纬度φ角，此时坐标系各轴与地面观测系平行，从地心坐标系到地面观测系的坐标变换矩阵为

图1嫦娥五号返回器回收过程

Fig.1Recovery process of the reentry capsule of CE-5

G_{E} = [\begin{matrix} c o s θ c o s ϕ & s i n θ c o s ϕ & s i n ϕ \\ - s i n θ & c o s θ & 0 \\ - c o s θ s i n ϕ & - s i n θ s i n ϕ & c o s ϕ \end{matrix}]

(1)

地面观测系与速度坐标系的变换如图2（b）所示，从速度坐标系到地面观测系的旋转顺序为：将地面观测系原点平移至飞行器质心o₁以便观察，首先绕z_v轴旋转γ角，使o₁y_v位于o₁xz平面内；然后绕新坐标系x^（1）_v轴的反向旋转ψ角，得到各坐标轴与地面观测系平行的坐标系。其中，γ角称为速度倾角，速度方向在水平面上方时为正；ψ角称为速度偏角，速度位于正东方向时，偏角为0，从东向北旋转为正。从速度坐标系到地面观测系的坐标变换矩阵为

G_{v} = [\begin{matrix} c o s γ & s i n γ & 0 \\ - c o s ψ s i n γ & - s i n ψ c o s γ & - s i n ψ \\ - s i n ψ s i n γ & s i n ψ c o s γ & c o s ψ \end{matrix}]

(2)

图2坐标系系统

Fig.2Coordinate system

1.2.2 动力学方程

为方便搜救指挥和地面观测，本文选取返回器质心地心距r、地心纬度φ和经度θ描述物伞系统的位置，以速度大小V、速度偏角ψ和速度倾角γ描述返回器质心速度矢量。返回器伞降阶段的动力学方程^[27]可以表示为

\{\begin{matrix} \dot{r} = V s i n γ \\ \dot{θ} = = \frac{V c o s γ c o s ψ}{r c o s ϕ} \\ \dot{ϕ} = \frac{V c o s γ s i n ψ}{r} \\ \dot{V} = L Λ_{V}^{L} + D Λ_{V}^{D} + g_{V} + R_{V} \\ \dot{γ} = \frac{L}{V} Λ_{γ}^{L} + \frac{D}{V} Λ_{γ}^{D} + g_{γ} + \frac{V}{r} c o s γ + R_{γ} \\ \dot{ψ} = \frac{L}{V c o s γ} Λ_{ψ}^{L} + \frac{D}{V c o s γ} Λ_{ψ}^{D} - \\ \frac{V}{r} c o s γ c o s ψ t a n ϕ + g_{ψ} + R_{ψ} \end{matrix}

(3)

其中，L和D分别为升力和阻力所产生的加速度，g_V、g_γ和g_ψ对应于重力所产生的加速度，R_V、R_γ和R_ψ对应于地球转动所产生的加速度，Λ^V_L、Λ^D_V、Λ^L_γ、Λ^D_γ、Λ^L_ψ和Λ^D_ψ为风干扰系数。式（3）中系数的详细表达可参考文献^[27]。

在返回器物伞系统的运动过程中，降落伞所产生的阻力远大于系统的升力，因此开伞后的动力学方程中的升力加速度项L可设为0。物伞系统所产生加速度可表示为

D = \frac{ρ C_{d} S V^{2}}{2 (m + m_{a})}

(4)

式中， ρ为大气密度，C_d为阻力系数，S为参考面积，m为系统当前质量，m_a为系统当前的附加质量。阻力面积C_dS和附加质量m_a可用以描述系统当前所处的物理状态。因此，采用合理的模型量化物伞系统在各个阶段的阻力系数及附加质量的变化规律即可对整个回收过程进行数值仿真。

1.2.3 降落伞充气模型

降落伞充气过程是衔接各个工作阶段的重要过程，也是物伞系统工作中最为复杂的过程。降落伞充气展开涉及柔性织物与气流的流固耦合作用、柔性体的几何非线性运动以及柔性体的摩擦接触，是一个具有强非线性的物理过程。依靠仿真计算模拟这一现象极为耗时，并且难以准确验证仿真结果的正确性。在工程实践中，利用大量试验数据总结出的经验公式已经在工程实践中取得了大量的应用，也被证明是一种简单可靠的方法^[28]。本文采用如式（5）所示近似模型描述降落伞在充气过程中的阻力面积的变化规律^[28]。

C_{d} S = {(C_{d} S)}_{1} + [{(C_{d} S)}_{2} - {(C_{d} S)}_{1}] {(\frac{t - t_{f l}}{t_{f} - t_{f l}})}^{n}

(5)

式中：（C_dS）₁为初始充气时刻t_f1时的阻力面积；（C_dS）₂为充气完成时刻t_f时的阻力面积；n为充气指数，一般可根据给定的降落伞型号查表得到。

不同于大气中飞行的一般飞行器，降落伞的附加质量会对物伞系统动力学特性产生显著的影响。工程中附加质量的简化计算公式^[29]为

m_{a} = K_{a} π ρ_{\infty} D_{0}^{3}

(6)

式中，K_a为附加质量系数，ρ_∞为伞质心处的大气密度，D₀为降落伞的名义直径。

2 返回器飞行管道计算方法

在着陆场返回器搜救任务中，搜救直升机飞行引导系统需要和指挥控制系统进行高频、及时的信息交互，为搜救直升机提供安全空域引导和前往区域坐标。返回器物伞系统的自身模型和环境参数具有诸多不确定性，要在指挥控制系统信息更新时间间隔内获取高维随机空间的概率分布信息是十分困难的。因此，本文采取的研究路线是采取高效的随机空间采样方法和期望弹道算法，在较少的采样计算次数下计算期望弹道和包络返回器及其分离物可能存在空间的飞行管道。尽管计算结果缺乏返回器的时空概率分布信息，但在搜救直升机安全空域引导和搜救方向引导方面是能够满足工程需求的。

2.1 基于Koopman算子理论计算期望弹道

设一动力学系统的状态空间为Ω，

Ω \subset R^{d}

，d为状态变量的维数。状态变量在空间Ω中的概率密度分布函数为f，f∈L¹且f≥0。考察测度空间（Ω，

A

，μ），其中

A

是σ代数^[21]，μ为概率度量，定义为

μ (A) = \int_{A} f (x) d x

(7)

式中，A为σ代数子集，A∈

A

。

系统还存在观测函数g:Ω→

R^{}

，g∈L^∞。系统从初始时刻到T时刻的动力学映射表示为S:Ω→Ω，系统初始时刻状态向量为x₀，对应在空间Ω中的概率密度为f₀，系统T时刻的状态向量为x_T，对应的空间Ω中的概率密度为f_T。定义x_T=S（x₀），Koopman算子U:L^∞→L^∞定义为

g (x_{T}) = U g (x_{0})

(8)

Frobenius-Perron算子P:L¹→L¹定义为

f_{T} = P f_{0}

(9)

其中，Koopman算子和Frobenius-Perron算子间的对偶关系^[21]表示为

⟨P f_{T}, g⟩ = ⟨f_{0}, U g⟩

(10)

式中，〈·〉表示内积运算，即

⟨ u, v ⟩ = \int_{Ω} u (x) v (x) d x

(11)

需要指出，式（10）成立要求观测函数g在状态空间中具有紧支性。

从初始时刻到T时刻，在系统动力学演化不加入任何新的随机因素的情况下，有

\int_{A} P f (x) d x = \int_{S^{- 1} (A)} f (x) d x, \forall A \in A

(12)

式中，S^-1（A）为σ代数子集A的逆映射集，若动力学映射S满足可逆和可微条件，则进一步有

P f (x) = f (S^{- 1} (x)) |\frac{d S^{- 1} (x)}{d x}|

(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] = ⟨f_{T}, g⟩ = ⟨P f_{0}, g⟩

(14)

利用Koopman算子和Frobenius-Perron算子间的对偶关系，观测量g在时刻T的期望值还可以表示为

E [g] = ⟨f_{0}, U g⟩

(15)

从式（14）右端来看，观测量g在T时刻的期望值为Frobenius-Perron算子作用下的概率密度函数与观测量的内积；从式（15）来看，观测量g在T时刻的期望值也等于系统T时刻的观测量与初始概率密度的内积，注意到在式（15）中，Koopman算子作用的观测量与初始状态的概率密度函数相关联，即将当前观测量拉回至与初始概率密度函数相关联，这一机制被称为Koopman算子的后拉机制。一般而言，观测量是离散形式的，因此，对观测量的期望计算涉及高维空间的离散数据的积分问题。对该类问题的处理可以采用Binning方法^[20]、Lobachevsky样条曲线^[30]、Radial Basis Function方法^[31]等进行计算。考虑到Halton采样的空间均匀性，为简化计算，本文中利用Binning方法计算观测量的期望值，即

E [g] = \frac{1}{\overset{N}{\underset{i = 1}{Σ}} f_{0}^{(i)}} \overset{N}{\underset{i = 1}{Σ}} g^{(i)} f_{0}^{(i)}

(16)

2.2 返回器飞行管道计算流程

返回器再入过程中包含多种不确定因素，具体可分为环境不确定因素和物-伞系统自身不确定因素。环境不确定因素包括风速、风向和大气密度，物-伞系统自身不确定因素包括初始位置（地心距、经度，地心纬度）、初始速度（速度大小、速度倾角、速度偏角）、阻力系数和返回器质量。本文中认为上述11个不确定因素服从正态分布且相互独立，联合概率密度分布函数f为上述11个不确定因素的函数。

在上述11个随机变量组成的高维随机空间中进行全面采样的计算成本是十分昂贵的，即使每个随机因素选取3个水平，采样数也会超过17万次。另一方面，根据搜救指挥总体任务需求，飞行管道的单次更新一般需要在3 s内完成，以便根据实时测量数据更新计算结果。本文中利用C++语言编写动力学计算程序，利用变步长龙格库塔法进行动力学积分^[32]，在一台CPU时钟频率为2.5 GHz的电脑上，一条弹道的平均计算时间约为2.9 ms，在保证任务需求前提下，采样计算的弹道总数最大不能超过1 020。由于计算时间的限制，必须对随机空间的采样策略进行合理的设计，以保证通过较少次数的采样计算得到精度较高的期望弹道和飞行管道。以返回器为例，对飞行管道的计算可分为如下3个步骤。

步骤1：期望弹道计算。将所有初始不确定状态变量和不确定参数的取值范围在[-3σ，3σ]内截断，用Halton采样方法^[33]在高维空间中取N-2个采样点，计算每个采样点概率密度值f^（ⁱ^）₀，i=1，2，···，N，随后计算所有采样点对应随机参数和初始条件下的弹道，存储相应数据，利用式（16）计算返回器或分离物的期望弹道。

步骤2：正交试验分析。利用正交设计方法获取使飞行管道半径和飞行时间最长的方案。采用标准正交表L₁₂（2¹¹）进行采样，将随机变量的第一水平取为-3σ，第二水平取为3σ，共得到12个采样点。计算每一个正交试验方案落点位置到期望落点的距离和落地时间，以每一试验方案到期望落点的距离和落地时间为试验指标进行极差分析，得到偏离期望落点最大的试验方案{x^（R）₀，p^（R）}和落地时间最晚的试验方案{x^（t）₀，p^（t）}，计算并存储相应的弹道数据。

步骤3：飞行管道半径和安全时间计算。利用Halton采样计算的N-2条弹道数据和正交试验分析得到的两条弹道数据，通过插值可得到任一高度状态量x^（ⁱ^）_h和最晚存在时刻t^（ⁱ^）_h，i=1，2，···，N。在最晚存在时刻t^（ⁱ^）_h之后，返回器被认为不会出现在h高度以上的空域，因此t^（ⁱ^）_h也被称作安全时间。以期望弹道在h高度的位置为圆心，作所有采样点在此高度上的最小包络圆，以此圆半径作为飞行管道半径。

重复以上步骤，可以计算出伞舱盖和减速伞等返回器分离物的期望弹道和飞行管道信息。在实际应用中，采样点个数N应根据硬件能力和系统更新频率要求进行限制，在本文中采样点个数设为500。

以返回器为例，判定搜救直升机是否位于远离返回器飞行管道的一般流程为：

1）输入当前时间和位置，根据当前高度分别计算该高度下返回器的安全时间以及飞行管道半径（安全半径）。

2）判定直升机是否在安全时间之后位于该高度，若是，直升机是安全的，否则转步骤3。步骤2称为时间安全判定。

3）计算得到直升机与返回器期望位置间的距离，判定该距离是否大于当前高度的管道半径，若是，则直升机是安全的，否则，当前位置存在风险，应沿返回器期望位置所在方向的相反向飞离。步骤3称为空间安全判定。

上述搜救直升机是否处于安全空域的判定流程如图3所示。重复上述过程可判定直升机是否位于伞舱盖和减速伞的安全空域。

3 返回器搜救任务不确定性分析

本节以CE-5返回器在不确定条件下的安全空域计算问题和搜救引导问题为例进行分析。安全空域计算结果确定了每一时刻的最低安全高度和返回器及其分离物的管道，为空中和地面搜救人员提供安全指引，同时给出返回器、伞舱盖以及减速伞的期望落点，引导搜救力量对返回器进行回收。CE-5返回器的主要标称参数如表1所示。利用探空气球测量的着陆场风场数据作为仿真风场输入计算返回器的飞行弹道，不同高度上的风速大小如图4所示。系统的不确定性参数均视为以名义值为期望值的正态分布，参考已有相关研究^[9]和CE-5返回器的自身特点，其3倍标准差的范围如表2所示。

图3安全空域的判定流程

Fig.3Determining process of the safe airspace

表1物伞系统参数名义值

Tab.1 Nominal value of the body-parachute system parameters

图4名义风场

Fig.4Nominal wind field

表2回收系统不确定参数及取值范围

Tab.2 Uncertain parameters and range of recovery system

不同采样方法得到的在地面处的飞行管道半径如图5所示，其中HS方法和MC方法的采样次数均为10万次，OD方法采用标准正交表L₁₂（2¹¹）进行采样，从图5中结果可见，HS方法得到的管道半径要明显大于MC方法得到的计算结果，正交设计采样得到的最大半径又要大于HS方法计算得到的管道半径，极差分析得到的正交设计最优方案（optimized plan of orthogonal design，OPOD）的管道半径又略有增大。通过对正交设计采样计算结果进行极差分析得到使飞行管道半径最大的随机参数取值方案，方案如表3所示，其中，ψ_w表示风向偏航角，V_w表示风速大小。在飞行管道半径的试验方案中，返回器的质量和速度都取最小值，其他因素均取最大值，这与返回器伞降过程的物理特性是一致的：返回器质量最小，初速度最小，风速和空气密度最大，能够使得返回器达到最大的风飘距离。

图5不同方法计算得到的飞行管道半径

Fig.5Fight pipeline radius calculated by different methods

期望弹道是返回器飞行管道计算的必要结果，由于更新时间的需要，期望弹道的计算量必须限制在一定的数量下。对比HS方法和MC方法计算期望落点东向坐标的收敛过程，如图6所示，从图中可见，HS方法的收敛速度要远优于MC方法，在采样次数小于1 000的限制下，MC方法计算的期望落点坐标仍处于不稳定的波动状态，而HS方法已经收敛至落点坐标25 m以内的状态。显然，HS方法能够更快地收敛至期望落点。

表3正交试验分析得到的偏离期望落点最大的方案

Tab.3 Scheme with the largest deviation from the expected landing point obtained by orthogonal test analysis

图6落点坐标期望值随采样次数增加的收敛曲线

Fig.6Convergence curve of the expected landing point coordinate with increasing of sampling times

在给定输入下计算得到的返回器、伞舱盖和减速伞的飞行管道分别如图7（a）~（c）所示。图中直角坐标系为原点在开伞点地面投影点的天东北坐标系，名义弹道指按名义参数计算得到的弹道。

图7返回器回收过程中的飞行管道

Fig.7Flight pipeline during the recovering of the reentry capsule

从飞行管道的计算结果可见，由于当地风场风向偏东，返回器、伞舱盖和减速伞均明显向东偏移。因为质量偏小，伞舱盖和减速伞受风场影响相比于返回器更明显，落地时刻伞舱盖和减速伞主伞包组合体的可能存在半径相比于返回器分别增大了25%和66%，计算结果反映了实际搜救任务中伞舱盖和减速伞难以找到的实际问题。通过对数据的处理，不同高度h下的返回器飞行管道半径R和安全时间T_s如表4所示，表中的安全时间从降落伞打开指令发出开始计时，飞行管道半径和安全时间的变化趋势如图8所示，伴随高度的增加，飞行管道半径和安全时间呈减小趋势，这与系统基本特性是相符的。

以搜救直升机在距离地面2 km高度工作为例，对返回器的搜救任务进行说明。一方面，根据计算结果，当搜救直升机在名义开伞指令下达后的489 s、551 s和987 s以后出现在距离地面2 km高度上，将几乎不可能与返回器、伞舱盖和减速伞发生空中碰撞。另一方面，在考虑随机因素扰动的情况下，返回器最晚会在名义开伞指令下达后的667 s落地。地面2 km高度上的飞行管道截面如图9所示，返回器期望落点包含于3个飞行管道之内，因此，要保证在返回器落地时，搜救直升机出现在落地点上空，搜救直升机必须承担一定与伞舱盖或减速伞相碰撞的风险。在保证搜救直升机不与返回器相撞的前提下，应尽可能合理地规划直升机飞行路线以降低搜救直升机与返回器分离物相撞的风险。

表4返回器在不同高度下的飞行管道半径和安全时间

Tab.4 Flight pipeline radius and safe time of the reentry capsule under different heights

图8不同高度下的飞行管道半径与安全时间

Fig.8Flight pipeline radius and safe time under different heights

图9地面2 km高度上的飞行管道截面

Fig.9Flight pipeline cross section at 2 km above ground

搜救直升机的搜索流程按如下步骤实施：①在开伞指令下达后的489 s之前，可由安全判定流程使搜救直升机位于2 km高度的飞行管道之外；②489 s以后，从返回器、伞舱盖和减速伞飞行管道最外边界选取一个最优出发点开始进行搜救。

为定量刻画搜救直升机与返回器分离物相撞的风险，定义碰撞风险指标为

\{\begin{matrix} f = \int_{s} (d_{2} + d_{3}) d s \\ d_{2} = \sqrt{{(y - y_{2})}^{2} + {(z - z_{2})}^{2}} \\ d_{3} = \sqrt{{(y - y_{3})}^{2} + {(z - z_{3})}^{2}} \end{matrix}

(17)

其中，s为飞行路线，y和z分别为飞行路线点在天东北坐标系下的东向和北向坐标，y₂和z₂为伞舱盖在2 km高度的期望坐标，y₃和z₃为减速伞在2 km高度的期望坐标。设搜救直升机的飞行速度为260 km/h，从出发点至期望落点，飞行路线为直线。通过计算得到从相对期望落点不同方向管道边界点进行搜救的风险指标数值和达到期望落点所需时间的曲线，如图10所示。最优进入点及搜索方向标注于图9中，飞行方向为北偏东1.6°，搜救出发点距离期望落点4.2 km，所需飞行时间为59 s。当搜救直升机达到期望落点时，名义状态下，返回器仍需24 s才能落地。相比于最优进入点，若搜救直升机从位于期望落点正东方向的管道边界点进行搜索，则对应的风险指标将增大116%，搜索所需时间将增大2.4倍，可见，搜索路线和搜索出发点的优化设计对于搜救任务的安全性和快速性均有较大的提升。

图10期望落点不同方向搜索出发点对应的风险指标和期望搜索时间

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

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

图1嫦娥五号返回器回收过程

Fig.1Recovery process of the reentry capsule of CE-5

下载: 全尺寸图片

图2坐标系系统

Fig.2Coordinate system

下载: 全尺寸图片

图3安全空域的判定流程

Fig.3Determining process of the safe airspace

下载: 全尺寸图片

图4名义风场

Fig.4Nominal wind field

下载: 全尺寸图片

图5不同方法计算得到的飞行管道半径

Fig.5Fight pipeline radius calculated by different methods

下载: 全尺寸图片

图6落点坐标期望值随采样次数增加的收敛曲线

Fig.6Convergence curve of the expected landing point coordinate with increasing of sampling times

下载: 全尺寸图片

图7返回器回收过程中的飞行管道

Fig.7Flight pipeline during the recovering of the reentry capsule

下载: 全尺寸图片

图8不同高度下的飞行管道半径与安全时间

Fig.8Flight pipeline radius and safe time under different heights

下载: 全尺寸图片

图9地面2 km高度上的飞行管道截面

Fig.9Flight pipeline cross section at 2 km above ground

下载: 全尺寸图片

图10期望落点不同方向搜索出发点对应的风险指标和期望搜索时间

Fig.10Risk index and expected search time corresponding to the departure point in different directions of the expected landing point

下载: 全尺寸图片

表1物伞系统参数名义值

下载: 全尺寸图片

表2回收系统不确定参数及取值范围

下载: 全尺寸图片

表3正交试验分析得到的偏离期望落点最大的方案

下载: 全尺寸图片

表4返回器在不同高度下的飞行管道半径和安全时间

下载: 全尺寸图片

图1嫦娥五号返回器回收过程

Fig.1Recovery process of the reentry capsule of CE-5

图2坐标系系统

Fig.2Coordinate system

图3安全空域的判定流程

Fig.3Determining process of the safe airspace

图4名义风场

Fig.4Nominal wind field

图5不同方法计算得到的飞行管道半径

Fig.5Fight pipeline radius calculated by different methods

图6落点坐标期望值随采样次数增加的收敛曲线

Fig.6Convergence curve of the expected landing point coordinate with increasing of sampling times

图7返回器回收过程中的飞行管道

Fig.7Flight pipeline during the recovering of the reentry capsule

图8不同高度下的飞行管道半径与安全时间

Fig.8Flight pipeline radius and safe time under different heights

图9地面2 km高度上的飞行管道截面

Fig.9Flight pipeline cross section at 2 km above ground

图10期望落点不同方向搜索出发点对应的风险指标和期望搜索时间

Fig.10Risk index and expected search time corresponding to the departure point in different directions of the expected landing point

表1物伞系统参数名义值

表2回收系统不确定参数及取值范围

表3正交试验分析得到的偏离期望落点最大的方案

表4返回器在不同高度下的飞行管道半径和安全时间

图1嫦娥五号返回器回收过程

Fig.1Recovery process of the reentry capsule of CE-5

图2坐标系系统

Fig.2Coordinate system

图3安全空域的判定流程

Fig.3Determining process of the safe airspace

图4名义风场

Fig.4Nominal wind field

图5不同方法计算得到的飞行管道半径

Fig.5Fight pipeline radius calculated by different methods

图6落点坐标期望值随采样次数增加的收敛曲线

Fig.6Convergence curve of the expected landing point coordinate with increasing of sampling times

图7返回器回收过程中的飞行管道

Fig.7Flight pipeline during the recovering of the reentry capsule

图8不同高度下的飞行管道半径与安全时间

Fig.8Flight pipeline radius and safe time under different heights

图9地面2 km高度上的飞行管道截面

Fig.9Flight pipeline cross section at 2 km above ground

图10期望落点不同方向搜索出发点对应的风险指标和期望搜索时间

Fig.10Risk index and expected search time corresponding to the departure point in different directions of the expected landing point

表1物伞系统参数名义值

表2回收系统不确定参数及取值范围

表3正交试验分析得到的偏离期望落点最大的方案

表4返回器在不同高度下的飞行管道半径和安全时间

刘伯阳, 唐学海, 孙中兴. 跳跃式返回航天器落点预报技术研究[J]. 南京航空航天大学学报,2021,53(增刊 1):45-50. LIU B Y, TANG X H, SUN Z X. Research on prediction of skip entry spacecraft′s landing site[J]. Journal of Nanjing University of Aeronautics & Astronautics,2021,53(Suppl 1):45-50.(in Chinese)

淡鹏, 姜宇, 李恒年. 航天器返回可达区域国土占比快速计算方法[J]. 航天返回与遥感,2022,43(3):25-32. DAN P, JIANG Y, LI H N. Rapid calculation method of the area proportion of spacecraft recovery footprint[J]. Spacecraft Recovery & Remote Sensing,2022,43(3):25-32.(in Chinese)

包进进, 刘大海, 荣伟, 等. 嫦娥五号探测器回收分系统研制与验证[J]. 航天器工程,2021,30(5):16-23. BAO J J, LIU D H, RONG W,et al. Development and validation of recovery system on Chang′e-5 spacecraft[J]. Spacecraft Engineering,2021,30(5):16-23.(in Chinese)

王立武, 高庆玉, 许望晶, 等. 嫦娥五号返回器降落伞回收动力学建模与分析[J]. 宇航学报,2021,42(8):1051-1056. WANG L W, GAO Q Y, XU W J,et al. Dynamic modeling and analysis of Chang′e-5 parachute deceleration system[J]. Journal of Astronautics,2021,42(8):1051-1056.(in Chinese)

MEREL M H, MULLIN F J. Analytic Monte Carlo error analysis[J]. Journal of Spacecraft and Rockets,1968,5(11):1304-1308.

ZHANG T F, SU H, GONG C L.hp-adaptive RPD based sequential convex programming for reentry trajectory optimization[J]. Aerospace Science and Technology,2022,130:107887.

张青斌, 唐乾刚, 彭勇, 等. 飞船返回舱降落伞系统动力学[M]. 北京: 国防工业出版社,2013. ZHANG Q B, TANG Q G, PENG Y,et al. Dynamics of parachute-capsule recovery system[M]. Beijing: National Defense Industry Press,2013.(in Chinese)

王慧娟. 火星进入过程的开伞控制方法研究[D]. 长沙: 国防科技大学,2015. WANG H J. Research on parachute opening triggering algorithms in Mars entry process[D]. Changsha: National University of Defense Technology,2015.(in Chinese)

SPENCER D A, BRAUN R D. Mars pathfinder atmospheric entry:trajectory design and dispersion analysis[J]. Journal of Spacecraft and Rockets,1996,33(5):670-676.

DESAI P N, BRAUN R D, POWELL R W,et al. Six-degree-of-freedom entry dispersion analysis for the METEOR recovery module[J]. Journal of Spacecraft and Rockets,1997,34(3):334-340.

BRAUN R D, MITCHELTREE R, CHEATWOOD F. Mars microprobe entry-to-impact analysis[J]. Journal of Spacecraft and Rockets,1999,36(3):412-420.

LEONARD A, ROGERS J, GERLACH A. Koopman operator approach to airdrop mission planning under uncertainty[J]. Journal of Guidance, Control,and Dynamics,2019,42(11):2382-2398.

李杰, 陈建兵. 随机动力系统中的概率密度演化方程及其研究进展[J]. 力学进展,2010,40(2):170-188. LI J, CHEN J B. Advances in the research on probability density evolution equations of stochastic dynamical systems[J]. Advances in Mechanics,2010,40(2):170-188.(in Chinese)

CHEN Z, YANG D D. Invariant measures and stochastic Liouville type theorem for non-autonomous stochastic reaction-diffusion equations[J]. Journal of Differential Equations,2023,353:225-267.

HALDER A, BHATTACHARYA R. Dispersion analysis in hypersonic flight during planetary entry using stochastic liouville equation[J]. Journal of Guidance, Control,and Dynamics,2011,34(2):459-474.

TRISOLINI M, COLOMBO C. Propagation and reconstruction of reentry uncertainties using continuity equation and simplicial interpolation[J]. Journal of Guidance, Control,and Dynamics,2021,44(4):793-811.

LEONARD A, KLEIN B, JUMONVILLE C,et al. Probabilistic algorithm for ballistic parachute transition altitude optimization[J]. Journal of Guidance, Control,and Dynamics,2017,40(12):3037-3049.

LEONARD A, ROGERS J, GERLACH A R,et al. A probabilistic approach to the optimal determination of a computed air release point[C]//Proceedings of the 24th AIAA Aerodynamic Decelerator Systems Technology Conference,2017.

KLEIN B, ROGERS J D. A probabilistic approach to unguided airdrop[C]//Proceedings of the 23rd AIAA Aerodynamic Decelerator Systems Technology Conference,2015.

LEONARD A. Probabilistic methods for decision making in precision airdrop[D]. Atlanta: Georgia Institute of Technology,2019.

LASOTA A, MACKEY M. Chaos,fractals and noise:stochastic aspects of dynamics[M].2nd ed. New York: Springer-Verlag,1994:37-49.

LEONARD A, ROGERS J, GERLACH A. Probabilistic release point optimization for airdrop with variable transition altitude[J]. Journal of Guidance, Control,and Dynamics,2020,43(8):1487-1497.

BLEDSOE K J, MORRIS A L. Development of the sasquatch drop test footprint tool[C]//Proceedings of the 21st AIAA Aerodynamic Decelerator Systems Technology Conference and Seminar,2011.

BLEDSOE K J, BERNATOVICH M A, GALVIN P J. Development and overview of CPAS sasquatch airdrop landing location predictor software[C]//Proceedings of the 23rd AIAA Aerodynamic Decelerator Systems Technology Conference,2015.

DONG G Y, LEMIEUX C. Dependence properties of scrambled Halton sequences[J]. Mathematics and Computers in Simulation,2022,200:240-262.

王静. 高强度正交表的构造及D-最优增广设计[D]. 上海: 上海师范大学,2021. WANG J. Construction of orthogonal arrays with high strength and D-optimal augmented designs[D]. Shanghai: Shanghai Normal University,2021.(in Chinese)

MANRIQUE J B. Advances in spacecraft atmospheric entry guidance[D]. Irvine: University of California, Irvine,2010.

贾贺, 荣伟, 包进进, 等. 探月返回器降落伞减速系统设计及试验验证[J]. 航天器工程,2020,29(4):26-33. JIA H, RONG W, BAO J J,et al. Design and verification of parachute deceleration system for lunar return capsule[J]. Spacecraft Engineering,2020,29(4):26-33.(in Chinese)

CRUZ J R, WAY D W, SHIDNER J D,et al. Parachute models used in the Mars science laboratory entry,descent,and landing simulation[C]//Proceedings of the AIAA Aerodynamic Decelerator Systems(ADS)Conference,2013.

ALLASIA G, CAVORETTO R, DE ROSSI A. Lobachevsky spline functions and interpolation to scattered data[J]. Computational and Applied Mathematics,2013,32(1):71-87.

KAZEM S, HATAM A. Scattered data interpolation:strictly positive definite radial basis/cardinal functions[J]. Journal of Computational and Applied Mathematics,2021,394:113580.

SHAMPINE L F, WATTS H A, DAVENPORT S M. Solving nonstiff ordinary differential equations:the state of the art[J]. SIAM Review,1976,18(3):376-411.

NIEDERREITER H. Random number generation and quasi-Monte Carlo methods[C]//Proceedings of CBMS-NSF Regional Conference Series in Applied Mathematics,1992.

首页

期刊介绍

投稿指南

编委会

出版声明

开放获取声明

联系我们

Email订阅

Rss

English

1 返回器降落伞回收阶段的动力学模型

2 返回器飞行管道计算方法

3 返回器搜救任务不确定性分析

4 结论