CN109752005B - 一种基于精确轨道模型的航天器初轨确定方法 - Google Patents

Abstract

本发明公开了一种基于精确轨道模型的航天器初轨确定方法，在精确轨道模型和拉普拉斯方程基础上建立了地基雷达测轨数据、天基光学测角类测轨数据、GNSS、北斗类测轨数据的加权最小二乘条件方程，以直接合成法或联邦合成法两类方法融合各种不同类型的测轨数据，迭代求解出初始位置、速度的改进值。本发明一种基于精确轨道模型的航天器初轨确定方法能够直接将完备或非完备外弹道观测数据、位置速度分量类观测数据、天基测角类数据融合起来进行星箭分离初轨计算，产生唯一的初轨输出。

Description

一种基于精确轨道模型的航天器初轨确定方法

技术领域

本发明属于航天导航技术领域，涉及一种基于精确轨道模型的航天器初轨确定方法，具体涉及一种基于精确轨道模型的、融合多类观测数据源的单位矢量法初轨确定方法。适用于航天器星箭分离后、及其它无动力飞行阶段的轨道的快速确定。

背景技术

航天器星箭分离后初轨的计算是航天器发射过程的最重要内容之一，是评价发射成功与否及进行航天器后续轨道控制的依据。星箭分离后初轨计算的实时性要求较高并且没有真实的参考轨道等特点，使得其与事后精密定轨在计算方法上存在较大差异。

陆本魁等公开了适用于航天器测轨计算的单位矢量法初轨计算方法。给出了对应常规地面观测数据的单位矢量法条件方程，利用传统f和g级数展开低阶项来近似轨道动力学模型。单位矢量法初轨计算方法由于避免了常规微分轨道改进方法要求复杂偏导数计算的问题，在初轨计算时具有计算效率高、方法稳定性好、收敛速度快等显著优点。茅永兴等人公开了航天发射任务入轨段遥外测数据融合定轨方法，此方法采用基于微分轨道改进基本原理，给出了适用于不同类型测轨数据联合定轨的动力学条件方程和方程权重计算方法，该方法中基于微分轨道改进原理建立各类测轨数据的动力学条件方程，需要推导和计算大量偏导数矩阵。

考虑到航天器无动力飞行轨道可进行精确建模，以及针对上述方法的不足，本发明以有摄单位矢量法和精确轨道模型为基础，融合多类测轨数据进行了初轨确定。

发明内容

本发明的目的是提供一种基于精确轨道模型的航天器初轨确定方法，直接将多类型的数据进行融合使用，产生唯一初轨输出。

本发明采用的技术方案是,一种基于精确轨道模型的航天器初轨确定方法，具体按照下述步骤进行：

步骤1，对航天器星箭分离后的轨道进行初始估计；

步骤2，融合不同的类型的测轨数据，在有摄单位矢量法和精确轨道模型基础上，基于迭代加权最小二乘法计算初始位置和速度的改进值；

步骤3，若初始位置速度前后两次迭代改进值满足收敛条件，则将初始位置、速度的迭代改进值；若不满足收敛条件，则重复步骤2至收敛为止。

本发明的特点还在于：

步骤2中根据直接合成法建立最小二乘条件方程具体按照下述步骤进行：

步骤2.1，对每一帧观测数据，以初值t₀、r₀、

进行动力学模型积分到当前数据帧时刻，得到其对应的位置速度估计值；

步骤2.2，建立条件方程矩阵；

建立最小二乘条件方程：

步骤2.3，分别计算地基雷达测轨数据、天基光学测角类测轨数据和GNSS、北斗类测轨数据计算得到条件方程各项的系数；

步骤2.4，将每个条件方程进行法化计算并将法化后的系数矩阵累加到系统条件方程中，并将累加了系数矩阵的条件方程组建条件方程组；

步骤2.5，判断测轨数据是否遍历完毕，若遍历完毕，得到条件方程矩阵；

若没有遍历完毕，则取下一帧观测数据，并重复步骤2.1-2.4，至每帧测轨数据遍历完毕。

步骤2中联邦合成法建立最小二乘条件方程组具体按照下述步骤进行：

