CN103940431A - 基于gnss精密定轨的圆轨道切向小推力在轨标定方法 - Google Patents

Abstract

本发明提供一种基于GNSS精密定轨的圆轨道切向小推力在轨标定方法，标定后的切向推力F用于航天器轨道控制。利用GNSS测量得到的航天器位置信息，采用Unscented卡尔曼滤波方法，得到J2000坐标系下航天器位置和速度信息的估计值。根据位置和速度信息的估计值，计算出航天器的瞬时轨道半长轴。针对每一个测量时刻，将该时刻前一个轨道交点周期内的轨道瞬时半长轴求平均值，得到该时刻的平均轨道半长轴。将圆轨道切向推力作用前后的平均轨道半长轴作差，得到轨道半长轴变化量Δa，根据Δa计算得到圆轨道切向推力标定值。本发明计算过程完全利用GNSS获得的实时轨道数据，无需地面测控站数据支持，方法标定结果准确、可靠，算法简便，容易实现。

Description

基于GNSS精密定轨的圆轨道切向小推力在轨标定方法

技术领域

本发明涉及一种基于GNSS（Global Navigation Satellite System，全球导航卫星***）精密定轨的圆轨道切向小推力在轨标定方法，属于航天器在轨标定技术领域。

背景技术

空间推进***通过推力作用实现对航天器的轨道与姿态控制。电推进作为一种先进的空间推进技术，除了具有高比冲、高效率、长寿命等突出特点之外，还具备的一个显著特点是能够产生小推力以实现轨道、姿态的高精度控制。采用电推进技术的推力发动机已经成功应用于空间***中，在地球同步轨道卫星位置保持、超静平台无拖曳控制、远深空探测等领域具有非常广泛的应用前景。

为提高电推进***在航天器轨道控制应用中的控制精度，需要对电推进***的推力大小、特别是小推力进行在轨标定。目前，电推进***的推力标定主要有两种方法：一种方法是依靠地面测试，即通过地面试验获得推力大小和电推进***工作时的电压、电流等信号的拟合曲线关系。电推进***在轨工作时，根据工作电压、电流的遥测信号，与地面试验结果比较，间接获得电推进***在轨产生的推力数值。由于地面环境和空间环境存在较大差异性，这种依靠地面试验数据对电推进***在轨工作的实际推力、特别是小推力进行标定的方法无法保证标定结果的可靠性和精确性。另一种方法是通过地面站对航天器轨道进行测控，根据测控信息判断电推进推力作用前后航天器轨道参数的变化，由此计算推力大小。受地面站位置及数量的制约，地面测控无法实现高精度、全实时的轨道参数测量，因此，地面测控精度很难实现对小推力的在轨标定。上述两种方法都无法保证标定结果的精确性，具有较大的局限性

发明内容

本发明的目的是为克服现有技术的不足，提供了一种基于GNSS测量数据的圆轨道切向小推力在轨标定方法，该方法克服了摄动因素和测量误差对微小推力标定精度的影响，提高了标定结果的可靠性和精确性。

本发明的技术解决方案是：

一种基于GNSS精密定轨的圆轨道切向小推力在轨标定方法，标定后的切向推力F用于航天器轨道控制，具体步骤如下：

步骤一、利用GNSS测量得到的测量时刻t_k（k＝0,1,...n）的WGS-84坐标系下航天器的位置信息，通过坐标转换得到tk时刻J2000坐标系下航天器的位置信息；

步骤二、以J2000坐标系下航天器的位置信息作为测量量，获得航天器的位置和速度信息的估计值；

步骤三、利用滤波处理后的航天器位置和速度信息估计值计算航天器tk时刻（k＝0,1,...n）的瞬时轨道半长轴

步骤四、针对每一个测量时刻t_k（k＝0,1,...n），将该时刻前一个轨道交点周期内的轨道瞬时半长轴求平均值，作为t_k时刻的平均轨道半长轴

步骤五、将圆轨道切向推力结束、并经过一个轨道交点周期之后的时刻的平均轨道半长轴和圆轨道切向推力开始前一个时刻的平均轨道半长轴作差，得到圆轨道切向推力引起的轨道半长轴的变化量Δa；

