CN107422342A

CN107422342A - Gnss卫星钟差实时估计质量控制方法

Info

Publication number: CN107422342A
Application number: CN201710658432.5A
Authority: CN
Inventors: 熊超; 崔红正; 周培源; 年劲飞; 杨赛男
Original assignee: Qianxun Position Network Co Ltd
Current assignee: Qianxun Position Network Co Ltd
Priority date: 2017-08-03
Filing date: 2017-08-03
Publication date: 2017-12-01
Also published as: WO2019024895A1

Abstract

本发明针对实时钟差估计对观测数据质量要求高、不能进行事后残差编辑的特点，提供了一种GNSS卫星钟差实时估计质量控制方法，包括两个方面：第一步是在参数估计之前，对原始观测数据的预处理，包括粗差剔除和周跳探测；第二步是在滤波观测更新之后，基于验后观测值残差对上一步中未探测出的粗差和周跳进行分析处理，具体包括以下步骤：建立平方根信息滤波算法模型，对噪声进行处理；进行平方根信息平滑；基于平方根信息滤波进行实时质量控制。本发明能实现实时钟差数据预处理与高效的质量控制，能高效地发现问题观测量，避免滤波器发散，从而实现实时钟差稳定估计，满足实时高精度定位的要求。

Description

GNSS卫星钟差实时估计质量控制方法

技术领域

本发明涉及质量控制技术领域，具体涉及一种GNSS卫星钟差实时估计质量控制方法。

背景技术

高质量的观测数据是实现高精度钟差确定的前提条件之一。尽管目前GNSS(Global Navigation SatelliteSystem，全球卫星导航系统)观测质量已大大提高，但在实际观测中，由于接收机的不稳定性、多路径效应影响及电离层闪烁等原因，常会造成卫星信号的失锁或观测值异常等数据质量问题，如果这些问题得不到较好的解决，就会直接影响GNSS实时钟差估计精度，甚至造成滤波发散。

GNSS实时钟差数据处理过程中，数学模型的有效性决定了参数估计的精度和正确性。观测量中的粗差会使平差模型偏离预先假设的数学模型，造成参数估计精度的下降，或是得到错误的平差结果。平差模型的有效性可以通过系统的验后残差进行检验，定位观测量中的粗差并加以消除。已有学者针对事后钟差估计的周跳和粗差问题进行了研究，并取得了较好的效果。但对实时钟差估计而言，基于验后残差的粗差探测与质量控制方法具有较大的局限性。实时钟差估计系统常常需要处理几十个甚至上百个参考站的大量观测数据，一方面相较于单站处理来说较难发现观测量中的粗差，另一方面，粗差探测的迭代适应过程相当耗时，系统的实时性难以保证。

平方根信息滤波(SRIF，The SquareRootInformationFilter)和平方根信息平滑(SRIS)算法，已成功应用于美国喷气动力实验室(JPL)开发的GPS(Global PositioningSystem，全球定位系统)数据处理软件GIPSY，其在处理GPS数据及其它卫星跟踪数据时能有效克服滤波器的发散，具有较高的数值稳健性和计算高效性。从最小二乘的基本原理出发，结合实时钟差估计的实际特征，导出了其在实时钟差估计软件中实现的详细公式。

尽管随着GPS数据预处理技术的不断发展和完善，部分周跳能被有效地探测，但仍然会有一定数量的周跳，特别是数值较小的周跳难以被发现。然而未被探测出来的周跳将极大地降低实时钟差估计的精度，因此自动地探测周跳是GPS数据处理阶段中最为关键的一环。针对平方根滤波的特点，提出了快速高效地探测和修复GPS观测数据中周跳的新方法。

因此，针对实时钟差估计系统的特点，研究高效地质量控制方法是非常有意义的，平方根信息滤波在采用Householder正交变换进行测量更新时，可以形成验后残差向量对观测值的敏感矩阵，能够保证输入系统的观测量的正确性，从而稳定地输出实时钟差满足实时定位用户的需求。

发明内容

本发明提出了一种GNSS卫星钟差实时估计质量控制方法，解决了实时钟差估计对观测数据质量要求高、不能进行事后残差编辑的技术问题。

本发明采用的技术方案是：

一种GNSS卫星钟差实时估计质量控制方法，其特征在于，包括以下步骤：

建立平方根信息滤波算法模型，对噪声进行处理；

进行平方根信息平滑；

