论文推荐 | 陈正生 : 电离层延迟变化自模型化的载波相位平滑伪距算法

《测绘学 》

构建与学术的桥梁 拉近与权威的距离

电离层延迟变化自模型化的载波相位平滑伪距算法

陈正生^1,2, 张清华^2,3

1. 火箭军工程大学, 陕西 西安 710025;

2. 地理信息工程国家重点实验室, 陕西 西安 710054;

3. 陆军工程大学, 江苏 南京 210007;

4. 信息工程大学, 河南 郑州 450052

收稿日期：2018-09-03；修回日期：2019-02-16

基金项目：国家重点研发计划(2016YFB0501701);国家自然科学基金(41604024;41674019;41804006);地理信息工程国家重点实验室开放研究基金(SKLGIE2016-M-2-3);陕西省自然科学基金(2018JQ4023)

关键词：电离层延迟变化 Hatch滤波 自模型化 载波相位平滑伪距 单频

An improved carrier phase smoothing pseudorange algorithm with self-modeling of ionospheric delay variation

CHEN Zhengsheng^1,2, ZHANG Qinghua^2,3, LI Linyang⁴, LI Xuerui^1,2, Lü Hao⁴

1. Rocket Force University of Engineering, Xi‘an 710025, China;

2. State Key Laboratory of Geo-Information Engineering, Xi‘an 710054, China;

3. Army Engineering University, Nanjing 210007, China;

4. Information Engineering University, Zhengzhou 450052, China

Foundation support: The National Key Research and Development Program of China (No. 2016YFB0501701); The National Natural Science Foundation of China (Nos. 41604024; 41674019; 41804006); The Open Research Fund of the State Key Laboratory of Geo-Information Engineering (No. SKLGIE2016-M-2-3); The National Natural Science Foundation of Shanxi Province of China (No. 2018JQ4023)

First author: CHEN Zhengsheng(1984—), male, PhD, lecturer, majors in GNSS precise data processing. E-mail:czsgeo@126.com.

Corresponding author: ZHANG Qinghua E-mail:sqmha@126.com.

Abstract: The traditional single-frequency carrier phase smoothing pseudo-range algorithm is prone to divergence and precision degradation due to the influence of ionospheric delay. However, the existing solutions have limited accuracy improvement or need external data support. In this paper, the regularity of ionospheric variation is studied and a regressive model is established. On this basis, a self-modelling algorithm for single-frequency carrier phase smoothing pseudorange with ionospheric delay variation is proposed. The algorithm uses the ionospheric delay information contained in the pseudo-range and carrier observations to model the ionospheric delay, and deducts the ionospheric delay variation between epochs from the smooth pseudo-range, thus effectively avoiding the divergence of the smooth pseudo-range. Carrier smoothing pseudo-range algorithm of ionospheric self-modeling is realized by using GNSSer software. Static and dynamic observation data are used to carry out positioning experiment and precision analysis. The example shows that:① the long period regular Hatch filtering is seriously affected by the ionosphere; ② the precision of the self-modelling ionospheric delay can reach centimeter level, and the linear moving window fitting method is the best in the 30-minute window; ③ self-modeling ionospheric correction can effectively eliminate the influence of smooth pseudo-range ionosphere. With the increase of time window, the precision does not decrease; ④ the proposed algorithm is used for epoch-by-epoch single-frequency smoothing pseudo-range single-point positioning, and the positioning accuracy reaches sub-decimeter level in both static and dynamic NEU directions. In the dynamic positioning test, the horizontal and elevation direction accuracy is 6.25 cm and 10.4 cm, which are 5.4 times and 3.3 times higher than the original pseudo-range respectively.

Key words: ionospheric delay variation Hatch filter self-modeling carrier phase smoothing pseudorange single frequency

