一种航天器飞行博弈中可交会的快速判别方法

本发明涉及一种航天器飞行博弈中可交会的快速判别方法，特别涉及一种在追逐航天器与逃避航天器进行博弈时，基于可达能力边界，对短时间内能否成功与目标交会进行快速判断的方法，属于航空航天领域。

背景技术：

随着人类空间技术的发展，空间对抗已经成为未来空间任务中的重要组成部分。在空间对抗任务设计时，根据博弈中的两航天器的初始位置和发动机参数限制，需要快速判断在逃避航天器逃避时，追逐航天器能否实现对目标的交会。针对一般轨道上航天器可达范围、博弈末端位置难以确定等问题，开展航天器可交会的快速判断问题研究。本发明结合给定的博弈初始时刻两航天器状态(位置、速度和质量)、燃料约束、推力幅值约束等，提出一种航天器飞行博弈中可交会的快速判别方法，可有效解决航天器可交会难以确定的问题。

在已发展的航天器可交会判断的研究中，在先技术[1](参见spacecraftrendezvousandpursuit/evasionanalysisusingreachablesets,venigalla,c.,andscheeres,d.j.,2018spaceflightmechanicsmeeting,kissimmee,florida,2018)，基于三维轨道根数可达边界(半长轴、偏心率和轨道倾角)，进行时间自由的航天器博弈交会性判断，但三维轨道根数无法与空间位置一一对应，且其研究难以考虑时间约束，从而无法解决此类航天器可交会的判断问题。

在先技术[2](参见pursuit-evasiongameforsatellitesbasedoncontinuousthrustreachabledomain,ieeetransactionsonaerospaceandelectronicsystems,gong,h.,gong,s.,andli,j.,2020)，基于航天器的位置可达性进行交会性判断。但其可达集仅限于连续推力的圆轨道，无法处理存在发动机关机情况和航天器处在一般轨道情况，且未给出合理的博弈末端近似时间，导致可达集的多次无效计算。所以，该方法无法解决一般航天器可交会性判断问题。

技术实现要素：

本发明的目的是提供一种航天器飞行博弈中可交会的快速判别方法，该方法能够在航天器博弈问题中，利用航天器的可达能力边界，快速判断目标可交会性。本发明具有计算速度快、收敛性好、适用性强的优点。

本发明的目的是通过下述技术方案实现的。

本发明公开的一种航天器飞行博弈中可交会的快速判别方法，在地球质心柱坐标惯性系下建立探测器动力学方程，并给定航天器动力学方程的初始条件和推力质量约束；给定时间上限，利用动力学方程进行无控积分，找到追逐航天器和逃避航天器相距最近的位置所对应的时刻td；分别求解td时刻追逐航天器和逃避航天器的可达能力范围；对两航天器可达能力范围边界的几何关系进行判断，若逃避航天器的可达能力范围完全包含于追逐航天器的可达能力范围，则逃避航天器一定可交会；否则，不一定能与逃避航天器交会。所述一定可交会指无论逃避航天器采用何种控制，一定存在一种控制方式使追逐航天器末端位置与逃避航天器重合。

本发明公开的一种航天器飞行博弈中可交会的快速判别方法，包括如下步骤：

步骤一、在地球质心柱坐标惯性系下建立探测器动力学方程，并给定航天器动力学方程的初始条件和推力质量约束；

追逐航天器已经发射到达大气层外后，开始与逃避航天器的博弈；由于整个过程历时较短，因此忽略轨道摄动；航天器的动力学方程如下：

其中r为航天器在柱坐标下的平面径向位置、θ为平面转角、z为高度方向位置、vr为平面径向速度、vθ为平面法向速度、vz为高度方向速度、m为航天器质量，isp为发动机比冲，g0为海平面重力加速度，μ为中心天体引力常数；tr,tθ,tz非别为三个方向推力，t为总推力大小，角标i＝p时代表追逐航天器；角标i＝e时代表逃避航天器；

对于推力，存在约束

在航天器博弈过程中，两航天器初始状态已知

xp(t0)＝xp0,xe(t0)＝xe0(3)

其中t0为初始时刻，xp为追逐航天器状态，xe为逃避航天器状态；

博弈结束时每个航天器的质量不能低于其最小质量

mp(tf)≥mpf,me(tf)≥mef(4)

其中tf为博弈结束时刻、mpf为追逐航天器的最小质量、mef为逃避航天器的最小质量；

另外，航天器在整个博弈对抗中，发动机推力存在上限，即

ti≤ci(5)

其中ci为发动机推力上限。

步骤二、给定时间上限tu，利用动力学方程进行无控积分，找到追逐航天器和逃避航天器相距最近的位置所对应的时刻td；

