CN107832528A - 一种空天飞行器三维空间可达域分析方法 - Google Patents

Abstract

本发明公开了一种空天飞行器的三维空间可达域分析方法，根据不同的推力方式和是否有时间约束条件可分为三种情况：其一，建立机动点位置、脉冲方向不同状态组合下的基于轨道六根数的三维空间可达域模型，分析脉冲变轨下的可达域；其二，建立具有时间约束的脉冲变轨三维空间可达域模型，分析不同脉冲大小、不通脉冲方向、不同时间长度下的可达域范围；其三，建立基于改进春分点轨道根数的有限推力变轨三维空间可达域模型，分析固定时间下是否考虑燃耗约束的可达域，分析有限推力对轨道六根数的可达域。

Description

一种空天飞行器三维空间可达域分析方法

技术领域

本发明涉及一种空天飞行器三维空间可达域分析方法，属于空天飞行器总体设计技术领域。

背景技术

空间机动目标的轨道预测和及时探测是天基空间目标探测系统的关键任务之一，更是空间态势感知的一个重要问题。随着我国航天事业的发展，对空间机动目标及时探测的需求更是越来越迫切。

天基空间目标探测系统在对空间重点目标进行观测时，常会由于地球遮挡、太阳光照、探测载荷能力等因素的影响而无法持续观测目标。若目标在系统不可观测时段内发生了非合作机动变轨，系统将无法在预报的下一次可观测到目标的位置捕获该目标，即目标丢失。因此，及时、快速地探测到该机动目标，并分析其变轨意图，便成为获取制空权和实施其他空间行动的重要基础。

空间机动目标的轨道预测的途径就是分析飞行器在空间的可达域，其设计结果的合理性直接关系到飞行器的安全与准确探测。考虑到空天飞行器空间任务的拓展，需要针对飞行器开展快速计算空间机动目标的可达域计算方法研究，以确保空间机动目标探测的可行，满足空天安全需求。

发明内容

本发明的目的是提出一种空天飞行器三维空间可达域分析方法，能够准确计算空天飞行器机动后所能到达的所有位置点的集合。

本发明的目的是通过以下技术方案实现的：

一种空天飞行器三维空间可达域分析方法，包括如下内容：

(1)建立机动点位置、脉冲方向不同状态组合下的基于轨道六根数的三维空间可达域模型，分析脉冲方式变轨下的三维空间可达域；具体为通过机动点位置f、脉冲方向对应的方位角α与高度角β、方向矢量参数κ,υ推导得到可达目标的位置公式。

(2)建立具有时间约束的脉冲变轨三维空间可达域模型，分析固定时间长度下脉冲方式变轨下的三维空间可达域；根据脉冲在初始轨道施加前后轨道根数之间的关系，得到脉冲机动后的位置矢量r和速度矢量v：

C₁为轨道坐标系到地心惯性坐标系的坐标转换矩阵，C₂为当地轨道坐标系到地心惯性坐标系的坐标转换矩阵，r₁为脉冲施加时刻的飞行器地心距，f为脉冲施加时刻的飞行器真近点角，e₀为飞行器初始轨道偏心率，p₀为半通径，μ为地球引力常数，Δv为施加脉冲的大小，(α,β)为脉冲施加方向角。

(3)建立基于改进春分点轨道根数的有限推力变轨三维空间可达域模型，分析有限推力方式变轨下的三维空间可达域；根据有限推力变轨三维空间可达域模型，分析时间固定条件下燃料约束对可达域的影响，得到有时间约束的有限推力变轨下三维空间可达域。

本发明与现有技术相比，具有如下有益效果：

(1)本发明考虑脉冲变轨方式和有限推力变轨方式，基于轨道六根数建立的三维空间可达域模型，可计算空间飞行器在任意机动点位置、任意脉冲方向和大小、任意时间长度下的飞行器可达空间最大包络，该包络可用于分析飞行器潜在的机动变轨能力；

(2)本发明提出的三维空间可达域分析方法可用于空间轨道典型任务分析，分析结果可为在轨服务和轨道拦截任务确定飞行器所达范围内的所有可服务目标或潜在目标，可为空间观测任务提供最佳观测位置的可选集合，可为非合作目标预估一定时间范围内机动后的区域范围提供参考依据。