GNSS载波相位观测量具有毫米级的精度，远高于伪距观测量，由文献[1]提出的Hatch滤波算法组合了载波相位和伪距两种观测量，以历元间电离层延迟变化较小为前提，通过多历元的叠加求出平均模糊度和电离层延迟量，并用于改进伪距观测值的精度。这种方法在短时段内可以显著提高伪距观测值精度，但是当历元间电离层变化较大或平滑时段较长时，平滑结果将产生严重的系统误差。基于此，国内外一些学者对算法进行改进。文献[2]利用差分基站提供的电离层物理信息，考虑到电离层在不同区域，日、季节的电离层变化和卫星仰角等因素，建立电离层模型并提出一种最优平滑伪距的DGPS单频接收机算法，比传统的Hatch滤波具有更高的定位精度和稳健性。文献[3]利用线性回归模型化电离层延迟的变化率，并综合考虑卫星高度截止角和多路径噪声的特性，提出一种自适应滤波算法，通过成本函数计算最优平滑时间，能够在一定程度上提高GPS静态定位和动态定位的精度。文献[4]利用Klobuchar模型补偿电离层延迟变化的变化率，提出了一种改进的单频载波相位平滑伪距方法，仿真算例表明其精度明显优于传统单频载波相位平滑伪距法。文献[5]利用Kalman滤波的新息向量，调整了噪声方差矩阵并得到最优平滑时间窗口的时间常数。此算法在仿真算例中有效改善了定位精度和稳定性。文献[6-7]指出在动态飞行器上搭载的单频WAAS和LAAS接收机，其载波相位平滑观测量经常受到电离层变化的影响，尤其在赤道和极地区域，并提出了减轻TEC影响的算法。文献[8]针对地基增强系统中载波相位平滑存在的误差问题，采用由卫星仰角、电离层变化及用户到参考站的距离坐标变量实时确定平滑窗口的宽度，其结果比传统的100 s Hatch滤波器的误差减少了68%。文献[9]为了应对传统Hatch滤波器的发散问题，提出利用电离层延迟变化量和卫星高度角自适应改变平滑窗口的宽度，其结果比经典的Hatch滤波器提高了约21%。为了应对电离层延迟的影响以及增强结果的稳定性，也有不少学者采用多频或差分数据来避免常规单频载波平滑伪距所带来的发散问题^[10-13]，但是单频或非差分用户无法使用。

由于电离层延迟变化是影响单频载波相位平滑伪距精度变差的主要原因，因此解决平滑伪距精度发散的关键在于准确计算并扣除历元间的电离层延迟变化值。目前存在的研究结果，部分是利用窗口减弱电离层影响，部分是利用外部的差分信息对电离层延迟进行改正，也有研究根据外部的物理和几何参数对电离层改正窗口的宽度进行估计，但大多需要依赖外部环境或数据，并且改正效果有限。

1 顾及电离层延迟变化的载波相位平滑伪距

忽略对伪距码和载波相位观测量影响基本相同的对流层延迟、接收机钟差、卫星钟差，以及在一定卫星高度截止角或多历元后均值后，可忽略不计的多路径效应等，k历元时某一频点的相位和伪距观测方程可以简化为

(1)

式中，Lk和Pk分别表示载波相位距离观测值和伪距观测值；ρk为站星几何距离；N为以距离为单位的模糊度参数；Ik为该频率的电离层距离延迟；εL、εP为载波相位距离和伪距观测的观测噪声。

令第1历元的电离层延迟为I₁，第k历元与第1历元电离层延迟差为ΔIk, 1=Ik－I₁，将式(1)中两方程相减，并将观测噪声吸收在观测量中，则

为方便表达，引入常数A=2I₁+N。在模糊度和电离层延迟未知的情况下，A需要由观测值来计算。考虑到各历元观测值存在噪声，在k历元内，可以求得A的估值Ak

(3)

式中，i为历元序号。由式(2)和式(3)可得到指定频点第k历元顾及电离层延迟变化的平滑伪距

相邻两历元相减，得到顾及电离层延迟变化的载波相位平滑伪距的递推公式

式中

(6)

将其代入式(5)，有

(7)

式(7)与式(5)等价，但是减少了Ak的迭代计算，具有更高的计算效率。如果整个时段中电离层变化为0，即对于任意历元i有ΔIi, 1=0，那么式(7)就变成了经典的Hatch滤波，也是众多文献中所指的载波平滑伪距加权平均公式^{[1, 4, 14]}，即经典Hatch滤波是顾及电离层延迟变化的载波平滑伪距的特殊形式。可以发现式(7)右边最后1项右括号内的内容与式(5)右边前3项是相同的，只包括前k-1历元的平滑值和当前历元载波变化改正值，称之为平滑伪距Pk的推估值Pk，即

式中，顾及电离层延迟变化的载波相位平滑伪距推估模型。若历元间电离层延迟变化为0，则是经典Hatch滤波的推估模型。式(7)可以写为

由式(3)和式(4)可知，历元间电离层延迟变化对于传统Hatch滤波在第k历元的平滑伪距，造成的偏差为

(10)

式中，

