实时精密定轨中轨道的星载快速多步积分方法 技术领域 本发明涉及导航卫星应用技术领域， 尤其涉及一种实时精密定轨中轨道的星载快 速多步积分方法。
背景技术 在精密定轨中为了获得高精度的卫星轨道和状态转移矩阵， 需要利用数值积分算 法求解卫星运动方程。当前在定轨中主要使用的积分算法为 Runge-Kutta 单步法和 Adams 多步法。Runge-Kutta 法间接引用泰勒展开式， 用区间 [ti ti+1] 上若干个右函数 f 的线性 组合来代替 f 的导数， 相应的组合系数由泰勒展开式确定。由于 Runge-Kutta 法在进行数 值积分时需要多次求解当前历元和积分历元间不同时刻的右函数值， 这样对于卫星精密定 轨中复杂右函数的数值积分将是较为耗时的工作 ( 主要时间消耗在重力场非球形摄动、 高 阶海潮、 大气潮汐计算上 )， 另外该方法还有截断误差难以估计的缺点， 因此在精度要求较 高的卫星运动方程数值积分中， 大多采用 Fehlberg 提出的 Runge-Kutta-Fehlberg 方法， 该 方法实质为嵌套的 Runge-Kutta 方法， 其同时给出 n 和 n+1 阶两组 Runge-Kutta 方程， 用 两组公式计算得到的积分历元卫星运动状态的差值来估计截断误差， 根据截断误差的大小 来控制步长。由于 n+1 阶与 n 阶 Runge-Kutta 计算公式相差较少， 只需多计算很少几次右 函数， 却能同时获得局部截断误差， 并且其稳定度较好， 能够保持积分所需精度。但是单步 法仅仅利用当前历元信息而无法使用之前历元信息， 因此其各步之间具有独立性。多步法 则可以有效利用已有历元信息， 综合求解下一积分历元卫星运动状态。其中使用最为广泛 的是 J.C.Adams 开发的 Adams 多步法。Adams 算法充分利用已有历元信息， 在数值积分下 一历元卫星运动状态时只需计算一次右函数， 从而大大减小计算量。为了控制积分精度， Bashforth 和 Moulton 等人对 Amdas 公式分别进行了改进， 得到了显示的 Adams-Bashforth 公式和隐式的 Adams-Moulton 公式， 两者的差别仅在于计算积分历元的右函数时是否使用 到该历元卫星的运动状态。 在卫星运动方程的数值积分计算时， 一般同时采用这两个公式， 先由显式公式计算出积分历元卫星运动状态近似值， 再由隐式公式校正该近似值， 在数学 上该过程称为 PECE 算法， 即： 预报 - 校正算法。相较于单步法， 多步法不能自行起算， 因此 其需要采用单步法推出足够的步点， 然后才能计算。 但是在同样阶次， 多步法相较于单步法 计算精度高， 运算速度快。
如何快速、 精确地求得卫星运动状态初值以及相应的状态转移矩阵用于参数估计 是实时精密定轨中关键问题。 Runge-Kutta 单步法由于需要计算多次右函数， 导致运算速度 过慢。 虽然多步法仅需多计算一次右函数， 但是由于在实施定轨时， 卫星运动状态会随观测 这 值而更新， 因此其需要使用观测更新的卫星状态重新计算右函数和相应状态转移矩阵， 样便使得右函数计算次数不再是一次而是多次， 同样会产生和单步法相同的问题， 即计算 速度过慢而不适合于卫星实时精密定轨。
发明内容 针对上述存在的技术问题， 本发明的目的是提供一种实时精密定轨中轨道的星载 快速多步积分方法， 以解决卫星实时定轨中轨道积分速度过慢的技术问题， 从而实现快速 轨道积分。
为达到上述目的， 本发明采用如下的技术方案 ：
①选择积分窗口 ；
②利用窗口末端历元处的观测值， 对整个窗口内所有历元卫星运动状态进行平 滑， 从而得到窗口内所有历元卫星状态平滑值 ；
③利用窗口内所有历元卫星运动状态和已有状态转移矩阵， 采用 Adams 预估 - 校 正算法预报窗口末端历元至下一历元的状态转移矩阵和该历元卫星的运动状态 ；
④更新窗口历元， 使窗口起始历元后移一个历元 ；
⑤判断当前历元是否存在观测值， 如果有观测值存在则执行步骤②， 否则执行步 骤⑥ ；
⑥结束。
所述步骤③采用 Adams 预估 - 校正方法 ：
Adams-Bashforth 公式 ：
Adams-Moulton 公式 ：
xn 表示在第 n 个积分历元卫星的运动状态， h 为积分步长， f 为积分函数， 而β和 γ 为多步法积分系数 ；
其中右函数 f 只需重新在校正时计算一次。
本发明具有以下优点和积极效果 ：
1) 历元更新时采用多步法， 相较采用单步法其数值精度高、 稳定性更好， 积分右函 数计算较少 ；
2) 利用已有状态转移矩阵信息， 避免了在进行状态更新时需要更新状态转移矩 阵， 从而减少积分中右函数计算次数， 有效提高计算速度 ；
3) 采用移动窗口方法， 可以根据不同计算精度要求变换窗口长度， 保证积分精度。
附图说明
图 1 是本发明提供的实时精密定轨中轨道的星载快速多步积分方法的流程图。 图 2 是本发明提供的实时精密定轨中轨道的星载快速多步积分方法的数据流图。具体实施方式
假设卫星运动状态在进行状态更新时其值变化不大， 那么由更新值计算得到的状 态转移矩阵变化也不大， 这样可以直接利用已有卫星状态转移矩阵， 而无需利用更新后的 卫星状态重新计算状态转移矩阵， 这样在采用多步法预报卫星在新的历元运动状态时， 仅 仅需要计算当前历元到下一历元的右函数， 这样右函数只需计算一次从而大大减小计算时 耗， 提高运算速度。
本发明提供的实时精密定轨中轨道星载快速多步积分方法， 其流程参见图 1， 具体 包括以下步骤 ：
步骤 S101 ： 选择积分窗口， 一般为多步法所需求的点数， 如果积分窗口选择为 1， 则退化为单步法 ；
假定所选择窗口的宽度为 N， 即窗口内有 N 个历元， 记为 ：
ti ti+1… ti+N-2 ti+N-1
步骤 S102 ： 利用窗口末端历元处的观测值， 对整个窗口内所有历元卫星运动状态 进行平滑， 从而得到窗口内所有历元卫星状态平滑值 ；
这一步采用均方根信息平滑算法完成， 其状态方程为 ：
xk ＝ Ф(tk， tk-1)xk-1+Γ(tk， tk-1)uk-1 (1)
式中 xk 和 xk-1 分别是卫星在 tk 和 tk-1 的运动状态， Ф(tk， tk-1) 是从 tk-1 时刻到 tk 时刻状态转移矩阵， Γ(tk， tk-1) 从 tk-1 时刻到 tk 时刻噪声转移矩阵， uk-1 为 tk-1 时刻状态噪 声。式中 xk-1 具有先验值 方程 ：
和先验方差将先验方差进行 Cholesky 分解， 构造虚拟观测式中ηk-1 为 卫 星 状 态 误 差， 其均值和先验方差分别为 而 式中 uk-1 先验值与真值关系可用下式来描述 ：
式中 αk-1 为 tk-1 时刻的状态噪声误差， 其均值和方差为 E[αk-1] ＝ 0， E[αk-1αk-1T] ＝ Q， 从而构建状态噪声 uk-1 的虚拟观测方程 ：
式中而而滤波的观测方程为 ： yk-1 ＝ Hk-1xk-1+εk-1 (5) 式中 yk-1 为观测向量， Hk-1 为设计矩阵， εk-1 为观测值误差， 其均值和方差分别为 根据最小方差准则， 也即使5
和 εk-1 平方和最小， 从而可以构建均方根信101853027 A CN 101853028说明书4/7 页息滤波算法观测更新性能函数
式中 ‖‖ 表示任意向量的 2 范数。 将 (6) 式写成矩阵形式可得 ：对 (7) 式进行正交变化可以得到 ：式中和 ek-1 分别是和 yk-1 进行正交变换的结果。同样可以根据最小方差准则， 可以构建均方根信息滤波算法状态更新性能函数
将 (9) 式写成矩阵形式 ：式中将 (10) 式做正交变换可以得到 ：其中和都为式 (10) 中对应值正交变换的结果。而求解卫星运动状态平滑解的递推公式为 ：式中为由 tj 时刻平滑 tj-1 时刻正交变换矩阵，与式 (11) 中含义相同， Γ(tj， tj-1) 与式 (1) 中含义相同， 是对
6进行正交变化的结果， 其起始值为上式右边 ()* 都表示为左边中对应值的正交变换值。 由此可得平滑解为 ： (13)其相应的方差和协方差矩阵为 ：101853027 A CN 101853028
说明书5/7 页步骤 S103 ： 利用窗口内所有历元卫星运动状态和已有状态转移矩阵， 采用 Adams 预估 - 校正算法预报窗口末端历元至下一历元的状态转移矩阵和该历元卫星的运动状态 ；
这一步采用 Adams 预估 - 校正算法， 其相应公式如下所述。
① Adams-Bashforth 公式 ( 隐式计算公式 ) ：
② Adams-Moulton 公式 ( 显示计算公式 ) ：
(15) 和 (16) 两式中， xn 表示在第 n 个积分历元卫星的运动状态， h 为积分步长， f 为积分函数， 而 β 和 γ 为多步法积分系数。其中右函数 f 只需在校正时重新计算一次。
步骤 S104 ： 更新窗口历元， 也即使窗口起始历元后移一个历元 ；
这一步也即将窗口由 ti ti+1… ti+N-2 ti+N-1 更新为 ti+1 ti+2… ti+N-1 ti+N
步骤 S105 ： 判断当前历元是否存在观测值， 如果有观测值存在则执行步骤 S102， 否则执行步骤 S106 ；
步骤 S106 ： 结束。
在上述的实时精密定轨中轨道的星载快速多步积分方法中， 假设卫星状态变化对 状态转移矩阵影响不大是核心内容， 只有这样才能有效地采用该算法。在实际中这种情况 较为容易满足， 也即滤波一旦收敛后卫星先验运动状态已经足够精确， 观测值对其改进已 经不太大， 因此状态转移矩阵变化也不大。另外积分窗口的选择可以根据积分精度变化而 定， 其并不影响积分中计算右函数的次数， 因此其并不影响积分时间消耗。
图 2 表示本发明的数据流， 主要展示了在实时精密定轨中窗口如何更新。其中黄 色表示当前窗口历元， 红色为预估历元， 灰色为过时历元， 该图横轴表示时间， 竖轴表示处 理进程。
表 1-3 对上述方法的假设条件 ( 即当卫星运动状态变化不大时， 相应状态转移矩 阵变化也不大 ) 进行了验证。表 1 中第二列 ( 即初值 1) 给出的是更新后的状态参数， 第三 列 ( 即初值 2) 给出的是更新前的状态参数 ( 前 6 行为卫星三维位置和速度， 后 6 行为力模 型参数， 其中位置单位为 km， 速度单位为 km/s， 力模型参数无量纲 )。表 2 中给出的是由表
1 中初值 1 积分得到的 8 个历元后状态转移矩阵。表 3 中给出的是由表 1 中初值 2 积分得 到的 8 个历元后状态转移矩阵 ( 表 2， 、 3 中第 2、 3、 4 列分别为当前卫星状态和力模型参数 对初始历元卫星位置和速度的偏导数 )。从表 1 中可以看出卫星轨道差值可达 cm 级， 而其 他力模型参数相差则更大， 但是利用两者分别计算得到的 8 个历元后的状态转移矩阵 ( 表 2 和表 3) 则相差不大， 这充分验证本发明假设的有效性。
表1： 卫星在某一时刻更新状态和初始状态
表2： 由卫星更新状态积分得到的状态转移矩阵
表3： 由卫星初始状态积分得到的状态转移矩阵