地球物理学报  2016, Vol. 59 Issue (1): 101-115   PDF    
GPS接收机仪器偏差的短期时变特征提取与建模
张宝成1,2, 袁运斌2, 欧吉坤2    
1. GNSS Research Centre, Department of Spatial Sciences, Curtin University, Perth 6845, Australia;
2. 中国科学院测量与地球物理研究所动力大地测量学国家重点实验室, 武汉 430077
摘要: 卫星和接收机仪器偏差(Differential Code Biases, DCB)是利用GPS(Global Positioning System)研究电离层的两类主要误差源.由于所处的空间环境恒定,且可被全球跟踪站连续观测,GPS卫星的DCB具备长期稳定性和较高的估计精度.但针对各类型接收机而言,受测站环境、硬件设施等影响,其DCB可能会呈现明显的短期变化.精确地模型化接收机DCB的短期变化特征,将有助于提高GPS电离层产品的可靠性,以及基于这些产品反演空间和地球科学现象的准确性.采用零/短基线GPS数据,本文改进了提取和分析接收机DCB变化的现有方案.随后,本文推导了一种能直接估计接收机DCB的函数模型.当检验出接收机DCB的短期变化服从随机游走时,通过对比接收机DCB的直接估值与间接提取值之间的符合程度,可"试探出"过程噪声标准差的最优经验值.实验分析选用4台双频接收机(共形成1条零基线和2条短基线,间距最大为15 m)多天的观测数据,主要结论包括:1)改进的接收机DCB提取方案能较好地克服低频伪距噪声和多路径效应的影响; 2)针对零基线,其接收机DCB在各天内的变化量级小于1 TECu,变化趋势则可采用过程噪声标准差为1.0~1.5 mm的随机游走加以描述; 3)对应于某短基线的接收机DCB在某天内的变化可达12 TECu,当采用随机游走描述其趋势时,过程噪声标准差的经验值超过2 mm.
关键词: GPS     精密单点定位     电离层     接收机仪器偏差     零/短基线    
Short-term temporal variability of GPS receiver's differential code biases(DCB):retrieving and modeling
ZHANG Bao-Cheng1,2, YUAN Yun-Bin2, OU Ji-Kun2    
1. GNSS Research Centre, Department of Spatial Sciences, Curtin University, Perth 6845, Australia;
2. State Key Laboratory of Dynamic Geodesy, Institute of Geodesy and Geophysics, Chinese Academy of Sciences, Wuhan 430077, China
Abstract: The satellite and receiver differential code biases(DCB) combined, account for the main error budget of GPS-based ionosphere investigations. As the space environment onboard the GPS satellites is quite constant, the long-term stability of GPS satellite DCB has been observed. At the same time, continuous GPS data collection from receivers of global coverage makes it possible to estimate GPS satellite DCB with high accuracy. These two facts, however, do not hold true for a variety of receivers' DCB. As a result of various operating environments as well as distinct firmware versions, receiver DCB may experience short-term variations over time. Precise modeling of receiver DCB's variation can raise the reliability of ionosphere products determined from GPS data, as well as ensure the correctness of conclusions drawn based on these products when investigating atmosphere/space effects and geodetic phenomena.Given zero/short-baseline GPS data, the customary scheme used to retrieve receiver DCB is further modified as follows:1) Precise point positioning(PPP) has been implemented, respectively, using GPS data collected by each of the receiver that forms the baseline. The slant ionosphere delays biased by satellite and receiver DCB, and the undifferenced, float-valued ambiguities can be estimated, among other parameters. 2) Those undifferenced ambiguities are then combined so as to form an independent set of double-differenced ambiguities that are fixable. 3) After taking these fixed ambiguities into consideration, the slant ionosphere delays determined by means of PPP can be further refined. The between-receiver, single-differenced values of these delays are then used to retrieve a time series of receiver DCB, the time resolution of which is equal to that of GPS observations in use. In addition, the ionosphere-fixed model with estimable receiver DCB has been derived. By characterizing the dynamic model of receiver DCB as random walk, the consistency between both the estimated and the formerly retrieved receiver DCB forms a basis to identify an optimal empirical value of the STD of the process noise.Numerical tests make use of multiple days' dual-frequency GPS data collected by 4 co-located receivers that form a total of one zero-baseline and two short-baselines, with a maximum separation of 15 m. The main conclusions include:1) The modified scheme outperforms the customary one, as being able to determine a time series of receiver DCB less affected by low-frequency code noise and multipath effects; 2) The intra-day variations of receiver DCB determined from the zero-baseline is less than 1 TECu, without apparent day-to-day repeatability. A random walk with STDs of process noise between 1.0 and 1.5 mm is sufficient to characterize different days' variation behaviors; 3) The size of receiver DCB variation corresponding to one of the short-baselines can exceed 12 TECu(roughly 2 m) during one day. To model it with random walk, the empirical STD of the process noise should be set no less than 2 mm. The proposed methods in this paper may serve as new ways to routinely monitor and calibrate receiver DCB.
Key words: GPS     PPP     Ionosphere     Receiver Differential Code Bias(DCB)     Zero/Short-baseline    
1 引言

与GPS(Global Positioning System)导航、定位等应用中尽量消除电离层延迟相反(Teunissen et al.,2010; Yang et al.,20112014; Yao et al.,2013; 高星伟等,2002; 刘经南和叶世榕,2002; 张宝成等,2010a2010b),利用双频GPS数据,可实时、准确地计算一系列的电离层参数,如总电子含量(TEC),电子密度、闪烁、行扰等(Dyrud et al.,2008; Hernández-Pajares et al.,2011; 袁运斌和欧吉坤,2005; Petrie et al.,2011; Øvstedal,2002; Yuan and Ou,2001; Li W et al.,2012Li Z S et al.,2012).这些参数从不同的角度反映了电离层的时空变化、局部异常等特性,可作为反演极端性空间天气事件、探究灾害性地学现象成因的依据和参考(毛田等,2009; 姚宜斌等,2102; 余涛等,2009; Yue et al.,2011; 张小红等,2013a2013b; Zhang et al.,2012; Li et al.,2015; Shen et al.,2015; Zheng and Xu,2015).

从双频GPS数据中获取TEC,可归纳为“提取-建模-分离”等步骤(Brunini and Azpilicueta,20092010; Ciraolo et al.,2007; Yuan et al.,2015),但易受以下几种误差的影响:

