CN110059361B - 一种基于抗差卡尔曼滤波算法的实时区域对流层建模方法 - Google Patents

Abstract

一种基于抗差卡尔曼滤波算法的实时区域对流层建模方法，包括：利用实时精密单点定位技术，实时估计GNSS参考站网内各参考站的天顶对流层延迟参数ZTD；选定数学函数模型、随机模型及动态模型的过程噪声；基于所述参考站网内各参考站的位置信息及实时ZTD信息，采用抗差卡尔曼滤波算法实时逐历元递归估计模型参数；对求解得到的模型参数进行编码，基于NTRIP协议通过互联网实时播发。不仅利用各GNSS参考站当前历元的对流层延迟信息，还可顾及区域对流层延迟的短时平稳变化特性；此外，引入抗差机制，有效避免了因个别GNSS参考站实时对流层延迟解算异常而引起的区域对流层延迟拟合的整体偏差，提升了建模的精度及稳健性。

Description

一种基于抗差卡尔曼滤波算法的实时区域对流层建模方法

技术领域

本发明涉及GNSS区域对流层建模技术领域，尤其涉及一种基于抗差卡尔曼滤波算法的实时区域对流层建模方法，主要适用于提升区域对流层建模的精度和稳健性。

背景技术

利用区域GNSS参考站网数据实时解算天顶对流层延迟，并将其进行空间拟合，生成区域对流层延迟模型，可用于增强区域内的GNSS定位、导航、授时等应用。目前实时区域对流层建模方法一般是基于各GNSS参考站当前历元解算的对流层延迟信息，采用合适的空间拟合函数，进行单历元最小二乘估计(本发明中简称其为LSETROP)。现有方法主要不足表现为以下两个方面：①建模过程将区域对流层模型参数随时间的变化视为白噪声过程，未顾及区域对流层延迟的短时平稳变化特性，时间域信息利用不充分；②建模过程对GNSS实时观测数据、卫星轨道与钟差产品的稳定度具有强依赖性，实时数据流传输的不稳定，会严重影响实时建模的精度。

发明内容

本发明的目的是克服现有技术中存在的建模精度低、建模稳健性差的缺陷与问题，提供一种建模精度高、建模稳健性好的基于抗差卡尔曼滤波算法的实时区域对流层建模方法。

为实现以上目的，本发明的技术解决方案是：一种基于抗差卡尔曼滤波算法的实时区域对流层建模方法，该方法包括以下步骤：

A、利用实时精密单点定位技术，实时估计GNSS参考站网内各参考站的天顶对流层延迟参数ZTD；

实时精密单点定位技术解算ZTD的公式如下：

式(1)中，s和r分别代表卫星和接收机号，f代表频率号，

和

为经精密卫星轨道与钟差产品改正后的伪距和相位观测值，

为测站至卫星方向的单位向量，i为历元标识，Δx为待估测站坐标，dt_r为待估接收机钟差，

为待估电离层延迟，

为待估模糊度参数，ε_P和ε_Φ分别表示伪距与相位观测噪声，T_r ^s为站星方向对流层总延迟；

T_r ^s进一步参数化为天顶方向对流层干延迟ZHD和湿延迟ZWD与相应干投影函数M_h、湿投影函数M_w乘积的形式，如下：

式(2)中，ZHD由Saastamoinen模型计算得出，ZWD作为未知参数进行求解；

B、选定对流层建模的数学函数模型、随机模型及动态模型的过程噪声；

数学函数模型为：

ZTD＝a₀₀+a₀₁·y+a₁₀·x+a₁₁·xy+b₁·h (3)

式(3)中，x、y和h表示参考站的位置信息参数，a₀₀为区域中心位置处天顶对流层延迟，a₀₁和a₁₀分别为ZTD沿x和y方向变化率，a₁₁为ZTD沿xy综合变化率，b₁为ZTD高度递减率；

随机模型采用等权方案，设定建模所用实时ZTD观测量的标准差为10mm；

模型参数随时间变化均模型化为随机游走过程，其过程噪声设置如下表：

C、基于所述参考站网内各参考站的位置信息及实时ZTD信息，采用抗差卡尔曼滤波算法实时逐历元递归估计区域对流层模型参数a₀₀、a₀₁、a₁₀、a₁₁和b₁；

卡尔曼滤波的状态方程和观测方程表示为：

式(4)中，X_k表示k历元时刻的状态向量，对应待估模型参数a₀₀、a₀₁、a₁₀、a₁₁和b₁，其协方差阵记为