(3)本发明建立的三维空间可达域模型可以计算飞行器不同约束条件下的三维可达域，这种可达域通过简单的坐标系变换可投影到地球表面上，为地面的测控布站或已有观测站点选取提供数据依据。

附图说明

图1为本发明RSW坐标系示意图；

图2为本发明变轨位置任意变轨方向任意情况下飞行器可达域；

图3为本发明任意机动点、固定脉冲方向的可达域；

图4为本发明任意机动点、固定脉冲方向的可达域在x-y平面投影图；

图5为本发明固定机动点、任意脉冲方向的可达域；

图6为本发明不同时间约束下的飞行器三维空间可达域；

图7为本发明有限推力方式的固定时间下不考虑燃料约束的飞行器三维空间可达域；

图8为本发明有限推力方式的固定时间下考虑燃料约束的飞行器三维空间可达域。

具体实施方式

下面结合附图，对本发明作进一步详细的说明，本发明的空天飞行器三维空间可达域分析方法，包括如下内容：

一、脉冲变轨方式的基于轨道六根数的三维空间可达域分析

1、机动点任意脉冲方向任意的三维空间可达域分析

(1)可达域目标位置方程

假设飞行器在一个半长轴为a₀，离心率为e₀的初始轨道上。单位矢量p在地心第二轨道坐标系S_xyz下κν可以表示成：

p＝[cosυcosκ sinυcosκ sinκ]^T (2-1)

这里κ为单位矢量高度角，υ为升交点幅角。

机动点r₀的选择完全依赖于真近点角f。为了覆盖初始轨道的所有范围，f的作用范围定义为[υ-π,υ+π)。定义RSW坐标系S₀(x₀y₀z₀)，具体的描述如图1，以飞行器质心为原点，x₀轴与r₀轴相一致，z₀轴与初始轨道角动量方向一致，则两个矢量r₀与v₀在该坐标系的投影为：

(r₀)₀＝[r₀ 0 0]^T (2-2)

下标0为矢量在坐标系S₀上的投影，r₀＝p₀/(1+e₀ cos f)，p₀指初始轨道的半通径。

由r₀,p决定的平面定义为M，如图1所示。图中还展示了一个新的坐标系S₁(x₁y₁z₁)，它是由坐标系S₀通过绕X轴旋转一个角度i得到的，转换矩阵如下：

由几何关系容易得到：

tan κ＝sin(υ-f)tan i，i∈(-π/2,π/2) (2-5)

反解出：

轨道机动的过程：

Δv＝Δvl(α,β) (2-7)

其中脉冲矢量大小Δv是由搭载的燃料决定的，而l(α,β)是由脉冲矢量的方向决定，其在坐标系S₁中表示为：l＝[cosαcosβ sinαcosβ sinβ]^T，其中α∈[0,2π)，

β∈[-π/2,π/2] (2-8)

α,β分别表示脉冲矢量Δv的方位角与高度角。转移轨道的初始状态量可通过下式获得：

方程(2-9)投影到坐标系S₁中可写为如下形式：

r₁＝[r₀ 0 0]^T (2-10)

为了确保通过矢量p到达目标位置，速度矢量v₁应该在平面M内，也即在坐标系S₁下z₁方向的分量应该为零，由方程(2-11)可得：

由于β的取值范围是[-π/2,π/2]，可以得出如下关系式：

把(2-12)，(2-13)代入方程(2-14)可以得到：

v₁＝v_M+Δv_Mj(α) (2-14)

其中：

j(α)＝[cosα sinα 0]^T (2-16)

从而得到v₁在水平方向与半径方向的分量分别为

在转移轨道面M内定义一个角度Δθ代表初始位置与目标位置的转角，由几何关系(三余弦公式)可以得到：

cosΔθ＝cos(υ-f)cos κ (2-19)

可以更进一步得出：

进一步定义θ₁为转移轨道上r₁对应的真近点角，则目标位置矢径为：

其中e₁,p₁分别指的是机动后的转移轨道偏心率与半通径。根据开普勒轨道理论v₁在水平方向与半径方向的分量可表示为：

这样就得到：

把上式代入(2-21)式，并结合飞行器动量矩h大小的计算公式可以得到：

式(2-24)即为可达域目标位置方程，它由控制变量f,α以及方向矢量参数κ,υ所决定。三个变量不是互相独立的，满足方程(2-13)。写为矢量形式为：