1)平滑误差:采用相位平滑伪距技术,可提取各卫星的电离层斜延迟,其中包含TEC、接收机和卫星仪器偏差(Differential Code Biases,DCB)(Conte et al.,2011; Nayir et al.,2007).其核心在于如何准确地“对齐”无几何影响的伪距和相位观测值,该“对齐”过程一般基于两类假设,即a)(加权)平均可充分消除伪距/相位多路径效应和观测噪声;b)在卫星连续弧段期间,接收机和卫星DCB不随时间变化.实际中,当两类假设不成立时,将引起所谓的平滑误差,其具体量级与观测环境有关(Ciraolo et al.,2007).其中,伪距多路径的贡献约为±1.4至±5.3 TECu,而接收机DCB单天内变化的贡献超过±1.4 TECu,最大达±8.8 TECu(Brunini and Azpilicueta,2009).

2)投影误差:采用薄层模型描述电离层时,通常借由某投影函数,将接收机至卫星方向的TEC转换为vTEC.与TEC相比,vTEC的变化规律与地磁和太阳活动、穿刺点处的时间和位置等因素的关联性更强,更易于模型化(Zhang et al.,2014).一般地,投影函数的准确性与所选取的电离层薄层高度(H)有关.H的取值与电离层F2层实际高度符合的越好,相应的投影函数将越准确.但在实际中,H常选用一个定值,如350、400 km或450 km,而电离层F2层高度则随时间变化,两者之间的差异将引起投影误差.Komjathy等(2005)分析表明,在2003年10月30日电离层异常活动期间,投影误差最大可达10 m(约61.7 TECu),而针对28日电离层平静活动时期,投影误差则降为0.8 m(约4.9 TECu).

3)模型误差:在地磁和太阳活动异常时期,描述vTEC时空变化的数学模型,如广义三角级数函数(袁运斌和欧吉坤,2005),双线性展开式(Brunini and Azpilicueta,2010)等,难以准确地刻画(中小尺度)行扰、闪烁、(垂向)梯度等电离层异常特性.另一方面,在联合估计vTEC和DCB等参数时,还通常假定:接收机和卫星DCB在一段时间,如若干小时(周东旭等,2011)、1—3天内(Brunini and Azpilicueta,20092010)等不存在明显的变化.但实际中,受测站环境显著变化(如温度,Zhang and Teunissen,2015)的影响,各类接收机DCB可能会在较短时间内发生变化(Coster et al.,2013; Kao et al.,2013; Zhang et al.,20102014).当采用单台GPS接收机实施电离层建模时,Brunini和Azpilicueta(2009)研究表明,在不同的电离层活动条件下,模型误差的范围约为2~15 TECu.

由以上论述可知,除投影误差外,另外两类误差的成因、大小均与接收机DCB是否稳定有密切的关系.因此,分析接收机DCB随时间的变化特性,揭示其短、中、长期变化的规律,发现影响其变化的因素,并探索适当的模型化方案,将有助于改善TEC计算结果的准确性.为实现该目的,现有研究采用了两种不同的技术方案:

1)从模型误差中间接分析接收机DCB的中、长期变化(周东旭等,2011; Zhang et al.,2014).其基本思路为:采用长期的双频GPS观测数据,提取电离层斜延迟.随后,基于单站电离层建模技术,固定精密已知的卫星DCB信息(来源于International GNSS Center,IGS),从电离层斜延迟中同时估计vTEC参数和接收机DCB.基于小波或傅里叶变换等方法,分析接收机DCB估值时间序列的趋势性变化规律,以建立其与外界因素(如硬件更替、观测异常等)之间的联系.

该方案存在一个典型的不足,即分析结果易受投影误差和vTEC建模误差的影响.对此,尽管可分别采用:a)强制约束相邻穿刺点处(如10 km范围内)vTEC完全相等(Zhang et al.,2014);和b)仅针对电离层活动较为平静时期(如夜间)实施电离层建模等措施(Kao et al.,2013),但仍难以充分消除它们对接收机DCB估值的影响(Zhang et al.,20102014).另一方面,受单站电离层建模技术的限制,同时估计vTEC参数和接收机DCB需要累积一定时长(若干小时或1天)的数据,并需假设接收机DCB在此期间恒定,以保证参数的可估性、可分离性和精度.这也意味着,该方案提供的接收机DCB估值时间分辨率较低,不便于研究接收机DCB的短期变化(周东旭等,2011).

2)从平滑误差中直接反映接收机DCB的短期变化,记为“Ciraolo方案”(Brunini and Azpilicueta,2009; Ciraolo et al.,2007).其基本原理可概括为:采用零/短基线GPS数据,逐测站地提取电离层斜延迟,并逐卫星地实施站间单差.理论上,该单差过程可完全消除TEC和卫星DCB,而仅余两台接收机DCB之间的差值.对比全部卫星的单差电离层斜延迟之间的符合程度,可描述接收机相对DCB随时间的整体变化趋势、量化其在某一时间范围内的变化程度.

可见,Ciraolo方案有效地避免了投影误差、vTEC建模误差对研究接收机DCB变化的影响,并显著地提高了接收机DCB估值的时间分辨率(Zhang et al.,2012).然而,当采用短基线实施分析时,将无法克服伪距多路径效应的影响.具体表现为,伪距多路径效应和接收机DCB两者的短期变化互相叠加,或在一定程度上相互抵消,从而降低了接收机DCB分析结果的可靠性(Brunini and Azpilicueta,2009).而当采用零基线分析时,尽管可显著削弱多路径效应,但伪距观测噪声经由站间单差过程被放大,仍将给接收机DCB估值带来不可忽略的误差.

为更加可靠地提取接收机DCB的短期变化,在Ciraolo方案的基础上,本文实施了系列的改进,采用非组合精密单点定位(PPP)算法提取电离层斜延迟(张宝成等,2011),并加入了零/短基线双差模糊度固定等约束条件(Teunissen,1995),以尽量克服平滑误差(记为“PPP方案”).随后,推导了一种能从零/短基线GPS数据中直接估计接收机DCB的数学模型.在检验出接收机DCB的短期变化可被近似为随机游走的前提下,对比分析由PPP方案提取和由模型直接估计的接收机DCB时间序列,可“试探出”最优的过程噪声标准差,进而实现接收机DCB的准确模型化.

2 算法描述