步骤六、根据轨道半长轴的变化量Δa，计算圆轨道切向推力的标定量其中M为推力作用期间的航天器质量，μ为引力常数，a为推力作用期间轨道瞬时半长轴的平均值，Δt为推力作用时间。

所述步骤二获得航天器的位置和速度信息的估计值是采用Unscented卡尔曼滤波方法实现，其中将t0时刻的滤波估计值及其对应的估计方差矩阵P_0|0作为Unscented卡尔曼滤波方法的初始化状态向量。

本发明与现有技术相比的优点在于：

本发明针对圆轨道切向小推力的在轨标定问题，通过对GNSS测量数据进行滤波处理，解算得到连续时刻的航天器平均轨道半长轴，由平均轨道半长轴的变化量计算得到推力大小。本发明的计算过程完全利用GNSS获得的实时轨道数据，无需地面测控站数据支持，方法标定结果准确、可靠，算法简便，容易实现。

附图说明

图1为本发明的流程框图；

图2为Unscented卡尔曼滤波算法流程图；

图3为切向小推力作用前后轨道平均半场轴计算结果图。

具体实施方式

如图1所示，为本发明的方法流程框图，主要步骤如下：

本发明一种基于GNSS测量数据的圆轨道切向小推力在轨标定方法，标定后的切向推力F用于航天器轨道控制，具体步骤：

步骤一、以GNSS测量的起点时间作为初始时刻t₀，以切向小推力作用结束、且又经过若干（大于等于1个）轨道交点周期后的时刻作为终止时刻t_n，将通过GNSS测量得到的t_k时刻（k＝0,1,...n）的航天器位置信息，通过坐标变换，由WGS-84坐标系转换到J2000坐标系，得到t_k时刻J2000坐标系下航天器的位置坐标分量测量量[x_ky_kz_k]；

本步骤中较佳的选取GPS测量得到的信息。

步骤二、以t_k时刻J2000坐标系下航天器的位置信息[x_ky_kz_k]作为测量量，采用Unscented卡尔曼滤波方法，得到J2000坐标系下位置和速度信息的估计值 ${\hat{X}}_{k|k}=\left[\begin{array}{cccccc}{\hat{x}}_k& {\hat{y}}_k& {\hat{z}}_k& {\hat{v}}_{x_k}& {\hat{v}}_{y_k}& {\hat{v}}_{z_k}\end{array}\right],(k=\mathrm{0,1},...n),$ 如图2所示，为Unscented卡尔曼滤波计算流程，步骤如下：

步骤201、定义状态向量：定义t_k时刻（k＝0,1,...n）J2000坐标系下航天器轨道状态向量的滤波估计状态向量及其估计方差矩阵P_k|k；由t_k时刻J2000坐标系下的航天器位置坐标分量和速度坐标分量的滤波估计值 $\left[\begin{array}{ccc}{\hat{x}}_k& {\hat{y}}_k& {\hat{z}}_k\end{array}\right]$ 和 $\left[\begin{array}{ccc}{\hat{v}}_{x_k}& {\hat{v}}_{y_k}& {\hat{v}}_{z_k}\end{array}\right]$ 组成，表示为 ${\hat{X}}_{k|k}=\left[\begin{array}{cccccc}{\hat{x}}_k& {\hat{y}}_k& {\hat{z}}_k& {\hat{v}}_{x_k}& {\hat{v}}_{y_k}& {\hat{v}}_{z_k}\end{array}\right]$ P_k|k为6×6的矩阵。

步骤202、初始化状态向量：确定初始时刻（k＝0），即t₀时刻的滤波估计值及其对应的估计方差矩阵P_0|0。由t₀时刻J2000坐标系下航天器的位置坐标分量测量值[x₀y₀z₀]以及估算出的速度坐标分量 $\left[\begin{array}{ccc}{v}_{x_0}& v_{y_0}& v_{z_0}\end{array}\right]$ 确定，即 ${\hat{X}}_{0|0}=\left[\begin{array}{cccccc}{x}_0& y_0& z_0& x_{x_0}& v_{y_0}& v_{z_0}\end{array}\right],$ P_0|0根据GPS定位误差及估算的初始时刻速度误差确定；

步骤203、计算采样点状态向量：对于t_k-1时刻（k＝1,2,...,n），根据计算13个采样点χ_i,k-1|k-1（i＝0,1，2,...12）；