由于常规的载波平滑伪距算法忽略了电离层延迟历元间变化的影响，当平滑时段变长，电离层延迟变化增加到一定程度后，式(10)的积累值变大，平滑伪距存在着发散的情况。

2 电离层延迟的演化规律与自模型化求解2.1 电离层延迟规律分析

2.1.1 采用双频载波计算电离层变化

对于L1和L2双频GNSS载波，考虑一阶电离层影响，其具有如下观测方程

(11)

式中，A为与传播路径电子总数相关的电离层参数A=－40.28c_TEC^[19]；T为对流层延迟；br、b_s为对应频点的接收机和卫星的未校准相位偏差，其他符号与式(1)相同。需要注意的是：这里的下标指的是对应的频点，而非历元编号。

令ΔL=L₁－L₂、Δb=bL₁－bL₂、ΔN=N₁－N₂、ΔB=(Δbr－Δbs)c+ΔN、ΔIA=A(1/f₁²－1/f₂²)，将式(11)中两个频率的载波观测值相减，可得

即频率间电离层延迟差分值可以通过载波距离观测值作差获得。由于短期内的硬件延迟通常变化很小^[21-23]，在不发生周跳的情况下，可以将式(12)中ΔB视为常量，与ΔI一起考虑。

将式(12)两边同时乘以ΔIA的系数，可得各频率电离层和相位偏差及模糊度的距离总和，以变量IB表示，这里只考虑L1载波的情况

(13)

历元作差，可以求得电离层历元间的变化值

(14)

通过两个载波测量值作差，可以消除大部分观测误差，只剩下电离层延迟、未校准相位偏差和模糊度。由误差传播定律可知，求出的电离层和相位偏差总和IB与载波相位观测量L属于同一数量级的测量精度，可达毫米级别。其中，如果不发生周跳，变化的只有电离层延迟。

2.1.2 电离层延迟变化计算

试验选取4个IGS观测站，其中两个赤道附近站(KIRI、NKLG)，一个高纬度站(LAMA)和一个中纬度站(SHAO)，采用全天观测数据，采样率为30 s，高度截止角为15°，计算GPS L1频率的电离层延迟变化率，其中NKLG的电离层延迟及其变化如图 1所示。

图 1 NKLG站全天电离层延迟及其变化率Fig. 1 All-day ionospheric delay and its rate of change at the NKLG station

图选项

4个站点电离层延迟变化的统计数据见表 1。

表 1 4个全球IGS站全天数据计算的电离层延迟率Tab. 1 Ionospheric delay rate calculated from all-day data of 4 global IGS stations

mm/s
站名	最大值	最小值	平均值(×10-8)	中误差	总历元数
KIRI	2.90	–3.04	3.603	0.423	25 422
NKLG	3.46	-2.66	1.760	0.401	25 121
LAMA	1.78	-1.99	1.519	0.265	22 515
SHAO	2.19	-2.06	-0.854	0.281	22 630

表选项

图 1(a)采用CODE球谐函数模型绘制^[24]，这里用于查看总体趋势。可以看出对于指定的站点，站星电离层具有明显的时间特性，在本地时间12—16时达到最大值，但是由于卫星的相对运动，各颗卫星信号的电离层延迟看上去连续，变化平缓，但是没有明显的运动规律；从实测数据的图 1(b)和从表 1可以看出电离层延迟变化量基本都在3 mm/s内与载波观测量噪声相当，且变化基本连续，可以认为各历元间的电离层延迟是连续变化的，因此在一定的时段内可以采用最小二乘拟合的方法求取其近似值。

2.2 电离层延迟变化的自模型化求解

(15)

在没有周跳的情况下，INk数值的变化直接体现为电离层延迟的变化。从式(15)中可以看出，采用伪距和载波作差获得的电离层延迟的精度比伪距精度高一倍。采用多个历元，通过最小二乘多项式拟合^[25]，可以精化其精度，从而求出站星电离层延迟变化的估值。在实时计算时，采用最近的时段窗口数据拟合；在事后或准实时处理时，拟合时段窗口中心历元。具体方法为：将各历元INk作为观测值，以时间作为自变量，连续多历元观测后，就可以组成误差方程进行平差计算，从而求得拟合系数，进而求得每个历元的拟合值及其变化。具体的拟合阶次和采用时间窗口大小，可以根据拟合残差大小进行判定。

2.3 电离层延迟自模型化效果分析