步骤2.a，对每一帧观测数据，以初值t₀、r₀、

进行动力学模型积分到当前数据帧时刻，得到其对应的位置速度估计值

步骤2.b，根据地基雷达测轨数据列写地基雷达条件方程组，并对地基雷达方程组求解得到地基雷达类初始位置和地基雷达类初始速度向量，并对地基雷达条件方程组进行方差统计；

根据天基光学测角类测轨数据列写天基光学测角条件方程组，并对天基光学测角方程组求解得到天基光学测角类初始位置和天基光学测角类初始速度向量，并对天基光学测角条件方程组进行方差统计；

根据GNSS、北斗类测轨数据列写地基雷达条件方程组，并对GNSS、北斗类方程组求解得到GNSS、北斗类初始位置和GNSS、北斗类初始速度向量，并对GNSS、北斗类方程组进行方差统计；

步骤2.c，根据地基雷达条件方程组的方差对地基雷达类初始位置和地基雷达类初始速度向量进行加权；

根据天基光学测角条件方程组的方差对天基光学测角类初始位置和天基光学测角类初始速度向量进行加权；

根据GNSS、北斗类方程组的方差对GNSS、北斗类初始位置和GNSS、北斗类初始速度向量进行加权；

将加权后的地基雷达类初始位置、天基光学测角类初始位置和GNSS、北斗类初始位置进行合成，得到新的初始位置

将加权后的地基雷达类初始速度、天基光学测角类初始速度和GNSS、北斗类初始速度向量进行合成，得到新的初始速度向量。

步骤2中在最小二乘迭代计算过程中，以精确轨道模型来计算位置速度估计值。

步骤3中具体按照下述方法判断新的初始位置和新的初始速度向量是否收敛：

给定收敛门限ξ，相邻两次得到的新的初始位置和新的初始速度向量分别为

和

若其满足

则得到的收敛，若不满足，则不收敛。

本发明的有益效果是

1)本发明一种基于精确轨道模型的航天器初轨确定方法能够直接将完备或非完备外弹道观测数据、位置速度分量类观测数据、天基测角类数据融合起来进行星箭分离初轨计算，产生唯一的初轨输出；

2)本发明一种基于精确轨道模型的航天器初轨确定方法中可直接使用单测角、单测距等非完备外弹道观测数据，有效利用了各类观测数据；

3)本发明一种基于精确轨道模型的航天器初轨确定方法中，算法采用高精度数值积分方法进行有摄状态外推，有效提高了计算精度；

4)本发明一种基于精确轨道模型的航天器初轨确定方法能够在在外测数据质量较差时能达到与较高精度GNSS、北斗类观测数据相当的结果精度。

5)本发明一种基于精确轨道模型的航天器初轨确定方法，具有一定的野值剔除及容错能力。

附图说明

图1是本发明一种基于精确轨道模型的航天器初轨确定方法中直接合成法的流程图；

图2是本发明一种基于精确轨道模型的航天器初轨确定方法中联邦合成法的流程图。

具体实施方式

下面结合附图和具体实施方式对本发明进行详细说明。

一种基于精确轨道模型的航天器初轨确定方法，如图1所示，具体按照下述步骤进行：

步骤1，系统状态初始化，J2000惯性系下进行航天器星箭分离后轨道的初始位置和航天器的初始速度向量估计；

对初始位置和初始速度向量的估计值可采用雷达测量数据确定初轨初始值方法得到，也可以直接使用位置速度类(GNSS、北斗等数据)观测数据中的一点得到；

步骤2，融合不同的类型的测轨数据，基于迭代加权最小二乘法计算初始位置和速度的改进值。

其中，步骤2可以采用直接合成法、联邦合成法进行计算。

当采用直接合成法时，如图1所示，具体按照下述步骤进行：

步骤2.1，对每一帧观测数据，以初值t₀、r₀、

进行动力学模型积分到当前数据帧时刻，得到其对应的位置速度估计值

动力学模型积分采用精确轨道模型，航天器在J2000地心惯性系下摄动运动加速度可近似为：

a＝g_E+g_M+g_S+a_S+a_g (1)