基于零/短基线GPS数据,本节简述了Ciraolo方案的模型算法和实施步骤.之后,通过改进电离层斜延迟提取技术,增强了该方案的可靠性.最后,推导了一种能估计接收机DCB的函数模型,用以定量描述接收机DCB的短期变化.

2.1 接收机DCB短期变化分析:Ciraolo方案

假定构成某零/短基线的两台接收机分别为M和N,且可同步观测P颗GPS卫星.仅以某台接收机为例,其对应于某历元的伪距和相位观测方程可表示为(Leick,2004)

其中,E(·)表示期望运算符;上标s代表第1,…,P颗卫星,i=1,…,n表示历元;下标r对应接收机MNj=1,2表示观测频率L1L2pr,js,iφr,js,i分别为伪距和相位观测值;ρrs,i表示所有与频率无关的未知参数总和;μj为电离层斜延迟ιrs,i的系数,且满足μj=λ2j/λ21λj为浮点模糊度zr,js的波长.

针对式(1),还需指出:

1)zr,js中包含了非差模糊度、伪距硬件延迟、相位初始偏差等几类参数,单位为周;剩余参数和两类观测值的单位均为m;

2)ιrs,i的具体形式为(张宝成等,2010a):

其中,Irs,i表示TEC;BriBs分别为接收机和卫星DCB.

基于式(1),分别对双频伪距和载波相位作差,消除ρrs,i,由此得到:

其中,Δprs,i=pr,2s,j-pr,1s,jΔφrs,ir,1s,j-φr,2s,j分别表示无几何影响的伪距和相位观测值;偏差参数Δzrs=λ1zr,1s-λ2zr,2s.

随后,对全部n个历元的Δφrs,iΔprs,i作差,并取平均,可得Δzrs的平均值

最后,从各历元的Δφrs,i中移除Δrs,并除以常数μ2-1,即可得到对应于ιrs,i的估计值rs,i

将同一卫星至两台接收机的 作差,并乘以常数μ2-1,可逐卫星、逐历元地计算接收机DCB,具体公式为

其中,分别表示卫星s至接收机M和N的电离层斜延迟估值;BMNi的物理含义为接收机M和N的DCB之差.

需要指出,平滑误差将导致同一时段、不同卫星的BMNi(或不同时段、对应同一卫星的BMNi)之间存在差异.对于短基线而言,伪距多路径和BMNi短期变化等两种效应相互叠加、彼此影响,因此难以单独地考察BMNi随时间变化的量级和趋势.这也意味着,充分、有效地消除平滑误差,将是提高Ciraolo方案可靠性的关键.

2.2 接收机DCB短期变化分析:PPP方案

线性化式(1)的伪距和相位观测方程,可得:

与式(1)相比,此处若干新符号的含义为 分别为“观测减计算”的伪距和相位观测值,其中,近似站星距、卫星钟差、若干>1 cm的系统误差均已被事先改正;xr表示测站近似位置增量,其系数crs,i表示站星单位方向向量;grs,i表示ZTD τr的投影函数;dtri为接收机钟差.显然,式(7)即为“非组合”PPP的满秩观测方程(张宝成等,2010a).

联合所有历元形如式(7)的观测方程,对两类观测值实施高度角加权,采用卡尔曼滤波算法,可递归地估计如下的状态向量Σ及其对称的协方差阵QΣ

由式(9)可知,两组估值(ιrs,i,zr,js)T(r=M,N)之间不存在相关性.此外,考虑到接收机M和N距离较近,各自观测方程的空间结构差异可以忽略,此时可合理地近似:

其中Cov()表示协方差函数.

经由双差运算,可从(zM,jszN,js)T中形成2(P-1)个独立整周模糊度ΔΔzj

其中 分别表示星间单差和站间单差运算,且分别将第一颗卫星(或第一台接收机)选取为基准星(或基准站);eI分别是元素均为1的列向量和单位矩阵,下标代表其维数.

根据误差传播定律,新向量 的协方差矩阵 计算公式为

其中,D=blg(IPIPD1D2),blg()表示块对角矩阵.

ΔΔzj的实数解被成功固定为一组整数(设为 )后,利用如下公式,可最终导出 的“固定解”

其中,的两类子矩阵,分别表示(ιs,iM,ιs,iN)TΔΔzj的协方差阵,以及ΔΔzj的协方差阵.

至此,类似于式(6)中的运算,通过对各历元、各卫星的 实施站间单差,并乘以常数μ2-1,可得

与式(6)类似,此处由PPP方案同样获得了接收机DCB的估值BMNi.

2.3 接收机DCB的估计和建模

考虑到零/短基线观测条件的特殊性,基于站间单差观测值,还可构建一种能直接估计接收机DCB的函数模型.针对接收机M和N,对应于式(7)的站间单差观测方程可表示为

其中,(·)MN=(·)M-(·)N表示站间单差运算,αjj/(μ2-1).分析可知:站间单差消除了全部与卫星相关的未知参数和两类大气延迟参数.对于零基线,两台接收机的测站位置相同,此时xMN=03×1.

为便于实施双差模糊度固定,可对式(16)中的dMNizsMN,j实施重新参数化.通过将某卫星(如q,作用类似于参考星)的双频模糊度zMN,jqdtMNi合并,可导出如下等价的相位观测方程:

其中,成为了与频率有关的接收机相位钟差; 则表示具备整周特性的双差模糊度.易知,对应于参考星q的 .

联合式(15)和式(17),可实现BMNi的直接估计,数据处理策略包括: 1)利用卡尔曼滤波实施参数估计; 2)对观测值的加权需同时考虑高度角相关和由站间单差引起的数学相关; 3)各类参数的动态模型分别为:xMN时不变参数(仅存在于短基线模型);dtMNi白噪声过程;BMNi随机游走; 时不变参数; 4)采用LAMBDA搜索并固定 ,并采用Ratio检验其可靠性(Teunissen and Verhagen,2009; Teunissen 1995).

上述数据处理可用于模型化接收机DCB的短期变化,其基本思路是:针对BMNi估计,通过“试探”不同的过程噪声标准差取值,将估值与PPP方案得到的结果相比较,当两者较好地符合时,即可确定BMNi过程噪声标准差的最优“经验值”.值得注意,该模型化方案的准确性取决于:相比Ciraolo方案,前述PPP方案能否更为可靠地提取BMNi;以及BMNi的实际动态模型是否可近似为随机游走.随后的实验分析将分别证明这两种前提的合理性.