Y_k表示k历元时刻的观测向量，对应各GNSS参考站实时解算的ZTD信息，A_k为设计矩阵；X_k-1表示k-1历元时刻的状态向量，其协方差阵记为

Φ_k，k-1表示由k-1历元时刻到k历元时刻的状态转移矩阵；Γ_k-1表示过程噪声，其协方差阵记为

Δ_k为k历元时刻观测噪声，其协方差阵记为

滤波值

按下列步骤进行求解：

首先，一步预测：

式(5)中，X_k，k-1表示k历元时刻的状态向量预报值，其协方差阵记为

Φ^T _k，k-₁表示Φ_k，k-1的转置矩阵；

然后，计算滤波值：

式(6)中，

为滤波增益矩阵，其中，

表示A_k的转置；I表示单位矩阵；

为k历元时刻的滤波值，其协方差阵记为

借鉴拟稳估计思想，采用如下IGGIII权函数调整上述卡尔曼滤波估计中的增益矩阵K_k，实现抗差估计；

式(7)中，

表示为经IGGIII权函数调整过后K_k矩阵中第p行、q列对应元素；K_pq为K_k矩阵中第p行、q列对应元素；k₀、k₁为经验值，取值分别为1.5和2.5；r_q为抗差因子，基于卡尔曼滤波残差及其协方差确定；

最后，给定系统状态向量初值X₀及其协方差阵

根据不断更新的区域参考站实时ZTD观测信息，利用抗差卡尔曼滤波算法递归估计得到当前观测历元的最优滤波值

及其协方差阵

D、对求解得到的

进行编码，基于NTRIP协议通过互联网实时播发。

步骤C中，所述抗差因子r_q的确定方法如下：

V_k，k-1＝Y_k-A_kX_k，k-1 (9)

式(8)、式(9)中，

为k历元时刻的滤波残差；V_k，k-1为新息向量，其协方差阵为：

将式(10)代入滤波增益矩阵K_k计算公式，则有：

根据式(6)、式(7)和式(9)，卡尔曼滤波残差

计算公式如下：

根据误差传播定律，卡尔曼滤波残差的协方差阵为：

式(13)中，

为

的转置；

故抗差因子r_q为：

式(14)中，v_q为第q个观测量的滤波残差，其相应的协方差为

与现有技术相比，本发明的有益效果为：

本发明一种基于抗差卡尔曼滤波算法的实时区域对流层建模方法采用抗差卡尔曼滤波算法实现区域对流层模型参数的实时逐历元递归估计，使得建模过程不仅利用各GNSS参考站当前历元解算的对流层延迟信息，还可充分顾及区域对流层延迟的短时平稳变化特性，充分利用了时间域信息；此外，在建模过程引入抗差机制，有效避免了因个别GNSS参考站实时对流层延迟解算异常而引起的区域对流层延迟拟合的整体偏差，提升了建模的精度及稳健性。因此，本发明提升了区域对流层建模的精度及稳健性。

附图说明

图1是一种基于抗差卡尔曼滤波算法的实时区域对流层建模方法的流程图。

图2是本发明的实施例中所用的香港SatRef跟踪站网的测站分布图。

图3是本发明的实施例中的分别利用图1所示EKFTROP建模方法与传统LSETROP建模方法实时计算的连续72小时的模型参数时间序列对比图。

图4是本发明的实施例中的每个建模历元所用的有效参考站数量序列图。

图5是本发明的实施例中的分别利用图1所示EKFTROP建模方法和传统LSETROP建模方法计算的两个测试站处连续72小时的实时ZTD序列与IGS最终ZTD产品序列比较图。

图3中，1号线表示传统LSETROP建模方法实时计算的连续72小时的模型参数时间序列，2号线表示EKFTROP建模方法实时计算的连续72小时的模型参数时间序列。

图5中，3号线表示IGS最终ZTD产品序列，4号线表示传统LSETROP建模方法计算的两个测试站处连续72小时的实时ZTD序列，5号线表示EKFTROP建模方法计算的两个测试站处连续72小时的实时ZTD序列。

具体实施方式

以下结合附图说明和具体实施方式对本发明作进一步详细的说明。

参见图1，一种基于抗差卡尔曼滤波算法的实时区域对流层建模方法，该方法包括以下步骤：

A、利用实时精密单点定位技术，实时估计GNSS参考站网内各参考站的天顶对流层延迟参数ZTD；

实时精密单点定位技术解算ZTD的公式如下：

式(1)中，s和r分别代表卫星和接收机号，f代表频率号，