利用初始条件(3)和动力学方程(1)，控制设定为ti＝0，在[t0,tu]时间内对两航天器进行积分，得到两航天器无控状态下的轨道xi(t)；对相同时间下两航天器的位置作差，得到两航天器相对位置

由航天器相对位置矢量得到相对位置最小的时刻td

步骤三、分别求解td时刻追逐航天器和逃避航天器的可达能力范围；

可达能力范围为航天器从固定初始状态x0，选取不同控制，末端时刻所能到达的所有位置的集合；由于可达集拥有平面对称性，因此在求解之前可将坐标系进行旋转，使航天器所在的速度面为z＝0平面；设航天器地心惯性系下初始位置为rs，速度为vs，则由地心惯性系到可达集计算惯性系的旋转矩阵为：

m＝[iajaka](8)

其中m为由地心惯性系到可达集计算惯性系的旋转矩阵，ia,ja,ka为矩阵m中分量；

则在可达集计算惯性系下，可达集关于xy平面对称，此时再将可达集计算惯性系转化为柱坐标进行求解，能够减少计算量；由笛卡尔坐标转换为柱坐标的方式为：

其中rcy,θcy,zcy,vrcy,vθcy,vzcy为柱坐标下六个状态分量，xca,yca,zca,vxca,vyca,vzca为笛卡尔坐标六个分量；

在短时间内，航天器可达集形状通常较规则，因此通过三步进行求解：

第一步：确定平面转角极值的性能指标，表示为：

minj＝±θ(tf)(11)

第二步：将转角范围进行离散，得到θs,s＝1,2,…,nθ，nθ为转角离散的节点数；分别求解确定转角的高度极值z(tf)

其中θs为离散后的可达的转角值；

第三步：在转角确定的情况下，将高度进行离散为zk,k＝1,2,…,nz，nz为高度离散的节点数；分别求解确定转角确定高度下的平面径向极值

其中zk为转角为θs时离散后的可达的高度值；

动力学约束为式(1)，初始位置约束为式(3)，此外还存在路径约束(2)、(5)和终端状态约束(4)；由于序列凸优化可以高效准确的求解最优控制问题，在此应用序列凸优化求解航天器可达集；需对上述约束进行凸化；为使方程更加简化，定义

其中，σ,ηr,ηθ,ηz,mz为中间变量；

则式(1)与质量和推力相关项转化为

则有

其中x＝[r,θ,z,vr,vθ,vz,mz]^t为状态量，u＝[ηr,ηθ,ηz,σ]^t为控制量，f(x,u)为式(1)变化后的动力学方程；

相应约束同时变为

mz(tf)≥ln(mf)(18)

质量相关的非线性项被转化为线性项，将求解难度降低；

将式(16)在(x^k,u^k)处进行线性化，其中(x^k,u^k)为给定初值，得

其中：

c(x^k,u^k)＝f(x^k,u^k)-a(x^k)x^k-bu^k

在约束中存在非凸的二次等式约束，现将式(17)凸化为

在可达集求取过程中，式(21)取得边界；

性能指标和所有约束皆转化为凸形式，则可将动力学与约束进行离散，并利用凸优化求解器进行求解；求得所有可达集计算惯性系下柱坐标的可达集边界点后，将边界点转换到笛卡尔坐标下，转换公式如下

然后利用式(8)中矩阵将笛卡尔坐标下的边界点转换到地心惯性系中，由可达集计算惯性系到地心惯性系的旋转矩阵为m^t；地心惯性系下所有可达集边界点构成td时刻追逐航天器和逃避航天器的可达能力范围边界。

步骤四：对两航天器可达能力范围边界的几何关系进行判断，若逃避航天器的可达能力范围完全包含于追逐航天器的可达能力范围，则逃避航天器一定可交会；否则，不一定能与逃避航天器交会；

由于可达集为紧集，因此如果目标可达范围边界包含于追逐航天器可达能力边界内，则认为逃避航天器的可达能力范围完全包含于追逐航天器的可达能力范围，因此需要对得到的每个逃避航天器可达能力边界点进行判定；对于地心惯性系下任一逃避航天器可达能力边界点xtec，将所述边界点xtec转到追逐航天器可达集计算惯性系下

xtkj＝mpxtec(23)