3 实验分析3.1 数据准备和处理

实验采用了4台GPS接收机所采集的双频伪距(C1+P2)和相位(L1+L2)观测值.3组接收机天线位于荷兰代尔伏特理工大学某建筑楼顶,其最大间距不超过15 m(见图 1).其中,接收机dlf4和dlf5共用一组天线,因此形成了一条零基线;接收机dlft和delf则被分别连接至另两组不同的天线.表 1汇总了全部接收机(天线)的硬件类型、观测时段、采样间隔等信息.需要指出,实验期间,荷兰采用夏时制,当地时(LT)与协调世界时(UTC)之间的关系为LT=UTC+2.

图 1 实验GPS接收机的空间分布及站间距信息
来源:http://gnss1.tudelft.nl/dpga/station/Delft.html # DELF.
Fig. 1 Spatial distribution and inter-site distances of experimental GPS receivers
Source: http://gnss1.tudelft.nl/dpga/station/Delft.html # DELF.

表 1 本文实验所采用的接收机特征 Table 1Characteristics of experimental receivers used in this paper

为验证本文提出的一系列模型与算法,将设计三类实验方案: 1)分别采用Ciraolo和PPP方案处理若干基线若干天的观测数据,对比验证PPP方案估计电离层斜延迟的可靠性、提取接收机DCB短期变化的准确性等; 2)基于相同基线部分天的观测数据,利用PPP方案提取接收机DCB时间序列.随后,验证该时间序列在多天内是否具有重复性,以及是否可采用随机游走描述; 3)采用与第2种方案同样的实验数据,直接估计接收机DCB,并将其动态模型选取为随机游走,通过“试探”不同的过程噪声标准差,寻求描述接收机DCB短期变化的最优经验模型.

与上述实验方案有关的数据处理设置为: 1)针对2.1节所介绍的Ciraolo方案,卫星截止高度角选为20°,以尽量避免低高度角伪距观测误差影响;仅处理超过1 h(即120个历元)的卫星连续弧段,以确保取平均能尽量消除各类观测噪声; 2)针对2.2节所介绍的PPP方案,卫星截止高度角选取为5°;伪距和相位的权比为1∶104,并考虑卫星高度角加权,其余设置前文已述; 3)针对2.3节所介绍的接收机DCB估计算法,由于不需要联合估计ZTD和位置参数,将卫星截止高度角选取为20°,伪距和相位的权比为1∶104,其余设置前文已述.特别地,采用卡尔曼滤波实施参数估计时,各类状态向量的初值由最小二乘平差首历元的全部观测值得到.

3.2 实验分析

实验处理了4台接收机所形成的6条基线共9天的观测数据.限于篇幅,在随后论述中,仅分析了4组代表性的接收机DCB时间序列,分别对应DOY 170和172的零基线dlf4-dlf5(其DCB变化较为平稳)和短基线dlft-delf(其DCB变化最为显著).

3.2.1 实验一

图 2表示从零基线GPS数据中提取的接收机C1-P2 DCB时间序列.分析DOY 172的结果可知:

图 2 基于Ciraolo方案,从零基线(dlf5-dlf4) GPS观测值中提取的接收机C1-P2 DCB时间序列
横轴:2010年积日170和172(当地时间);纵轴:不同的颜色代表不同的卫星.
Fig. 2 The C1-P2 receiver DCB time series retrieved from GPS observations collected by a zero-baseline (dlf5-dlf4) using Ciraolo scheme
X-axis: days 170 and 172 of year 2010 (local time); Y-axis: different colors correspond to different satellites.

1)对应不同卫星的接收机DCB时间序列互不重合,某些时段(如14∶00 LT前后)内的差异可超过0.5 TECu(约8 cm).造成该差异的主因是低频伪距观测噪声,取平均运算一般难以有效地处理这类噪声.

2)对应不同卫星、同一时段的接收机DCB时间序列变化趋势一致,证实了接收机DCB的确存在短期变化,其范围约为0.5~1 TECu.该变化量级较小,一种可能的原因是:构成该零基线的两台接收机硬件类型、观测环境一致,各自接收机绝对DCB的量级、变化具备一定程度的相似性,并已通过站间单差大大抵消.

该接收机DCB在2天内的变化既有类似性,又存在差异.类似性体现为,针对各天的某些特定时刻,如06∶00 LT和22∶00 LT附近,全部卫星的接收机DCB提取值存在“跃迁”现象,其原因可能是环境温度的上升或下降(注:06∶00 LT和22∶00 LT分别对应当地近似的日出和日落时刻1)).差异性则体现为:日边界处的接收机DCB结果并不连续,表明除短期变化外,接收机DCB还存在长时间尺度,如至少1天的变化,这是由不同天内的观测数据质量不同,以及实施逐天数据处理等两方面因素造成的.

1)来源见:http://www.time and date.com/worldclock/astronomy.html?n=16&month=6&year=2010&obj=sun&afl=-11&day=1

相应地,图 3对应于短基线dlft-delf(间距约为15 m)的接收机C1-P2 DCB时间序列.仍以DOY172为例,与图 2比较可发现:受多路径效应的影响,针对同一观测时段,不同接收机DCB时间序列之间的差异更为明显.与图 2类似,同一天内的接收机DCB变化受气温影响,不同天内接收机DCB变化程度不同.

图 3 基于Ciraolo方案,从短基线(dlft-delf) GPS观测值中提取的接收机C1-P2 DCB时间序列
横轴:2010年积日170和172(当地时间);纵轴:不同的颜色代表不同的卫星.
Fig. 3 The C1-P2 receiver DCB time series retrieved from GPS observations collected by a short-baseline (dlft-delf) using Ciraolo scheme
X-axis: days 170 and 172 of year 2010 (local time); Y-axis: different colors correspond to different satellites.

综合上述讨论可知,Ciraolo方案的典型不足可概括为:1)针对零基线分析,当接收机DCB在某一时段内的变化<0.5 TECu时,分析结果的可靠性易受低频观测噪声的不利影响;2)针对短基线分析,受多路径效应的影响,对应相同时段、不同卫星的接收机DCB时间序列存在若干TECu的差异,该差异在一定程度上“淹没”了接收机DCB的真实变化特性,降低了分析的准确性.

