航天器轨道追逃态势分析的水平集方法

杨傅云翔; 杨乐平; 朱彦伟; 张乘铭; YANG Fuyunxiang; YANG Leping; ZHU Yanwei; ZHANG Chengming

航天器轨道追逃态势分析的水平集方法

doi: 10.11887/j.cn.202403004

杨傅云翔，杨乐平，朱彦伟，张乘铭

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

基金项目: 国防科技大学自主创新科学基金资助项目(22-ZZCX-083)

详细信息

作者简介

杨傅云翔（1996—），男，江苏镇江人，博士研究生，E-mail: yangfuyunxiang@alumni.nudt.edu.cn

通讯作者

朱彦伟（1981—），男，河南平舆人，教授，博士，硕士生导师，E-mail: zywnudt@163.com

中图分类号: V448

文献标识码: A

文章编号: 1001-2486(2024)03-030-09

Situation analysis method based on level set for spacecraft pursuit-evasion game

YANG Fuyunxiang ， YANG Leping ， ZHU Yanwei ， ZHANG Chengming

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

摘要

航天器轨道追逃是当前航天动力学与控制领域的研究热点。针对轨道追逃定性问题开展研究，提出了一种综合运用降维动力学模型和后向可达集的航天器近距离轨道追逃态势分析方法，以支撑任务可行性分析；通过在视线旋转坐标系下推导博弈系统降维动力学模型，建立近距离追逃定性问题模型，减少了状态空间维度；使用目标集的后向可达集描述捕获区并划分追逃状态空间，基于水平集方法建立可达集在降维动力学模型中演化的动态HJI(Hamilton-Jacobi-Isaacs)偏微分方程，并设计WENO-TVD求解器数值计算HJI方程终值问题粘性解，完成了追逃目标集的准确描述并避免了可能出现的终端奇异现象。通过不同推力构型的追逃场景数值仿真验证了方法的有效性，展现了一次计算批量化处理初始态势的功能。

关键词

航天器追逃博弈 / 可行性分析 / 水平集方法

Abstract

Spacecraft pursuit-evasion game is currently a research hotspot in the field of aerospace dynamics and control. Qualitative spacecraft pursuit-evasion game was studied in order to provide feasibility support for strategy design, and a situation analysis method for scenarios in close range was proposed by comprehensively using dimension-reduction dynamics and backwards reachability set. A dimension-reduction dynamic model was derived in the line-of-sight rotation coordinate system, and the pursuit-evasion qualitative model of the game system was determined to reduce state space dimension. The backwards reachable set of the target set was used to divide the pursuit-evasion state space and describe the capture zone. A time-dependent HJI (Hamilton-Jacobi-Isaacs) PDE (partial differential equation) was established to describe the evolution of backwards reachable set in the dimension-reduction dynamics, based on level set method, and a WENO-TVD(weighted essentially non-oscillatory-total variation diminishing) solver was designed to numerically calculate the viscous solution of the final value problem of the HJI PDE. These measures achieve the accurate description of pursuit-evasion target set and avoid the possible terminal singularity. The effectiveness of the method was demonstrated by numerical simulations of several pursuit-evasion scenarios with different thrust configurations, and the function of batch processing of initial situations in a single calculation was demonstrated.

Keywords

spacecraft pursuit-evasion game / feasibility analysis / level set method

1 问题建模 1.1 视线运动方程 1.2 降维动力学模型 1.3 定性微分对策 2 求解方法 2.1 问题重构 2.2 方程求解 2.2.1 网格生成 2.2.2 空间离散 2.2.3 时间离散 2.3 态势评估 3 数值仿真 3.1 推力构型Ⅰ 3.2 推力构型Ⅱ 4 结论

近年来，航天器接近与交会操作技术的发展使得轨道追逃问题成为航天动力学与控制领域的研究热点。该问题实质上是一个包含双边控制的连续动态对抗问题，其中被称为追踪方（pursuer，P）与逃逸方（evader，E）的局中人具有相互冲突的目标^[1]。

整个博弈过程可以用零和微分对策问题进行描述，并根据有无具体的支付泛函，可以进一步细分为定量问题、定性问题两类。其中，定量问题关注于最优化具体支付泛函，求解均衡点，规划局中人最优轨迹；而定性问题只研究可行性，即判断某种博弈结局能否实现。当前航天器轨道追逃的研究成果主要集中于定量问题领域，形成了诸如配点法、进化算法等一系列计算方法^[1]，并借助近期活跃的人工智能技术，产生了一批机器学习方法^[2]，为具体策略的设计提供了辅助。而作为策略设计的前提，需要判定当前博弈态势导向某种特定结局是否可行，这便需要研究定性问题。