利用式(10)将得到的追逐航天器可达集计算惯性系下逃避航天器可达能力边界点转换到柱坐标下，记为(rtkj,θtkj,ztkj)，与追逐航天器可达范围进行比较；首先比较平面转角方向，若θpmin≤θtkj≤θpmax，则设定pθ＝1，否则pθ＝0，其中θpmin,θpmax分别为追逐航天器可达范围的平面转角最小值和最大值，pθ为转角方向判断参数；若pθ＝1，则比较高度方向，若ztkj≤zpmax(θtkj)，则pz＝1，否则pz＝0，其中zpmax(θtkj)为追逐航天器平面转角为θtkj时对应的高度方向可达范围的最大值，pz为高度方向判断参数；若pθ＝pz＝1，则比较平面径向方向，若rpmin(θtkj,ztkj)≤rtkj≤rpmax(θtkj,ztkj)，则设pr＝1，否则pr＝0，其中rpmin(θtkj,ztkj)和rpmax(θtkj,ztkj)分别为追逐航天器平面转角为θtkj，高度为ztkj时对应的平面径向方向可达范围的最小值和最大值，pr为平面径向方向判断参数；若pθ＝pz＝pr＝1，则(rtkj,θtkj,ztkj)在追逐航天器可达集中，否则不在；若(rtkj,θtkj,ztkj)在追逐航天器可达集中，则验证下一个逃避航天器可达能力边界点，否则可直接确定不一定能与逃避航天器交会；遍历所有逃避航天器可达能力边界点后，若全部在追逐航天器可达集中，则逃避航天器一定可交会。

步骤五：根据步骤四得到的航天器飞行博弈中可交会判别结果，采用相应航天器飞行博弈控制策略，进而解决航天器飞行博弈领域相关技术问题。

有益效果：

1、本发明公开的一种航天器飞行博弈中可交会的快速判别方法，未对可达集的形状、航天器所在轨道作出限制，因此适用于一般轨道上航天器的可达集求解。

2、本发明公开的一种航天器飞行博弈中可交会的快速判别方法，通过考虑对航天器的燃料约束，因此适用于考虑航天器燃料约束的目标可交会的判断。

3、本发明公开的一种航天器飞行博弈中可交会的快速判别方法，通过判断目标可交会的准确时刻，因此可能够避免不同时刻的可达集的重复计算，提高收敛性，为目标可交会的判断提高效率，所述提高效率为相较于先技术[2]提高效率。

附图说明

图1是本发明公开的一种航天器飞行博弈中可交会的快速判别方法的流程图；

图2是本发明公开的一种航天器飞行博弈中可交会的快速判别方法的地心惯性系下追逐航天器可达能力范围图像；

图3是本发明公开的一种航天器飞行博弈中可交会的快速判别方法的地心惯性系下逃避航天器可达能力范围图像；

图4是本发明公开的一种航天器飞行博弈中可交会的快速判别方法的地心惯性系下两航天器可达能力范围位置关系图像。

具体实施方式

为了更好的说明本发明的目的和优点，下面结合附图和实例对发明内容做进一步说明。

实施例1：

刚发射升空的飞行器与低轨目标卫星进行博弈，利用所提出的方法进行目标可交会性的判断。设定追逐航天器发动机提供推力3000n，逃避航天器推力为800n，两航天器发动机比冲均为280s，质量均为200kg。其中逃避航天器可用质量为40kg，追逐航天器可用质量为50kg。两航天器初始状态参数如表1所示。

表1两航天器轨道参数

如图1所示，本实施例公开的一种面向航天器博弈中质量约束的处理方法，具体实现步骤如下：

步骤一：在地球质心柱坐标惯性系下建立探测器动力学方程，并给定航天器动力学方程的初始条件和推力质量约束。

追逐航天器已经发射到达大气层外后，开始与逃避航天器的博弈。由于整个过程历时较短，因此忽略轨道摄动。航天器的动力学方程如下：

其中r,θ,z,vr,vθ,vz,m分别为航天器在柱坐标下的平面径向、法向和高度方向的位置、速度和航天器质量，isp为发动机比冲，g0为海平面重力加速度。tr,tθ,tz分别为三个方向推力，t为总推力大小，角标i＝p,e代表追逐航天器或逃避航天器。

对于推力，存在约束

在航天器博弈过程中，两航天器初始状态已知

xp(t0)＝xp0,xe(t0)＝xe0

其中t0为初始时刻，xp为追逐航天器状态，xe为逃避航天器状态。

博弈结束时每个航天器的质量不能低于其最小质量

mp(tf)≥mpf,me(tf)≥mef

其中mpf＝150kg,mef＝160kg

另外，航天器在整个博弈对抗中，发动机推力存在上限，即

ti≤ci

其中cp＝3000n,ce＝800n

步骤二：给定时间上限tu，利用动力学方程进行无控积分，找到追逐航天器和逃避航天器相距最近的位置所对应的时间td。

在此，给定时间上限tu＝1800s，利用初始条件(3)和动力学方程(1)，控制设定为ti＝0，在[t0,tu]时间内对两航天器进行积分，得到两航天器无控状态下的轨道xi(t)。对相同时间下两航天器的位置作差，可得两航天器相对位置

由航天器相对位置矢量得到相对位置最小的时刻td

在此，航天器相对位置矢量得到相对位置最小的时刻td＝685.19s。

步骤三：分别求解td时刻追逐航天器和逃避航天器的可达能力范围。