分别与图 2图 3相对应,图 4图 5中给出了利用PPP方案计算的接收机DCB时间序列.其中,受卡尔曼滤波收敛性的影响,各天内前30 min的结果已事先剔除.

图 4 与图2相对应,基于PPP方案,从零基线(dlf5-dlf4) GPS观测值中提取的接收机C1-P2 DCB时间序列 为避免滤波初始化的影响,此处剔除了各天起始30 min内的结果 Fig. 4 Corresponding to Fig.2, it shows here the C1-P2 receiver DCB time series retrieved from GPS observations collected by a zero-baseline (dlf5-dlf4) using PPP scheme. To get rid of possible impacts due to Kalman-filter initialization, the first 30 min′s solution of receiver DCB has been excluded for each experimental day

图 5 与图3相对应,基于PPP方案,从同样的短基线(dlft-delf)GPS数据中 提取的接收机C1-P2 DCB时间序列,单位为TECu Fig. 5 Corresponding to Fig.3, use of the PPP scheme has been made here to retrieve the C1-P2 receiver DCB time series from GPS observations collected by a short-baseline (dlft-delf)

对比图 4图 2的结果可发现:1)对应于不同卫星的接收机DCB时间序列互相重合,相互之间的差异可以忽略,表明改进方案可以较好地克服低频观测噪声的影响,而进一步突出了接收机DCB真实的变化特性; 2)尽管两图所反映的接收机DCB变化行为大体上类似,但图 4可揭示更多的细节:以DOY 172为例,接收机DCB在每天的凌晨(02∶00 LT)附近最小,约为3.5 TECu.而对应于接收机DCB最大的两个时刻分别为08∶00和20∶00 LT,接收机DCB的取值均约为4.5 TECu.这表明,该接收机DCB变化中可能还存在一个周期约为12 h的趋势项,即半日周期项,其幅度约为0.8 TECu.

同时,对比分析图 5图 3中的结果可知: 1)实施PPP方案后,同一时段、不同接收机DCB时间序列之间的差异由先前的4~6 TECu,显著地减少至不超过0.5 TECu.联合前述分析可知,PPP方案能同时克服低频观测噪声和多路径效应的影响,以准确、可靠地还原接收机DCB的本质变化特征; 2)此外,图 5所反映的接收机DCB变化趋势中,不存在类似于图 4中的半日周期项.这说明,该周期项可能与零基线接收机的内部电子元器件温度变化有关,而非由外部环境因素所造成.

总之,与Ciraolo方案相比,利用PPP方案分析接收机DCB时,能有效地处理低频观测噪声、多路径效应的综合影响,可靠地还原接收机DCB不同时间尺度的短期变化特性.基于PPP方案的结果,下一节将验证能否采用随机游走描述接收机DCB在一天内的变化行为.

3.2.2 实验二

首先,基于图 4中DOY 170的结果,针对每一个历元时刻,选取不同接收机DCB时间序列的中位数作为“采样值”,以形成一组接收机DCB中位数时间序列.对该时间序列实施历元间差分,相应结果的频率分布直方图如图 6所示.同时还绘出了两组零均值正态分布的概率密度函数曲线,其标准差分别为0.004和0.005,以便于揭示该直方图的总体分布规律.显见,该直方图整体上接近正态分布,进而表明接收机DCB时间序列确实可被近似为随机游走,其过程噪声标准差的经验值介于0.004,0.005 TECu之间.

图 6 历元间差分的接收机DCB中位数时间序列:频率分布直方图(该时间序列是由图4中年积日170的结果所导出) Fig. 6 The normalized histogram of time-differenced time series of median receiver DCB (The time series has been derived based on the DOY 170 results given in Fig.4)

基于图 4中DOY 172的结果,可形成另一组历元间差分的接收机DCB中位数时间序列,图 7中绘出了与其对应的频率分布直方图.类似地,该直方图也接近正态分布,可同样选取随机游走描述该组接收机DCB随时间的短期变化.

图 7 历元间差分的接收机DCB中位数时间序列:频率分布直方图(该时间序列是由图4中年积日172的结果所导出) Fig. 7 The normalized histogram of time-differenced time series of median receiver DCB (The time series has been derived based on the DOY 172 results given in Fig.4)

紧接着,基于图 5中DOY 170和172的结果,分别计算了另两组历元间差分的接收机DCB中位数时间序列,各自的频率分布直方图分别绘于图 8图 9.对比图 8图 6的结果可知:1)图 8中采样值的跨度变大,绝大部分处于-0.03,0.03 TECu之间,但部分采样值出现在-0.06 TECu附近.这意味着,与图 6相比,该组接收机DCB在相邻历元间的差异更为明显;2)图 8中直方图的分布规律仍然接近正态分布,但其标准差介于0.012和0.014 TECu之间.当采用随机游走描述该组接收机DCB时,相应的过程噪声标准差可经验地限定为略小于0.012 TECu,以确保能在一定程度上“隔离”残余多路径效应的影响.

图 8 历元间差分的接收机DCB中位数时间序列:频率分布直方图(该时间序列是由图5中年积日170的结果所导出) Fig. 8 The normalized histogram of time-differenced time series of median receiver DCB (The time series has been derived based on the DOY 170 results given in Fig.5)

图 9 历元间差分的接收机DCB中位数时间序列:频率分布直方图(该时间序列是由图5中年积日172的结果所导出) Fig. 9 The normalized histogram of time-differenced time series of median receiver DCB (The time series has been derived based on the DOY 172 results given in Fig.5)

另一方面,图 9中采样值的零均值特性不如图 8明显.直观上,采样值出现负值的频率略高于正值.通过分析图 5中相应的结果,可归纳出导致上述差异的原因,如接收机DCB在两天内的变化幅度、趋势明显不同,且在DOY 172内出现了量级较大的显著变化等.尽管如此,当忽略个别区间(如0.03,0.04)内的若干异常采样时,图 9中的频率分布直方图仍近似满足正态分布,相应的标准差也略低于0.012 TECu.

至此,综合上述可知:图 6—9中的频率分布直方图均可以被合理地近似为正态分布,由此表明,对应于图 4图 5中的接收机DCB时间序列均可采用随机游走加以描述,且各自过程噪声的标准差可被经验地限定于0.004,0.005(约0.6~0.8 mm,对应dlf5-dlf4零基线)和0.01,0.012(约1.6~1.9 mm,对应dlft-delf短基线)TECu两个区间范围附近.

3.2.3 实验三