${\mathrm{\χ}}_{0,k-1|k-1}={\hat{X}}_{k-1|k-1},$

${\mathrm{\χ}}_{i,k-1|k-1}={\hat{X}}_{k-1|k-1}+{\left(\sqrt{(6+\mathrm{\λ})P_{k-1|k-1}}\right)}_i(i=1,...,6),$

${\mathrm{\χ}}_{i,k-1|k-1}={\hat{X}}_{k-1|k-1}+{\left(\sqrt{(6+\mathrm{\λ})P_{k-1|k-1}}\right)}_{i-6}(i=7,...,12),$

其中λ为采样常数，表示矩阵平方根的第i列；

步骤204、根据采样点χ_i,k-1|k-1（i＝0,1，2,...12）计算t_k时刻的状态向量一步预测值分别以χ_i,k-1|k-1（i＝0,1，2,...12）为初值，采用4阶Runge-Kutta数值解法，求解航天器轨道动力学微分方程，得到χ_i,k|k-1（i＝0,1，2,...12），其中轨道动力学微分方程考虑地球非球形引力摄动J₂、J₃、J₄项。由χ_i,k|k-1（i＝0,1,...12）得到 ${\hat{X}}_{k-k-1}=\underset{i=0}{\overset{12}{\mathrm{\Σ}}}{w}_i^m{\mathrm{\χ}}_{i,k|k-1}$ 其中W_i ^m为向量权重系数， $W_0^m=\mathrm{\λ}/(6+\mathrm{\λ}),W_i^m=1/\left[2(6+\mathrm{\λ})\right]$ $(i=1,...,12);$

步骤205、根据和计算状态估计方差和测量量的一步预测值Pk|k-1和计算公式为 $P_{k|k-1}=\underset{i=0}{\overset{12}{\mathrm{\Σ}}}{W}_i^c({\mathrm{\χ}}_{i,k|k-1}-{\hat{X}}_{k|k-1}){({\mathrm{\χ}}_{i,k|k-1}-{\hat{X}}_{k|k-1})}^T+Q,z_{i,k|k-1}={\mathrm{H\χ}}_{i,k|k-1}$ 其中Q为***噪声协方差常值矩阵，H为观测方程常值矩阵，Wic为矩阵权重系数， $W_0^c=\mathrm{\λ}/(6+\mathrm{\λ})+3,W_i^c=1/\left[2(6+\mathrm{\λ})\right](i=1,...,12);$

步骤206、计算测量量方差矩阵的一步预测P_z,k|k-1和相关方差矩阵P_xz,k|k-1： $P_{z,k|k-1}=\underset{i=0}{\overset{12}{\mathrm{\Σ}}}{W}_i^c{(z_{i,k\backslash k-1}-{\hat{Z}}_{k|k-1})}^T+R,P_{\mathrm{xz},k|k-1}=\underset{i=0}{\overset{12}{\mathrm{\Σ}}}{W}_i^c({\mathrm{\χ}}_{i,k|k-1}-{\hat{X}}_{k|k-1}){(z_{i,k|k-1}-{\hat{Z}}_{k|k-1})}^T,$ 其中R为测量噪声协方差常值矩阵；

步骤207、计算增益矩阵K_k，

步骤208、计算t_k时刻的状态向量估计值及其估计方差矩阵P_k|k，计算公式分别为 ${\hat{X}}_{k|k}={\hat{X}}_{k|k-1}+K_k(Z_k-{\hat{Z}}_{k|k-1}),P_{k|k}=P_{k|k-1}-K_k{P}_{z,k|k-1}{K}_k^T,$ 其中Z_k为t_k时刻J2000坐标系下航天器的位置坐标分量测量值，即Z_k＝[x_ky_kz_k]。

在第一次执行完步骤208后得到t_k时刻（k＝1）状态向量估计值及其估计方差矩阵P_1|1，基于此时得到的和P_1|1，返回步骤203，计算t_k时刻（k＝2）状态向量估计值及其估计方差矩阵P_2|2，然后反复执行，直至t_k时刻（k＝n）的状态向量估计值及其估计方差矩阵P_n|n计算完为止。

