专利名称：实时精密定轨中轨道的星载快速多步积分方法
技术领域：
本发明涉及导航卫星应用技术领域，尤其涉及一种实时精密定轨中轨道的星载快 速多步积分方法。
背景技术：
在精密定轨中为了获得高精度的卫星轨道和状态转移矩阵，需要利用数值积分算 法求解卫星运动方程。当前在定轨中主要使用的积分算法为Rimge-Kutta单步法和Adams 多步法。Runge-Kutta法间接引用泰勒展开式，用区间[、ti+1]上若干个右函数f的线性 组合来代替f的导数，相应的组合系数由泰勒展开式确定。由于Rimge-Kutta法在进行数 值积分时需要多次求解当前历元和积分历元间不同时刻的右函数值，这样对于卫星精密定 轨中复杂右函数的数值积分将是较为耗时的工作(主要时间消耗在重力场非球形摄动、高 阶海潮、大气潮汐计算上)，另外该方法还有截断误差难以估计的缺点，因此在精度要求较 高的卫星运动方程数值积分中，大多采用Fehlberg提出的Runge-Kutta-Fehlberg方法，该 方法实质为嵌套的Rimge-Kutta方法，其同时给出n和n+1阶两组Runge-Kutta方程，用 两组公式计算得到的积分历元卫星运动状态的差值来估计截断误差，根据截断误差的大小 来控制步长。由于n+1阶与n阶Rimge-Kutta计算公式相差较少，只需多计算很少几次右 函数，却能同时获得局部截断误差，并且其稳定度较好，能够保持积分所需精度。但是单步 法仅仅利用当前历元信息而无法使用之前历元信息，因此其各步之间具有独立性。多步法 则可以有效利用已有历元信息，综合求解下一积分历元卫星运动状态。其中使用最为广泛 的是J. C. Adams开发的Adams多步法。Adams算法充分利用已有历元信息，在数值积分下 一历元卫星运动状态时只需计算一次右函数，从而大大减小计算量。为了控制积分精度， Bashforth和Moulton等人对Amdas公式分别进行了改进，得到了显示的Adams-Bashforth 公式和隐式的Adams-Moulton公式，两者的差别仅在于计算积分历元的右函数时是否使用 到该历元卫星的运动状态。在卫星运动方程的数值积分计算时，一般同时采用这两个公式， 先由显式公式计算出积分历元卫星运动状态近似值，再由隐式公式校正该近似值，在数学 上该过程称为PECE算法，即预报-校正算法。相较于单步法，多步法不能自行起算，因此 其需要采用单步法推出足够的步点，然后才能计算。但是在同样阶次，多步法相较于单步法 计算精度高，运算速度快。如何快速、精确地求得卫星运动状态初值以及相应的状态转移矩阵用于参数估计 是实时精密定轨中关键问题。Rimge-Kutta单步法由于需要计算多次右函数，导致运算速度 过慢。虽然多步法仅需多计算一次右函数，但是由于在实施定轨时，卫星运动状态会随观测 值而更新，因此其需要使用观测更新的卫星状态重新计算右函数和相应状态转移矩阵，这 样便使得右函数计算次数不再是一次而是多次，同样会产生和单步法相同的问题，即计算 速度过慢而不适合于卫星实时精密定轨。
发明内容
针对上述存在的技术问题，本发明的目的是提供一种实时精密定轨中轨道的星载 快速多步积分方法，以解决卫星实时定轨中轨道积分速度过慢的技术问题，从而实现快速 轨道积分。为达到上述目的，本发明采用如下的技术方案①选择积分窗口 ；②利用窗口末端历元处的观测值，对整个窗口内所有历元卫星运动状态进行平 滑，从而得到窗口内所有历元卫星状态平滑值；③利用窗口内所有历元卫星运动状态和已有状态转移矩阵，采用Adams预估-校 正算法预报窗口末端历元至下一历元的状态转移矩阵和该历元卫星的运动状态；
骤⑥；
④更新窗口历元，使窗口起始历元后移一个历元；
⑤判断当前历元是否存在观测值，如果有观测值存在则执行步骤②，否则执行步
)结束■
所述步骤③采用Adams预估-校正方法 Adams-Bashforth 公式 Adams-Moulton 公式 、表示在第n个积分历元卫星的运动状态，h为积分步长，f为积分函数，而日和 Y为多步法积分系数；其中右函数f只需重新在校正时计算一次。本发明具有以下优点和积极效果1)历元更新时采用多步法，相较采用单步法其数值精度高、稳定性更好，积分右函 数计算较少；2)利用已有状态转移矩阵信息，避免了在进行状态更新时需要更新状态转移矩 阵，从而减少积分中右函数计算次数，有效提高计算速度；3)采用移动窗口方法，可以根据不同计算精度要求变换窗口长度，保证积分精度。
图1是本发明提供的实时精密定轨中轨道的星载快速多步积分方法的流程图。图2是本发明提供的实时精密定轨中轨道的星载快速多步积分方法的数据流图。
具体实施例方式假设卫星运动状态在进行状态更新时其值变化不大，那么由更新值计算得到的状 态转移矩阵变化也不大，这样可以直接利用已有卫星状态转移矩阵，而无需利用更新后的 卫星状态重新计算状态转移矩阵，这样在采用多步法预报卫星在新的历元运动状态时，仅 仅需要计算当前历元到下一历元的右函数，这样右函数只需计算一次从而大大减小计算时 耗，提高运算速度。本发明提供的实时精密定轨中轨道星载快速多步积分方法，其流程参见图1，具体 包括以下步骤步骤S101 选择积分窗口，一般为多步法所需求的点数，如果积分窗口选择为1， 则退化为单步法；假定所选择窗口的宽度为N，即窗口内有N个历元，记为 步骤S102 利用窗口末端历元处的观测值，对整个窗口内所有历元卫星运动状态 进行平滑，从而得到窗口内所有历元卫星状态平滑值；这一步采用均方根信息平滑算法完成，其状态方程为 式中xk和x^分别是卫星在tk和t^的运动状态，①(tk，t^)是从tn时刻到tk 时刻状态转移矩阵，r (tk，从k时刻到tk时刻噪声转移矩阵，Uh为、时刻状态噪 声。式中Xh具有先验值和先验方差巧，将先验方差进行Cholesky分解，构造虚拟观测 方程 式中巧，！！㈣为卫星状态误差，其均值和先验方差分别为 = o，￡[nrl = = RtR，而巧=h ’k-、=。⑴式中 uh 先验值
与真值关系可用下式来描述 式中ciH为^时刻的状态噪声误差，其均值和方差为E[ak_J 二。』^。^^] =Q，从而构建状态噪声Uh的虚拟观测方程 而滤波的观测方程为 式中yn为观测向量，L为设计矩阵，e H为观测值误差，其均值和方差分别为 五 根据最小方差准则，也即使％_,、^！^平方和最小，从而可以构建均方根信
5
息滤波算法观测更新性能函_改人-丨
式中1II表示任意向量的2范
将(6)式写成矩阵形式可得
(6)
对(7)式进行正交变化可以得到
(7)
(8)式中民_,、和知分别是、和yH进行正交变换的结果。同样可以根据最小方差准则，可以构建均方根信息滤波算法状态更新性能函数
将(9)式写成矩阵形式 R,.
(9)
"t-i
(10)
其中&、K
而求解卫星运动状态平滑解的递推公式为
、和都为式(10)中对应值正交变换的结果。 T,
(12)式中为由tj时刻平滑时刻正交变换矩阵，RUj、与式(11)中含义相 同，r (tj, tjM)与式(i)中含义相同，<是对<+,r(/y+1々)进行正交变化的结果，其起始值为上式右边(广都表示为左边中对应值的正交变换值。 由此可得平滑解为
(13)
其相应的方差和协方差矩阵为[oo67] 步骤S103 利用窗口内所有历元卫星运动状态和已有状态转移矩阵，采用Adams 预估-校正算法预报窗口末端历元至下一历元的状态转移矩阵和该历元卫星的运动状态；这一步采用Adams预估-校正算法，其相应公式如下所述。①Adams-Bashforth公式(隐式计算公式) ‘ k-l ②Adams-Moulton公式(显示计算公式)
(15)和(16)两式中，、表示在第n个积分历元卫星的运动状态，h为积分步长，f 为积分函数，而0和Y为多步法积分系数。其中右函数f只需在校正时重新计算一次。步骤S104 更新窗口历元，也即使窗口起始历元后移一个历元；这一步也即将窗口由、ti+1".ti+N_2ti+N_i 更新为 ti+1 ti+2…‘h 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 由卫星更新状态积分得到的状态转移矩阵
Kbsl BERN
4.115577330647144E-8 -9.878391316129633E-5 -9.517202236014446E-8 表3:由卫j
L初始状态积分得到的状态转移矩阵
参数 PXSAT
权利要求
一种实时精密定轨中轨道的星载快速多步积分方法，其特征在于，包括①选择积分窗口；②利用窗口末端历元处的观测值，对整个窗口内所有历元卫星运动状态进行平滑，从而得到窗口内所有历元卫星状态平滑值；③利用窗口内所有历元卫星运动状态和已有状态转移矩阵，采用Adams预估-校正算法预报窗口末端历元至下一历元的状态转移矩阵和该历元卫星的运动状态；④更新窗口历元，使窗口起始历元后移一个历元；⑤判断当前历元是否存在观测值，如果有观测值存在则执行步骤②，否则执行步骤⑥；⑥结束。
2.根据权利要求1中所述的实时精密定轨中轨道的星载快速多步积分方法，其特征在于所述步骤③采用Adams预估-校正方法 Adams-Bashforth 公式 Adams-Moulton 公式 xn表示在第η个积分历元卫星的运动状态，h为积分步长，f为积分函数，而β和γ为 多步法积分系数；其中右函数f只需重新在校正时计算一次。
全文摘要
本发明涉及导航卫星应用技术领域，尤其涉及一种实时精密定轨中轨道的星载快速多步积分方法。本发明选择积分窗口，利用窗口末端历元处的观测值对整个窗口内所有历元卫星运动状态进行平滑，利用窗口内所有历元卫星更新后运动状态和已有状态转移矩阵积分计算新的历元卫星运动状态，然后更新窗口历元，也即使窗口起始历元后移一个历元，利用更新窗口末端历元的观测信息，平滑窗口内所有历元卫星运动状态。本发明历元更新时采用多步法，相交采用单步法其数值精度高、稳定性更好，积分右函数计算较少，利用已有状态转移矩阵信息，避免了在进行状态更新时需要更新状态转移矩阵，从而减少积分中右函数计算次数，有效提高计算速度，可以根据不同计算精度要求变换窗口长度，保证积分精度。
文档编号G05D1/00GK101853027SQ201010184669
公开日2010年10月6日 申请日期2010年5月21日 优先权日2010年5月21日
发明者叶世榕, 施闯, 赵齐乐, 郭靖 申请人:武汉大学