基于零基线dlf5-dlf4 DOY 170观测数据,滤波估计了一系列接收机DCB时间序列,具体结果见图 10.其中,接收机DCB的动态模型为随机游走,并分别选取了3种不同的过程噪声标准差:0.8,1.0和1.5 mm,各自的滤波解对应于图中的黄色、绿色和红色线条.为便于对比分析,图 10中还绘出了图 4中DOY 170的结果,即基于PPP方案提取的接收机DCB时间序列(蓝色线),用于衡量各滤波时间序列的可靠性.

图 10 零基线(dlf5-dlf4)的接收机C1-P2 DCB时间序列(年积日170),和一系列采用随机游走描述其动态变化得到的滤波值时间序列 Fig. 10 The time series of receiver C1-P2 DCB corresponding to zero-baseline (dlf5-dlf4) & DOY170, and its filtered solutions obtained by characterizing its dynamic behavior as random walk

分析图 10可知:

1)当接收机DCB呈近似线性的平稳变化时,如首尾两端的若干小时内,对应于3种过程噪声标准差的滤波时间序列相互重合,其差异较小.这表明,当接收机DCB不存在短期的“突变”时,可考虑选取较小的过程噪声标准差(如0.8 mm),以增强滤波解的可靠性;

2)然而,当接收机DCB存在明显的短期变化时,对应于0.8 mm标准差的滤波时间序列(黄色线)过于平滑,与对应的提取值之间存在明显的差异.而对应于另外两个较大标准差的滤波时间序列均可以较好地描述这类变化;

3)进一步对比对应于标准差1.0 mm和1.5 mm的两个滤波时间序列(分别为绿线和红线)可发现,红色线与接收机DCB提取值之间的整体吻合程度略优于绿色线,原因在于其更有效地刻画了更多接收机DCB的“细节”变化行为,如短期波动等.

类似地,图 11绘出了从零基线dlf5-dlf4 DOY172观测数据中滤波估计的3组接收机DCB时间序列.通过与图 10对比分析,可得出类似的结论.但另一方面,图 11还进一步表明:尽管接收机DCB在DOY 170和172的提取值存在明显区别,但它们的变化行为都被准确地模型化为过程噪声约为1.5 mm的随机游走.该经验标准差的选取一方面确保了接收机DCB滤波值与提取值时间序列在一天内的整体最优符合,同时又能可靠、有效地反映接收机DCB在短期内的不规则变化.

图 11 零基线(dlf5-dlf4)的接收机C1-P2 DCB时间序列(年积日172),和一系列采用随机游走描述其动态变化得到的滤波值时间序列 Fig. 11 The time series of receiver C1-P2 DCB corresponding to zero-baseline (dlf5-dlf4) & DOY 172, and its filtered solutions obtained by characterizing its dynamic behavior as random walk

针对短基线dlft-delf,图 12绘出了DOY 170内3组不同的接收机DCB滤波解,分别对应过程噪声标准差1.0、1.5 mm和2.0 mm.分析可知:当标准差取值为1.0 mm时,相应滤波解(黄线)在某些特定时段(如14∶00—18∶00 LT)与提取值的差异最大约为2 TECu,因此其整体可靠性最差;而当将标准差增大至1.5 mm(该取值对应图 10图 11中的“最优”经验标准差)或2.0 mm时,可明显地改善滤波解与提取值的符合程度.需要指出,在若干时段,如10∶00 LT后出现的第一个U形变化期间,以及22∶00 LT之后的两小时内,3种滤波解与提取值之间均存在一定的差异,量级甚至超过1 TECu.造成这种现象的原因是接收机DCB提取值受一定量的误差,如多路径效应等的影响.

图 13描述了对应短基线dlft-delf在DOY 172内的3组接收机DCB滤波时间序列.与图 12相类似,对应于2.0 mm标准差的滤波解在个别时段仍与提取值存在略低于0.5 TECu的偏差,且仍可归结为多路径效应的影响.

图 12 短基线(dlft-delf)的接收机C1-P2 DCB时间序列(年积日170),和一系列采用随机游走描述其动态变化得到的滤波值时间序列 Fig. 12 The time series of receiver C1-P2 DCB corresponding to short-baseline (dlft-delf) & DOY 170, and its filtered solutions obtained by characterizing its dynamic behavior as random walk

图 13 短基线(dlft-delf)的接收机C1-P2 DCB时间序列(年积日172),和一系列采用随机游走描述其动态变化得到的滤波值时间序列 Fig. 13 The time series of receiver C1-P2 DCB corresponding to short-baseline (dlft-delf) & DOY172, and its filtered solutions obtained by characterizing its dynamic behavior as random walk

综合分析图 1013,可得两个主要结论: 1)针对零基线dlf5-dlf4,其接收机DCB在DOY 170和172两天内的变化均可采用过程噪声标准差为1.0~1.5 mm的随机游走加以描述. 2)针对短基线dlft-delf,针对每天起始的若干小时内,选取较小的标准差如1.0 mm仍能得到较为可靠的滤波解,但在另外两个典型时段,对应于2.0 mm标准差的滤波解仍与提取值存在0.5~1 TECu的差异,该差异应主要归因于多路径效应的影响.

上述结论还表明,采用单独的过程噪声标准差可能不足以准确地刻画接收机DCB在一天内的变化.建议引入一种能探测接收机DCB变化的统计检验量,以自适应地确定不同时段内合适的过程噪声标准差,最终实现接收机DCB滤波值的准确、可靠估计.

4 总结与讨论

接收机DCB的量级和变化与测站环境、硬件设施等密切相关,是利用GPS研究电离层的主要误差源之一.特别地,当接收机DCB存在明显的短期变化时,将会降低各种电离层参数如vTEC的准确性.现有分析接收机DCB变化规律的各种方案均存在不足,其结果难以避免地受平滑误差或(和)模型误差的影响.基于零/短基线GPS数据,本文改进了Ciraolo方案,通过精化电离层延迟估计技术(利用PPP取代相位平滑伪距,并引入双差整周模糊度约束),削弱了低频观测噪声、多路径效应的影响,增强了接收机DCB提取值的可靠性.在明确接收机DCB可被近似为随机游走后,为确定适当的过程噪声标准差,本文导出了从站间单差观测值中直接估计接收机DCB的算法,通过考察接收机DCB提取值和估计值相互之间的符合程度,经验地模型化了接收机DCB的短期变化.

