CN101872021A - Gps双频实时星载数据处理方法 - Google Patents

Abstract

本发明涉及星载数据处理技术领域，尤其涉及一种GPS双频实时星载数据处理方法。本发明包括利用窗口末端历元处的观测值，对整个窗口内所有历元卫星运动状态进行平滑、利用窗口内所有历元卫星运动状态和已有状态转移矩阵，采用多步法积分预报窗口末端历元至下一历元的状态转移矩阵和该历元卫星的运动状态、更新窗口历元，也即向后移动窗口至新的历元、利用预报历元新的观测信息平滑窗口内所有历元卫星运动状态等步骤。本发明历元更新时采用多步法，相交采用单步法其数值精度高、稳定性更好，积分右函数计算次数较少，可以有效提高计算速度，采用移动窗口方法，可以根据不同计算精度要求变换窗口长度，灵活多变且保证积分精度。

Description

GPS双频实时星载数据处理方法

技术领域

本发明涉及星载数据处理技术领域，尤其涉及一种GPS双频实时星载数据处理方法。

背景技术

在轨道确定中，参数估计主要包括两种方法：批处理算法和序贯处理算法。

所谓批处理就是将所有需要处理的观测资料一起解算，参考轨道可以循环迭代，数据处理精度相对较高，数值稳定性较好，因此该方法适合于事后精密定轨，但对估计参数初值、计算机的内存和计算速度要求较高。

所谓序贯处理算法是一种递推算法，对观测资料逐条处理，参考轨道不能循环迭代，对计算机的内存要求不高，可适用于实时定轨。

Kalman于1960年提出的滤波算法是最为常用的序贯处理算法，称为经典Kalman滤波。由于经典Kalman滤波存在如下三个缺陷：

1)如果参考轨道和真实轨道相差较大，使得线性化残差较大，导致滤波发散；

2)随着观测量增多状态参数的方差-协方差矩阵趋近于0，使得观测更新时先验状态权重过大，导致解对观测值不敏感；

3)在未知参数方差-协方差矩阵更新时，由于受到计算机位数限制，导致方差-协方差矩阵失去对称、正定性质。

针对以上Kalman滤波算法的缺点，许多学者分别提出相应算法加以解决。为了减少线性化带来的截断误差，通常采用推广Kalman滤波算法，该算法不使用事先预报的标称轨道作为线性化的初值，而是使用滤波估计的轨道进行线性化。为了避免经典Kalman滤波中未知参数方差-协方差矩阵趋近于0，在滤波状态方程中引入均值为0而协方差不为0阵的过程噪声，通过引入过程噪声在保持卫星运动状态估计值不变的情况下，可以通过状态噪声的方差-协方差有效改善卫星状态协方差矩阵特性，避免其趋于0阵。为了克服由于计算机内存限制所引起的滤波发散问题，美国喷气动力实验室(JPL)的C.L.Thornton于1967年提出均方根滤波算法，该算法将法阵的逆阵——信息矩阵根据Cholesky分解为一个上三角矩阵与其转置阵的乘积，这样避免了卫星状态的方差-协方差矩阵负定所引起的滤波发散，均方根信息滤波算法已成功应用于美国喷气动力实验室(JPL)开发的GPS数据处理软件GIPSY。

以上滤波算法，主要根据当前历元观测值来更新其预报值，但是对之前历元的卫星运动状态不做更新，此时采用更新的当前历元卫星运动状态积分预报下一历元卫星运动状态，则需要采用单步法积分求解，但是单步法需要多次计算右函数值，这样会增加计算时间，因此其不适合于快速实时精密定轨。如果采用多步法，由于之前历元卫星运动状态未进一步更新，如此会使得积分历元卫星运动状态的精度较低。

发明内容

针对上述存在的技术问题，本发明的目的是提供一种GPS双频实时星载数据处理方法，以解决实时精密定轨中快速估计卫星轨道的技术问题，提高卫星实时定轨精度和速度。

为达到上述目的，本发明采用如下的技术方案：

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

②利用窗口内所有历元卫星运动状态和已有状态转移矩阵，采用多步法积分预报窗口末端历元至下一历元的状态转移矩阵和该历元卫星的运动状态；