r＝r(f,α；κ,υ)p(κ,υ) (2-25)

(2)可达域边界方程

假设r_in，r_out表示可达域在p方向上的边界点，故也是函数r(f,α)的极限点，因此，边界点可以通过求r(f,α)的边界值来找到。定义x＝[α f]^T，函数r(x)在边界点处的梯度为零，即满足

从而边界方程的推导转化为以上两个偏微分的推导。

首先，由式(2-24)可得：

根据则式(2-27)化简为：

由关系式(2-18)可以得到上式中的两个偏微分为：

带入式(2-28)，则式(2-26)中第一个偏微分最终结果为：

其中，

其次，由方程(2-24)可以得到：

由关系式h＝r₀v_u可以得到：

带入到式(2-32)中可以得到：

其中：

由Q(f,α)的表达式可以看出，仍然存在一些偏导数是未知的。

由方程(2-6)可以推算出新的微分关系式为：

因此可以得到：

结合方程(2-18)可以得到：

而通过式(2-17)可以得到：

将式(2-37)代入(2-40)得到：

初始矢径r₀可以通过r₀＝p₀/(1+e₀ cos f)求出，所以

对式(2-19)关于f求偏导数有

因此可以知道

因此通过以上式子就可以求解出Q(f,α)。

通过引入变量H(x)＝[P Q]^T，则r(x)可以被写为

因此位置矢径的边界点r_in和r_out可以通过求解下式

H(x)＝[P Q]^T＝0 (2-46)

得到x^*＝[f^* α^*]，从而求得r(f,α)。

3)机动点任意脉冲方向任意的三维空间可达域分析

考虑的约束条件分别为：

脉冲方向：κ_min≤κ≤κ_max,l∈[0,2π]

变轨点真近点角：f∈[0,2π]

轨道半长轴为a₀＝10000km，偏心率为e₀＝0.2，脉冲大小为0.5km/s。

给定速度脉冲约束后得到飞行器可达域，其中初始轨道在可达域的内部。所得到的可达域如图2所示。

2、机动点任意脉冲方向固定的三维空间可达域分析

1)机动点任意脉冲方向固定的三维空间可达域模型

对于脉冲方向固定，机动点任意的情形，固定的脉冲方向指脉冲矢量l在坐标系S₀中是一个给定的常矢量，因此式(2-9)中的矢量投影到坐标系S₀有：

(r₁)₀＝[r₀ 0 0]^T (2-47)

其中和ψ是描述脉冲方向相对于坐标系S₀的两个给定的方向角。相似地，将公式(2-48)转换到坐标系S₁中并令其z₁分量为0，得到约束方程如下：

Acos(f+φ)+B＝0，f∈[υ₀-π,υ₀+π] (2-49)

其中，

方程(2-50)中只有一个未知变量f，解存在的条件是|A|≥|B|。当|A|>|B|时，方程存在两个解。由于脉冲方向给定，p方向的可达目标位置完全由f确定，因此这个条件意味着p方向上只有两个单独的点，而不再有一条可达线段。在三维空间中旋转p，每个方向的两个可达点连接起来将形成一个二维的可达曲面，而不再是三维可达域。

因此，类似(2-12)至(2-24)，可以得到机动点任意脉冲方向固定时的可达域位置方程为：

r＝r(κ,υ)p(κ,υ)

其中

p(κ,υ)＝[cosυcosκ sinυcosκ sinκ]^T

这里A、B和φ的定义见式(2-50)。

可达域边界方程：

r＝r(f^*；κ,υ)p(κ,υ)

其中f^*由以下方程求取：

2)机动点任意脉冲方向固定的三维空间可达域分析

轨道半长轴a₀＝10000km，偏心率e₀＝0.3。

脉冲大小设为Δv＝0.5km/s，两个方向角分别取ψ＝70°。

所得到的可达域如图3所示。由于脉冲方向总是指向初始轨道外侧，所以可达域总是在初始轨道外侧，如图4所示。可达域只有图3中环形管状包络面，其内部空间不是。

3、机动点固定脉冲方向任意的三维空间可达域分析

1)机动点固定脉冲方向任意的三维空间可达域模型

