[发明专利]电推进转移轨道航天器星载GNSS在轨实时定轨方法

权利要求书

1.一种电推进转移轨道航天器星载GNSS在轨实时定轨方法，其特征在于，包括以下具体步骤：

步骤S1：实时获取GNSS观测数据、卫星姿态以及电推力广播数据；

步骤S2：当判断先验轨道信息有效标记有效时，计算并统计验前伪距观测值残差的均值和标准差，当判断先验轨道信息有效标记无效时，跳转到步骤S4；

步骤S3：根据步骤S2中得到的伪距观测值残差的均值和标准差信息判断先验轨道信息是否准确和接收机各通道是否含有伪距粗差，当先验轨道信息不准确时，调整补偿加速度的过程噪声；

步骤S4：完成顾及伪距粗差的伪距单点定轨,

星载GPS/BDS伪距观测方程表示为：

式(1)中，分别为经过电离层改正后的星载GPS/BDS的伪距观测值；δt^G和δt^C分别为GPS/BDS接收机钟差；c为真空中的光速；和分别为第i颗GPS卫星和第j颗BDS卫星的钟差，分别由GPS/BDS广播星历计算得到；和分别为单频GPS/BDS电离层延迟改正值，和分别为GPS/BDS伪距观测值噪声；为相应伪距观测值对应的站星几何距离，用下式表示：

式(2)中，(x,y,z)分别为GNSS接收机在地心地固系内三维位置坐标向量，分别为第i颗GPS卫星和第j颗BDS卫星在地心地固系内三维位置坐标向量，分别由GPS/BDS广播星历计算得到；

假定GNSS接收机的位置坐标(x,y,z)和接收机钟差(cδt_G,cδt_C)组成的向量为X＝(x,y,z,cδt^G,cδt^C)^T，其初值和改正值分别为和ΔX＝(dx,dy,dz,dcδt^G,dcδt^C)^T，则GPS/BDS伪距观测值关于X的偏导数为：

如果t时刻星载GPS/BDS接收机共观测到m颗GPS和n颗BDS卫星，假设步骤S3中第k颗GPS卫星和第l颗BDS卫星被标识含有伪距粗差，则这两颗卫星不参与计算观测值矩阵A和残差向量b，其计算公式如下：

式(4)中，和分别为第i颗GPS卫星和第j颗BDS卫星伪距观测值关于X的偏导数，初始计算时，是将X＝X₀带入式(3)；分别为GPS卫星实测的伪距观测值和伪距观测值的计算值；分别为BDS卫星实测的伪距观测值和伪距观测值的计算值；

步骤S5：完成顾及电推力影响的定轨时间更新过程，在地心惯性系中，卫星运动方程用一阶微分方程组来表示：

式(5)中，分别为卫星的位置、速度；为卫星受到的各种摄动加速度之和，可以表示为：

式(6)中，为近地航天器所受的惯性系总加速度；为由保守力引起的加速度；为由非保守力引起的加速度；为电推力引起的加速度；为人为引入的经验补偿加速度，用于补偿无法模型化或错误模型的微小摄动力的影响，采用一阶高斯-马尔可夫随机模型对径向、切向、法向3个方向进行动力学模型补偿；T是RTN坐标系到地心惯性系的转换矩阵；

由电推力引起的加速度计算如下式所示：

式(7)中，为卫星本体系到地心惯性系的转换矩阵，由步骤S1实时获取的姿态广播数据计算；为电推力参考坐标系到卫星本体系的转换矩阵，由离子电推进器在卫星本体系具体安装方向计算得到；为离子电推进器在三个方向上的电推力矢量，单位为N，由步骤S1实时获取的电推力广播数据得到；在精密定轨中，除了要考虑位置参数和速度参数之外，还要考虑动力学参数且有因此将式(5)进行扩展，并令则有

当设定一初始状态后，通过运动方程积分得到卫星的参考状态向量该初值从步骤S4中得到星载GNSS接收机位置作为粗略轨道获得，由于上面的状态方程是非线性的，所以对在上展开取一次项，可得：

令则有：其解写为：

则有：

式(7)称为变分方程，其中I为单位矩阵，Φ(t,t₀)被称为状态转移矩阵；

动力学定轨中，通过数值积分方法，按式(8)和式(11)，同时求解卫星的运动方程和变分方程；

伪距定轨卡尔曼滤波的状态方程为：

X_k＝Φ_k,k-1X_k-1+W_k (12)

式(5)中,X_k＝(r,v,cδt^G,cδt^C,C_d,C_r,ω)^T为伪距定轨卡尔曼滤波状态量，Φ_k,k-1为离散形式的状态转移矩阵，W_k为系统噪声矩阵；

式(13)中，Φ_rr,Φ_rv,Φ_vr,Φ_vv分别位置和速度分量的状态转移矩阵；位置和速度分别关于大气阻力和太阳光压的状态转移矩阵；Φ_rw,Φ_vw分别位置和速度关于RTN方向的补偿加速度的状态转移矩阵，Φ_ww为补偿加速度关于其自身的状态转移矩阵；

伪距定轨卡尔曼滤波的观测方程为：

Z_k＝H_kX_k-1+V_k (14)

式(14)中，Z_k为观测矢量，H_k为观测矩阵，V_k为观测噪声矢量，且系统噪声W_k和观测噪声V_k，为零均值白噪声序列，Q_k和R_k分别为系统动态噪声协方差阵和观测噪声协方差阵；

伪距定轨时间更新具体过程是，通过数值积分方法，同时求解卫星的运动方程和变分方程，得到Φ_k,k-1，再按式(12)完成对滤波状态量的时间更新，再完成对状态误差协方差矩阵的时间更新；

式(15)中，为k-1时刻的状态误差协方差矩阵，为k时刻的状态误差协方差矩阵；

步骤S6：完成顾及伪距粗差的定轨测量更新过程，式(3)中伪距观测值关于伪距定轨滤波器状态量的偏导数H可以表示为：

式(16)中，U^T为J2000惯性系到地固系的转换矩阵，为航天器相对于GPS卫星i和BDS卫星j的视线向量；

顾及伪距粗差的定轨测量更新过程为：按式(17)，依次对于步骤S3中没有被标识含有伪距粗差的GPS/BDS伪距观测值进行处理，更新绝对定轨卡尔曼滤波的状态量及状态协方差阵；

式(16)中，为当前k时刻第i颗GPS/BDS对应的卡尔曼增益矩阵；为式(16)中的观测矩阵；为观测噪声协方差阵，为第i颗GPS/BDS测量更新后的滤波状态量；为步骤S5中当前k时刻时间更新后的滤波状态量，第i颗GPS/BDS实测的伪距观测值；为式(15)中当前k时刻时间更新后的状态误差协方差矩阵；为当前k时刻测量更新后的状态误差协方差矩阵；

步骤S7：完成先验轨道的预报及先验轨道精度计算，标识先验轨道信息有效，在步骤S6中定轨测量更新最后得到的滤波状态量的基础上，通过数值积分方法，按式(10)求解卫星运动方法，得到下一个时刻的先验轨道PreOrb，先验轨道精度Accu由在步骤S6中状态误差协方差矩阵计算；当步骤S6中参与测量更新的卫星数大于0时，标识先验轨道信息有效，否则，标识先验轨道信息无效；

步骤S8：输出定轨滤波结果和先验轨道预报结果。