CN109094818A - 一种空间飞行器的远程交会制导方法 - Google Patents

Abstract

本发明公开一种空间飞行器的远程交会制导方法，包括：根据飞行器的飞行参数和交会位置，计算得到初始的增益速度；根据飞行器的飞行参数、增益速度和交会位置，外推计算得到交会偏差；根据飞行器的飞行参数、终端位置和速度计算当前时刻位置和速度的偏导数，得到偏导数矩阵；利用偏导数矩阵和交会偏差计算出修正的增益速度。本发明提供一种空间飞行器的远程交会制导方法，该方法对于空间飞行器的远程交会制导可兼顾计算精度和计算速度，可以确保远程(数小时级)交会制导计算精度为百米级，同时计算时间与二体解析算法相当，可以在线使用，方法简便可行。

Description

一种空间飞行器的远程交会制导方法

技术领域

本发明涉及空间远程制导技术领域。更具体地，涉及一种空间飞行器的远程交会制导方法。

背景技术

空间飞行器的交会制导方法受终端约束的影响很大，其中一类终端约束是飞行时间固定和终端位置固定。针对该类约束，现有的制导方法是基于二体轨道理论的兰伯特算法，该算法在近程或中程交会制导时精度尚可，但应用到远程交会制导时则由于交会误差较大而无法满足要求，因此在远程交会制导时需要对兰伯特算法进行修正。

目前已有算法中，轨道外推的精确计算一般采用积分算法，而且为了保证计算精度，常需要选用较小的计算步长，这导致远程交会制导的计算周期过长而无法实现飞行器的在线使用；偏导数矩阵的计算一般采用求差法，即先计算初始时刻位置、速度的轨道外推结果，再计算位置矢量的三个分量、速度矢量的3个分量各自发生微小改变后的轨道外推结果，对外推结果求差，得到偏导数矩阵，这种方法需要进行多达7次轨道外推，计算周期较长，无法满足要求。

因此，目前现有的计算方法仍然无法满足交会制导的计算需求。

发明内容

为了解决上述问题至少之一，本发明的一个目的在于提供一种空间飞行器的远程交会制导方法。

为达到上述目的，本发明采用下述技术方案：

一种空间飞行器的远程交会制导方法，包括：

根据飞行器的飞行参数和交会位置，计算得到初始的增益速度；

根据飞行器的飞行参数、增益速度和交会位置，计算得到交会偏差；

根据飞行器的飞行参数、终端位置和速度计算当前时刻位置和速度的偏导数，得到偏导数矩阵；

利用偏导数矩阵和交会偏差计算出修正的增益速度。

优选地，所述飞行参数包括飞行器的飞行时刻、位置、速度、交会时刻以及交会位置。

优选地，所述初始的增益速度根据二体轨道理论、采用兰伯特算法计算。

优选地，所述交会偏差为采用中间轨道法，进行外推计算得到交会时刻飞行器的位置，并与已知的交会位置进行求差计算而得到。

优选地，所述偏导数矩阵用于计算起始时刻位置或速度的微小摄动对外推时刻位置、速度的影响。

优选地，所述初始的增量速度采用第一公式计算，所述第一公式为：

其中，所述z为普适变量，x(z)、S(z)和y(z)为z的函数，A为根据飞行器起始点和终点的位置得出的常值，μ是地球引力常数。

优选地，所述x(z)、S(z)和y(z)由第一方程组和第二方程组计算得到，所述第一方程组为：

所述第二方程组为：

优选地，所述外推算法的计算包括：

将当前时刻的J2000坐标系下的位置矢量转化为为椭球坐标系下的位置矢量；

计算Jacobi前三个常数，其中，第一常数为：

第二常数为：

第三常数为：

计算第一函数和第二函数的四次式，所述四次式包括：

F(ρ)＝μ[c²p₀(1-S₀)+(ρ²+c²)(γ₀ρ²+2ρ-p₀)]，

F(ρ)＝μγ₁(γρ²+2ρ-p)(ρ²+2A₁ρ+B₁)，

G(η)＝μ[-p₀(1-S₀)+(1-η²)(p₀+2δη+c²γ₀η²)]，

G(η)＝μS₁p₀(S+2Pη-η²)(1+P₁η-Q₁η²)；

初始化六个积分常数，其中，所述六个积分常数包括第一至第六积分常数，其计算公式为：

N₂＝D₁[u+k₁T₂/2+k₁T₂/2+3k₁ ²T₄/8+5k₁ ³T₆/16]，

其中，T₀＝u，T₁＝1-cosu，T_k＝[(k-1)T_k-2-cosusin^k-1u]/k,k＝2,...,6；