基于平方根信息滤波进行实时质量控制。

进一步地，建立平方根信息滤波算法模型包括滤波器的状态方程和观测方程，分别如下式所示：

z_j＝A_P(j)P_j+A_X(j)X_j+A_Y(j)Y+v_j (12)

其中，X_j+1为第j+1步的随时间变化的状态参数，P_j+1为第j+1步的相关随机过程噪声参数，V_X(j)为第j步随时间变化的状态参数的系数矩阵，V_P(j)为第j步相关随机过程噪声参数的系数矩阵，V_Y(j)为第j步不随时间变化的状态参数的系数矩阵，I为单位矩阵，z_j为第j步的观测量，A_P(j)为第j步相关随机过程噪声参数的系数矩阵，A_X(j)为第j步随时间变化的状态参数的系数矩阵，A_Y(j)为第j步不随时间变化的状态参数的系数矩阵，v_j为第j步的误差矩阵，P_j为相关随机过程噪声参数；X_j为随时间变化的状态参数，但并不直接作为随机过程噪声处理的参数；Y为不随时间变化的状态参数；w_j为滤波器的随机过程噪声；M_j为一对角矩阵，取决于相关噪声类型，如果为一阶马尔科夫过程则可表示为：

M_j＝Diag[exp(-Δt_j/τ₁)，…，exp(-Δt_j/τ_np)] (13)

其中Δt_j＝t_j+1-t_j，为t_j+1时刻与t_j时刻的时刻差，τ_k为用户定义的第k个随机参数的相关时间，np为随机参数的总个数。

进一步地，对噪声进行处理具体包括以下步骤：

根据状态参数的类型，构建第j步平方根信息矩阵等式：

其中，为第j步的随时间变化的状态参数的验后权方差矩阵的平方根矩阵，为第j步的随时间变化的状态参数和相关随机过程噪声参数的验后权方差矩阵的平方根矩阵，为第j步的随时间变化的状态参数和不随时间变化的状态参数的验后权方差矩阵的平方根矩阵，为第j步的相关随机过程噪声参数的验后权方差矩阵的平方根矩阵，为第j步的相关随机过程噪声参数和不随时间变化的状态参数的验后权方差矩阵的平方根矩阵，为第j步的不随时间变化的状态参数的验后权方差矩阵的平方根矩阵，为第j步的验后随时间变化的状态参数，为第j步的相关随机过程噪声参数，为验后的随时间变化的状态参数，为验后的相关随机过程噪声参数，为验后的不随时间变化的状态参数；

根据状态方程公式(11)得到状态向量X从t_j到t_j+1时刻的投影方程：

X_j+1＝V_X(j)X_j+V_P(j)P_j+V_Y(j)Y (15)

将公式(15)代入到公式(14)得：

其中，V_X为第j步随时间变化的状态参数的系数矩阵，V_P为第j步相关随机过程噪声参数的系数矩阵，V_Y为第j步不随时间变化的状态参数的系数矩阵，为第j+1步先验随时间变化的状态参数；

考虑到公式(14)中：

P_j+1＝M_jP_j+w_j (17)

其中，M_j为第j过程噪声参数状态转移矩阵；

w_j是独立的随机过程噪声，并服从N(0，σ²)分布，用下式描述：

r_ww_j＝z_w (18)

其中，r_w为相应随机过程噪声的标准差，z_w为随机过程噪声的虚拟观测量；

综合公式(17)和公式(18)并代入公式(16)，并在等式两边同乘以正交矩阵算子可得：

其中，R_X为第j步变换后的随时间变化的状态参数的验后权方差矩阵的平方根矩阵，R_XP为第j步变换后的随时间变化的状态参数和相关随机过程噪声参数的验后权方差矩阵的平方根矩阵，R_XY为第j步变换后的随时间变化的状态参数和不随时间变化的状态参数的验后权方差矩阵的平方根矩阵，R_P为第j步变换后的相关随机过程噪声参数的验后权方差矩阵的平方根矩阵，R_PY为第j步变换后的相关随机过程噪声参数和不随时间变化的状态参数的验后权方差矩阵的平方根矩阵，R_Y为第j步的不随时间变化的状态参数的验后权方差矩阵的平方根矩阵，为第j+1步的先验随时间变化的状态参数，为第j+1步的先验相关随机过程噪声参数，z_X为验后的随时间变化的状态参数，z_P为验后的相关随机过程噪声参数，z_Y为验后的不随时间变化的状态参数；