但相比于定量问题，定性问题的研究要缺乏得多。求解定性微分对策问题的关键在于划分状态空间为捕获区与逃逸区，这通常意味着需要构造复杂界栅，轨道动力学模型的加入进一步增加了计算难度。目前的研究成果基本都是针对简易场景使用Isaacs建立的界栅理论^[3]分析获得的：Anderson等^[4]对低轨航天器轨道追逃问题进行了研究，通过将运动方程在圆参考轨道处线性化，得到了平面追逃界栅的解析解；张秋华等^[5-6]等对近地共面轨道上的航天器追逃问题进行了分析，研究了视线坐标系下局中人的最优控制策略和界栅构造方法，以及径向连续可变小推力假设下以轨道根数为状态变量的界栅解析表达；Hafer等^[7]通过转化定量问题数值求解三维界栅并明确捕获区、逃逸区；祝海^[8]基于CW（Clohessy-Witshire）方程进行分析，将航天器控制量近似为剩余捕获时间的高阶多项式，从而得到界栅的近似解析表达。这类分析方法统称Isaacs方法，基于状态回转的思想，以目标集边界的可用部分边界作为初始曲线，积分倒退轨迹方程（retrograde path equation，RPE）求解界栅以划分状态空间为捕获区与逃逸区，从而确定博弈的捕获（逃逸）条件。

然而，随着追逃问题的复杂度不断提升，Isaacs方法也越来越难以奏效，主要源于以下三点：一是状态空间维度不断升高导致RPE越来越难以解析求解；二是计算过程中会出现终端协态奇异现象；三是博弈结局分析时大量测试迭代会带来计算负担。

由于定性微分对策等价于可达性问题^[9]，这为研究定性问题提供了新的研究思路。Mitchell^[10]从可达分析的角度出发，使用水平集方法（level set method，LSM），将定性微分对策重构为求解动态哈密顿-雅可比-埃萨克斯（Hamilton-Jacobi-Isaacs，HJI）偏微分方程（partial differential equation，PDE），成功处理一类飞行器碰撞避免问题。目前Venigalla等^[11-12]尝试借鉴该研究路线对航天器追逃态势进行分析，其针对脉冲推进场景，以轨道根数为自变量，顺时间流向计算个体航天器可达集，通过两方前向可达集的交集分析当前态势并为最优轨迹规划提供依据。但该方法存在以下两个问题：第一，多数情况下局中人个体可达集交集并不能描述微分对策的目标集，结论参考价值有限；第二，初值敏感，不同的初始状态就需要重新计算，无法批量处理场景。

针对上述问题，本论文综合运用降维动力学模型、可达分析理论和水平集方法，创造性提出了一种高效的航天器轨道追逃态势分析框架。对追逃场景而非个体航天器进行建模，基于目标集生成计算结果；同时考虑到虽然场景初始状态易变，但目标集由任务需求确定一般固定，故从目标集出发逆时间流向分析问题，计算后向可达集（backwards reachable set，BRS）来划分初始状态空间，做到批量分析。

1 问题建模

在研究航天器近距离追逃问题时，通常会以虚拟的参考航天器为原点建立主星轨道LVLH（local vertical，local horizontal）坐标系，使用CW方程描述轨道追逃的动力学特征。一方面线性特性增加了求解方法的选择范围，但另一方面也导致追逃变成一个至少6维的高维问题，增加了求解的难度。为了降低问题的维度，引入视线旋转坐标系。

1.1 视线运动方程

本小节将推导视线旋转坐标系下的视线运动方程，为建立描述近距离追逃的降维动力学模型提供基础。惯性空间中，从P质心指向E质心的单位矢量称为视线e_r。追踪航天器与目标航天器在地心惯性（earth centrical inertial， ECI）坐标系中的相对位置关系如图1所示。

图1追逃双方相对位置关系

Fig.1Relative position relationship between pursuer and evader

图1中：

r = r_{E} - r_{P} = r e_{r}

(1)

式中，r为追逃航天器之间相对距离。

航天器之间相对位置的变化会导致视线在惯性空间中发生旋转，定义视线的瞬时旋转角速度为：

ω_{s} = ω_{s} e_{ω}

(2)

式中，ω_s为视线的瞬时转率。以追踪航天器质心为原点，e_r、 e_ω以及

e_{θ} ≜ e_{ω} \times e_{r}

组成视线旋转坐标系。

任意时刻，视线的瞬时转率ω_s与方向e_ω都可能发生变化。在定义e_ω的旋转角速度为

Ω_{s} = Ω_{s} e_{r}

(3)

之后，可以获得（e_r，e_ω，e_θ）的运动方程如下：

\frac{d e_{r}}{d t} = ω_{s} e_{θ}

(4)

\frac{d e_{θ}}{d t} = - ω_{s} e_{r} + Ω_{s} e_{ω}

(5)

\frac{d e_{ω}}{d t} = - Ω_{s} e_{θ}

(6)

1.2 降维动力学模型

将视线运动方程与航天器相对运动方程结合起来，可以得到视线旋转坐标系下的航天器近距离追逃动力学模型。

对式（1）求导：

\dot{r} = \dot{r} e_{r} + r \frac{d e_{r}}{d t} = v = v_{E} - v_{P}

(7)

其中，v_P、 v_E分别表示追逃两方的速度矢量，v表示双方的相对速度矢量。将式（4）代入式（7）可得：

\dot{r} e_{r} + r ω_{s} e_{θ} = v = v_{E} - v_{P}

(8)

对上式求导，并将式（4）与式（5）代入可得：