计算Jacobi后三个常数，其中，所述后三个常数的计算公式为：

t_i+β₁＝R₁(ρ_i)+c²N₁(η_i)，

β₂＝-α₂R₂(ρ_i)+α₂N₂(η_i)，

β₃＝φ_i+c²α₃R₃(ρ_i)-α₃N₃(η_i)，

其中，

其中，

将椭球坐标系下的位置矢量通过第一转化公式转换为J2000坐标系的位置矢量，

所述第一转化公式为：

其中，

优选地，所述偏导数矩阵采用解析方程组计算，所述解析方程组为：

式中，

其中，R₀和V₀为t₀时刻位置、速度的大小，R和V为t_f时刻位置、速度的大小。

本发明的有益效果如下：

本发明提供一种空间飞行器的远程交会制导方法，该方法对于空间飞行器的远程交会制导可兼顾计算精度和计算速度，可以确保远程(数小时级)交会制导计算精度为百米级，同时计算时间与二体解析算法相当，可以在线使用，方法简便可行。

附图说明

下面结合附图对本发明的具体实施方式作进一步详细的说明。

图1示出本发明的一个实施方式提供的一种空间飞行器的远程交会制导方法流程示意图。

图2示出现有技术中采用兰伯特算法全程制导精度变化曲线示意图。

图3示出现有技术中采用兰伯特算法最后阶段制导精度变化曲线示意图。

图4示出本发明的一个实施方式提供的一种空间飞行器的远程交会制导方法全程制导精度变化曲线示意图。

图5示出本发明的一个实施方式提供的一种空间飞行器的远程交会制导方法最后阶段制导精度变化曲线示意图。

具体实施方式

为了更清楚地说明本发明，下面结合优选实施例和附图对本发明做进一步的说明。附图中相似的部件以相同的附图标记进行表示。本领域技术人员应当理解，下面所具体描述的内容是说明性的而非限制性的，不应以此限制本发明的保护范围。

在本发明的一个实施方式中，提供一种空间飞行器的远程交会制导方法，包括：

S1：根据飞行器的飞行参数和交会位置，计算得到初始的增益速度。

具体的，以表示起始点和终点的飞行器位置、速度，根

据二体理论，这四个矢量之间的相互关系见公式(1)：

其中f,g是普适变量z的函数；z＝(ΔE)²，即航天器轨道偏近点角变化的平方。由公式(1)，可以推导得到航天器在处的速度

由公式(2)分析可以得到，只要能计算出f,g的值，就可以求出 的值，从而解决问题。而f,g是普适变量z的函数，f(z),g(z)的具体表达公式由下列一组方程式求出：

r₁、r₂为已知参数，μ为地球引力常数，A可以由已知条件直接求出，求解公式如公式(4)所示：

经过以上分析以后，问题就转化为求解y(z)函数的值。下面给出了y(z) 函数的求解公式，如公式(5)所示：

联立公式(3)、(4)、(5)，就可以根据普适变量z的值,已知条件的值求解f(z),g(z)，随即就能求出飞行器在处的速度因此兰伯特问题求解的关键就归结为对普适变量z值的求解。

下面给出了兰伯特问题中飞行时间t与普适变量z的函数关系式，如公式(6)所示：

由关系式(6)，满足兰伯特问题飞行时间ΔT的航天器运动方程满足下列方程式

其中x(z)、S(z)、y(z)为普适变量z的函数，参见关系式(5)、(6)；A是由已知条件得出的常值，参见关系式(4)；μ是地球引力常数。因此兰伯特问题求解的关键就转化为(7)式对z的求解。

S2：根据飞行器的飞行参数、增益速度和交会位置，计算得到交会偏差。

具体的，采用中间轨道法，外推得到交会时刻飞行器位置，与事先已知的交会位置进行求差计算，得到交会偏差。

下面给出解析外推算法的步骤：

1.将t_i时给定的J2000坐标系的位置矢量x(t_i)，转换为椭球坐标系

其中x(t_i)由x_i，y_i，z_i，和组成，并且

d＝(r_i ²-c²)+δ(2z_i+δ) (9)

r_i ²＝x_i ²+y_i ²+z_i ² (10)

D²＝(ρ_i ²+c²)(1-η_i ²) (14)

2.计算Jacobi前三个常数(α₁，α₃，α₂)

3.F(ρ)和G(η)的四次式

F(ρ)＝μ[c²p₀(1-S₀)+(ρ²+c²)(γ₀ρ²+2ρ-p₀)] (18)

F(ρ)＝μγ₁(γρ²+2ρ-p)(ρ²+2A₁ρ+B₁) (19)