式中，g_E、g_M、g_S、a_S、a_g分别为地球引力、月球引力、太阳引力、光压摄动、大气阻力摄动等摄动力产生的加速度。对位置速度r、

的外推计算可通过对有摄加速度方程的数值法积分进行计算。而非传统的利用F、G级数展开低阶项来近似轨道动力学模型的模型迭代技术。

步骤2.2，以各类测轨数据驱动列写加权最小二乘条件方程。

在轨道计算中主要采用地基雷达测轨数据、天基光学测角类测轨数据和GNSS、北斗类测轨数据计算得到条件方程各项的系数；

为建立各类测轨数据中各类测量元素的量测模型，引入空间飞行器摄动运动状态转移方程(拉普拉斯方程)：

其中位置与速度矢量分别表示为r、

初始位置及速度向量分别为r₀、

f_r、f_v、f_n、

为拉普拉斯方程系数，可由r₀、

计算得到，具体算法参见拉普拉斯方程定义；

根据地基雷达测轨数据计算得到条件方程各项系数的过程具体按照下述步骤进行：

步骤2.3.1，在地心惯性系中，由空间飞行器、地心和观测点的几何关系可得到：

r＝M^T·ρ_S+R (4)

其中ρ_S为观测矢量在测站跟踪坐标系的投影，ρ_S＝[ρ,0,0]^T，ρ为测距值，M为从J2000地心惯性坐标系到测站跟踪坐标系的转换矩阵；R为测站原点在地心惯性坐标系的位置矢量；

将M^T的三个列向量分别标记为L_ρ、L_A、L_E，则：

M＝[L_ρ L_A L_E]^T (5)

计算测速数据：

其中，

为测站原点惯性系下的速度矢量；

根据测速数据得到：

步骤2.3.2，根据每一帧观测数据的位置和速度向量和空间飞行器、地心和观测点的几何关系得到地基雷达测轨观测数据的加权最小二乘条件方程：

其中，w_ρ为测距的权系数，w_A为方位角的权系数，w_E为仰角的权系数，

为测速数据的权系数；

步骤2.3.3，根据公式(8)、(9)、(10)和(11)计算得到最小二乘条件方程各项的系数。

根据天基光学测角类测轨数据计算得到条件方程各项的系数具体按照下述步骤进行：

步骤2.3.a，计算天基跟踪坐标系三轴：

其中，α为飞行器相对于天基平台的目标赤经，β为飞行器相对于天基平台的目标纬经，r_e＝[0,0,1]^T；

步骤2.3.b，根据式(4)、(12)、(13)和(14)得到：

L_x、L_y、L_z表示M^T矩阵的各行向量。

步骤2.3.c，向式(15)中带入位置矢量得到天基光学测角类测轨数据的小二乘条件方程：

其中，w_α为飞行器相对于天基平台的目标赤经的测量权重系数，w_β为飞行器相对于天基平台的目标纬经的测量权重系数；

步骤2.3.d，根据式(16)和式(17)计算得到最小二乘条件方程各项的系数。

根据GNSS、北斗类测轨数据计算得到条件方程各项的系数具体按照下述步骤进行：

步骤2.3.α，将GNSS、北斗类测轨数据先转换到地固系(协议地球参考系)下。

步骤2.3.β，计算地心固连系下位置的观测方程：

r_G＝M_G·r (18)

地心固连系下速度观测方程：

其中

(用L_x、L_y、L_z分别记录M_G各行向量)；

同样，将M_G的导数记为：

其中，M_G与

的值可根据地固系与J2000惯性系转换关系计算得出。

步骤2.3.γ，根据式(18)、(19)(20)得到GNSS、北斗类测轨数据下加权最小二乘条件方程：

其中，w_x为位置在x轴方向分量的权系数，w_y为位置在y轴方向分量的权系数，w_z为位置在z轴方向分量的权系数，

为速度在x轴方向分量的权系数，

为速度在y轴方向分量的权系数；

步骤2.3.δ，根据式(21)、(22)、(23)、(24)、(25)和式(26)计算得到最小二乘条件方程各项的系数。

步骤2.4，将每个条件方程进行法化计算并将法化后的系数矩阵累加到系统条件方程中；