a_{E} - a_{P} = (\ddot{r} - r ω_{s}^{2}) e_{r} + (r \dot{ω_{s}} + 2 \dot{r} ω_{s}) e_{θ} + r ω_{s} Ω_{s} e_{ω}

(9)

其中，加速度a_E、 a_P分为引力项与推力项两部分。考虑到近距离追逃情况下，相对距离r远小于||r_E||、||r_P||，故忽略引力差的影响，同时定义

v_{r} ≜ \dot{r} ， v_{θ} ≜ r ω_{s}

，可得：

u_{E} - u_{P} = ({\dot{v}}_{r} - \frac{v_{θ}^{2}}{r}) e_{r} + ({\dot{v}}_{θ} + \frac{v_{r} v_{θ}}{r}) e_{θ} + v_{θ} Ω_{s} e_{ω}

(10)

其中，u_E、 u_P表示局中人推力控制量。

观察式（10）可知，e_ω轴向的加速度只改变Ω_s与v_θ两个速度量，对相对距离r没有影响。同时考虑到追逃问题一般关注相对距离变化，目标集约束只包含终端相对距离，故追逃问题的状态空间维度可以减少至3维。降维动力学模型如下：

\{\begin{matrix} \dot{r} = v_{r} \\ {\dot{v}}_{r} - \frac{v_{θ}^{2}}{r} = u_{E}^{r} - u_{P}^{r} \\ {\dot{v}}_{θ} + \frac{v_{r} v_{θ}}{r} = u_{E}^{θ} - u_{P}^{θ} \end{matrix}

(11)

再引入r_ref、v_ref对上式进行去量纲化处理，记：

\{\begin{matrix} x_{r} ≜ \frac{r}{r_{r e f}} \\ x_{v_{r}} ≜ \frac{v_{r}}{v_{r e f}} \\ x_{v_{θ}} ≜ \frac{v_{θ}}{v_{r e f}} \\ \hat{t} ≜ \frac{v_{r e f}}{r_{r e f}} t \end{matrix}

(12)

可以得到去量纲化的降维动力学模型，即：

\{\begin{matrix} \frac{d x_{r}}{d \hat{t}} = x_{v_{r}} \\ \frac{d x_{v_{r}}}{d \hat{t}} = \frac{x_{v_{θ}}^{2}}{x_{r}} + \frac{r_{r e f}}{v_{r e f}^{2}} (u_{E}^{r} - u_{P}^{r}) \\ \frac{d x_{v_{θ}}}{d \hat{t}} = - \frac{x_{v_{r}} x_{v_{θ}}}{x_{r}} + \frac{r_{r e f}}{v_{r e f}^{2}} (u_{E}^{θ} - u_{P}^{θ}) \end{matrix}

(13)

1.3 定性微分对策

描述航天器追逃的定性微分对策问题通常分为动力学模型和目标集两部分。动力学模型一般使用常微分方程描述如下：

\{\begin{matrix} \frac{d x}{d t} = \dot{x} = f (x, u_{E}, u_{P}) \\ x (t_{0}) = x_{0} \end{matrix}

(14)

其中，t₀≤0表示初始时刻，

x \in R^{n}

表示状态量，

u_{E} \in U_{E} \subseteq R^{n_{E}}

是逃逸方控制量，

u_{P} \in U_{P} \subseteq R^{n_{P}}

为追踪方控制量。追逃动力学模型f:

R^{n} \times U_{E} \times U_{P} \to R^{n}

有界且Lipschitz连续。

目标集

G_{0} \subseteq R^{n}

一般为闭集，并且可以表示成一个有界且Lipschitz连续的函数g:

R^{n}

→

R^{}

的零水平子集，即：

G_{0} = \{x (0) \in R^{n} ∣ g (x (0)) ⩽ 0\}

(15)

其中，x（0）表示终端状态，追踪方通过选择

u_{p} \in U_{p} \subseteq R^{n_{p}}

力图实现g（x（0））≤0，从而达到捕获逃逸方的目的。但是，逃逸方则选择

u_{E} \in U_{E} \subseteq R^{n_{E}}

，使得g（x（0））＞0，从而躲避追踪方的追捕。

2 求解方法

为了弥补传统方法中局中人个体可达集交集并不能准确描述博弈目标集以及无法批量处理场景的缺陷，本文从场景而非局中人出发，建立以博弈目标集逆时间演化的BRS，通过初始状态与可达集的关系进行博弈态势分析。该方法主要分为问题重构、方程求解以及态势评估三部分，整体流程如图2所示，其中箭头表示数据流向，箭头上方标签表示数据名称，方框表示执行操作。

图2方法流程

Fig.2Method process

2.1 问题重构

在明确局中人双方控制量的情况下，总存在一条特定轨迹

ξ_{f} (s; x_{0}, t_{0}, u_{E}, u_{P}) : [t_{0}, 0] \to R^{n}

(16)

作为式（14）的解。s∈[t₀，0]表示时间变量，分号区分变量与参数，其余为表示轨迹特征的参数。

定性分析过程中需要使用BRS^[9]描述捕获区，从而划分状态空间，如图3所示。BRS G（τ）定义如下：

\begin{matrix} G (τ) ≜ \{x \in R^{n} ∣ \exists u_{P} \in U_{P}, \forall u_{E} \in U_{E}, \\ ξ_{f} (0; x, t_{0}, u_{E}, u_{P}) \in G_{0}\} \end{matrix}

(17)

其中，τ=-t₀表示倒退时间。

图3后向可达集

Fig.3Backwards reachable set

使用LSM^[13]，定义水平集函数φ（x，t）:

R^{n} \times [t_{0} ， 0] \to R

为：

ϕ (x, t) = \underset{u_{P} \in U_{P}}{i n f} \underset{u_{E} \in U_{E}}{s u p} g (ξ_{f} (0; x_{0}, t_{0}, u_{E}, u_{P}))

(18)

式中，inf和sup分别表示一个集合最大的下界和最小的上界。

G（τ）可以表示为φ（x，t）的零水平子集：

G (τ) ≜ \{x \in R^{n} ∣ ϕ (x, t) ⩽ 0\}

(19)

且φ（x，t）自身随时间的演化规律可以通过求解如下HJI PDE的粘性解获得：

\{\begin{matrix} \frac{\partial ϕ (x, t)}{\partial t} + m i n [0, H (x, p)] = 0 \\ ϕ (x, 0) = g (x) \end{matrix}

(20)

其中，

H (x, p) = \underset{u_{E} \in U_{E}}{m a x} \underset{u_{P} \in U_{P}}{m i n} p^{T} f (x, u_{E}, u_{P})

(21)

表示哈密顿函数，

p = \nabla_{x} ϕ

表示空间梯度。

2.2 方程求解

2.2.1 网格生成

在综合考虑降维动力学模型的适用范围以及计算消耗之后，将空间计算域限制在一个有限的区域内，并在此区域内生成均匀笛卡儿网格^[14]，各空间维度的步长大小一致，均为Δx。另外在设置时间步长Δt时，为了保证计算稳定性，需要满足柯朗-弗里德里希斯-列维（Courant-Friedrichs-Lewy，CFL）条件:

| Δ t | (\frac{m a x {| f |}}{Δ x}) = η

(22)

其中，CFL数η需要在（0，1）之间进行选取，以保证方程收敛^[15]。

2.2.2 空间离散

水平集函数空间梯度p使用五阶加权本质无振荡（weighted essentially non-oscillatory，WENO）格式近似。以第j维偏导

ϕ_{x_{j}} ≜ \frac{\partial ϕ}{\partial x_{j}}

为例进行说明，定义：

\{\begin{matrix} D_{i}^{-} ϕ_{x_{j}} = D_{i - 1 / 2}^{1} ϕ_{x_{j}} = \frac{ϕ_{i} - ϕ_{i - 1}}{Δ x} \\ D_{i}^{+} ϕ_{x_{j}} = D_{i + 1 / 2}^{1} ϕ_{x_{j}} = \frac{ϕ_{i + 1} - ϕ_{i}}{Δ x} \end{matrix}

(23)

取{φ_i_-3，φ_i_-2，φ_i_-1，φ_i，φ_i₊₁，φ_i₊₂}，计算

\{D_{i - 2}^{-} ϕ_{x_{j}} ， D_{i - 1}^{-} ϕ_{x_{j}} ， D_{i}^{-} ϕ_{x_{j}} ， D_{i + 1}^{-} ϕ_{x_{j}} ， D_{i + 2}^{-} ϕ_{x_{j}}\}

作为插值模板，根据插值格式可得式（24）。

\{\begin{matrix} {(D_{i}^{1} ϕ_{x_{j}})}^{0} = \frac{1}{3} D_{i - 2}^{-} ϕ_{x_{j}} - \frac{7}{6} D_{i - 1}^{-} ϕ_{x_{j}} + \frac{11}{6} D_{i}^{-} ϕ_{x_{j}}, I_{0} = \{D_{i - 2}^{-} ϕ_{x_{j}}, D_{i - 1}^{-} ϕ_{x_{j}}, D_{i}^{-} ϕ_{x_{j}}\} \\ {(D_{i}^{1} ϕ_{x_{j}})}^{1} = - \frac{1}{6} D_{i - 1}^{-} ϕ_{x_{j}} + \frac{5}{6} D_{i}^{-} ϕ_{x_{j}} + \frac{1}{3} D_{i + 1}^{-} ϕ_{x_{j}}, I_{1} = \{D_{i - 1}^{-} ϕ_{x_{j}}, D_{i}^{-} ϕ_{x_{j}}, D_{i + 1}^{-} ϕ_{x_{j}}\} \\ {(D_{i}^{1} ϕ_{x_{j}})}^{2} = \frac{1}{3} D_{i}^{-} ϕ_{x_{j}} + \frac{5}{6} D_{i + 1}^{-} ϕ_{x_{j}} - \frac{1}{6} D_{i + 2}^{-} ϕ_{x_{j}}, I_{2} = \{D_{i}^{-} ϕ_{x_{j}}, D_{i + 1}^{-} ϕ_{x_{j}}, D_{i + 2}^{-} ϕ_{x_{j}}\} \end{matrix}

(24)

其中，I₀、I₁、I₂表示模板。

ϕ_{x_{j}}

在i处的左近似

{(ϕ_{x_{j}}^{-})}_{i}

可以写成

{(D_{i}^{1} ϕ_{x_{j}})}^{0}

、

{(D_{i}^{1} ϕ_{x_{j}})}^{1}

、

{(D_{i}^{1} ϕ_{x_{j}})}^{2}

三者的凸组合：

{(ϕ_{x_{j}}^{-})}_{i} = ω_{0} {(D_{i}^{1} ϕ_{x_{j}})}^{0} + ω_{1} {(D_{i}^{1} ϕ_{x_{j}})}^{1} + ω_{2} {(D_{i}^{1} ϕ_{x_{j}})}^{2}

(25)

非线性权重ω_k（k=0，1，2）为：

\{\begin{matrix} ω_{k} = \frac{α_{k}}{\sum_{l = 0}^{2} α_{l}} \\ α_{k} = \frac{c_{k}}{{(ε + β_{k})}^{2}} \end{matrix}

(26)

其中：ε为灵敏度参数，是一个小量，用于避免分母为零，其取值一定程度上影响收敛精度，通常取为ε=10^-6；^[16]c_k为理想权重，为确保在光滑区域能够获得最优近似，取值分别为

c_{0} = \frac{1}{10} ， c_{1} = \frac{3}{5} ， c_{2} = \frac{3}{10}

；β_k为光滑因子，表达式为

\{\begin{matrix} β_{0} = \frac{13}{12} {(D_{i - 2}^{-} ϕ_{x_{j}} - 2 D_{i - 1}^{-} ϕ_{x_{j}} + D_{i}^{-} ϕ_{x_{j}})}^{2} + \\ \frac{1}{4} {(D_{i - 2}^{-} ϕ_{x_{j}} - 4 D_{i - 1}^{-} ϕ_{x_{j}} + 3 D_{i}^{-} ϕ_{x_{j}})}^{2} \\ β_{1} = \frac{13}{12} {(D_{i - 1}^{-} ϕ_{x_{j}} - 2 D_{i}^{-} ϕ_{x_{j}} + D_{i + 1}^{-} ϕ_{x_{j}})}^{2} + \\ \frac{1}{4} {(D_{i - 1}^{-} ϕ_{x_{j}} - D_{i + 1}^{-} ϕ_{x_{j}})}^{2} \\ β_{2} = \frac{13}{12} {(D_{i}^{-} ϕ_{x_{j}} - 2 D_{i + 1}^{-} ϕ_{x_{j}} + D_{i + 2}^{-} ϕ_{x_{j}})}^{2} + \\ \frac{1}{4} {(3 D_{i}^{-} ϕ_{x_{j}} - 4 D_{i + 1}^{-} ϕ_{x_{j}} + D_{i + 2}^{-} ϕ_{x_{j}})}^{2} \end{matrix}

(27)

取{φ_i_-2，φ_i_-1，φ_i，φ_i₊₁，φ_i₊₂，φ_i₊₃}，以

{D_{i - 2}^{+} ϕ_{x_{j}} {， D}_{i - 1}^{+} ϕ_{x_{j}} ， D_{i}^{+} ϕ_{x_{j}} ， D_{i + 1}^{+} ϕ_{x_{j}} ， D_{i + 2}^{+} ϕ_{x_{j}}}

作为插值模板，重复上述计算流程便可以获得

{(ϕ_{x_{j}}^{+})}_{i}

。之后将左右导数组成空间梯度的左右近似p^-、 p⁺：

p^{-} = [ϕ_{x_{1}}^{-}; ϕ_{x_{2}}^{-}; \dots; ϕ_{x_{j}}^{-}; \dots; ϕ_{x_{n}}^{-}]

(28)

p^{+} = [ϕ_{x_{1}}^{+}; ϕ_{x_{2}}^{+}; \dots; ϕ_{x_{j}}^{+}; \dots; ϕ_{x_{n}}^{+}]

(29)

为了避免解的振荡，使用Lax-Friedrichs格式^[17]近似哈密顿量，即

H (x, p^{+}, p^{-}) = H (x, \frac{p^{+} + p^{-}}{2}) - \frac{1}{2} λ^{T} (p^{+} - p^{-})

(30)

其中，第j维耗散系数λ_j为

λ_{j} = \underset{p_{j} \in [p_{jmin}, p_{jmax}]}{m a x} |\frac{\partial H}{\partial p_{j}}|

(31)

且

p_{j m i n} = m i n (p_{j}^{+} ， p_{j}^{-}) ， p_{j m a x} = m a x (p_{j}^{+} ， p_{j}^{-})

。

在完成空间梯度p以及哈密顿量的近似之后，动态HJI PDE就退化为常微分方程：

\frac{d ϕ (x, t)}{d t} = L (ϕ (x, t))

(32)

其中，L（·）为空间离散算子。

2.2.3 时间离散

时间离散使用三阶显式全局总变差减少（total variation diminishing， TVD）Runge-Kutta格式^[18]，假设第l个时间层上的φ（x，t）值φ^l已知，则：

ϕ^{l + 1} = \frac{1}{3} ϕ^{l} + \frac{2}{3} h_{1} + \frac{2}{3} Δ t \cdot L (h_{2})

(33)

其中，

\{\begin{matrix} h_{1} = ϕ^{l} + Δ t \cdot L (ϕ^{l}) \\ h_{2} = \frac{3}{4} ϕ^{l} + \frac{1}{4} h_{1} + \frac{1}{4} Δ t \cdot L (h_{1}) \end{matrix}

(34)

2.3 态势评估

在水平集函数φ（x，t）计算完成后，根据式（19）确定最大任务时间t_max的后向可达集G（τ_max），完成对状态空间的划分。根据初始态势x₀在状态空间中的位置进行判断：若x₀∈G（τ_max），则说明初始态势x₀位于捕获区内，追踪方可以在t_max内完成拦截任务，该场景追踪方占优；若

x_{0} \notin G (τ_{m a x})

，则说明初始态势x₀位于逃逸区内，逃逸方可以在t_max内完成逃逸任务，该场景逃逸方占优。在完成态势评估、确定博弈导向特定结局的可行性之后，问题从定性转向定量，博弈双方便进入根据具体指标设计优化策略的阶段。

3 数值仿真

本节将通过分析不同推力构型下的多组航天器追逃场景验证方法的可行性。

设置最大任务时间t_max=3 h，去量纲参数 r_ref与v_ref分别为1 000 m以及100 m/s后，配置10组追逃场景，初始相对状态如表1所示。

表1追逃场景初始相对状态

Tab.1 Initial relative state of pursuit-evasion scenarios

追逃目标集为：

G_{0} = \{x \in R^{n} ∣ x_{r} ⩽ 1\}

(35)

即：

g (x) = x_{r} - 1

(36)

场景计算区域网格信息如表2所示，CFL数η=0.5。

表2视线旋转坐标系下追逃系统状态空间网格

Tab.2 State space grid of the pursuit-evasion system in the line-of-sight rotation coordinate system

3.1 推力构型Ⅰ

此构型中，推力满足如下约束：

\{\begin{matrix} u_{i}^{r} = T_{i} c o s α_{i} \\ u_{i}^{θ} = T_{i} s i n α_{i} \end{matrix} i = E, P

(37)

其中：α_i∈[0，2π]表示推力角；T_E=0.000 4g、T_P=0.001g分别表示逃逸方、追踪方最大推进加速度，g=9.78m/s²表示水平面重力加速度。

将式（13）、式（37）代入式（21）可知最优控制量：

\{\begin{matrix} c o s α_{E}^{*} = \frac{p_{v_{r}}}{\sqrt{p_{v_{r}}^{2} + p_{v_{θ}}^{2}}} = c o s α_{P}^{*} \\ s i n α_{E}^{*} = \frac{p_{v_{θ}}}{\sqrt{p_{v_{r}}^{2} + p_{v_{θ}}^{2}}} = s i n α_{P}^{*} \end{matrix}

(38)

哈密顿函数：

\begin{matrix} H (x, p) = p_{r} x_{v_{r}} + p_{v_{r}} \frac{x_{v_{θ}}^{2}}{x_{r}} - p_{v_{θ}} \frac{x_{v_{r}} x_{v_{θ}}}{x_{r}} + \\ 0.00006 g \sqrt{p_{v_{r}}^{2} + p_{v_{θ}}^{2}} \end{matrix}

(39)

图4展示了计算过程中可达集逆向演化的整个过程，其中红色曲面表示捕获区边界。当τ=5 400 s时，可达集开始出现第二个区域，并在之后的演化过程中不断扩大。随着逆向演化的进行，计算域中有越来越多的区域水平集函数值由正变负，即从逃逸区转变为捕获区。对应到现实场景则表示随着最大任务时间的增加，追踪方的可选初始态势增加，拦截任务的条件逐渐放宽，任务难度降低；对于逃逸方而言则正好相反，随着逃逸区范围的缩小，逃逸任务的条件逐渐严格，任务难度上升。图5给出了G（τ_max）对初始状态空间的划分结果。10组初始场景中，场景1和8位于捕获区内，捕获任务具有可行性，其他场景均位于逃逸区。说明场景1和8均为追踪方占优，而其余场景均为逃逸方占优。

图4后向可达集演化过程（推力构型Ⅰ）

Fig.4Evolutionary process of backwards reachable set (thrust configuration Ⅰ)

图5场景分析结果（推力构型Ⅰ）

Fig.5Analysis results of scenarios (thrust configuration Ⅰ)

3.2 推力构型Ⅱ

此构型中，推力满足如下约束：

\{\begin{matrix} u_{i}^{r} \in [- T_{i}^{r}, T_{i}^{r}] \\ u_{i}^{θ} \in [- T_{i}^{θ}, T_{i}^{θ}] \end{matrix} i = E, P

(40)

其中：T^r_E=0.000 4g、T^θ_E=0.000 3g，T^r_P=0.001g、T^θ_P=0.000 9g分别表示逃逸方、追踪方相应方向的最大推进加速度。

则将式（13）、式（40）代入式（21）可知最优控制量为：

\begin{matrix} {(u_{i}^{r})}^{*} = T_{i}^{r} s g n (p_{v_{r}}) \\ {(u_{i}^{θ})}^{*} = T_{i}^{θ} s g n (p_{v_{θ}}) \end{matrix}

(41)

哈密顿函数为：

\begin{matrix} H (x, p) = p_{r} x_{v_{r}} + p_{v_{r}} \frac{x_{v_{θ}}^{2}}{x_{r}} - p_{v_{θ}} \frac{x_{v_{r}} x_{v_{θ}}}{x_{r}} + \\ 0.00006 g |p_{v_{r}}| - 0.00006 g |p_{v_{θ}}| \end{matrix}

(42)

图6展示了计算过程中可达集逆向演化的整个过程，其中红色曲面表示捕获区边界。当 τ=4 050 s时，可达集开始出现第二个区域，并在之后的演化过程中不断扩大；在τ=9 450 s时两个可达区域又合并成一个连通区域。与推力构型Ⅰ场景得到的结论一致，随着逆向演化的进行，拦截任务的任务难度降低，逃逸任务难度上升。图7给出了G（τ_max）对初始状态空间的划分结果以及对10组初始场景的分析结果，其中场景1、3、4、8、10位于捕获区内，捕获任务具有可行性，其他场景均位于逃逸区。

图6后向可达集演化过程（推力构型Ⅱ）

Fig.6Evolutionary process of backwards reachable set (thrust configuration Ⅱ)

图7场景分析结果（推力构型Ⅱ）

Fig.7Analysis results of scenarios (thrust configuration Ⅱ)

4 结论

本文关注于解决航天器追逃定性问题，针对同步轨道近距离博弈场景，提出了基于水平集方法计算后向可达集，分析博弈态势以及确定特定追逃结局能否实现的研究方法。主要研究结果如下：

1）通过状态与可达集的关系，给出了一种直观判断空间博弈态势优劣的处理方法，可用于态势分析、任务可行性判断，为具体轨迹规划任务提供支撑；

2）基于视线旋转坐标系推导降维动力学模型，实现了追逃系统状态空间的降维，有效降低了问题的复杂度；

3）将定性微分对策问题转化为可达分析问题，再通过水平集方法重构为HJI PDE终值问题，避免了传统方法中出现的终端协态奇异问题，并可以实现状态空间划分演化过程的可视化；

4）基于目标集逆时间流向分析，可以做到批量处理初始态势，不再需要具体场景具体计算。

但该方法还存在以下局限：

1）降维动力学的推导是基于追逃目标集不包含Ω_s进行的，当目标集中包含Ω_s约束时，无法实现降维；

2）因为该方法需要划分网格，故分析区域不宜过大，且由于使用的是均匀笛卡儿网格，目前只适用于规则区域。

今后将针对方法局限性开展以下研究：

1）针对复杂目标集与高维动力学模型进行研究，尝试总结航天器追逃过程中的一般性规律；

2）为计算捕获区长时间的演化过程寻找稳定可靠的无网格方法；

3）与深度学习相结合，尝试使用物理信息神经网络（physics informed neural network，PINN）计算HJI PDE终值问题。

图1追逃双方相对位置关系

Fig.1Relative position relationship between pursuer and evader

下载: 全尺寸图片

图2方法流程

Fig.2Method process

下载: 全尺寸图片

图3后向可达集

Fig.3Backwards reachable set

下载: 全尺寸图片

图4后向可达集演化过程（推力构型Ⅰ）

Fig.4Evolutionary process of backwards reachable set (thrust configuration Ⅰ)

下载: 全尺寸图片

图5场景分析结果（推力构型Ⅰ）

Fig.5Analysis results of scenarios (thrust configuration Ⅰ)

下载: 全尺寸图片

图6后向可达集演化过程（推力构型Ⅱ）

Fig.6Evolutionary process of backwards reachable set (thrust configuration Ⅱ)

下载: 全尺寸图片

图7场景分析结果（推力构型Ⅱ）

Fig.7Analysis results of scenarios (thrust configuration Ⅱ)

下载: 全尺寸图片

表1追逃场景初始相对状态

下载: 全尺寸图片

表2视线旋转坐标系下追逃系统状态空间网格

下载: 全尺寸图片

图1追逃双方相对位置关系

Fig.1Relative position relationship between pursuer and evader

图2方法流程

Fig.2Method process

图3后向可达集

Fig.3Backwards reachable set

图4后向可达集演化过程（推力构型Ⅰ）

Fig.4Evolutionary process of backwards reachable set (thrust configuration Ⅰ)

图5场景分析结果（推力构型Ⅰ）

Fig.5Analysis results of scenarios (thrust configuration Ⅰ)

图6后向可达集演化过程（推力构型Ⅱ）

Fig.6Evolutionary process of backwards reachable set (thrust configuration Ⅱ)

图7场景分析结果（推力构型Ⅱ）

Fig.7Analysis results of scenarios (thrust configuration Ⅱ)

表1追逃场景初始相对状态

表2视线旋转坐标系下追逃系统状态空间网格

图1追逃双方相对位置关系

Fig.1Relative position relationship between pursuer and evader

图2方法流程

Fig.2Method process

图3后向可达集

Fig.3Backwards reachable set

图4后向可达集演化过程（推力构型Ⅰ）

Fig.4Evolutionary process of backwards reachable set (thrust configuration Ⅰ)

图5场景分析结果（推力构型Ⅰ）

Fig.5Analysis results of scenarios (thrust configuration Ⅰ)

图6后向可达集演化过程（推力构型Ⅱ）

Fig.6Evolutionary process of backwards reachable set (thrust configuration Ⅱ)

图7场景分析结果（推力构型Ⅱ）

Fig.7Analysis results of scenarios (thrust configuration Ⅱ)

表1追逃场景初始相对状态

表2视线旋转坐标系下追逃系统状态空间网格

罗亚中, 李振瑜, 祝海. 航天器轨道追逃微分对策研究综述[J]. 中国科学: 技术科学,2020,50(12):1533-1545. LUO Y Z, LI Z Y, ZHU H. Survey on spacecraft orbital pursuit-evasion differential games[J]. Scientia Sinica(Technologica),2020,50(12):1533-1545.(in Chinese)

YANG F, YANG L P, ZHU Y W,et al. A DNN based trajectory optimization method for intercepting non-cooperative maneuvering spacecraft[J]. Journal of Systems Engineering and Electronics,2022,33(2):438-446.

ISAACS R. Differential games:a mathematical theory with applications to warfare and pursuit,control and optimization[J]. Journal of the Royal Statistical Society: Series A,1966,129(3):474-475.

ANDERSON G M, GRAZIER V W. Barrier in pursuit-evasion problems between two low-thrust orbital spacecraft[J]. AIAA Journal,1976,14(2):158-163.

张秋华, 赵小津, 孙毅. 空间飞行器在视线坐标系中的追逃界栅[J]. 航天控制,2007,25(1):26-30. ZHANG Q H, ZHAO X J, SUN Y. Pursuit-evasion barrier of two spacecrafts based on the sightline coordinate system[J]. Aerospace Control,2007,25(1):26-30.(in Chinese)

张秋华, 孙毅, 黄明明, 等. 近地共面轨道上两飞行器在径向连续小推力下的追逃界栅[J]. 控制与决策,2007,22(5):530-534. ZHANG Q H, SUN Y, HUANG M M,et al. Pursuit-evasion barrier of two spacecrafts under minute continuous radial thrust in coplanar orbit[J]. Control and Decision,2007,22(5):530-534.(in Chinese)

HAFER W T, REED H L, TURNER J D,et al. Sensitivity methods applied to orbital pursuit evasion[J]. Journal of Guidance, Control,and Dynamics,2015,38(6):1118-1126.

祝海. 基于微分对策的航天器轨道追逃最优控制策略[D]. 长沙: 国防科技大学,2017. ZHU H. Optimal control of spacecraft orbital pursuit-evasion based on differential game[D]. Changsha: National University of Defense Technology,2017.(in Chinese)

BANSAL S, CHEN M, HERBERT S,et al. Hamilton-Jacobi reachability:a brief overview and recent advances[C]//Proceedings of the IEEE 56th Annual Conference on Decision and Control(CDC),2017:2242-2253.

MITCHELL I M. The flexible,extensible and efficient toolbox of level set methods[J]. Journal of Scientific Computing,2008,35(2/3):300-329.

VENIGALLA C, SCHEERES D J. Spacecraft rendezvous and pursuit/evasion analysis using reachable sets[C]//Proceedings of the 2018 Space Flight Mechanics Meeting,2018.

VENIGALLA C, SCHEERES D J. Delta-V-based analysis of spacecraft pursuit-evasion games[J]. Journal of Guidance, Control,and Dynamics,2021,44(11):1961-1971.

ZHANG F, LIU T G, LIU M B. A high-order maximum-principle-satisfying discontinuous Galerkin method for the level set problem[J]. Journal of Scientific Computing,2021,87(2):45.

BECKERS D, ELDREDGE J D. Planar potential flow on Cartesian grids[J]. Journal of Fluid Mechanics,2022,941: A19.

STRIKWERDA J C. Finite difference schemes and partial differential equations[M].2nd ed. Philadelphia: Society for Industrial and Applied Mathematics,2007.

JIANG G S, SHU C W. Efficient implementation of weighted ENO schemes[J]. Journal of Computational Physics,1996,126(1):202-228.

JIANG Y Q, ZHOU S G, ZHANG X,et al. High-order weighted compact nonlinear scheme for one-and two-dimensional Hamilton-Jacobi equations[J]. Applied Numerical Mathematics,2022,171:353-368.

AHMAT M, QIU J X. Direct WENO scheme for dispersion-type equations[J]. Mathematics and Computers in Simulation,2023,204:216-229.

首页

期刊介绍

投稿指南

编委会

出版声明

开放获取声明

联系我们

期刊订阅

Rss

English

1 问题建模

2 求解方法

3 数值仿真

4 结论