G(η)＝μ[-p₀(1-S₀)+(1-η²)(p₀+2δη+c²γ₀η²)] (20)

G(η)＝μS₁p₀(S+2Pη-η²)(1+P₁η-Q₁η²) (21)

4.初始化六个积分系数(R₁，R₂，R₃，N₁，N₂，N₃)

计算R₁，R₂，R₃，需要A_k和W_k，其中A₀＝1，通过因数分解得到A₁及其它，W₀＝V₀＝W为真近地点角，W₁＝(W+eV₁)/p，V₁＝sinW。

N₂＝D₁[u+k₁T₂/2+k₁T₂/2+3k₁ ²T₄/8+5k₁ ³T₆/16] (26)

计算N₁，N₂，N₃，需要C_1k，C_2k和T_k，其中

T₀＝u (32)

T₁＝1-cosu (33)

T_k＝[(k-1)T_k-2-cosusin^k-1u]/k,k＝2,...,6 (34)

注意，在t_i时，需要精确的确定椭圆积分振幅u和通用变量其中是 t_i时Vinti的解。

5.计算Jacobi后三个常数(β₁，β₂，β₃)

t_i+β₁＝R₁(ρ_i)+c²N₁(η_i) (35)

β₂＝-α₂R₂(ρ_i)+α₂N₂(η_i) (36)

β₃＝φ_i+c²α₃R₃(ρ_i)-α₃N₃(η_i) (37)

使用初始环境和已初始化的系数。在程序中，

时间β₁＝-T

(38)

近地点幅角β₂＝ω

(39)

升交点赤经β₃＝Ω

(40)

6.使用Jacobi常数(α₁，α₃，α₂，β₁，β₂，β₃)替换运动方程，求解给定时间t_f时的ρ_f，η_f和φ_f

t_f+β₁＝R₁(ρ_f)+c²N₁(η_f) (41)

β₂＝-α₂R₂(ρ_f)+α₂N₂(η_f) (42)

β₃＝φ_i+c²α₃R₃(ρ_f)-α₃N₃(η_f) (43)

第一个运动方程是开普勒方程的通式。只要计算使用椭球坐标系，相应的值由下式计算：

其中

7.将时间t_f时椭球坐标系下的位置矢量X(t_f)，转换为J2000坐标系的位置矢量x(t_f)

其中

计算得到外推时刻的位置后，与事先已知的交会位置进行求差计算，得到交会偏差。

S3：根据飞行器的飞行参数、终端位置和速度计算当前时刻位置和速度的偏导数，得到偏导数矩阵。

具体的，该矩阵用于计算起始时刻位置或速度的微小摄动对外推时刻位置、速度的影响。

下面给出偏导数矩阵的解析表达式：

设t₀时刻飞行器位置和速度矢量为r₀和v₀，则二体轨道理论中t_f时刻飞行器的位置和速度矢量为：

式中，f,g,称为拉格朗日系数。由式(50)可得

式中，I为3×3的单位矩阵；P_rv为t₀时刻速度变化对t_f时刻位置影响的偏导数矩阵。同理可得偏导数矩阵P_rr，P_vr和P_vv。

偏导数矩阵具体解析表达式如下：

式中，

R₀和V₀为t₀时刻位置、速度的大小，R和V为t_f时刻位置、速度的大小。根据R₀、V₀、t_f-t₀可以计算得到半长轴a、偏近点角变化值E-E₀，进而计算得到α、χ，代入U_n(χ,α)公式得到U₅、U₄、U₂和C。最后代入方程中可以得到偏导数矩阵

S4：利用偏导数矩阵和交会偏差计算出修正的增益速度。

具体的，得到偏导数矩阵后，需要根据交会偏差修正兰伯特算法计算出的增益速度。

计算步骤如下：

根据初始时刻位置速度偏差计算外推时刻位置、速度变化量的方程为：

式中，δr_f、δv_f分别为t₀时刻飞行器位置和速度矢量分别改变δr₀、δv₀所导致t_f时刻飞行器位置和速度矢量的改变量。其中，仅当δr₀、δv₀为小量时式(57)成立。

若兰伯特算法计算出的增益速度为Δv，考虑摄动影响后的零控交会偏差为δr_f(该量为小量)，假设此时在Δv基础上令飞行器产生速度增量δΔv，δΔv正好产生δr_f的交会偏差，则最终交会偏差将会变为零。

由式(57)可得δr_f＝P_rrδr₀+P_rvδv₀，在速度矢量改变的瞬间，δr₀＝0，则δr_f＝P_rvδv₀，因此故实际增益速度应为