机动点固定，即真近角f为常值，由关系式(2-13)可直接计算得到脉冲方向角β的值。而使可达目标矢径r达到极致的条件只剩下公式(2-31)。此种情况下，可达域位置方程为：

r＝r(α；κ,υ)p(κ,υ)

其中

p(κ,υ)＝[cosυcosκ sinυcosκ sinκ]^T

h＝r₀v_u，

可达域边界方程：

r＝r(α^*；κ,υ)p(κ,υ)

其中α^*由以下方程求取：

2)机动点固定脉冲方向任意的三维空间可达域分析

轨道半长轴a₀＝10000km，偏心率e₀＝0.3。

脉冲大小设为Δv＝0.5km/s，固定机动点的真近角取为f＝120°。

所得到的可达域如图5所示。

4、有时间约束的脉冲三维空间可达域分析

1)有时间约束的脉冲三维空间可达域模型

令t₀为飞行器在初始轨道上的初始时刻；t_d为可能机动的最大时间；脉冲施加时间为t₁，则有t₀≤t₁≤t_d；脉冲施加方向角为(α,β)。定义Δt_d＝t_d-t₀、Δt₁＝t₁-t₀，则有0≤Δt₁≤Δt_d。根据脉冲推力的特点，空间目标在t₁时刻施加方向角为(α,β)的脉冲推力后到达t_d时刻的可达域可按照以下流程计算。

设空间目标的初始轨道六根数为：半长轴a₀、偏心率e₀、轨道倾角i₀、升交点赤经Ω₀、近地点角距ω₀和真近点角f₀。

首先列出求平近点角、偏近点角、真近点角的公式：

M＝nt (4-1)

E-esin E＝M (4-2)

其中t为时间，n为轨道平均角速度，M为轨道平近点角，E为轨道偏近点角，e为轨道偏心率，f为轨道真近点角。

初始轨道平均角速度为初始真近点角f₀已知，因此由公式(4-3)和公式(4-2)可以求得初始平近点角M₀，再由公式(4-1)可以求得机动时刻平近点角M₁＝M₀+n₀Δt₁，然后通过公式(4-2)和公式(4-3)求得机动时刻真近点角f₁。

则t₁时刻空间目标地心距为

式中

进而得到目标在地心惯性坐标系下的位置矢量和速度矢量为：

r_E＝C₁[r₁ cos f₁ r₁ sin f₁ 0]^T (4-5)

式中C₁为轨道坐标系到地心惯性坐标系的坐标转换矩阵，具体形式如下

C₁＝M₃(-Ω₀)M₁(-i₀)M₃(-ω₀) (4-7)

式中，

此时，施加的脉冲推力的速度增量在地心惯性坐标系下表示为

Δv_E＝C₂[Δv cosβcosα Δv cosβsinα Δv sinβ]^T (4-8)

式中，C₂为当地轨道坐标系到地心惯性坐标系的坐标转换矩阵，具体形式如下

C₂＝M₃(-Ω₀)M₁(-i₀)M₃(-θ) (4-9)

式中，θ＝ω₀+f₁。

因此脉冲机动后的位置矢量和速度矢量为

然后通过求解开普勒问题求得Δt＝t_d-t₁后的位置矢量r_d和速度矢量v_d。

为了求得飞行器所有可能的机动情况，即考虑脉冲施加时间和脉冲施加方向的任意性，令脉冲施加时间t₁在区间(t₀,t_d)内均匀分布(即Δt₁在(0,Δt_d)内均匀分布)；在对应的每个t₁时刻，令脉冲施加方向角α在区间[0,2π)内均匀分布，β在[-π/2,π/2]内均匀分布。利用以上计算流程，通过对Δt₁、α和β不同选取，即可得到一组目标在t_d时刻的位置矢量。这些位置矢量组成的集合就是飞行器在有限时间Δt_d后的可达域，该可达域包含了飞行器所有可能运动到的位置。

2)有时间约束的脉冲的三维空间可达域分析

取初始轨道为轨道高度H₁＝500km的近圆轨道，轨道的六根数如表1所示。速度脉冲大小Δv取0.1km/s。考虑的Δt_d分别取1000s、2000s、3000s、4000s。

表1飞行器初始轨道六根数