通过正交变换得到公式(20)和公式(21)：

其中，为第j步正交矩阵算子变换后的相关随机过程噪声参数的验后权方差矩阵的平方根矩阵，为第j步正交矩阵算子变换后的相关随机过程噪声参数和随时间变化的状态参数的验后权方差矩阵的平方根矩阵，为第j步正交矩阵算子变换后的相关随机过程噪声参数和相关随机过程噪声参数的验后权方差矩阵的平方根矩阵，为第j步正交矩阵算子变换后的相关随机过程噪声参数和不随时间变化的状态参数的验后权方差矩阵的平方根矩阵，为第j步正交矩阵算子变换后的随机过程噪声的虚拟观测量；

其中，为第j步先验的随时间变化的状态参数的验后权方差矩阵的平方根矩阵，为第j步先验的随时间变化的状态参数和相关随机过程噪声参数的验后权方差矩阵的平方根矩阵，为第j步先验的随时间变化的状态参数和不随时间变化的状态参数的验后权方差矩阵的平方根矩阵，为第j步先验的相关随机过程噪声参数的验后权方差矩阵的平方根矩阵，为第j步先验的相关随机过程噪声参数和不随时间变化的状态参数的验后权方差矩阵的平方根矩阵，为第j步的不随时间变化的状态参数的先验权方差矩阵的平方根矩阵，为第j+1步的先验随时间变化的状态参数，为第j+1步的先验相关随机过程噪声参数，为第j+1步的先验的随时间变化的状态参数，为第j+1步的先验的相关随机过程噪声参数，为第j+1步的先验的不随时间变化的状态参数；

根据公式(21)为实现考虑相关噪声的状态更新步骤，公式(20)用于后向平方根信息平滑。

进一步地，如果在公式(21)底部增加第j+1步的观测量：

z_j+1＝A_P(j+1)P_j+1+A_X(j+1)X_j+1+A_Y(j+1)Y+v_j+1 (22)

就可在软件中实现平方根信息逐次滤波，其中，z_j+1为第j+1步的观测量，A_P(j+1)为第j+1步相关随机过程噪声参数的系数矩阵，P_j+1为第j+1步相关随机过程噪声参数，A_X(j+1)为第j+1步随时间变化的状态参数的系数矩阵，X_(j+1)为第j+1步随时间变化的状态参数，A_Y(j+1)为第j+1步不随时间变化的状态参数的系数矩阵，v_(j+1)为第j+1步的误差矩阵。

讲一步地，所述讲行平方根信息平滑具体包括以下步骤：

已知公式(15)中的V_X(j)、V_P(j)和V_Y(j)；

在每一步平方根信息滤波过程中保存公式(20)中的：和

最终平方根信息滤波结果；

假设已知第j+1步的平方根信息平滑结果[X_j+1 P_j+1 Y]^T，根据公式(20)式得：

根据公式(15)得：

X_j＝(V_X(j))^-1[X_j+1-V_P(j)P_j-V_Y(j)Y] (24)

由公式(23)和公式(24)得到第j步的平滑值[X_j P_j Y]^T，逐次递推，即可实现平方根信息平滑。

进一步地，基于平方根信息滤波进行实时质量控制具体包括以下步骤：

2)进行正交变换：

和

其中T₀为正交矩阵，A₀为误差方差系数矩阵，为先验权方差阵的平方根矩阵，为变换后的权方差阵的平方根矩阵，为虚拟观测量，z为实际观测量，为变换后的虚拟观测量；

得到历元j的正交变换：

和

其中，T_i为第i历元的一个正交矩阵，为第i历元的先验权方差的平方根矩阵，A_i为第i历元的实际观测值的系数矩阵，为第i历元的验后权方差的平方根矩阵，T_i的确定与有关，将公式(26)中T_i矩阵分解为验后残差e对先验值的敏感矩阵和验后残差对应观测量的敏感矩阵S得敏感方程：

目

设验后残差向量为v，对应的权阵为P，则有：e^Te＝v^TPv；

2)利用验后残差服从正态分布的统计特性构造统计检验量来探测模型误差，对验后残差e进行分析，构造检验统计量，进行χ²检验：

假设统计量：

H₀：：H₁：

其中，n为第i历元的观测值个数，为验后残差理论单位权中误差；若备选假设H₁成立则说明模型误差存在异常，需对异常误差历元进一步定位，可以实时地定位有问题的历元；