在求解之前可将坐标系进行旋转，使航天器所在的速度面为z＝0平面。由地心惯性系到可达集计算惯性系的旋转矩阵为：

则在可达集计算惯性系下，可达集关于xy平面对称，所得两航天器在各自可达集计算惯性系下状态为：

其中xca,yca,zca,vxca,vyca,vzca为笛卡尔坐标六个分量。

再将其转化为柱坐标进行求解。由笛卡尔坐标转换为柱坐标的方式为：

其中rcy,θcy,zcy,vrcy,vθcy,vzcy为柱坐标下六个状态分量，所得结果为

在短时间内，航天器可达集形状通常较规则，因此可通过三步进行求解：

第一步为确定平面转角的极值，其性能指标为

minj＝±θf

第二步为将转角范围进行离散为θs,s＝1,2,…,nθ，nθ为转角离散的节点数。

分别求解确定转角的高度极值

第三步为在转角确定的情况下，将高度进行离散为zk,k＝1,2,…,nz，nz为高度离散的节点数。分别求解确定转角确定高度下的平面径向极值

动力学约束为式(1)，初始位置约束为式(3)，此外还存在路径约束(2)、(5)和终端状态约束(4)。由于序列凸优化可以高效准确的求解最优控制问题，在此应用序列凸优化求解航天器可达集。需对上述约束进行凸化。为使方程更加简化，定义

则式(1)与质量和推力相关项转化为

相应约束同时变为

mz(tf)≥mzmin

看到，质量相关的非线性项被转化为线性项，将求解难度降低。

此时状态量为x＝[r,θ,z,vr,vθ,vz,mz]^t，控制量为u＝[ηr,ηθ,ηz,σ]^t。将上式在(x^k,u^k)处进行线性化，可得

x'＝a(x^k)x+b(x^k)u+c(x^k)

在约束中存在非凸的二次等式约束，现将其凸化为

可以证明，在可达集求取过程中，该松弛项可取得边界，即满足原问题。

可见，性能指标和所有约束皆转化为凸形式，则可将动力学与约束进行离散，并利用凸优化求解器进行求解。求得所有可达集计算惯性系下柱坐标的可达能力范围边界点后，转换到笛卡尔坐标下转换公式如下

然后利用式(8)中矩阵将其转换到地心惯性系中，由可达集计算惯性系到地心惯性系的旋转矩阵为m^t。所得在地心惯性系的追逐航天器可达能力范围图像如图2所示，逃避航天器图像如图3所示。

步骤四：对两航天器可达能力边界的几何关系进行判断，若逃避航天器的可达能力范围完全包含于追逐航天器的可达能力范围，则逃避航天器一定可交会；否则，不一定能与逃避航天器交会。

由于工程上一般认为可达集为紧集，因此如果目标可达范围边界包含于追逐航天器可达能力边界内，则认为逃避航天器的可达能力范围完全包含于追逐航天器的可达能力范围，因此需要对得到的每个逃避航天器可达能力边界点进行判定。对于地心惯性系下某个逃避航天器可达能力边界点xtec，将其转到追逐航天器可达集计算惯性系下

xtkj＝mpxtec

利用式(10)将得到的追逐航天器可达集计算惯性系下逃避航天器可达能力边界点转换到柱坐标下，记为(rtkj,θtkj,ztkj)，与追逐航天器可达范围进行比较。首先比较平面转角方向，θpmin≤θtkj≤θpmax，则pθ＝1，其中θpmin,θpmax分别为追逐航天器可达范围的平面转角最小值和最大值。由于pθ＝1，则比较高度方向，ztkj≤zpmax(θtkj)，则设pz＝1，其中zpmax(θtkj)为追逐航天器平面转角为θtkj时对应的高度方向可达范围的最大值。由于pθ＝pz＝1，则比较平面径向方向，rpmin(θtkj,ztkj)≤rtkj≤rpmax(θtkj,ztkj)，则设pr＝1，其中rpmin(θtkj,ztkj)和rpmax(θtkj,ztkj)分别为追逐航天器平面转角为θtkj，高度为ztkj时对应的平面径向方向可达范围的最小值和最大值。由于pθ＝pz＝pr＝1，则(rtkj,θtkj,ztkj)在追逐航天器可达集中。由于(rtkj,θtkj,ztkj)在追逐航天器可达集中，则验证下一个逃避航天器可达能力边界点。遍历所有逃避航天器可达能力边界点后，全部在追逐航天器可达集中，则可知逃避航天器一定可交会。两个可达能力范围的位置关系图如图4所示。

步骤五：根据步骤四得到的航天器飞行博弈中可交会判别结果，采用相应航天器飞行博弈控制策略，进而解决航天器飞行博弈领域相关技术问题。

以上所述的具体描述，对发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。