本文研究表明: 1)与相位平滑伪距相比,PPP提取的电离层延迟受低频观测噪声和多路径效应影响更小.针对零/短基线GPS数据处理而言,附加测站间双差独立整周模糊度约束,可进一步增强电离层延迟估值的可靠性; 2)相比Ciraolo方案,PPP方案可以获取更为准确的接收机DCB时间序列,进而能分析出更为细节的接收机DCB短期变化趋势; 3)针对本文实验所采用的零、短基线,各自对应的接收机DCB在一天内的最大变化约为1.5和12 TECu,但均可以被模型化为随机游走,相应的过程噪声标准差约为1.5和2 mm; 4)特别地,针对短基线而言,单一的过程标准差难以准确地模型化其接收机DCB在一个试验天内的变化,建议引入能自动检验接收机DCB“突变”的统计指标,以便于自适应地确定最优的过程噪声标准差.

最后,需要指出,本文提取和估计的接收机DCB实际上是两台接收机的绝对DCB之差,即接收机相对DCB.针对利用GPS数据估计电离层参数而言,可估的接收机DCB同样是一个相对量,因此本文的相关算法和结论对于GPS电离层研究而言具有直接的参考意义.但针对精密授时等应用而言,则需要标定接收机绝对DCB,以获取无偏的时频信息.此时,本文建议将其中一台接收机连接至跨接器(jumper),进而模拟GPS双频信号在接收机内部各通道的传播路径,即可准确地确定该接收机DCB的绝对量级(Wilson and Mannucci 1993).此后,基于该已知的接收机绝对DCB,结合本文估计的接收机相对DCB时间序列,可实现任意一台接收机绝对DCB的检校和标定.

致谢 中科院测量与地球物理研究所iGMAS(International GNSS Monitoring and Assessment System)分析中心提供软件平台.