3)对有问题的历元的验后残差e进行假设检验，观测值z_m被探测为异常观测值，z_m为观测向量z中第m个观测值；假设异常观测值z_m与真实观测值具有偏差Δz_m，则利用正交变换的同时得到的异常观测值与验后残差的敏感矩阵S_m，S_m为S的第m列，获得观测偏差与残差的对应关系：

其中，为异常值的验后偏差，P_m为第i历元的观测值权矩阵，则引起的验后残差偏差

5)附加条件最小，可以获得偏差Δz_m，如果偏差Δz_m不小于一个周跳的量值，而且明显小于||e||²，则认为在z_m处发生了周跳，将滤波器状态向量扩维并增加一个整周模糊度参数，后续观测值将估计新的整周模糊度；

6)将e-Δe_m作为新的残差向量e_new，新的残差向量e_new与原验后残差具有相同的分布，对e_new进行假设检验，重新计算假设统计量：

H₀：H₁：如果备选假设H₁成立，则表明仍然有观测值被检验为异常观测值，则重复4)、5)步骤直到所有的观测值都被完全接受；

7)如果偏差已经确认，可以使用异常观测与验后残差之间的敏感方程，将(28)式作为新的观测值加入滤波观测方程更新中，滤波器中实现校正步骤。

进一步地，所属有问题的历元包括残差的方差达到观测值噪声的1.2倍或者单个残差平方达到观测值噪声方差的6倍。

本发明使用的是验后残差，而非预测残差，有益效果在于，能实现实时钟差数据预处理与高效的质量控制，能高效地发现问题观测量，避免滤波器发散，探测周跳准确可靠，从而实现实时钟差稳定估计，提高终端定位服务精度，提高播发的SSR(State-SpaceRepresentation，状态空间表示信息)改正信息中钟差的稳定性，满足时间同步和实时高精度定位对于钟差的需求。

具体实施方式

本发明从GNSS卫星实时钟差估计的观测值层面出发，解决了实时解算中的数据质量控制问题，包括两个方面：第一步是在参数估计之前，对原始观测数据的预处理，包括粗差剔除和周跳探测；第二步是在滤波观测更新之后，基于验后观测值残差对上一步中未探测出的粗差和周跳进行分析处理。

1)数据预处理

数据预处理是实时钟差解算质量控制的重要一环，尤其对于实时钟差估计系统而言，数据预处理的好坏将直接影响实时解算的效率和精度，预处理不当甚至会引起处理系统的崩溃。本发明的实时钟差估计是基于双频GNSS观测数据，因此解决的是基于双频观测数据的数据预处理问题。

2)滤波质量控制

数据预处理一般能准确探测出较大的粗差和周跳，但很难将一些小的周跳全部探测出来，故在滤波测量更新时需要对数据质量做进一步分析处理，以保证滤波解算的稳定可靠。

下文中结合实施例对本发明作进一步阐述。

1、平方根信息滤波基本原理

为了描述平方根信息滤波方法，不失一般性，假设t0时刻拥有一先验信息，并将其当成具有先验权方差阵的虚拟观测值由于阵的正定性可以用CholeskV分解将其分解成两上三角阵的乘积，表示如下：

式中为先验权方差阵的平方根矩阵，是一个上三角阵，则虚拟观测值可表示为：

式中为虚拟零均值随机向量。现如果得到实际观测：

z₀＝A₀x+v₀ (3)

式中，A₀为误差方差系数矩阵，x为估计参数矩阵，z₀为实际观测值，v₀为实际观测零均值随机向量，与虚拟观测组成正态方程：

在公式(4)两边同乘一正交矩阵T₀，进行QR分解(正交三角分解)得：

T₀矩阵可以利用Householder变换得到，则公式(5)可以化为：

式中为变换后的权方差阵的平方根矩阵，为变换后的虚拟观测量，e₀为变换后的实际观测量，为变换后的虚拟观测量的残差，v_e0为变换后的实际观测量的残差。

考虑到正交矩阵的特性，误差方差函数J(x)可以表示为：

其中由公式(6)知e₀＝v_e0，故||e₀||²为残差平方和，要使J(x)最小，则要求：

式中是一个上三角矩阵，则x的解可表示为下式：

x的解的协方差阵为：

式中，Γ为解的协方差阵，R₀为解的协方差阵的平方根矩阵。

