一种解析构造航天器追逃界栅和判断捕获逃逸区域的方法与流程

本发明涉及航天器追逃博弈路径规划方法，尤其涉及一种航天器追逃界栅的构造及捕获逃逸区域判断方法。

背景技术：

空间交会规避机动是空间安全的一项关键技术。随着空间资源的日益紧张，在未来可能的空间对抗中，航天器将面临被追踪航天器捕获或杀伤的威胁。当对抗两航天器均能自主决策和机动时，该对抗便成为一个空间连续动态博弈的问题。

微分对策是一种用微分方法研究连续时间无限动态博弈的理论，非常适合对这类问题进行建模和分析。在应用该理论定性分析追逃的结局能否实现时，涉及到对捕获区、逃逸区以及二者的分界面——界栅的构造求解。而航天器追逃界栅和捕获逃逸区域分析方法则是将微分对策理论应用于航天器追逃博弈问题里的成果。对微分对策及界栅概念的详细介绍可以参看李登峰所著的《微分对策及其应用》一书。

构造求解界栅的方法通常可分为数值方法和解析方法。而目前对于航天器追逃界栅的求解基本上采用的都是数值方法。其原因之一是对策模型复杂：对抗在三维空间进行，同时进行追逃博弈的两航天器运动需要满足轨道动力学规律；二是对策求解困难：基于微分对策理论求解航天器轨道追逃的最优控制策略时，通常需要求解复杂的微分方程组。

但数值方法相比解析方法不便于快速简便地分析和计算，更不利于直接反映追逃博弈的运动规律。此外，采用数值方法求得的界栅应用起来也有一定局限性，对捕获区和逃逸区的分析更多止于两种区域的划分，不适合针对特定对策状态判断航天器当前所处区域是捕获区还是逃逸区。

技术实现要素：

为解决上述问题，本发明特提出一种用解析方法构造求解航天器追逃界栅和判断捕获逃逸区域的方法。该种解析方法能够快速计算求解追逃界栅，直观地反映追逃博弈的运动规律，进而明确判别航天器在追逃博弈中处于捕获区还是逃逸区，其结果可为空间威胁预警和博弈路径规划提供有效的参考。

本发明的技术方案和实施步骤如下：

一种解析构造航天器追逃界栅的方法，包括以下步骤：

s1，用微分对策方法构造初始模型；

进一步的，所述s1包括以下步骤：

s101，基于c-w方程构造微分对策的hamilton函数，

1)建立lvlh(localverticallocalhorizontal)坐标系，构造基于c-w方程的运动状态方程；

进一步的，以追逃两航天器附近的一颗虚拟航天器作为参考航天器建立当地轨道坐标系，该坐标系的原点位于参考航天器的质心o，ox轴沿参考航天器的径向，oz轴沿参考航天器轨道面的法向，oy轴沿参考航天器运动的轨迹切向，并与ox、oz轴构成右手坐标系；使用下式所示的c-w方程描述逃逸航天器或追踪航天器相对虚拟航天器的运动：

式中，ω为虚拟航天器做圆周运动的角速度，ax、ay、az为逃逸或追踪航天器分别在径向、迹向和法向推力加速度分量，若设追踪航天器和逃逸航天器的相对运动状态分别为(本文中的下标p和e分别表示追踪航天器和逃逸航天器)且均满足式给出的c-w方程；令为航天器追逃微分对策的状态变量，ae.x、ae.y、ae.z、ap.x、ap.y、ap.z分别为两航天器三个方向的推力加速度分量，则由式(1)可得式：

设两航天器均采用连续推力控制，最大加速度大小分别为tp、te，推力加速度方向用α和β表示，分别为航天器推力加速度方向的偏航角(ti在xoy平面投影和x轴夹角)和俯仰角(ti和xoy平面夹角)，αi∈[0,2π]、(下标i＝pore)，则可得微分对策的运动状态方程：

2)构造微分对策的hamilton函数；

将式(3)写成状态空间的表达式：