和

为经精密卫星轨道与钟差产品改正后的伪距和相位观测值，

为测站至卫星方向的单位向量，i为历元标识，Δx为待估测站坐标，dt_r为待估接收机钟差，

为待估电离层延迟，

为待估模糊度参数，ε_P和ε_Φ分别表示伪距与相位观测噪声，T_r ^s

为站星方向对流层总延迟；

T_r ^s进一步参数化为天顶方向对流层干延迟ZHD和湿延迟ZWD与相应干投影函数M_h、湿投影函数M_w乘积的形式，如下：

式(2)中，ZHD由Saastamoinen模型计算得出，ZWD作为未知参数进行求解；

B、选定数学函数模型、随机模型及动态模型的过程噪声；

数学函数模型为：

ZTD＝a₀₀+a₀₁·y+a₁₀·x+a₁₁·xy+b₁·h (3)

式(3)中，x、y和h表示参考站的位置信息参数，a₀₀为区域中心位置处天顶对流层延迟，a₀₁和a₁₀分别为ZTD沿x和y方向变化率，a₁₁为ZTD沿xy综合变化率，b₁为ZTD高度递减率；

随机模型采用等权方案，设定建模所用ZTD观测量的标准差为10mm；

模型参数随时间变化均模型化为随机游走过程，其过程噪声设置如下表：

C、基于所述参考站网内各参考站的位置信息及实时ZTD信息，采用抗差卡尔曼滤波算法实时逐历元递归估计区域对流层模型参数a₀₀、a₀₁、a₁₀、a₁₁和b₁；

卡尔曼滤波的状态方程和观测方程表示为：

式(4)中，X_k表示k历元时刻的状态向量，对应待估模型参数a₀₀、a₀₁、a₁₀、a₁₁和b₁，其协方差阵记为

Y_k表示k历元时刻的观测向量，对应各GNSS参考站实时解算的ZTD信息，A_k为设计矩阵；X_k-1表示k-1历元时刻的状态向量，其协方差阵记为

Φ_k，k-1表示由k-1历元时刻到k历元时刻的状态转移矩阵；Γ_k-1表示过程噪声，其协方差阵记为

Δ_k为k历元时刻观测噪声，其协方差阵记为

滤波值

按下列步骤进行求解：

首先，一步预测：

式(5)中，X_k，k-1表示k历元时刻的状态向量预报值，其协方差阵记为

Φ^T _k，k-₁表示Φ_k，k-1的转置矩阵；

然后，计算滤波值：

式(6)中，

为滤波增益矩阵，其中，

表示A_k的转置；I表示单位矩阵；

为k历元时刻的滤波值，其协方差阵记为

借鉴拟稳估计思想，采用如下IGGIII权函数调整上述卡尔曼滤波估计中的增益矩阵K_k，实现抗差估计；

式(7)中，

表示为经IGGIII权函数调整过后K_k矩阵中第p行、q列对应元素；K_pq为K_k矩阵中第p行、q列对应元素；k₀、k₁为经验值，取值分别为1.5和2.5；r_q为抗差因子，基于卡尔曼滤波残差及其协方差确定；

最后，给定系统状态向量初值X₀及其协方差阵

根据不断更新的区域参考站ZTD观测信息，利用抗差卡尔曼滤波算法递归估计得到当前观测历元的最优滤波值

及其协方差阵

D、对求解得到的

进行编码，基于NTRIP协议通过互联网实时播发。

步骤C中，所述抗差因子r_q的确定方法如下：

V_k，k-1＝Y_k-A_kX_k，k-1 (9)

式(8)、式(9)中，

为k历元时刻的滤波残差；V_k，k-1为新息向量，其协方差阵为：

将式(10)代入滤波增益矩阵K_k计算公式，则有：

根据式(6)、式(7)和式(9)，卡尔曼滤波残差

计算公式如下：

根据误差传播定律，卡尔曼滤波残差的协方差阵为：

式(13)中，

为

的转置；

故抗差因子r_q为：

式(14)中，v_q为第q个观测量的滤波残差，其相应的协方差为

本发明的原理说明如下：

ZHD通过Saastamoinen模型计算得出，其输入参数为标准气象参数。

等权方案是指各观测量(本设计中具体指GNSS参考站网内各参考站解算的ZTD)在平差过程中赋予的权重相同。

拟稳估计思想指的是测量数据处理中相关平差理论的公知。

卡尔曼滤波基于观测信息序列及动力学模型信息，通过不断预测与修正过程，实现函数模型中未知参数X的递归估计。