步骤2.5，判断测轨数据是否遍历完毕，若遍历完毕，求解条件方程矩阵，得到初始位置、速度向量的一次迭代改进值；

若没有遍历完毕，则取下一帧观测数据，并重复步骤2.1-2.4，至每帧测轨数据遍历完毕。

其中，步骤2也可以使用联邦融合法完成，如图2所示，具体按照下述步骤进行：

步骤2.a，对每一帧观测数据，以初值t₀、r₀、

进行动力学模型积分到当前数据帧时刻，得到其对应的位置速度估计值

步骤2.b，对地基雷达测轨数据列写地基雷达条件方程组，求解该方程组得到地基雷达类观测数据估算的初始位置、速度向量改进值，并对地基雷达观测数据进行方差统计；

具体按照下面方法列写地基雷达条件方程组：

步骤2.b.1，建立地基雷达测轨数据加权最小二乘条件方程：

其中，w_ρ为测距的权系数，w_A为方位角的权系数，w_E为仰角的权系数，

为测速数据的权系数；

步骤2.b.3，根据公式(8)、(9)、(10)和(11)计算得到最小二乘条件方程各项的系数；

步骤2.b.4，将得到的各个最小二乘条件方程进行法化计算，并将法化系数矩阵累加到地基雷达观测数据系统条件方程组中；

判断地基雷达测轨数据是否遍历完毕，若遍历完毕，求解条件方程矩阵，得到地基雷达测轨数据计算的初始位置和速度向量改进值；

若没有遍历完毕，则取下一帧观测数据，重复步骤2.b1-2.b.4，至每帧地基雷达测轨数据遍历完毕。

根据天基光学测角类测轨数据列写天基光学测角数据条件方程组，求解该方程组得到天基光学测角类数据计算的初始位置、速度向量改进值，并对天基光学测角类数据进行方差统计；

具体按照下面方法列写天基光学测角条件方程组：

步骤2.b.a，建立天基光学测角数据的加权最小二乘条件方程：

其中，w_α为赤经数据的权重系数，w_β为赤纬数据的权重系数；

步骤2.3.c，根据式(16)和式(17)计算得到天基光学测角数据的最小二乘条件方程各项系数；

步骤2.b.d，将得到的各最小二乘条件方程进行法化计算并将法化矩阵累加到天基光学测角数据的最小二乘条件方程中；

判断天基光学测轨数据是否遍历完毕，若遍历完毕，求解条件方程矩阵，得到天基光学测角数据计算的初始位置、速度向量改进值；

若没有遍历完毕，则取下一帧观测数据，重复步骤2.b.a-2.b.d，至每帧测轨数据遍历完毕。

根据GNSS、北斗类测轨数据列写加权最小二乘条件方程组，求解该方程组得到GNSS、北斗类数据计算的初始位置、速度向量改进值，并对GNSS、北斗类测轨数据进行方差统计；

具体按照下面方法列写GNSS、北斗类测轨数据的条件方程：

步骤2.b.α，列写GNSS、北斗类测轨数据下的加权最小二乘条件方程：

其中，w_x为位置在x轴方向分量的权系数，w_y为位置在y轴方向分量的权系数，w_z为位置在z轴方向分量的权系数，

为速度在x轴方向分量的权系数，

为速度在y轴方向分量的权系数；

步骤2.3.γ，根据式(21)、(22)、(23)、(24)、(25)和式(26)计算得到GNSS、北斗类测轨数据计算的最小二乘条件方程各项的系数；

步骤2.b.δ，将得到的各最小二乘条件方程进行法化计算，并将法化矩阵累加到GNSS、北斗类测轨数据的系统最小二乘条件方程中；

判断GNSS、北斗类测轨数据是否遍历完毕，若遍历完毕，求解条件方程矩阵，得到GNSS、北斗类测轨数据计算的初始位置、速度向量改进值；

若没有遍历完毕，则取下一帧观测数据，重复步骤2.b.a-2.b.d，至每帧测轨数据遍历完毕。

步骤3，对地基雷达测轨数据进行的初始位置、速度改进值、天基光学测角数据计算的改进值、GNSS及北斗类数据计算的改进值进行加权合成得到总的初始位置、速度改进值。