如果不考虑随机过程噪声，可将t₀时刻变换后的信息矩阵作为下一步t₁时刻观测量的先验信息加入新的观测值可得到下一步的观测信息矩阵：重复由公式(4)-(6)进行下一步正交变换，即可实现逐次滤波过程。可以发现：如果不输出中间结果，就不需要求解方程(8)，整个处理过程只涉及到扩展的观测信息矩阵：...等，称之为平方根信息滤波。

其中，为t₁时刻的先验权方差阵的平方根矩阵，为t₁时刻的虚拟观测量，A₁为t₁时刻的误差方差系数矩阵，z₁为t₁时刻的实际观测量。

2、平方根信息滤波算法在实时钟差估计软件中的实现

1)建立平方根信息滤波算法模型

基于平方根信息滤波的基本原理，考虑到实现GNSS实时钟差估计的复杂性，可将待估状态参数分为以下三类进行处理：

(1)P_j为相关随机过程噪声参数，j为第j步(第j个历元)，例如：卫星及接收机钟差和对流层参数等可以作为此类参数处理。

(2)X_j为随时间变化的状态参数，但并不直接作为随机过程噪声处理的参数，例如：当加速度被认为随机量(动态定位，简化动力学定轨)时的位置和速度等。

(3)Y为不随时间变化的状态参数，包括站坐标，(动力学定轨时的)卫星星历参数，模糊度参数等。

此时，滤波器的状态方程和观测方程可表示为：

z_j＝A_P(j)P_j+A_X(j)X_j+A_Y(j)Y+v_j (12)

式中，X_j+1为第j+1步的随时间变化的状态参数，P_j+1为第j+1步的相关随机过程噪声参数，V_X(j)为第j步随时间变化的状态参数的系数矩阵，V_P(j)为第j步相关随机过程噪声参数的系数矩阵，V_Y(j)为第j步不随时间变化的状态参数的系数矩阵，I为单位矩阵，z_j为第j步的观测量，A_P(j)为第j步相关随机过程噪声参数的系数矩阵，A_X(j)为第j步随时间变化的状态参数的系数矩阵，A_Y(j)为第j步不随时间变化的状态参数的系数矩阵，v_j为第j步的误差矩阵，w_j为滤波器的随机过程噪声，M_j为一对角矩阵，取决于相关噪声类型，如果为一阶马尔科夫过程则可表示为：

M_j＝Diag[exp(-Δt_j/τ₁)，…，exp(-Δt_j/τ_np)] (13)

其中Δt_j＝t_j+1-t_j，为t_j+1时刻与t_j时刻的时刻差，到τ_k为用户定义的第k个随机参数的相关时间，np为随机参数的总个数。

2)相关噪声处理

考虑到上述不同类型的状态参数，根据公式(8)可得到第j步信息矩阵等式：

式中，为第j步的随时间变化的状态参数的验后权方差矩阵的平方根矩阵，为第j步的随时间变化的状态参数和相关随机过程噪声参数的验后权方差矩阵的平方根矩阵，为第j步的随时间变化的状态参数和不随时间变化的状态参数的验后权方差矩阵的平方根矩阵，为第j步的相关随机过程噪声参数的验后权方差矩阵的平方根矩阵，为第j步的相关随机过程噪声参数和不随时间变化的状态参数的验后权方差矩阵的平方根矩阵，为第j步的不随时间变化的状态参数的验后权方差矩阵的平方根矩阵，为第j步的验后随时间变化的状态参数，为第j步的相关随机过程噪声参数，为验后的随时间变化的状态参数，为验后的相关随机过程噪声参数，为验后的不随时间变化的状态参数。

根据状态方程公式(11)，可知状态向量X从t_j到t_j+1时刻的投影方程：

X_j+1＝V_X(j)X_j+V_P(j)P_j+V_Y(j)Y (15)

代入到公式(14)得：

式中，V_X为第j步随时间变化的状态参数的系数矩阵，V_P为第j步相关随机过程噪声参数的系数矩阵，V_Y为第j步不随时间变化的状态参数的系数矩阵，为第j+1步先验随时间变化的状态参数。

考虑到公式(14)中：

P_j+1＝M_jP_j+w_j (17)

式中，M_j为第j过程噪声参数状态转移矩阵。

由于w_j可以认为是一独立的随机过程噪声，并服从N(0，σ²)分布(正太分布)，可以用下式描述：