实施例：

参见图1，一种基于抗差卡尔曼滤波算法的实时区域对流层建模方法，该方法包括以下步骤：

A、利用实时精密单点定位技术(PPP)，实时估计GNSS参考站网内各参考站的天顶对流层延迟参数ZTD；

实时精密单点定位技术解算ZTD的公式如下：

式(1)中，s和r分别代表卫星和接收机号，f代表频率号，

和

为经精密卫星轨道与钟差产品改正后的伪距和相位观测值，

为测站至卫星方向的单位向量，i为历元标识，Δx为待估测站坐标，dt_r为待估接收机钟差，

为待估电离层延迟，

为待估模糊度参数，ε_P和ε_Φ分别表示伪距与相位观测噪声，T_r ^s为站星方向对流层总延迟；

T_r ^s进一步参数化为天顶方向对流层干延迟ZHD和湿延迟ZWD与相应干投影函数M_h、湿投影函数M_w乘积的形式，如下：

式(2)中，由于ZHD较稳定，ZHD由Saastamoinen模型计算得出，ZWD作为未知参数进行求解；

利用式(1)和式(2)可实现ZTD参数的精确解算，ZTD参数解算过程中所用数据及主要参数处理策略如表1所示，数据时间为2018年年积日第195天至197天；

表1.ZTD参数估计过程中所用数据及主要参数处理策略

B、选定对流层建模的数学函数模型、随机模型及动态模型的过程噪声，为后续抗差卡尔曼滤波估计做准备；

数学函数模型为：

ZTD＝a₀₀+a₀₁·y+a₁₀·x+a₁₁·xy+b₁·h (3)

式(3)中，x、y和h表示参考站的位置信息参数，a₀₀为区域中心位置处天顶对流层延迟，a₀₁和a₁₀分别为ZTD沿x和y方向变化率，a₁₁为ZTD沿xy综合变化率，b₁为ZTD高度递减率；可以看出该公式中待估的模型参数共有5个，即a₀₀、a₀₁、a₁₀、a₁₁和b₁，因此，应用时要求区域参考站网内有效GNSS参考站数量最少为5个；

随机模型采用等权方案，设定建模所用ZTD观测量的标准差为1Omm；

模型参数随时间变化均模型化为随机游走过程，其过程噪声设置如表2：

表2.模型参数过程噪声

C、基于所述参考站网内各参考站的位置信息及实时ZTD信息，采用抗差卡尔曼滤波算法实时逐历元递归估计区域对流层模型参数a₀₀、a₀₁、a₁₀、a₁₁和b₁；

与单历元最小二乘估计不同，卡尔曼滤波算法可实现区域对流层模型参数的实时逐历元递归估计，卡尔曼滤波的状态方程和观测方程表示为：

式(4)中，X_k表示k历元时刻的状态向量，对应待估模型参数a₀₀、a₀₁、a₁₀、a₁₁和b₁，其协方差阵记为

Y_k表示k历元时刻的观测向量，对应各GNSS参考站实时解算的ZTD信息，A_k为设计矩阵；X_k-1表示k-1历元时刻的状态向量，其协方差阵记为

Φ_k，k-1表示由k-1历元时刻到k历元时刻的状态转移矩阵；Γ_k-1表示过程噪声，其协方差阵记为

Δ_k为k历元时刻观测噪声，其协方差阵记为

滤波值

按下列步骤进行求解：

首先，一步预测：

式(5)中，X_k，k-1表示k历元时刻的状态向量预报值，其协方差阵记为

Φ^T _k，k-1表示Φ_k，k-1的转置矩阵；

然后，计算滤波值：

式(6)中，

为滤波增益矩阵，其中，

表示A_k的转置；I表示单位矩阵；

为k历元时刻的滤波值，其协方差阵记为

借鉴拟稳估计思想，采用如下IGGIII权函数调整上述卡尔曼滤波估计中的增益矩阵K_k，实现抗差估计；

式(7)中，

表示为经IGGIII权函数调整过后K_k矩阵中第p行、q列对应元素；K_pq为Kk矩阵中第p行、q列对应元素；k₀、k₁为经验值，取值分别为1.5和2.5；r_q为抗差因子，基于卡尔曼滤波残差及其协方差确定；

最后，给定系统状态向量初值X₀及其协方差阵

根据不断更新的区域参考站ZTD观测信息，利用抗差卡尔曼滤波算法递归估计得到当前观测历元的最优滤波值

及其协方差阵