本发明提供一种空间飞行器的远程交会制导方法，该方法对于空间飞行器的远程交会制导可兼顾计算精度和计算速度，可以确保远程(数小时级)交会制导计算精度为百米级，同时计算时间与二体解析算法相当，可以在线使用，方法简便可行。

下面以某空间飞行器为例，进行说明。

该空间飞行器的当前飞行参数和目标参数见表1所示，飞行时间固定为 9700s。

表1空间飞行器当前飞行参数和目标参数

分别采用兰伯特算法和本算法进行了全程制导计算，利用STK对制导结束时的制导误差进行了评估，见表2所示，飞行过程中制导精度的变化情况见图2至图5所示，两种算法的计算速度见表3所示。

表2制导结束时的误差(与STK比对)

表3制导计算时间

由上所述，可见该方法的计算时间与兰伯特算法相当，但计算精度大大提高。

显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定，对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动，这里无法对所有的实施方式予以穷举，凡是属于本发明的技术方案所引伸出的显而易见的变化或变动仍处于本发明的保护范围之列。

Claims (9)

1.一种空间飞行器的远程交会制导方法，其特征在于，包括：

根据飞行器的飞行参数和交会位置，计算得到初始的增益速度；

根据飞行器的飞行参数、增益速度和交会位置，外推计算得到交会偏差；

根据飞行器的飞行参数、终端位置和速度计算当前时刻位置和速度的偏导数，得到偏导数矩阵；

利用偏导数矩阵和交会偏差计算出修正的增益速度。

2.根据权利要求1所述方法，其特征在于，所述飞行参数包括飞行器的飞行时刻、位置、速度、交会时刻以及交会位置。

3.根据权利要求1所述方法，其特征在于，所述初始的增益速度根据二体轨道理论、采用兰伯特算法计算。

4.根据权利要求1所述方法，其特征在于，所述交会偏差为采用中间轨道法，进行外推计算得到交会时刻飞行器的位置，并与已知的交会位置进行求差计算而得到。

5.根据权利要求1所述方法，其特征在于，所述偏导数矩阵用于计算起始时刻位置或速度的微小摄动对外推时刻位置、速度的影响。

6.根据权利要求4所述方法，其特征在于，所述初始的增益速度采用第一公式计算，所述第一公式为：

其中，所述z为普适变量，x(z)、S(z)和y(z)为z的函数，A为根据飞行器起始点和终点的位置得出的常值，μ是地球引力常数。

7.根据权利要求6所述方法，其特征在于，所述x(z)、S(z)和y(z)由第一方程组和第二方程组计算得到，所述第一方程组为：

所述第二方程组为：

8.根据权利要求1所述方法，其特征在于，所述外推算法的计算包括：

将当前时刻的J2000坐标系下的位置矢量转化为为椭球坐标系下的位置矢量；

计算Jacobi前三个常数，其中，第一常数为：

第二常数为：

第三常数为：

计算第一函数和第二函数的四次式，所述四次式包括：

F(ρ)＝μ[c²p₀(1-S₀)+(ρ²+c²)(γ₀ρ²+2ρ-p₀)]，

F(ρ)＝μγ₁(γρ²+2ρ-p)(ρ²+2A₁ρ+B₁)，

G(η)＝μ[-p₀(1-S₀)+(1-η²)(p₀+2δη+c²γ₀η²)]，

G(η)＝μS₁p₀(S+2Pη-η²)(1+P₁η-Q₁η²)；

初始化六个积分常数，其中，所述六个积分常数包括第一至第六积分常数，其计算公式为：

N₂＝D₁[u+k₁T₂/2+k₁T₂/2+3k₁ ²T₄/8+5k₁ ³T₆/16]，

其中，T₁＝1-cosu，T_k＝[(k-1)T_k-2-cosusin^k-1u]/k,k＝2,...,6；

计算Jacobi后三个常数，其中，所述后三个常数的计算公式为：

t_i+β₁＝R₁(ρ_i)+c²N₁(η_i)，

β₂＝-α₂R₂(ρ_i)+α₂N₂(η_i)，

β₃＝φ_i+c²α₃R₃(ρ_i)-α₃N₃(η_i)，

其中，

其中，

将椭球坐标系下的位置矢量通过第一转化公式转换为J2000坐标系的位置矢量，

所述第一转化公式为：

其中，

9.根据权利要求1所述方法，其特征在于，所述偏导数矩阵采用解析方程组计算，所述解析方程组为：

式中，

其中，R₀和V₀为t₀时刻位置、速度的大小，R和V为t_f时刻位置、速度的大小。