式中a、b为状态方程的系数矩阵，u＝te-tp；

构造微分对策的hamilton函数：

式中λ＝[λ1,λ2,λ3,λ4,λ5,λ6]^t为协态向量；

s102，求解微分对策的协态方程和最优控制方程；

1)求解微分对策的协态方程；

协态方程可以写为：

其解析表达式为：

λ(t)＝φλ(t,t0)λ(t0)(7)

式中状态转移矩阵φλ(t,t0)满足：

设tf为对策结束时刻，记τ＝tf-t为剩余捕获时间，可得：

假设追踪航天器的捕获半径为r0，则该微分对策的终端目标集边界为：

则终端时刻的横截条件为：

式中μ为乘子变量，为实数，在该边界下，优选的0.5；由上式可得协态量的末值条件：

式中为目标集边界上一点相对坐标系的偏航角与俯仰角；

(2)求解微分对策的最优控制方程；

最优控制方程即为该微分对策的鞍点，所谓鞍点，就是泛函的极值点(除掉最大值点及最小值点)，当微分对策的hamilton函数为追逃双方控制量的连续可导函数，鞍点可用公式表达为：

将hamilton函数代入上式得最优控制方程如下：

s2，求解初始模型，得到界栅解析解；

s201，据协态方程和末值条件，求解协态量和最优控制量的解析式；

进一步的，求解协态量的解析表达式为：

将协态量的解析式代入并在目标集边界τ＝0附近进行泰勒展开，综合考虑计算的工作量和精度的，优选的，保留至二阶项，可得追逃两航天器最优控制量与剩余捕获时间τ的解析表达式：

s202，将最优控制量的解析式代回c-w方程积分，求解航天器追逃界栅解析表达式；

具体的，将最优控制量的解析式代回式进行积分求解，可得界栅关于剩余捕获时间τ的解析表达式：

式中各系数如下：

其中：

特别的，当两航天器在共面圆轨道追逃博弈时，界栅解析解可表达为：

式中各系数如下：

其中：

至此，完成界栅的构造求解。

本发明同时公布了一种基于界栅解析解判断航天器捕获逃逸区域方法，具体方案如下：

s1，用微分对策方法构造初始模型；

s2，求解初始模型，得到界栅解析解；

进一步的，给定两航天器博弈情景中的参数，初始化界栅的解析解模型，即将参数带入下式：

s3，用流程方法判断航天器所处区域；

进一步的，所述s3包括以下步骤：

s301：给定对策的当前状态和最长博弈时间τs，令θ＝0°、τ＝τs；

s302：由式根据x＝xs、y＝ys、z＝zs求解出另外三个参数r、s5和s4；

s303：由式根据上步求解得到的一组参数，计算得到和

s304：保持其余参数不变，对每一个τ∈[0,τs]，分别计算出对应的和并从中找到使：

的τ，若存在这样的τ且此时r＞0，则进入s306，否则进入s307；

s305：返回s301，逐步增大θ和并重复s302～s304，直至两个方位角达到360°；

s306：对得到的解，若存在多组满足式的解，取其中r最小的那组解；

s307：根据得到的解对捕获逃逸区域进行判断：如果r＞r0，则对策当前状态处于逃逸区内；如果r＜r0，则对策当前状态处于捕获区内；如果r＝r0，则对策当前状态恰好处于界栅上；若不存在这样的一组解，则对策当前状态处于逃逸区内；

至此，完成了航天器处于捕获逃逸区域的判断。

本发明一种解析构造航天器追逃界栅和判断捕获逃逸区域的方法的优点如下：

1.本发明采用解析方法构造求解航天器的追逃界栅，并基于解析法求得的界栅明确判断航天器当前状态所在区域是捕获区还是逃逸区。经验证，这种解析方法在计算结果上精度较高，与数值方法结果吻合度较好，却更为简便快捷，更好地反映博弈运动规律。这种对航天器所在区域的判断方法，相比于对捕获逃逸区域的划分方法，结果更加明确、实用。