所述抗差因子r_q的确定方法如下：

V_k，k-1＝Y_k-A_kX_k，k-1 (9)

式(8)、式(9)中，

为k历元时刻的滤波残差；V_k，k-1为新息向量，其协方差阵为：

将式(10)代入滤波增益矩阵K_k计算公式，则有：

根据式(6)、式(7)和式(9)，卡尔曼滤波残差

计算公式如下：

根据误差传播定律，卡尔曼滤波残差的协方差阵为：

式(13)中，

为

的转置；

故抗差因子r_q为：

式(14)中，v_q为第q个观测量的滤波残差，其相应的协方差为

结合式(12)和式(14)可以看出，当第q个观测量存在粗差时，则其滤波残差v_q变大，对应的抗差因子r_q增大，K_pq则相应减小，从而减弱了观测粗差对滤波估值的影响，实现了抗差估计；

D、对求解得到的

(对应对流层模型参数)进行编码，基于NTRIP协议通过互联网实时播发。

本实施例选取中国香港SatRef GNSS观测网络部分站的实时GPS观测数据，测站分布如图2所示(图片源于http：//www.geodetic.gov.hk/en/satref/satref.htm)，测站总数为15个(HKTK、HKFN T430、HKKT、HKLT、HKCL、HKNP、HKMW、HKPC、HKSC、HKST、HKSS、HKKS、HKQT、HKOH、HKLM)，其中参与建模的测站数量为13个(HKTK、HKFN T430、HKKT、HKLT、HKCL、HKNP、HKMW、HKPC、HKSC、HKST、HKSS、HKOH、HKLM)，另外2个选作测试站(HKSL和HKWS，均为IGS站)，测区面积约为60km×30km，测站间最大高差约为345m。

图3显示了在上述测区分别利用图1所示对流层建模方法(EKFTROP)和传统LSETROP建模方法实时计算的连续72小时的模型参数时间序列对比图。从图中可明显看出，传统的LSETROP模型系数具有较大噪声，这主要是由于LSETROP建模方法仅利用当前历元的实时ZTD信息进行空间建模，而实时GNSS轨道与钟差产品估计的实时ZTD通常具有较大噪声，加之，实时数据流传输的不稳定性，更增加了LSETROP方法实时建模结果的不确定性。

图4给出了整个测试段内每个建模历元所用的有效参考站数量序列图。从图中可看出，在第8个小时附近，建模可用有效测站数由13个急剧减少为5个。由于参与建模的参考站数量及其位置分布在此历元的急剧变化，因此，造成了图3中LSETROP模型系数(b₁：ZTD的高度递减)在此历元出现了不合理的急剧减小。而图1所示对流层建模方法不仅利用各参考站当前历元的ZTD信息，还充分顾及了对流层延迟的短时平稳变化特性，因此，图1所示对流层建模方法系数随时间变化较为平稳，还可有效避免个别参考站实时数据传输的不稳定对建模过程的不利影响，使得实时区域对流层建模更加稳健。

为进一步比较分析图1所示EKFTROP建模方法与传统LSETROP建模方法的应用效果，本设计选取香港地区两个未参与建模的IGS站(HKSL和HKWS)作为测试站，并以IGS提供的最终ZTD产品(IGS-ZTD)为参考，评估分析了应用上述两种方法建立的实时区域对流层模型的应用精度。图5显示了分别利用图1所示EKFTROP建模方法和传统LSETROP建模方法计算的两个测试站处连续72小时的实时ZTD序列(时间分辨率5min)与IGS-ZTD序列比较图。从图中可明显看出，利用传统LSETROP建模方法计算的ZTD序列具有较大噪声，且个别历元存在不合理的短时急剧变化，主要原因同上述的分析内容。而图1所示EKFTROP建模方法充分顾及了对流层延迟的短时平稳变化特性，在时间域对建模过程增加了合理的约束，因此，利用图1所示EKFTROP建模方法计算的ZTD序列随时间变化较为平稳，有效避免了不合理的短时急剧变化。此外，以IGS-ZTD产品为参考，统计上述两种建模方法的应用精度，结果见表3。

表3.EKFTROP与LSETROP方法在HKWS和HWSL两个测试站处的应用精度统计结果(mm)

结合图5和表3可以看出，在模型稳健性及精度方面，图1所示EKFTROP建模方法(5号线)显示出了一定的优越性，最大瞬时精度提升约达40mm。

综上，图1所示EKFTROP建模方法相对于现有方法可有效提升现有实时区域对流层建模的精度及稳定度。

Claims (1)

1.一种基于抗差卡尔曼滤波算法的实时区域对流层建模方法，其特征在于，该方法包括以下步骤：

A、利用实时精密单点定位技术，实时估计GNSS参考站网内各参考站的天顶对流层延迟参数ZTD；

实时精密单点定位技术解算ZTD的公式如下：

式(1)中，s和r分别代表卫星和接收机号，f代表频率号，

和

为经精密卫星轨道与钟差产品改正后的伪距和相位观测值，

为测站至卫星方向的单位向量，i为历元标识，Δx为待估测站坐标，dt_r为待估接收机钟差，

为待估电离层延迟，

为待估模糊度参数，ε_P和ε_Φ分别表示伪距与相位观测噪声，T_r ^s为站星方向对流层总延迟；

T_r ^s进一步参数化为天顶方向对流层干延迟ZHD和湿延迟ZWD与相应干投影函数M_h、湿投影函数M_w乘积的形式，如下：

式(2)中，ZHD由Saastamoinen模型计算得出，ZWD作为未知参数进行求解；

B、选定对流层建模的数学函数模型、随机模型及动态模型的过程噪声；

数学函数模型为：

ZTD＝a₀₀+a₀₁·y+a₁₀·x+a₁₁·xy+b₁·h (3)

式(3)中，x、y和h表示参考站的位置信息参数，a₀₀为区域中心位置处天顶对流层延迟，a₀₁和a₁₀分别为ZTD沿x和y方向变化率，a₁₁为ZTD沿xy综合变化率，b₁ 为ZTD高度递减率；

随机模型采用等权方案，设定建模所用实时ZTD观测量的标准差为10mm；

模型参数随时间变化均模型化为随机游走过程，其过程噪声设置如下表：

C、基于所述参考站网内各参考站的位置信息及实时ZTD信息，采用抗差卡尔曼滤波算法实时逐历元递归估计区域对流层模型参数a₀₀、a₀₁、a₁₀、a₁₁和b₁；

卡尔曼滤波的状态方程和观测方程表示为：

式(4)中，X_k表示k历元时刻的状态向量，对应待估模型参数a₀₀、a₀₁、a₁₀、a₁₁和b₁，其协方差阵记为

Y_k表示k历元时刻的观测向量，对应各GNSS参考站实时解算的ZTD信息，A_k为设计矩阵；X_k-1表示k-1历元时刻的状态向量，其协方差阵记为

Φ_k，k-1表示由k-1历元时刻到k历元时刻的状态转移矩阵；Γ_k-1表示过程噪声，其协方差阵记为

Δ_k为k历元时刻观测噪声，其协方差阵记为

滤波值

按下列步骤进行求解：

首先，一步预测：

式(5)中，X_k，k-1表示k历元时刻的状态向量预报值，其协方差阵记为

Φ^T _k，k-1表示Φ_k，k-1的转置矩阵；

然后，计算滤波值：

式(6)中，

为滤波增益矩阵，其中，

表示A_k的转置；I表示单位矩阵；

为k历元时刻的滤波值，其协方差阵记为

借鉴拟稳估计思想，采用如下IGGIII权函数调整上述卡尔曼滤波估计中的增益矩阵K_k，实现抗差估计；

式(7)中，

表示为经IGGIII权函数调整过后K_k矩阵中第p行、q列对应元素；K_pq为K_k矩阵中第p行、q列对应元素；k₀、k₁为经验值，取值分别为1.5和2.5；r_q为抗差因子，基于卡尔曼滤波残差及其协方差确定；

所述抗差因子r_q的确定方法如下：

V_k，k-1＝Y_k-A_kX_k，k-1 (9)

式(8)、式(9)中，

为k历元时刻的滤波残差；V_k，k-1为新息向量，其协方差阵为：

将式(10)代入滤波增益矩阵K_k计算公式，则有：

根据式(6)、式(7)和式(9)，卡尔曼滤波残差

计算公式如下：

根据误差传播定律，卡尔曼滤波残差的协方差阵为：

式(13)中，

为

的转置；

故抗差因子r_q为：

式(14)中，v_q为第q个观测量的滤波残差，其相应的协方差为

最后，给定系统状态向量初值X₀及其协方差阵

根据不断更新的区域参考站实时ZTD观测信息，利用抗差卡尔曼滤波算法递归估计得到当前观测历元的最优滤波值

及其协方差阵

D、对求解得到的

进行编码，基于NTRIP协议通过互联网实时播发。

Images