③更新窗口历元，向后移动窗口至新的历元；

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

⑤结束。

所述步骤②采用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{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=\{\begin{array}{cc}1,(m=0)& \\ 0,(m\≠0)& \end{array}\end{array}\right.$

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

所述步骤①具体采用如下子步骤：

首先引入平滑的简记符号：

${\hat{R}}_{i+N-1}{\hat{x}}_{i+N-1}^{i+N-1}={\hat{b}}_{i+N-1}$

${{\stackrel{\‾}{R}}_u}_j{\hat{u}}_{j-1}={\tilde{b}}_{u_j}-{{\stackrel{\‾}{R}}_{\mathrm{ux}}}_j{\hat{x}}_j^{i+N-1}$ j＝i+N-1，…，i

$x_{j-1}^{i+N-1}={\mathrm{\Φ}}^{-1}(t_i,t_{i-1})({\hat{x}}_i^m-\mathrm{\Γ}(t_i,t_{i-1}){\hat{u}}_{i-1}^{i+N-1})$

式中

表示利用t_i+N-1处的观测值来平滑t_j时刻的观测值，上标^表示估计值，

表示t_j时刻卫星状态向量x_k的估计值，

表示t_j时刻状态噪声向量u_k的估计值；

利用最小方差准则可以求得平滑的递推公式为：

$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]$

式中

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

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

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

整个平滑解为：

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

$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}=R_{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}.$

$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}$

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

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

2、采用移动窗口方法，可以根据不同计算精度要求变换窗口长度，灵活多变且保证积分精度；

3、在观测更新时同时平滑窗口内所有历元卫星运动状态，可以提高积分历元卫星运动状态预报精度。

附图说明

图1是本发明提出的GPS双频实时星载数据处理方法的流程图。

图2是本发明中数据流图。

具体实施方式

本发明提供的GPS双频实时星载数据处理方法，首先选择具有固定个数的历元(具有观测值或者预估值的时刻)构成一个数据窗口，然后利用窗口末端历元卫星观测值平滑窗口内所有历元卫星运动状态，从而在基于状态转移矩阵变化不大的假设下采用多步法数值积分卫星运动方程从而预报卫星在下一历元的运动状态，参见图1所示的流程图，具体包括以下步骤：

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

步骤S2：利用窗口内所有历元卫星运动状态和已有状态转移矩阵，采用多步法积分预报窗口末端历元至下一历元的状态转移矩阵和该历元卫星的运动状态；

步骤S3：更新窗口历元，也即向后移动窗口至新的历元；

步骤S4：判断当前历元是否存在观测值，如果有观测值存在则执行步骤1，否则执行步骤5；

步骤S5：结束。

在上述的GPS双频实时星载数据处理方法中，其不同于已有序贯方法之处在于，该方法采用移动窗口模式，在进行观测更新时不仅仅是更新观测历元处卫星状态，其同样对窗口内所有历元卫星运动状态进行平滑。从而在假设卫星状态变化对状态转移矩阵影响不大的情况下，采用多步法预报窗口末端下一历元卫星的运动状态，这样做不仅能有效保证积分精度，同样可以减少由于计算右函数所带来的时间损耗。

下面以具体实施例结合附图对本发明作进一步说明：

假定所选择窗口宽度为N，即该窗口内有N个历元，记为：t_i t_i+1…t_i+N-2 t_i+N-1。为了便于叙述，首先给出均方根信息滤波算法，其状态方程为：

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为卫星状态误差，其均值和先验方差分别为E[η_k-1]＝0，

而