为了验证电离层延迟自模型化的效果，确定拟合阶次和窗口大小，采用两个IGS参考站：赤道站MBAR和中纬度站SHAO，在2018年001日，分别采用电离层延迟最活跃的赤道正午前后和最不活跃的中纬度夜间时段，采用不同拟合阶次(1~3)，采样率分别为10 s和1 s的实测卫星观测数据，进行逐历元移动开窗最小二乘多项式事后拟合计算，将双频载波计算数值作为真值，统计并比较不同参数下计算结果的中误差。

2.3.1 赤道MBAR早上到午后G23

MBAR站G23卫星，有效观测时段从2018-01-01 05:31到2018-01-01 14:27，跨度达9 h，由于地处赤道且包含正午时段，该段数据受电离层延迟影响大且变化剧烈。图 2是G23卫星电离层延迟偏差的拟合情况，其中黄色载波计算值作为真值，蓝色是伪距计算值，红色曲线是依据伪距值拟合的结果，深绿色是卫星高度角对应的图(a)的次竖轴。1 s采样率和10 s采样率在不同拟合阶次和不同时段窗口下的拟合统计结果见表 2，其中加粗字体表示当前数值为最小值。

图 2 G23 30 min窗口线性拟合电离层延迟Fig. 2 Ionospheric delay of G23 with a linear fit of 10 s sampling rate using a 30-minute window

图选项

表 2 MBAR G23星电离层拟合精度比较Tab. 2 Comparison of ionospheric fitting of MBAR G23

时段窗口/m	1 s采样率(32 141历元)				10 s采样率(3214历元)
	历元数	线性	2次	3次	历元数	线性	2次	3次
10	600	0.071	0.072	0.208	600	0.059	0.070	0.163
20	1200	0.069	0.069	0.070	1200	0.050	0.056	0.062
30	1800	0.068	0.069	0.070	1800	0.045	0.057	0.059
40	2400	0.070	0.069	0.069	2400	0.044	0.044	0.054
60	3600	0.076	0.067	0.067	3600	0.047	0.041	0.046
90	5400	0.095	0.072	0.068	5400	0.066	0.047	0.045
120	7200	0.134	0.079	0.070	7200	0.085	0.043	0.039

表选项

2.3.2 中纬度SHAO夜间到早上G17

SHAO站的G17卫星的时段为:UTC 16:30—23:00，对应本地时间为夜间到早上，22:30—7:00，共7 h，由于地处中纬度在夜间时段，该段数据电离层延迟较小且变化缓慢。表 3显示了该星在1 s采样率和10 s采样率在不同拟合阶次和不同时段窗口下的拟合统计结果。

表 3 SHAO G17星电离层拟合精度比较Tab. 3 Comparison of ionospheric fitting of SHAO G17

时段窗口/m	1 s采样率(32 141历元)				10 s采样率(3214历元)
	历元数	线性	2次	3次	历元数	线性	2次	3次
10	600	0.059	0.070	0.163	600	0.076	0.101	9.313
20	1200	0.050	0.056	0.062	1200	0.058	0.071	0.081
30	1800	0.045	0.057	0.059	1800	0.050	0.060	0.068
40	2400	0.044	0.044	0.054	2400	0.046	0.055	0.061
60	3600	0.047	0.041	0.046	3600	0.039	0.047	0.056
90	5400	0.066	0.047	0.045	5400	0.045	0.038	0.043
120	7200	0.085	0.043	0.039	7200	0.055	0.038	0.039

表选项

通过以上试验，可以发现由于短期内站星电离层延迟变化缓慢，采用移动开窗最小二乘多项式对单频伪距和载波组合计算的电离层延迟进行拟合，都取得了厘米级别的精度，但是不同的时段和拟合阶次具有不同结果，概括如下：

(2) 由单频载波和伪距计算的电离层延迟，在多项式拟合后整体精度可达4~7 cm，精度高于伪距噪声。

(3) 在相同时段内，采用不同的采样率数据，拟合结果有所区别，但是并不明显，因而可以认为电离层拟合与时段长度相关，与采样率相关不大。

3 自模型化平滑伪距的实现与精度分析3.1 自模型化平滑伪距的实现

图 3 电离层自模型化的平滑伪距算法流程Fig. 3 Flow chart of the ionospheric self-modeling smooth pseudorange algorithm

图选项

3.2 电离层变化改正的平滑伪距残差分析

图 4 常规Hatch滤波加电离层改正前后的平滑伪距残差Fig. 4 Residual of smooth pseudorange before and after ionospheric correction of Hatch filter

图选项