参考文献
[1] Brunini C, Azpilicueta F J. 2009. Accuracy assessment of the GPS-based slant total electron content. J. Geod., 83(8):773-785.
[2] Brunini C, Azpilicueta F. 2010. GPS slant total electron content accuracy using the single layer model under different geomagnetic regions and ionospheric conditions. J. Geod., 84(5):293-304.
[3] Ciraolo L, Azpilicueta F, Brunini C, et al. 2007. Calibration errors on experimental slant total electron content(TEC) determined with GPS. J. Geod., 81(2):111-120.
[4] Conte J F, Azpilicueta F, Brunini C. 2011. Accuracy assessment of the GPS-TEC calibration constants by means of a simulation technique. J. Geod., 85(10):707-714.
[5] Coster A, Williams J, Weatherwax A, et al. 2013. Accuracy of GPS total electron content:GPS receiver bias temperature dependence. Radio Sci., 48(2):190-196.
[6] Dyrud L, Jovancevic A, Brown A, et al. 2008. Ionospheric measurement with GPS:Receiver techniques and methods. Radio Sci., 43(6), doi:10.1029/2007RS003770.
[7] Gao X W, Liu J N, Ge M R. 2002. An ambiguity searching method for network RTK baselines between base stations at single epoch. Acta Geodaetica et Cartographica Sinica(in Chinese), 31(4):305-309.
[8] Hernández-Pajares M, Juan J M, Sanz J, et al. 2011. The ionosphere:effects, GPS modeling and the benefits for space geodetic techniques. J. Geod., 85(12):887-907.
[9] Kao S, Tu Y, Chen W, et al. 2013. Factors affecting the estimation of GPS receiver instrumental biases. Surv. Rev., 45(328):59-67.
[10] Komjathy A, Sparks L, Mannucci A J, et al. 2005. The ionospheric impact of the October 2003 storm event on wide area augmentation system. GPS Solut., 9(1):41-50.
[11] Leick A. 2004. GPS Satellite Surveying. 3rd ed.Hoboken, NJ:Wiley.
[12] Li W, Cheng P F, Bei J Z, et al. 2012. Calibration of regional ionospheric delay with uncombined precise point positioning and accuracy assessment. J. Earth Syst. Sci., 121(4):989-999.
[13] Li Z S, Yuan Y B, Li H, et al. 2012. Two-step method for the determination of the differential code biases of COMPASS satellites. J. Geod., 86(11):1059-1076.
[14] Li J Y, You X Z, Zhang R, et al. 2015. Ionospheric total electron content disturbance associated with May 12, 2008, Wenchuan earthquake. Geod Geodyn, 6(2):126-134.
[15] Liu J N, Ye S R. 2002. GPS Precise point positioning using undifferenced phase observation. Geomat. Inf. Sci. Wuhan Univ.(in Chinese), 27(3):234-240.
[16] Mao T, Wang J S, Yang G L, et al. 2010. Effects of typhoon Matsa on ionospheric TEC. Chinese Sci. Bull., 55(8):712-717.
[17] Nayir H, Arikan F, Arikan O, et al. 2007. Total electron content estimation with Reg-Est. J. Geophys. Res., 112:A11313, doi:10.1029/2007JA012459.
[18] Øvstedal O. 2002. Absolute positioning with single-frequency GPS receivers. GPS Solut., 5(4):33-44.
[19] Petrie E J, Hernández-Pajares M, Spalla P, et al. 2011. A review of higher order ionospheric refraction effects on dual frequency GPS. Surv. Geophys., 32(3):197-253.
[20] Shen Y Z, Chen Q J, Xu H Z. 2015. Monthly gravity field solution from GRACE range measurements using modified short arc approach. Geod Geodyn, 6(4):261-266.
[21] Teunissen P J G. 1995. The least-squares ambiguity decorrelation adjustment:a method for fast GPS integer ambiguity estimation. J. Geod., 70(1-2):65-82.
[22] Teunissen P J G, Verhagen S. 2009. The GNSS ambiguity ratio-test revisited:a better way of using it. Surv. Rev., 41(312):138-151.
[23] Teunissen P J G, Odijk D, Zhang B C. 2010. PPP-RTK:results of CORS network-based PPP with integer ambiguity resolution. J. Aero. Astron. Aviat., 42(4):223-230.
[24] Wilson B D, Mannucci A J. 1993. Instrumental biases in ionospheric measurement derived from GPS data.//Proceedings of the 6th International Technical Meeting of the Satellite Division of the Institute of Navigation. Salt Lake City, UT, 1343-1351.
[25] Yang Y X, Li J L, Xu J Y, et al. 2011. Contribution of the compass satellite navigation system to global PNT users. Chinese Sci. Bull., 56(26):2813-2819.
[26] Yang Y X, Li J L, Wang A B, et al. 2014. Preliminary assessment of the navigation and positioning performance of BeiDou regional navigation satellite system. Sci. China Earth Sci., 57(1):144-152.
[27] Yao Y B, Chen P, Wu H, et al. 2012. Analysis of ionospheric anomalies before the 2011Mw9.0 Japan earthquake. Chinese Sci. Bull., 57(5):500-510.
[28] Yao Y B, Zhang R, Song W W, et al. 2013. An improved approach to model regional ionosphere and accelerate convergence for precise point positioning. Adv. Space Res., 52(8):1406-1415.
[29] Yu T, Mao T, Wang Y G, et al. 2009. Study of the ionospheric anomaly before the Wenchuan earthquake. Chinese Sci. Bull., 54(6):1080-1086.
[30] Yuan Y B, Ou J K. 2001. An improvement to ionospheric delay correction for single-frequency GPS user-the APR-I scheme. J. Geod., 75(5-6):331-336.
[31] Yuan Y B, Ou J K. 2004. A generalized trigonometric series function model for determining ionospheric delay. Prog. Nat. Sci., 14(11):1010-1014.
[32] Yuan Y B, Li Z S, Wang N B, et al. 2015. Monitoring the ionosphere based on the Crustal Movement Observation Network of China. Geod Geodyn, 6(2):73-80.
[33] Yue X, Schreiner W S, Hunt D C, et al. 2011. Quantitative evaluation of the low Earth orbit satellite based slant total electron content determination. Space Weather, 9:S09001, doi:10.1029/2011SW000687.
[34] Zhang B C, Ou J K, Yuan Y B, et al. 2010a. Precise point positioning algorithm based on original dual-frequency GPS code and carrier-phase observations and its application. Acta Geodaetica et Cartographica Sinica(in Chinese), 39(5):478-483.
[35] Zhang B C, Ou J K, Yuan Y B, et al. 2010b. Yaw attitude of eclipsing GPS satellites and its impact on solutions from precise point positioning. Chinese Sci. Bull., 55(32):3687-3693.
[36] Zhang B C, Ou J K, Li Z S, et al. 2011. Determination of ionospheric observables with precise point positioning. Chinese J. Geophys.(in Chinese), 54(4):950-957, doi:10.3969/j.issn. 0001-5733. 2011.04.009.
[37] Zhang B C, Ou J K, Yuan Y B, et al. 2012. Extraction of line-of-sight ionospheric observables from GPS data using precise point positioning. Sci. China Earth Sci., 55(11):1919-1928.
[38] Zhang B C, Teunissen P J G. 2015. Characterization of multi-GNSS between-receiver differential code biases using zero and short baselines. Sci. Bull., 60(21):1840-1849.
[39] Zhang D H, Zhang W, Li Q, et al. 2010. Accuracy analysis of the GPS instrumental bias estimated from observations in middle and low latitudes. Ann. Geophys., 28(8):1571-1580.
[40] Zhang D H, Shi H, Jin Y Q, et al. 2014. The variation of the estimated GPS instrumental bias and its possible connection with ionospheric variability. Sci. China Technol. Sci., 57(1):67-79.
[41] Zhang X H, Ren X D, Wu F B, et al. 2013a. A new method for detection of pre-earthquake ionospheric anomalies. Chinese J. Geophys.(in Chinese), 56(2):441-449, doi:10.6038/cjg20130208.
[42] Zhang X H, Tang L, Guo B F. 2013b. Research on medium-scale traveling ionospheric disturbances using a modified SRTI method. Chinese J. Geophys.(in Chinese), 56(12):3953-3959, doi:10.6038/cjg20131201.
[43] Zheng W, Xu H Z. 2015. Progress in satellite gravity recovery from implemented CHAMP, GRACE and GOCE and future GRACE Follow-On missions. Geod Geodyn, 6(4):241-247.
[44] Zhou D X, Yuan Y B, Li Z S, et al. 2011. Analysis of long-term variations of GPS receivers' differential code bias. J. Geod. Geodyn.(in Chinese), 31(5):114-118.
[45] 高星伟, 刘经南, 葛茂荣. 2002. 网络RTK基准站间基线单历元模糊度搜索方法. 测绘学报, 31(4):305-309.
[46] 刘经南, 叶世榕. 2002. GPS非差相位精密单点定位技术探讨. 武汉大学学报:信息科学版, 27(3):234-240.
[47] 毛田, 王劲松, 杨光林等. 2009. 台风"麦莎"对电离层TEC的影响. 科学通报, 54(24):3858-3863.
[48] 姚宜斌, 陈鹏, 吴寒等. 2012. 2011年3月11日日本地震震前电离层异常变化分析. 科学通报, 57(5):355-365.
[49] 余涛, 毛田, 王云冈等. 2009. 汶川特大地震前电离层主要参量变化. 科学通报, 54(4):493-499.
[50] 袁运斌, 欧吉坤. 2005. 广义三角级数函数电离层延迟模型. 自然科学进展, 15(8):1015-1019.
[51] 张宝成, 欧吉坤, 袁运斌等. 2010a. 基于GPS双频原始观测值的精密单点定位算法及应用. 测绘学报, 39(5):478-483.
[52] 张宝成, 欧吉坤, 袁运斌等. 2010b. GPS卫星姿态异常及其对精密单点定位估值的影响. 科学通报, 55(27-28):2712-2718.
[53] 张宝成, 欧吉坤, 李子申等. 2011. 利用精密单点定位求解电离层延迟. 地球物理学报, 54(4):950-957, doi:10.3969/j.issn.0001-5733.2011.04.009.
[54] 张小红, 任晓东, 吴风波等. 2013a. 震前电离层TEC异常探测新方法. 地球物理学报, 56(2):441-449, doi:10.6038/cjg20130208.
[55] 张小红, 唐龙, 郭博峰. 2013b. 利用改进的SRTI法研究中尺度电离层行波扰动. 地球物理学报, 56(12):3953-3959, doi:10.6038/cjg20131201.
[56] 周东旭, 袁运斌, 李子申等. 2011. GPS接收机仪器偏差的长期变化特性分析. 大地测量与地球动力学, 31(5):114-118.