2.特别地，航天器的追逃博弈更可能是在轨实时进行，由于星载计算机的计算能力相对较弱，因此解析计算追逃界栅和判断捕获逃逸区域的方法，更适合星载计算机快速准确地预警和求解博弈机动策略，对于实际工程任务具有更好的适用性。

附图说明

图1为本发明方法中构造求解界栅解析解的技术方案图。

图2为本发明方法中判断捕获逃逸区域的技术方案图。

图3为本发明实施例一的精度对比图。

图4为本发明实施例二的精度对比图。

图5为本发明实施例三的精度对比图。

图6为本发明实施例四的精度对比图。

图7为本发明实施例五的基本流程图。

图8为本发明实施例五的追逃验证图。

具体实施方式

实施例一：

如图1所示，本实施例一种解析构造航天器追逃界栅方法的步骤包括：

s1，用微分对策方法构造初始模型；

s101，基于c-w方程构造微分对策的hamilton函数，

1)建立lvlh(localverticallocalhorizontal)坐标系，构造基于c-w方程的运动状态方程；

以追逃两航天器附近的一颗虚拟航天器作为参考航天器建立当地轨道坐标系，该坐标系的原点位于参考航天器的质心o，ox轴沿参考航天器的径向，oz轴沿参考航天器轨道面的法向，oy轴沿参考航天器运动的轨迹切向，并与ox、oz轴构成右手坐标系；使用下式所示的c-w方程描述逃逸航天器或追踪航天器相对虚拟航天器的运动：

式中，ω为虚拟航天器做圆周运动的角速度，ax、ay、az为逃逸或追踪航天器分别在径向、迹向和法向推力加速度分量，若设追踪航天器和逃逸航天器的相对运动状态分别为且均满足式给出的c-w方程；令为航天器追逃微分对策的状态变量，ae.x、ae.y、ae.z、ap.x、ap.y、ap.z分别为两航天器三个方向的推力加速度分量，则由式(1)可得：

设两航天器均采用连续推力控制，最大加速度大小分别为tp、te，推力加速度方向用α和β表示，分别为航天器推力加速度方向的偏航角(ti在xoy平面投影和x轴夹角)和俯仰角(ti和xoy平面夹角)，αi∈[0,2π]、(下标i＝pore)，则可得微分对策的运动状态方程：

2)构造微分对策的hamilton函数；

将式(3)写成状态空间的表达式：

式中a、b为状态方程的系数矩阵，u＝te-tp；

构造微分对策的hamilton函数：

式中λ＝[λ1,λ2,λ3,λ4,λ5,λ6]^t为协态向量；

s102，求解微分对策的协态方程和最优控制方程；

1)求解微分对策的协态方程；

协态方程可以写为：

其解析表达式为：

λ(t)＝φλ(t,t0)λ(t0)(7)

式中状态转移矩阵φλ(t,t0)满足：

设tf为对策结束时刻，记τ＝tf-t为剩余捕获时间，可得：

假设追踪航天器的捕获半径为r0，则该微分对策的终端目标集边界为：

则终端时刻的横截条件为：

式中μ＝0.5；由上式可得协态量的末值条件：

式中为目标集边界上一点相对坐标系的偏航角与俯仰角；

(2)求解微分对策的最优控制方程；

鞍点可用公式表达为：

将hamilton函数代入上式得最优控制方程如下：

s2，求解初始模型，得到界栅解析解；

s201，据协态方程和末值条件，求解协态量和最优控制量的解析式；

本实施例中，求解协态量的解析表达式为：

将协态量的解析式代入并在目标集边界τ＝0附近进行泰勒展开，保留至二阶项，可得追逃两航天器最优控制量与剩余捕获时间τ的解析表达式：

s202，将最优控制量的解析式代回c-w方程积分，求解航天器追逃界栅解析表达式：

式中各系数如下：

其中：

特别的，当两航天器在共面圆轨道追逃博弈时，界栅解析解可表达为：