地球惯性坐标系下的三维空间可达域如图6所示。图中椭圆形轨道为飞行器初始轨道，Ⅰ处的点为飞行器的初始位置(仅表示位置，没有大小)，Ⅱ区域为1000s后的可达域，Ⅲ区域为2000s后的可达域，Ⅳ区域为3000s后的可达域，Ⅴ区域为4000s后的可达域。

二、有限推力变轨方式的基于轨道六根数的三维空间可达域分析

1、有限推力变轨方式的三维空间可达域模型

由于经典轨道要素高斯摄动方程在轨道偏心率e＝0或e＝1和轨道倾角i＝0或π时存在奇异，而有限推力转移过程中轨道根数连续变化，出现上述奇异情况可能性很大，因此这里采用消除奇异的改进春分点轨道根数。这组轨道根数只有在轨道倾角为180°时才会产生奇异，而这种情况在有限推力轨道转移问题中几乎不存在。改进春分点轨道根数定义如下：

这里a为轨道半长轴，e为轨道偏心率，ω为近地点幅角，Ω为升交点赤经，i为轨道倾角，θ为升交点幅角。

对于一般的空间转移轨道，轨道倾角i≠0，因此上式定义的改进春分点轨道根数h和k一般也不等于零。从而得到相应的动力学模型为：

式中w＝1+f cos L+g sin L，s²＝1+h²+k²；μ为地球引力常数；f_r,f_t,f_n分别代表轨道平面上径向、切向和法向的推力加速度大小，I_sp发动机比冲，T为发动机推力，g₀为重力加速度。

以经典轨道根数和改进春分点轨道根数描述的动力学方程在形式上有一定的相似之处，它们的矩阵形式都可表示为：

式中f_p为摄动加速度，m为质量。采用改进春分点轨道根数时x＝[p,f,g,h,k,L]^T，矩阵M的形式如下：

指向推力方向的单位矢量α满足约束α^Tα＝1，定义哈密尔顿函数为：

式中λ与λ_m为协态向量。利用控制方程可得到最优推力方向为：

基于庞特里亚金极大值原理，对于推力大小可调的推进器而言，其发动机控制策略的开关函数为：

将哈密尔顿函数对状态求偏导数，得到协状态方程，形式如下：

则基于改进春分点根数的有限推力轨道的协状态方程为：

质量的协态方程为：

通过式(5-9)可见，基于改进春分点根数的有限推力轨道协状态方程由两部分组成，第一个括号内的各项由与二体动力学模型相关的量组成，而第二个括号内的各项则由摄动模型确定，在摄动模型复杂或是考虑的摄动因素较多的情况下这部分的形式十分复杂。由于M是一个6×3的矩阵，而要分别对6个状态变量求偏导，因此会产生108个偏导数表达式，由于数量较多就不在此列出。

动力学方程(5-2)、协状态方程(5-9)、边界条件组成一个两点边值问题。求解两点边值问题可得到基于轨道根数的有限推力飞行器最优飞行轨迹，同时获得协状态的变化曲线。通过控制方程(5-7)可确定相应的最优控制曲线。

2、有限推力变轨方式的三维空间可达域分析

1)固定时间不考虑燃耗约束的可达域分析

该情况不考虑转移过程中的燃耗约束，即给定时间范围内，发动机一直处于开机的状态。采用改进春分点轨道根数动力学模型(5-2)中f_r,f_t,f_n分别是推力在径向、切向和法向所产生的加速度。推力方向也可以通过俯仰和偏航控制角来描述，具体形式如式(5-12)：

其中，俯仰控制角α为推力矢量在轨道平面内的投影与当地水平面之间的夹角；偏航控制角β为推力矢量与密切轨道平面之间的夹角。

因此，该问题对应的可达范围则由控制角α和β共同决定，穷举这两个参数进行打靶，可以得到整个有限推力的空间可达范围。

假设飞行器初始质量为m₀＝4500kg，飞行时间t_tof＝1525s，发动机最大推力T_max＝1000N，发动机比冲I_sp＝311.2s。初始轨道的经典轨道根数如表2，打靶结果边界参数如表3所示。

表2初始轨道的经典轨道根数

表3打靶结果边界参数

有限推力转移轨道可达域空间示意图如图7。所有示意图的绘制均选择惯性坐标系，下同。

2)固定时间考虑燃耗约束的可达域分析