r_ww_j＝z_w (18)

其中，z_w为随机过程噪声的虚拟观测量，r_w为相应过程噪声的标准差，而r_w考虑到w_j为零均值，可以设为零。

此时综合公式(17)和公式(18)并代入公式(16)，并在等式两边同乘以正交矩阵算子可得：

式中，R_X为第j步变换后的随时间变化的状态参数的验后权方差矩阵的平方根矩阵，R_XP为为第j步变换后的随时间变化的状态参数和相关随机过程噪声参数的验后权方差矩阵的平方根矩阵，R_XY为第j步变换后的随时间变化的状态参数和不随时间变化的状态参数的验后权方差矩阵的平方根矩阵，R_P为第j步变换后的相关随机过程噪声参数的验后权方差矩阵的平方根矩阵，R_PY为第j步变换后的相关随机过程噪声参数和不随时间变化的状态参数的验后权方差矩阵的平方根矩阵，R_Y为第j步的不随时间变化的状态参数的验后权方差矩阵的平方根矩阵，为第j+1步的先验随时间变化的状态参数，为第j+1步的先验相关随机过程噪声参数，z_X为验后的随时间变化的状态参数，z_P为验后的相关随机过程噪声参数，z_Y为验后的不随时间变化的状态参数。

完成公式(19)，正交变换可得到：

式中，为第j步正交矩阵算子变换后的相关随机过程噪声参数的验后权方差矩阵的平方根矩阵，为第j步正交矩阵算子变换后的相关随机过程噪声参数和随时间变化的状态参数的验后权方差矩阵的平方根矩阵，为第j步正交矩阵算子变换后的相关随机过程噪声参数和相关随机过程噪声参数的验后权方差矩阵的平方根矩阵，为第j步正交矩阵算子变换后的相关随机过程噪声参数和不随时间变化的状态参数的验后权方差矩阵的平方根矩阵，为第j步正交矩阵算子变换后的随机过程噪声的虚拟观测量。

式中，为第j步先验的随时间变化的状态参数的验后权方差矩阵的平方根矩阵，为第j步先验的随时间变化的状态参数和相关随机过程噪声参数的验后权方差矩阵的平方根矩阵，为第j步先验的随时间变化的状态参数和不随时间变化的状态参数的验后权方差矩阵的平方根矩阵，为第j步先验的相关随机过程噪声参数的验后权方差矩阵的平方根矩阵，为第j步先验的相关随机过程噪声参数和不随时间变化的状态参数的验后权方差矩阵的平方根矩阵，为第j步的不随时间变化的状态参数的先验权方差矩阵的平方根矩阵，为第j+1步的先验随时间变化的状态参数，为第j+1步的先验相关随机过程噪声参数，为第j+1步的先验的随时间变化的状态参数，为第j+1步的先验的相关随机过程噪声参数，为第j+1步的先验的不随时间变化的状态参数。

根据公式(21)可实现考虑相关噪声的状态更新步骤，而公式(20)与后续滤波过程无关，但可存储于一文件中，用于后向平方根信息平滑。如果在(21)式底部增加第j+1步的观测量：

z_j+1＝A_P(j+1)P_j+1+A_X(j+1)X_j+1+A_Y(j+1)Y+v_j+1 (22)

式中，z_j+1为第j+1步的观测量，A_P(j+1)为第j+1步相关随机过程噪声参数的系数矩阵，P_j+1为第j+1步相关随机过程噪声参数，A_X(j+1)为第j+1步随时间变化的状态参数的系数矩阵，X_(j+1)为第j+1步随时间变化的状态参数，A_Y(j+1)为第j+1步不随时间变化的状态参数的系数矩阵，v_(j+1)为第j+1步的误差矩阵。

类似公式(4)-(6)就可在软件中实现平方根信息逐次滤波。

观察公式(14)和公式(21)，可以发现等式左边的信息矩阵具有相同的维数，在程序设计过程中可以使用相同数组存储，这样既可节省计算机存储空间又可改善计算效率。

3、平方根信息平滑的实现

软件中实现平方根信息平滑需要在滤波的同时准备如下必备信息：

1)已知公式(15)中的：V_X(j)、V_P(j)和V_Y(j)。

2)在每一步平方根信息滤波过程中保存公式(20)中的：和

3)最终平方根信息滤波结果。

现假设已知第j+1步的平滑结果[X_j+1 P_j+1 Y]^T，根据公式(20)得：

