CN101853027A - 实时精密定轨中轨道的星载快速多步积分方法 - Google Patents

Abstract

本发明涉及导航卫星应用技术领域，尤其涉及一种实时精密定轨中轨道的星载快速多步积分方法。本发明选择积分窗口，利用窗口末端历元处的观测值对整个窗口内所有历元卫星运动状态进行平滑，利用窗口内所有历元卫星更新后运动状态和已有状态转移矩阵积分计算新的历元卫星运动状态，然后更新窗口历元，也即使窗口起始历元后移一个历元，利用更新窗口末端历元的观测信息，平滑窗口内所有历元卫星运动状态。本发明历元更新时采用多步法，相交采用单步法其数值精度高、稳定性更好，积分右函数计算较少，利用已有状态转移矩阵信息，避免了在进行状态更新时需要更新状态转移矩阵，从而减少积分中右函数计算次数，有效提高计算速度，可以根据不同计算精度要求变换窗口长度，保证积分精度。

Description

实时精密定轨中轨道的星载快速多步积分方法

技术领域

本发明涉及导航卫星应用技术领域，尤其涉及一种实时精密定轨中轨道的星载快速多步积分方法。

背景技术

在精密定轨中为了获得高精度的卫星轨道和状态转移矩阵，需要利用数值积分算法求解卫星运动方程。当前在定轨中主要使用的积分算法为Runge-Kutta单步法和Adams多步法。Runge-Kutta法间接引用泰勒展开式，用区间[t_i t_i+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公式：

$\left\{\begin{array}{c}{x}_{n+1}=x_n+h\underset{i=0}{\overset{k-1}{\mathrm{\Σ}}}{\mathrm{\β}}_{\mathrm{ki}}{f}_{n-i}\\ {\mathrm{\β}}_{\mathrm{ki}}={(-1)}^i\underset{\text{m=1}}{\overset{k-1}{\mathrm{\Σ}}}\left(\begin{array}{c}m\\ i\end{array}\right){\mathrm{\γ}}_m,(k=\mathrm{1,2,3},...)\\ {\mathrm{\γ}}_m+\frac{1}{2}{\mathrm{\γ}}_{m-1}+\frac{1}{3}{\mathrm{\γ}}_{m-2}+...+\frac{1}{m+1}{\mathrm{\γ}}_0=1\end{array}\right.$

Adams-Moulton公式：

$\left\{\begin{array}{c}{x}_{n+1}=x_n+h\underset{i=0}{\overset{k-1}{\mathrm{\Σ}}}{\mathrm{\β}}_{\mathrm{ki}}{f}_{n+1-i}\\ {\mathrm{\β}}_{\mathrm{ki}}={(-1)}^i\underset{m=1}{\overset{k-1}{\mathrm{\Σ}}}\left(\begin{array}{c}m\\ i\end{array}\right){\mathrm{\γ}}_m,(k=\mathrm{1,2,3},...)\\ {\mathrm{\γ}}_m+\frac{1}{2}{\mathrm{\γ}}_{m-1}+\frac{1}{3}{\mathrm{\γ}}_{m-2}+...+\frac{1}{m+1}{\mathrm{\γ}}_0=\left\{\begin{array}{c}1,(m=0)\\ 0,(m\≠0)\end{array}\right.\end{array}\right.$

x_n表示在第n个积分历元卫星的运动状态，h为积分步长，f为积分函数，而β和γ为多步法积分系数；

其中右函数f只需重新在校正时计算一次。

本发明具有以下优点和积极效果：

1)历元更新时采用多步法，相较采用单步法其数值精度高、稳定性更好，积分右函数计算较少；

2)利用已有状态转移矩阵信息，避免了在进行状态更新时需要更新状态转移矩阵，从而减少积分中右函数计算次数，有效提高计算速度；

3)采用移动窗口方法，可以根据不同计算精度要求变换窗口长度，保证积分精度。

附图说明

图1是本发明提供的实时精密定轨中轨道的星载快速多步积分方法的流程图。

图2是本发明提供的实时精密定轨中轨道的星载快速多步积分方法的数据流图。

具体实施方式

假设卫星运动状态在进行状态更新时其值变化不大，那么由更新值计算得到的状态转移矩阵变化也不大，这样可以直接利用已有卫星状态转移矩阵，而无需利用更新后的卫星状态重新计算状态转移矩阵，这样在采用多步法预报卫星在新的历元运动状态时，仅仅需要计算当前历元到下一历元的右函数，这样右函数只需计算一次从而大大减小计算时耗，提高运算速度。

本发明提供的实时精密定轨中轨道星载快速多步积分方法，其流程参见图1，具体包括以下步骤：

步骤S101：选择积分窗口，一般为多步法所需求的点数，如果积分窗口选择为1，则退化为单步法；

假定所选择窗口的宽度为N，即窗口内有N个历元，记为：

t_i t_i+1…t_i+N-2 t_i+N-1

步骤S102：利用窗口末端历元处的观测值，对整个窗口内所有历元卫星运动状态进行平滑，从而得到窗口内所有历元卫星状态平滑值；

这一步采用均方根信息平滑算法完成，其状态方程为：

x_k＝Ф(t_k，t_k-1)x_k-1+Γ(t_k，t_k-1)u_k-1    (1)

式中x_k和x_k-1分别是卫星在t_k和t_k-1的运动状态，Ф(t_k，t_k-1)是从t_k-1时刻到t_k时刻状态转移矩阵，Γ(t_k，t_k-1)从t_k-1时刻到t_k时刻噪声转移矩阵，u_k-1为t_k-1时刻状态噪声。式中x_k-1具有先验值

和先验方差

将先验方差进行Cholesky分解，构造虚拟观测方程：

${\stackrel{\‾}{b}}_{k-1}={\stackrel{\‾}{R}}_{k-1}{x}_{k-1}+{\stackrel{\‾}{\mathrm{\η}}}_{k-1}---\left(2\right)$

式中η_k-1为卫星状态误差，其均值和先验方差分别为

而

式中u_k-1先验值与真值关系可用下式来描述：

${\stackrel{\‾}{u}}_{k-1}=u_{k-1}+{\mathrm{\α}}_{k-1}---\left(3\right)$

式中α_k-1为t_k-1时刻的状态噪声误差，其均值和方差为E[α_k-1]＝0，E[α_k-1α_k-1 ^T]＝Q，从而构建状态噪声u_k-1的虚拟观测方程：

${\stackrel{\‾}{b}}_{u_{k-1}}=R_u{\stackrel{\‾}{u}}_{k-1}=R_u{u}_{k-1}+{\stackrel{\‾}{\mathrm{\α}}}_{k-1}---\left(4\right)$

式中

而

而滤波的观测方程为：

y_k-1＝H_k-1x_k-1+ε_k-1    (5)

式中y_k-1为观测向量，H_k-1为设计矩阵，ε_k-1为观测值误差，其均值和方差分别为 $E\left[{\mathrm{\ϵ}}_{k-1}\right]=0,E\left({\mathrm{\ϵ}}_{k-1}{\mathrm{\ϵ}}_{k-1}^T\right)=I.$

根据最小方差准则，也即使和ε_k-1平方和最小，从而可以构建均方根信息滤波算法观测更新性能函数

${\hat{J}}_{k-1}={\left|\right|{\stackrel{\‾}{R}}_{k-1}{x}_{k-1}-{\stackrel{\‾}{b}}_{k-1}\left|\right|}^2+{(H_{k-1}{x}_{k-1}-y_{k-1})}^2+{\left|\right|R_u{u}_{k-1}-{\stackrel{\‾}{b}}_{u_{k-1}}\left|\right|}^2---\left(6\right)$

式中‖‖表示任意向量的2范数。

将(6)式写成矩阵形式可得：

${\hat{J}}_{k-1}={\left|\right|\left[\begin{array}{c}{\stackrel{\‾}{R}}_{k-1}\\ H_{k-1}\end{array}\right]x_{k-1}-\left[\begin{array}{c}{\stackrel{\‾}{b}}_{k-1}\\ y_{k-1}\end{array}\right]\left|\right|}^2+{\left|\right|R_u{u}_{k-1}-{\stackrel{\‾}{b}}_{u_{k-1}}\left|\right|}^2---\left(7\right)$

对(7)式进行正交变化可以得到：

${\hat{J}}_{k-1}={\left|\right|\left[\begin{array}{c}{\hat{R}}_{k-1}\\ 0\end{array}\right]x_{k-1}-\left[\begin{array}{c}{\hat{b}}_{k-1}\\ e_{k-1}\end{array}\right]\left|\right|}^2+{\left|\right|R_u{u}_{k-1}-{\stackrel{\‾}{b}}_{u_{k-1}}\left|\right|}^2---\left(8\right)$

式中

和e_k-1分别是

和y_k-1进行正交变换的结果。

同样可以根据最小方差准则，可以构建均方根信息滤波算法状态更新性能函数

${\stackrel{\‾}{J}}_k={\left(e_{k-1}\right)}^2+{\left|\right|{\hat{R}}_{k-1}{\mathrm{\Φ}}^{-1}(t_k,t_{k-1})(x_k-\mathrm{\Γ}(t_k,t_{k-1})u_{k-1})-{\hat{b}}_{k-1}\left|\right|}^2+{\left|\right|R_u{u}_{k-1}-{\stackrel{\‾}{b}}_{u_{k-1}}\left|\right|}^2---\left(9\right)$

将(9)式写成矩阵形式：

${\stackrel{\‾}{J}}_k={\left(e_{k-1}\right)}^2+\left|\right|\left[\begin{array}{cc}{R}_u& 0\\ -{\tilde{R}}_k\mathrm{\Γ}(t_k,t_{k-1})& {\tilde{R}}_k\end{array}\right]\left[\begin{array}{c}{u}_{k-1}\\ x_k\end{array}\right]-\left[\begin{array}{c}{\stackrel{\‾}{b}}_{u_{k-1}}\\ {\hat{b}}_{k-1}\end{array}\right]\left|\right|---\left(10\right)$

式中

将(10)式做正交变换可以得到：

${\stackrel{\‾}{J}}_k={\left(e_{k-1}\right)}^2+\left|\right|\left[\begin{array}{cc}{\stackrel{\‾}{R}}_{u_k}& {\stackrel{\‾}{R}}_{{\mathrm{ux}}_k}\\ 0& {\stackrel{\‾}{R}}_k\end{array}\right]\left[\begin{array}{c}{u}_{k-1}\\ x_k\end{array}\right]-\left[\begin{array}{c}{\tilde{b}}_{u_{k-1}}\\ {\stackrel{\‾}{b}}_{k-1}\end{array}\right]\left|\right|---\left(11\right)$

其中

和

都为式(10)中对应值正交变换的结果。

而求解卫星运动状态平滑解的递推公式为：

$T_{j-1}^{*}\left[\begin{array}{ccc}{\stackrel{\‾}{R}}_{u_j}+{\stackrel{\‾}{R}}_{{\mathrm{ux}}_j}\mathrm{\Γ}(t_j,t_{j-1})& {\stackrel{\‾}{R}}_{{\mathrm{ux}}_j}\mathrm{\Γ}(t_j,t_{j-1})& {\tilde{b}}_{u_j}\\ R_j^{*}\mathrm{\Γ}(t_j,t_{j-1})& R_j^{*}\mathrm{\Γ}(t_j,t_{j-1})& b_j^{*}\end{array}\right]=\left[\begin{array}{ccc}{R}_{u_{j-1}}^{*}& R_{{\mathrm{ux}}_{j-1}}^{*}& b_{u_{j-1}}^{*}\\ 0& R_{j-1}^{*}& b_{j-1}^{*}\end{array}\right]---\left(12\right)$

式中为由t_j时刻平滑t_j-1时刻正交变换矩阵，

与式(11)中含义相同，Γ(t_j，t_j-1)与式(1)中含义相同，

是对进行正交变化的结果，其起始值为

上式右边()^*都表示为左边中对应值的正交变换值。

由此可得平滑解为：

$R_{j-1}^{*}{\hat{x}}_{j-1}^{i+N-1}=b_{j-1}^{*}$ (13)

$R_{u_{j-1}}^{*}{\hat{u}}_{j-1}^{i+N-1}=b_{u_{j-1}}^{*}-R_{{\mathrm{ux}}_{j-1}}^{*}{\hat{x}}_{j-1}^m$

其相应的方差和协方差矩阵为：

$P_{x_{j-1}}^{t+N-1}=R_{j-1}^{*-1}{R}_{j-1}^{*-T}$

$P_{u_{j-1}}^{t+N-1}=P_{u_{j-1}}^{*-1}{R}_{u_{j-1}}^{*-T}+R_{u_{j-1}}^{*-1}{R}_{{\mathrm{ux}}_{j-1}}^{*}{R}_{j-1}^{*-1}{R}_{j-1}^{*-T}{R}_{{\mathrm{ux}}_{j-1}}^{*T}{R}_{u_{j-1}}^{*-T}---\left(14\right)$

$P_{{\mathrm{ux}}_{j-1}}^{i+N-1}=-R_{u_{j-1}}^{*-1}{R}_{{\mathrm{ux}}_{j-1}}^{*}{R}_{j-1}^{*-1}{R}_{j-1}^{*-T}$

步骤S103：利用窗口内所有历元卫星运动状态和已有状态转移矩阵，采用Adams预估-校正算法预报窗口末端历元至下一历元的状态转移矩阵和该历元卫星的运动状态；

这一步采用Adams预估-校正算法，其相应公式如下所述。

①Adams-Bashforth公式(隐式计算公式)：

$\left\{\begin{array}{c}{x}_{n+1}=x_n+h\underset{i=0}{\overset{k-1}{\mathrm{\Σ}}}{\mathrm{\β}}_{\mathrm{ki}}{f}_{n-i}\\ {\mathrm{\β}}_{\mathrm{ki}}={(-1)}^i\underset{\text{m=1}}{\overset{k-1}{\mathrm{\Σ}}}\left(\begin{array}{c}m\\ i\end{array}\right){\mathrm{\γ}}_m,(k=\mathrm{1,2,3},...)\\ {\mathrm{\γ}}_m+\frac{1}{2}{\mathrm{\γ}}_{m-1}+\frac{1}{3}{\mathrm{\γ}}_{m-2}+...+\frac{1}{m+1}{\mathrm{\γ}}_0=1\end{array}\right.---\left(15\right)$

②Adams-Moulton公式(显示计算公式)：

$\left\{\begin{array}{c}{x}_{n+1}=x_n+h\underset{i=0}{\overset{k-1}{\mathrm{\Σ}}}{\mathrm{\β}}_{\mathrm{ki}}{f}_{n+1-i}\\ {\mathrm{\β}}_{\mathrm{ki}}={(-1)}^i\underset{m=1}{\overset{k-1}{\mathrm{\Σ}}}\left(\begin{array}{c}m\\ i\end{array}\right){\mathrm{\γ}}_m,(k=\mathrm{1,2,3},...)\\ {\mathrm{\γ}}_m+\frac{1}{2}{\mathrm{\γ}}_{m-1}+\frac{1}{3}{\mathrm{\γ}}_{m-2}+...+\frac{1}{m+1}{\mathrm{\γ}}_0=\left\{\begin{array}{c}1,(m=0)\\ 0,(m\≠0)\end{array}\right.\end{array}\right.---\left(16\right)$

(15)和(16)两式中，x_n表示在第n个积分历元卫星的运动状态，h为积分步长，f为积分函数，而β和γ为多步法积分系数。其中右函数f只需在校正时重新计算一次。

步骤S104：更新窗口历元，也即使窗口起始历元后移一个历元；

这一步也即将窗口由t_i t_i+1…t_i+N-2 t_i+N-1更新为t_i+1 t_i+2…t_i+N-1 t_i+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：由卫星初始状态积分得到的状态转移矩阵

Claims (2)

1.一种实时精密定轨中轨道的星载快速多步积分方法，其特征在于，包括：

①选择积分窗口；

②利用窗口末端历元处的观测值，对整个窗口内所有历元卫星运动状态进行平滑，从而得到窗口内所有历元卫星状态平滑值；

③利用窗口内所有历元卫星运动状态和已有状态转移矩阵，采用Adams预估-校正算法预报窗口末端历元至下一历元的状态转移矩阵和该历元卫星的运动状态；

④更新窗口历元，使窗口起始历元后移一个历元；

⑤判断当前历元是否存在观测值，如果有观测值存在则执行步骤②，否则执行步骤⑥；

⑥结束。

2.根据权利要求1中所述的实时精密定轨中轨道的星载快速多步积分方法，其特征在于：

所述步骤③采用Adams预估-校正方法：

Adams-Bashforth公式：

$\left\{\begin{array}{cc}{x}_{n+1}=x_n+h\underset{i=0}{\overset{k-1}{\mathrm{\Σ}}}{\mathrm{\β}}_{\mathrm{ki}}{f}_{n-i}& \\ {\mathrm{\β}}_{\mathrm{ki}}={(-1)}^i\underset{m=1}{\overset{k-1}{\mathrm{\Σ}}}\left(\begin{array}{c}m\\ i\end{array}\right){\mathrm{\γ}}_m& (k=\mathrm{1,2,3},\·\·\·)\\ {\mathrm{\γ}}_m+\frac{1}{2}{\mathrm{\γ}}_{m-1}+\frac{1}{3}{\mathrm{\γ}}_{m-2}+\·\·\·+\frac{1}{m+1}{\mathrm{\γ}}_0=1& \end{array}\right.$

Adams-Moulton公式：

$\left\{\begin{array}{cc}{x}_{n+1}=x_n+h\underset{i=0}{\overset{k-1}{\mathrm{\Σ}}}{\mathrm{\β}}_{\mathrm{ki}}{f}_{n+1-i}& \\ {\mathrm{\β}}_{\mathrm{ki}}={(-1)}^i\underset{m=1}{\overset{k-1}{\mathrm{\Σ}}}\left(\begin{array}{c}m\\ i\end{array}\right){\mathrm{\γ}}_m& (k=\mathrm{1,2,3},\·\·\·)\\ {\mathrm{\γ}}_m+\frac{1}{2}{\mathrm{\γ}}_{m-1}+\frac{1}{3}{\mathrm{\γ}}_{m-2}+\·\·\·+\frac{1}{m+1}{\mathrm{\γ}}_0=\left\{\begin{array}{c}1,(m=0)\\ 0,(m\≠0)\end{array}\right.& \end{array}\right.$

x_n表示在第n个积分历元卫星的运动状态，h为积分步长，f为积分函数，而β和γ为多步法积分系数；

其中右函数f只需重新在校正时计算一次。

Images