表 4 不同算法在电离层改正前后平滑伪距残差中误差Tab. 4 Errors in smoothing pseudo-range residuals before and after ionospheric correction

m
名称/时段窗口/min	历元数/权	固定权Hatch	电离层改正(窗口20 min)
常规Hatch	—	4.977 8	0.400 5
120	720	3.082 4	0.400 3
60	360	1.896 5	0.400 1
40	240	1.390 5	0.399 7
30	180	1.111 9	0.399 5
20	120	0.820 8	0.398 8
10	60	0.541 5	0.396 7
5	30	0.432 8	0.391 8

表选项

由表 4可以看出，随着时间的积累，受电离层变化影响，常规Hatch滤波产生了较大的系统性偏差。采用加权限制平滑伪距权值后，偏差得到改善，但是在长时段中的改善效果有限。考虑到平滑伪距需要一定的历元数量才能达到较好的平滑效果，因此采用加权窗口的方法，在降低系统误差的同时也降低了平滑精度。采用自模型化数据对电离层变化进行改正后，无论是常规的Hatch滤波，还是开窗加权Hatch滤波，系统性偏差基本得到了消除。本试验中，残差稳定在0.4 m左右。由于此处采用原始伪距作为比较基准，其本身具有较大噪声，因此本数值并不能代表平滑精度。下面通过定位计算对算法进行验证。

表 5 模型设置与参数估计Tab. 5 Model setting and parameter estimation

类型	处理策略和方法
观测值	GPS L1伪距及其载波相位平滑伪距
卫星定权方法	高度角定权，截止角10°
星历与钟差	IGS发布的精密产品
电离层延迟	采用双频非差非组合计算结果[27]
对流层延迟	GPT2模型+VMF1映射函数[28]
估计参数	4参数(3个位置参数和一个钟差参数)
平差方法	逐历元独立参数平差

表选项

4.1 MBAR站静态定位试验

采用MBAR站中午11—12点、采样率为1 s、总计3600历元的数据，进行伪距逐历元4参数平差定位计算，分别采用原始伪距，常规Hatch滤波，加权Hatch滤波以及电离层自模型化改正后的平滑伪距进行定位计算，最后，以IGS发布日解的坐标为真值，将定位结果进行作差比较，每个历元的定位偏差见图 5、图 6和表 6。其中, 符号De、Dn、Du分别是东北天(ENU)3方向的差值。

图 5 MBAR站电离层改正前后定位偏差Fig. 5 Positioning deviation before and after ionospheric correction of MBAR station

图选项

图 6 不同算法在ENU 3方向定位中误差Fig. 6 RMS of ENU in positioning of different algorithms

图选项

表 6 电离层改正前后的定位残差中误差Tab. 6 RMS of the positioning residual before and after ionospheric correction

m
算法	历元数/权	电离层改正前				电离层改正后
		E方向	N方向	U方向	水平	E方向	N方向	U方向	水平
原始伪距	–	0.156	0.201	0.481	0.254	–	–	–	–
常规Hatch	–	0.466	0.425	1.132	0.631	0.07	0.096	0.113	0.119
30 m加权	1800	0.241	0.262	0.402	0.356	0.069	0.097	0.103	0.119
20 m加权	1200	0.127	0.155	0.228	0.200	0.068	0.097	0.099	0.118
10 m加权	600	0.084	0.114	0.17	0.142	0.066	0.097	0.099	0.117
5 m加权	300	0.067	0.106	0.144	0.125	0.066	0.097	0.101	0.117
2 m加权	120	0.065	0.083	0.119	0.105	0.066	0.096	0.108	0.116
1 m加权	60	0.068	0.095	0.142	0.117	0.068	0.097	0.134	0.118

表选项

图 5为MBAR站不同伪距定位方法ENU 3方向的定位偏差。可以发现：①载波相位平滑伪距起到了很好的滤波作用，使得历元间的定位结果连续，且窗口越长，定位结果越平滑；②载波相位平滑伪距可以在很短的时间内(10 s)，达到分米级别的定位精度；③随着历元窗口的增加，常规Hatch滤波受电离层影响而发散(图 5(a))，高程方程最大偏移达3 m，而采用窗口加权后，系统偏差得到了抑制，但是历元间的定位结果波动变大；④采用自模型化电离层改正后的平滑伪距，无论是常规Hatch滤波，还是加权窗口滤波的电离层系统偏差都得到了很好的消除或减弱，并且定位结果基本稳定。电离层改正前后的定位残差中误差