该情况考虑转移过程中的燃耗约束，即给定时间范围内，发动机不能一直处于开机的状态。因此，该问题对应的可达范围不止需要考虑控制角α和β，仍需考虑发动机开机时间。同理，可达范围边界必定消耗最大燃耗，发动机最长开机时间可由最大燃耗求得，即：

这里m_p为最大燃料可消耗量，T_max为发动机最大推力。穷举α、β和不同的开机时间段(总开机时间为t_max)进行打靶，可以得到整个有限推力的空间可达范围。

假设飞行器初始质量为m₀＝4500kg，飞行时间t_tof＝3000s，最大燃耗m_p＝500kg，发动机最大推力T_max＝1000N，发动机比冲I_sp＝311.2s。初始轨道的经典轨道根数如表4，打靶结果边界参数如表5所示。

表4初始轨道的经典轨道根数

表5打靶结果边界参数

有限推力转移轨道可达域空间示意图如图8所示。

本发明未详细说明部分属本领域技术人员公知常识。

Claims (4)

1.一种空天飞行器三维空间可达域分析方法，其特征在于：包括如下内容：

(1)确定空天飞行器的变轨方式，如果为脉冲方式变轨则进入步骤(2)，如果为有限推力变轨方式则进入步骤(5)；

(2)按照要求获得可达域；判断是否考虑时间约束，如果不考虑时间约束则进入步骤(3)，如果考虑时间约束则进入步骤(4)；

(3)建立机动点位置、脉冲方向不同状态组合下的基于轨道六根数的三维空间可达域模型，获得脉冲方式变轨下的三维空间可达域；

(4)建立具有时间约束的脉冲变轨三维空间可达域模型，获得固定时间长度下脉冲方式变轨下的三维空间可达域；

(5)建立基于改进春分点轨道根数的有限推力变轨三维空间可达域模型，获得有限推力方式变轨下的三维空间可达域。

2.根据权利要求1或2所述的一种空天飞行器三维空间可达域分析方法，其特征在于：步骤(3)中基于轨道六根数的三维空间可达域模型由可达域位置方程与边界方程组成，具体为：

可达域位置方程：

r＝r(f,α；κ,υ)p(κ,υ)

其中

p(κ,υ)＝[cosυcosκ sinυcosκ sinκ]^T

可达域边界方程：

r＝r(f^*,α^*；κ,υ)p(κ,υ)

其中边界上的真近点角f^*，能达到边界的推力方向α^*，由以下方程求取f的值作为f^*，求取α值作为α^*：

其中κ为单位矢量高度角，υ为升交点幅角，f为真近点角，(α,β)为脉冲施加方向角，Δv的方位角与高度角，h为飞行器动量矩，Δθ代表初始位置与目标位置的转角，v_u和v_r分别为速度在水平方向与径向的分量，e₀为初始轨道的离心率，p₀为初始轨道的半通径，μ为地心引力常数。

3.根据权利要求1或2所述的一种空天飞行器三维空间可达域分析方法，其特征在于：步骤(4)中具有时间约束的脉冲变轨三维空间可达域模型为：

其中C₁为轨道坐标系到地心惯性坐标系的坐标转换矩阵，C₂为当地轨道坐标系到地心惯性坐标系的坐标转换矩阵，r₁为脉冲施加时刻的飞行器地心距，f₁为脉冲施加时刻的飞行器真近点角，e₀为飞行器初始轨道偏心率，p₀为半通径，μ为地球引力常数，Δv为施加脉冲的大小，(α,β)为脉冲施加方向角。

4.根据权利要求1或2所述的一种空天飞行器三维空间可达域分析方法，其特征在于：步骤(5)中改进春分点轨道根数的有限推力变轨三维空间可达域模型为：

其中a为轨道半长轴，e为轨道偏心率，ω为近地点幅角，Ω为升交点赤经，i为轨道倾角，θ为升交点幅角；w＝1+fcosL+gsinL，s²＝1+h²+k²；μ为地球引力常数；f_r,f_t,f_n分别代表轨道平面上径向、切向和法向的推力加速度大小，I_sp发动机比冲，T为发动机推力，g₀为重力加速度，α为推力矢量俯仰控制角，β为推力矢量偏航控制角，m_p为最大燃料可消耗量，T_max为发动机最大推力。