(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[ε_k-1]＝0，

根据最小方差准则，也即使和ε_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}{\stackrel{\‾}{R}}_{k-1}\\ 0\end{array}\right]x_{k-1}-\left[\begin{array}{c}{\stackrel{\‾}{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)中对应值正交变换的结果。

GPS双频实时星载数据处理方法的数据处理步骤如下：

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

首先引入平滑的简记符号：

${\hat{R}}_{i+N-1}{\hat{x}}_{i+N-1}^{i+N-1}={\hat{b}}_{i+N-1}$

${\stackrel{\‾}{R}}_{u_j}{\hat{u}}_{j-1}={\tilde{b}}_{u_j}-{\stackrel{\‾}{R}}_{{\mathrm{ux}}_j}{\hat{x}}_j^{i+N-1}$ j＝i+N-1，…，i    (12)

$x_{j-1}^{i+N-1}={\mathrm{\Φ}}^{-1}(t_i,t_{i-1})({\hat{x}}_i^{i+N-1}-\mathrm{\Γ}(t_i,t_{i-1}){\hat{u}}_{i-1}^{i+N-1})$

式中

表示利用t_i+N-1处的观测值来平滑t_j时刻的观测值，上标^表示估计值，对于卫星状态向量x_k其估计值从式(8)中求得，状态噪声向量u_k估计值

由(11)式求得。

同样利用最小方差准则可以求得平滑的递推公式为：

$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}=R_{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}$

步骤2，利用窗口内所有历元卫星运动状态和已有状态转移矩阵，采用多步法积分预报窗口末端历元至下一历元的状态转移矩阵和该历元卫星的运动状态；

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

1)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{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)$

2)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}{}\begin{array}{c}1,(m=0)\\ 0,(m\≠0)\end{array}\end{array}\right.\end{array}\right.---\left(16\right)$

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

步骤3，更新窗口历元，也即向后移动窗口至新的历元；

这一步也即将窗口由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

步骤4，判断当前历元是否存在观测值，如果有观测值存在则执行步骤1，否则执行步骤5；

步骤5：结束。

附图2给出了本发明的数据流图，图中黄色为当前数据窗口其历元总数为N，而红色表示新增历元，灰色表示失效历元。横轴表示时间轴，而纵轴表示数据处理过程。整个数据处理算法由两个更新组成：状态(时间)更新和观测更新。状态更新也即利用窗口内卫星运动状态预估新增历元卫星运动状态，图中由窗口下方曲线表示；观测更新也即利用新增历元卫星的观测值，图中由窗口上方曲线表示。

上述实例用来解释说明本发明，而不是对本发明进行限制，在本发明的精神和权利要求的保护范围内，对本发明做出任何的修改和改变，都落入本发明的保护范围。

Claims (3)

1.一种GPS双频实时星载数据处理方法，其特征在于，包括以下步骤：

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

②利用窗口内所有历元卫星运动状态和已有状态转移矩阵，采用多步法积分预报窗口末端历元至下一历元的状态转移矩阵和该历元卫星的运动状态；

③更新窗口历元，向后移动窗口至新的历元；

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

⑤结束。

2.根据权利要求1所述的GPS双频实时星载数据处理方法，其特征在于：

所述步骤②采用Adams预估-校正方法来实现，其中Adams-Bashforth公式为：

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

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

3.根据权利要求1所述的GPS双频实时星载数据处理方法，其特征在于：

所述步骤①具体采用如下子步骤：

首先引入平滑的简记符号：

${\hat{R}}_{i+N-1}{\hat{x}}_{i+N-1}^{i+N-1}={\hat{b}}_{i+N-1}$

${\stackrel{\‾}{R}}_{u_j}{\stackrel{\‾}{u}}_{j-1}={\stackrel{\‾}{b}}_{u_j}-{\stackrel{\‾}{R}}_{{\mathrm{ux}}_j}{\hat{x}}_j^{i+N-1}$ j＝i+N-1，…，i

$x_{j-1}^{i+N-1}={\mathrm{\Φ}}^{-1}(t_i,t_{i-1})({\hat{x}}_i^m-\mathrm{\Γ}(t_i,t_{i-1}){\hat{u}}_{i-1}^{i+N-1})$

式中

表示利用t_i+N-1处的观测值来平滑t_j时刻的观测值，上标^表示估计值，

表示卫星状态向量x_k的估计值，

表示状态噪声向量u_k的估计值；

利用最小方差准则可以求得平滑的递推公式为：

$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]$

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

是对

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

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

整个平滑解为：

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

$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}=R_{u_{j-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}.$

$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}$

Images