根据公式(15)得：

X_j＝(V_X(j))^-1[X_j+1-V_P(j)P_j-V_Y(j)Y] (24)

由(23)式和(24)式可以得到第j步的平滑值[X_j P_j Y]^T逐次递推，即可在软件中实现平方根信息平滑。

4、基于平方根信息滤波的实时质量控制

平方根滤波器在进行正交变换的同时，可构成验后残差对应观测量的敏感方程(即公式(27))，充分利用这一特点可在滤波器中构建一个高效的质量控制模块。

具体实施步骤如下：

1)首先类似于公式(5)进行正交变换：

和

式中，为虚拟观测量，z为实际观测量，为变换后的虚拟观测量。

可得到历元i的正交变换：

和

其中T_i为第i历元的一个正交矩阵，为第i历元的先验权方差的平方根矩阵，A_i为第i历元的实际观测值的系数矩阵，为第i历元的验后权方差的平方根矩阵，T_i的确定与有关，而与无关，将公式(26)式中T_i矩阵分解为验后残差e对先验值的敏感矩阵和验后残差对应观测量的敏感矩阵S得敏感方程：

日

2)对验后残差e进行分析，可以实时地定位有问题的历元，比如残差的方差达到观测值噪声的1.2倍或者单个残差平方达到观测噪声方差的6倍，则认为此历元可能出现粗差观测。

3)对有问题的历元的验后残差e进行假设检验，某些观测值z_m就被探测为异常观测值。

4)假设异常观测z_m(观测向量z中第m个观测值)与真实观测具有相应的偏差Δz_m，则利用正交变换的同时得到的异常观测值与验后残差的敏感矩阵S_m(敏感矩阵S的第m列)，可以获得观测偏差与残差的对应关系：

式中，为异常值的验后偏差，Δe_m为残差向量。

5)附加条件最小，可以获得异常观测值的偏差Δz_m，如果异常观测值的偏差Δz_m不小于一个跳周的量值，而且明显小于||e||²，则认为在z_m处发生了周跳，将滤波器状态向量扩维并增加一整周模糊度参数，后续观测值将估计新的整周模糊度。

6)将e-Δe_m作为新的残差向量e_new，新残差向量与原验后残差应具有相同的分布，对e_new进行假设检验，如果仍然有观测值被检验为异常观测，则重复4)、5)两步操作直到所有的观测值都被完全接受。

7)如果偏差已经确认，使用相同的敏感方程，在滤波器中就很容易实现校正步骤。

上述操作充分顾及了平方根信息滤波的特点，利用问题观测值及残差之间的敏感矩阵，通过增加观测值的偏差量，对问题观测值进行松弛。该方法的实际效果相当于在该历元增加了一新的模糊度参数，但却不需要利用观测量重新求解，极大地减少了计算量，具有较高的计算效率。

本发明虽然已以较佳实施例公开如上，但其并不是用来限定本发明，任何本领域技术人员在不脱离本发明的精神和范围内，都可以利用上述揭示的方法和技术内容对本发明技术方案做出可能的变动和修改，因此，凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化及修饰，均属于本发明技术方案的保护范围。

Claims

1.一种GNSS卫星钟差实时估计质量控制方法，其特征在于，包括以下步骤：

建立平方根信息滤波算法模型，对噪声进行处理；

进行平方根信息平滑；

基于平方根信息滤波进行实时质量控制。

2.如权利要求1所述的一种GNSS卫星钟差实时估计质量控制方法，其特征在于，建立平方根信息滤波算法模型包括滤波器的状态方程和观测方程，分别如下式所示：

z_j＝A_P(j)P_j+A_X(j)X_j+A_Y(j)Y+v_j (12)

M_j＝Diag[exp(-Δt_j/τ₁)，…，exp(-Δt_j/τ_np)] (13)

3.如权利要求2所述的一种GNSS卫星钟差实时估计质量控制方法，其特征在于，对噪声进行处理具体包括以下步骤：

根据状态参数的类型，构建第j步平方根信息矩阵等式：

X_j+1＝V_X(j)X_j+V_P(j)P_j+V_Y(j)Y (15)

将公式(15)代入到公式(14)得：

考虑到公式(14)中：

P_j+1＝M_jP_j+w_j (17)

其中，M_j为第j过程噪声参数状态转移矩阵；

r_ww_j＝z_w (18)