并进行收敛判断，具体按照下述方法进行：

给定收敛门限ξ，相邻两次得到的初始位置和速度向量改进值分别为

和

若其满足

则认为收敛(此时直接将新的改进值输出)，若不满足，则不收敛(未收敛时重复步骤2进行迭代改进)。

Claims (3)

1.一种基于精确轨道模型的航天器初轨确定方法，其特征在于，具体按照下述步骤进行：

步骤1，对航天器星箭分离后的轨道进行初始估计；

步骤2，融合不同的类型的测轨数据，在有摄单位矢量法和精确轨道模型基础上，基于迭代加权最小二乘法计算初始位置和速度的改进值，最小二乘条件方程根据直接合成法或联邦合成法建立；

步骤3，若初始位置速度前后两次迭代改进值满足收敛条件，则将初始位置、速度的迭代改进值输出；若不满足收敛条件，则重复步骤2至收敛为止；

所述步骤2中根据直接合成法建立最小二乘条件方程具体按照下述步骤进行：

步骤2.1，对每一帧观测数据，以初值t₀、r₀、

进行动力学模型积分到当前数据帧时刻，得到其对应的位置速度估计值，初始位置及速度向量分别为r₀、

步骤2.2，建立条件方程矩阵；

建立最小二乘条件方程：

步骤2.3，分别计算地基雷达测轨数据、天基光学测角类测轨数据和GNSS、北斗类测轨数据计算得到条件方程各项的系数；

步骤2.4，将每个条件方程进行法化计算并将法化后的系数矩阵累加到系统条件方程中，并将累加了系数矩阵的条件方程组建条件方程组；

步骤2.5，判断测轨数据是否遍历完毕，若遍历完毕，得到条件方程矩阵；

若没有遍历完毕，则取下一帧观测数据，并重复步骤2.1-2.4，至每帧测轨数据遍历完毕；

所述步骤2中根据联邦合成法建立最小二乘条件方程具体按照下述步骤进行：

步骤2.a，对每一帧观测数据，以初值t₀、r₀、

进行动力学模型积分到当前数据帧时刻，得到其对应的位置速度估计值

步骤2.b，根据地基雷达测轨数据列写地基雷达条件方程组，并对地基雷达方程组求解得到地基雷达类初始位置和地基雷达类初始速度向量，并对地基雷达条件方程组进行方差统计；

根据天基光学测角类测轨数据列写天基光学测角条件方程组，并对天基光学测角方程组求解得到天基光学测角类初始位置和天基光学测角类初始速度向量，并对天基光学测角条件方程组进行方差统计；

根据GNSS、北斗类测轨数据列写地基雷达条件方程组，并对GNSS、北斗类方程组求解得到GNSS、北斗类初始位置和GNSS、北斗类初始速度向量，并对GNSS、北斗类方程组进行方差统计；

步骤2.c，根据地基雷达条件方程组的方差对地基雷达类初始位置和地基雷达类初始速度向量进行加权；

根据天基光学测角条件方程组的方差对天基光学测角类初始位置和天基光学测角类初始速度向量进行加权；

根据GNSS、北斗类方程组的方差对GNSS、北斗类初始位置和GNSS、北斗类初始速度向量进行加权；

将加权后的地基雷达类初始位置、天基光学测角类初始位置和GNSS、北斗类初始位置进行合成，得到新的初始位置

将加权后的地基雷达类初始速度、天基光学测角类初始速度和GNSS、北斗类初始速度向量进行合成，得到新的初始速度向量。

2.根据权利要求1所述的一种基于精确轨道模型的航天器初轨确定方法，其特征在于，所述步骤2中在最小二乘迭代计算过程中，以精确轨道模型来计算位置速度估计值。

3.根据权利要求1所述的一种基于精确轨道模型的航天器初轨确定方法，其特征在于，所述步骤3中具体按照下述方法判断新的初始位置和新的初始速度向量是否收敛：

给定收敛门限ξ，相邻两次得到的新的初始位置和新的初始速度向量分别为

和

若其满足

x、y、z代表直角坐标的三个方向的分量，

则得到的收敛，若不满足，则不收敛。

Images