步骤三、根据Unscented卡尔曼滤波得到的t_k时刻（k＝0,1,...n）航天器位置、速度估计状态向量 ${\hat{X}}_{k|k}=\left[\begin{array}{cccccc}{\hat{x}}_k& {\hat{y}}_k& {\hat{z}}_k& {\hat{v}}_{x_k}& {\hat{v}}_{y_k}& {\hat{v}}_{z_k}\end{array}\right],$ 计算对应时刻的的瞬时轨道半长轴计算公式为 ${\hat{a}}_k=-\mathrm{\μ}/({\hat{v}}_k^2-2\mathrm{\μ}/{\hat{r}}_k),$ 其中μ为引力常数， ${\hat{r}}_k=\sqrt{{\hat{x}}_k^2+{\hat{y}}_k^2+{\hat{z}}_k^2},$ ${\hat{v}}_k=\sqrt{{\hat{v}}_{x_k}^2+{\hat{v}}_{y_k}^2+{\hat{v}}_{z_k}^2},k=\mathrm{0,1},...,n;$

步骤四、计算t_k时刻（k＝0,1;...n）的平均轨道半长轴计算公式为其中m为一个轨道交点周期内的的采样个数，表示从t_k时刻向前递推j个时刻的瞬时轨道半长轴；采用上述公式计算得到平均轨道半长轴，可以有效消除非球形引力摄动等保守力引起的轨道半长轴变化量，也可以减少随机测量噪声对半长轴测量精度的影响，从而最大程度反映切向推力对半长轴的作用结果。

步骤五、计算切向推力引起的轨道半长轴的变化量Δa，计算公式为其中表示切向推力结束、并经过一个轨道交点周期之后的时刻的平均轨道半长轴，表示切向推力开始前一个时刻的平均轨道半长轴；由于非球形引力摄动等保守力不会改变平均轨道半长轴，因此，采用上述公式计算得到的平均轨道半长轴变化量反映了切向推力引起的平均半长轴的改变程度。

步骤六、计算切向推力F，计算公式为其中M为推力作用期间航天器的质量，Δt为推力作用时间，a为推力作用期间轨道瞬时半长轴的平均值，d为推力作用时间Δt期间的采样个数，表示推力作用时间Δt期间各采样时刻的瞬时轨道半长轴。

通常情况轨道控制需要基于标定的切向推力来实现，因此在需要对航天器轨道进行控制时，即可基于该切向推力F实现。

实施例

以某航天器的小推力在轨标定结果为例，应用本发明方法解算得到的航天器轨道平均半长轴结果如图3所示。应用本发明方法计算得到的小推力结果和推力器的推力地面测试标称值以及传统电流电压测量法结果比较：

本发明说明书中未作详细描述的内容属本领域技术人员的公知技术。

综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (2)

1.基于GNSS精密定轨的圆轨道切向小推力在轨标定方法，标定后的切向推力F用于航天器轨道控制，其特征在于步骤如下：

步骤一、利用GNSS测量得到的测量时刻t_k(k＝0，1，…n)的WGS-84坐标系下航天器的位置信息，通过坐标转换得到t_k时刻J2000坐标系下航天器的位置信息：

步骤二、以J2000坐标系下航天器的位置信息作为测量量，获得航天器的位置和速度信息的估计值；

步骤三、利用滤波处理后的航天器位置和速度信息估计值计算航天器t_k时刻的瞬时轨道半长轴

步骤四、针对每一个测量时刻，将该时刻前一个轨道交点周期内的轨道瞬时半长轴求平均值，作为t_k时刻的平均轨道半长轴

步骤五、将圆轨道切向推力结束、并经过一个轨道交点周期之后的时刻的平均轨道半长轴和圆轨道切向推力开始前一个时刻的平均轨道半长轴作差，得到圆轨道切向推力引起的轨道半长轴的变化量Δa；

步骤六、根据轨道半长轴的变化量Δa，计算圆轨道切向推力的标定量其中M为推力作用期间的航天器质量，μ为引力常数，a为推力作用期间轨道瞬时半长轴的平均值，Δt为推力作用时间。

2.根据权利要求1所述基于GNSS精密定轨的圆轨道切向小推力在轨标定方法，其特征在于，所述步骤二获得航天器的位置和速度信息的估计值是采用unscented卡尔曼滤波方法实现，其中将t₀时刻的滤波估计值及其对应的估计方差矩阵P_0|0作为Unscented卡尔曼滤波方法的初始化状态向量。