通过正交变换得到公式(20)和公式(21)：

根据公式(21)可实现考虑相关噪声的状态更新步骤，公式(20)用于后向平方根信息平滑。

4.如权利要求3所述的一种GNSS卫星钟差实时估计质量控制方法，其特征在于，如果在公式(21)底部增加第j+1步的观测量：

z_j+1＋A_P(j+1)P_j+1+A_X(j+1)X_j+1+A_Y(j+1)Y+v_j+1 (22)

5.如权利要求3所述的一种GNSS卫星钟差实时估计质量控制方法，其特征在于，所述进行平方根信息平滑具体包括以下步骤：

已知公式(15)中的V_X(j)、V_P(j)和V_Y(j)；

在每一步平方根信息滤波过程中保存公式(20)中的：和

最终平方根信息滤波结果；

根据公式(15)得：

X_j＝(V_X(j))^-1[X_j+1-V_P(j)P_j-V_Y(j)Y] (24)

6.如权利要求5所述的一种GNSS卫星钟差实时估计质量控制方法，其特征在于，基于平方根信息滤波进行实时质量控制具体包括以下步骤：

1)进行正交变换：

和

得到历元i的正交变换：

和

且

设验后残差向量为v，对应的权阵为P，则有：e^Te＝v^TPv.

假设统计量：

<mrow> <msup> <mover> <mi>&sigma;</mi> <mo>^</mo> </mover> <mn>2</mn> </msup> <mo>=</mo> <mfrac> <mrow> <msup> <mi>v</mi> <mi>T</mi> </msup> <mi>P</mi> <mi>v</mi> </mrow> <mi>n</mi> </mfrac> <mo>=</mo> <mfrac> <mrow> <msup> <mi>e</mi> <mi>T</mi> </msup> <mi>e</mi> </mrow> <mi>n</mi> </mfrac> </mrow>

<mfenced open = "" close = ""> <mtable> <mtr> <mtd> <mrow> <msub> <mi>H</mi> <mn>0</mn> </msub> <mo>:</mo> <mi>E</mi> <mrow> <mo>(</mo> <msup> <mover> <mi>&sigma;</mi> <mo>^</mo> </mover> <mn>2</mn> </msup> <mo>)</mo> </mrow> <mo>=</mo> <msubsup> <mi>&sigma;</mi> <mn>0</mn> <mn>2</mn> </msubsup> <mo>,</mo> </mrow> </mtd> <mtd> <mrow> <msub> <mi>H</mi> <mn>1</mn> </msub> <mo>:</mo> <mi>E</mi> <mrow> <mo>(</mo> <msup> <mover> <mi>&sigma;</mi> <mo>^</mo> </mover> <mn>2</mn> </msup> <mo>)</mo> </mrow> <mo>></mo> <msubsup> <mi>&sigma;</mi> <mn>0</mn> <mn>2</mn> </msubsup> </mrow> </mtd> </mtr> </mtable> </mfenced>

<mrow> <msub> <mi>S</mi> <mi>m</mi> </msub> <mo>&times;</mo> <msub> <mi>&Delta;z</mi> <mi>m</mi> </msub> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mi>&Delta;</mi> <msub> <mover> <mi>z</mi> <mo>^</mo> </mover> <mi>m</mi> </msub> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>&Delta;e</mi> <mi>m</mi> </msub> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>28</mn> <mo>)</mo> </mrow> </mrow> 4

<mrow> <msup> <msub> <mover> <mi>&sigma;</mi> <mo>^</mo> </mover> <mi>b</mi> </msub> <mn>2</mn> </msup> <mo>=</mo> <mfrac> <mrow> <msup> <msub> <mi>e</mi> <mrow> <mi>n</mi> <mi>e</mi> <mi>w</mi> </mrow> </msub> <mi>T</mi> </msup> <msub> <mi>e</mi> <mrow> <mi>n</mi> <mi>e</mi> <mi>w</mi> </mrow> </msub> </mrow> <mi>n</mi> </mfrac> </mrow>

如果备选假设H₁成立，则表明仍然有观测值被检验为异常观测值，则重复4)、5)步骤直到所有的观测值都被完全接受；

7.如权利要求5所述的一种GNSS卫星钟差实时估计质量控制方法，其特征在于，所属有问题的历元包括残差的方差达到观测值噪声的1.2倍或者单个残差平方达到观测值噪声方差的6倍。