式中各系数如下：

其中：

输入情景参数：

得到界栅分量随剩余捕获时间变化曲线，见图2中解析求解结果。

为验证解析求解结果计算的可靠性与精度，本实施例中也将解析方法求解得到的界栅结果与数值方法进行比较，见图3中数值求解结果。从图中可以看出，解析公式求解的结果和数值积分得到的结果吻合得较好，这表明本发明所提出的界栅求解方法精度较高。另外，从图中还可以看出，解析求解方法在方向上得到的位置和速度结果的精度比x方向上的要高，解析求解结果中精度最低的是方向x上的速度随着剩余捕获时间τ的增大，的计算结果发散较快。

实施例二：

本实施例与实施例一的步骤基本相同，其主要不同点为：本实施例中，捕获点方位角不再是给定值，而是作为自变量。因而通过计算可以得到不同捕获点方位角下微分对策在τ＝1200s时刻的界栅。

再比较数值积分的结果和解析求解的结果，见图4。从图中可以看出，方位角仅对界栅y方向上位置的解析求解精度有一定影响，当或附近时误差明显，而在其它界栅状态量的求解上，方位角影响较小。

实施例三：

本实施例与实施例二的步骤基本相同，其主要不同点为：本实施例中，自变量是捕获点相对速度

图5中分别给出了不同捕获点相对速度下计算得到的对策在τ时刻的界栅，并比较了数值积分结果和解析求解结果。从图中可以看出，对不同的解析计算的结果与数值积分得到的结果均吻合较好。同时从图中还可以看出，在给定的范围内，界栅在x和y两个方向上的位置均随着的增大而减小，而在两个方向上的速度则随之增大。

实施例四：

本实施例与实施例二的步骤基本相同，其主要不同点为：本实施例中，自变量是捕获半径r0。

图6中分别给出了不同捕获半径下计算得到的对策在τ时刻的界栅，并比较了数值积分结果和解析求解结果。从图中可以看出，捕获半径r0对界栅的求解影响较小，且解析求解带来的误差随r0的变化也较小。

实施例五：

本实施例中，基于解析求解的界栅判断捕获逃逸区域方法的步骤包括：

s1，用微分对策方法构造初始模型，此步骤同实施例一；

s2，求解初始模型，得到界栅解析解：

给定两航天器博弈情景中的参数，初始化界栅的解析解模型，即将参数带入下式：

输入情景参数：

然后初始化界栅的解析解；

s3，用流程方法判断航天器所处区域，具体步骤的流程图见图7：

s301：给定对策的当前状态和最长博弈时间τs，令θ＝0°、τ＝τs；

s302：由式根据x＝xs、y＝ys、z＝zs求解出另外三个参数r、s5和s4；

s303：由式根据上步求解得到的一组参数，计算得到和

s304：保持其余参数不变，对每一个τ∈[0,τs]，分别计算出对应的和并从中找到使：

的τ，若存在这样的τ且此时r＞0，则进入s206，否则进入s207；

s305：返回s301，逐步增大θ和并重复s302～s304，直至两个方位角达到360°；

s306：对得到的解，若存在多组满足式的解，取其中r最小的那组解；

s307：根据得到的解对捕获逃逸区域进行判断：如果r＞r0，则对策当前状态处于逃逸区内；如果r＜r0，则对策当前状态处于捕获区内；如果r＝r0，则对策当前状态恰好处于界栅上；若不存在这样的一组解，则对策当前状态处于逃逸区内；

至此，完成了航天器处于捕获逃逸区域的判断。

下表给出了共面轨道上四种不同初始对策状态下捕获逃逸区域的计算结果，相对轨迹仿真结果见图8。

表1航天器追逃博弈捕获逃逸区域分析

从表中数据可以看出，第1、3种情况下对策状态处于捕获区内，而另外两种情况下的对策状态处于逃逸区内，其中第4种情况中在给定时间范围内不存在最小微分对策的解。