2. 大地测量与地球动力学国家重点实验室, 武汉 430077;
3. 中国科学院光电研究院, 北京 100094
2. State Key Laboratory of Geodesy and Earth's Dynamics, Wuhan 430077, China;
3. Academy of Opto-Electronics., Chinese Academy of Sciences, Beijing 100094, China
相比于差分定位,PPP(Precise Point Positioning)技术仅利用一个接收机的观测数据以及卫星精密轨道和钟差产品即可实现高精度定位,其工作方式更加灵活,已广泛应用于大地测量等领域(Bisnath and Gao, 2008;Kouba and Héroux,2001;Zumberge et al., 1997). 在实时PPP中,模糊度参数经过一段时间完成收敛,可获得高精度定位结果(Li et al., 2011),即使通过实现PPP模糊度固定,也需要大约10分钟的时间(Geng et al., 2011),其初始化速度比差分模式明显偏慢(Giorgi et al., 2012).当观测环境不好,信号出现遮挡或中断时,模糊度参数需要重新初始化,将导致实时PPP长时间无法获得高精度定位结果,直接影响PPP的实际应用效果,因而必须探求有效的解决措施.
通常情况下,卫星硬件延迟在短时间内的变化可忽略不计(Zhang and Li, 2012),卫星信号中断等因素造成的模糊度参数变化(周跳)可认为是一个整数;如果能够确定该变化量,结合中断前的模糊度参数估计结果,便可精确确定中断后的模糊度,快速获得高精度定位结果(Li et al., 2013;Shi and Gao, 2011);因此,要实现实时PPP快速重新初始化,其关键在于检测和确定相位观测值的周跳,实现周跳修复.
在仅利用单个接收机的观测数据时,周跳检测修复的方法主要可以分为两类:一类方式是采用多种线性组合观测值进行周跳检测与修复(Dai et al., 2008;de Lacy et al., 2008;Xu and Kou, 2011);另一类方式是将周跳或模糊度参数作为未知参数进行估计,进而尝试固定,其中,Banville and Langley(2009)提出分别固定宽巷组合和无几何组合的周跳;Geng等(2010)通过引入预报的电离层延迟信息实现宽巷模糊度快速固定,从而加速PPP模糊度固定;Zhang和Li(2012)利用电离层延迟变化预报结果修正历元间差分观测值中的电离层延迟 变化,进而估计得到周跳参数的浮点解,并尝试固定.
由此可见,在第二类方法中,通过引入电离层延迟约束信息,提高周跳/模糊度参数估计精度和修复成功率,是目前解决模糊度重新初始化问题较好的思路之一.然而,就该思路而言,电离层延迟预报结果的精度将直接影响其效果.由于电离层延迟与测站位置、太阳活动等因素存在密切关系(Shi et al., 2012),不同情形下,电离层延迟预报精度存在显著差异.因此,引入电离层延迟预报信息时必须合理考虑该信息的精度,保证其可靠性,避免其对参数估计结果产生不利的影响.
针对此问题,本文提出了一种充分顾及电离层预报精度的实时PPP快速重新初始化算法,首先讨论了在不依赖其他外部改正信息情形下的周跳修复观测模型及数据处理流程;然后设计了一种有效的电离层预报模型和预报信息定权方法;最后,利用不同区域和不同时期的观测数据全面分析了上述方法在不同情形下的应用效果. 2 周跳修复的观测模型
GPS原始观测方程可表示为
其中,p、m分别表示接收机和卫星,fi(i=1,2)为信号频率; Pmpi、Lmpi 分别为 Li 频率的伪距和相位观测值,ρmp 为卫星与接收机之间的几何距离,Imp1 为 L1 频率的电离层延迟,Tmp 为对流层延迟,dtp、dtm 分别为接收机和卫星的钟差; wpi为Li 频率的接收机相位缠绕误差; λi为Li 频率的波长,Nmpi为Li 频率的模糊度,εmpi、δmpi 为观测值的随机误差.其他系统误差,如潮汐改正、相位中心偏差改正等直接采用模型改正.在动态PPP中,首先采用星间差分方式消除接收机相位缠绕误差,生成的星间差分观测方程可以表示为
其中,m,n为进行星间差分的卫星.假设从 n0 至 nt 历元,卫星信号保持连续,此时模糊度参数 Nmnpi 始终为一常数;由于信号中断或观测数据发生周跳,nt+j 时刻的相位观测数据发生周跳,对历元 nt、nt+j 的观测数据进行历元间差分;并将卫星轨道和钟差固定为已知,同时考虑到对流层延迟在短时间内变化较小,将天顶对流层湿延迟参数固定为实时PPP中的最新估计结果,则生成的观测方程可以表示为
其中,Δ 表示历元间差分; CSmnpi为Li 频率上卫星m,n的周跳之差,为一整数,若该卫星对未发生周跳,则对应的周跳参数为0.可以看到,获取周跳参数的固定解,即可避免模糊度参数重新初始化.依据公式(3)进行参数估计时,宽巷组合周跳参数的浮点解等于采用历元间差分的MW组合的计算结果,仅利用少数几个历元的观测数据时,浮点解精度偏低,难以准确固定为整数;并且,即使宽巷组合固定成功,窄巷组合波长较短,而其浮点解由无电离层组合的伪距观测值决定,因而仍然难以准确固定,所以仅依据公式(3)是无法可靠的实现周跳修复的.针对于此,相关学者提出引入其他类型的观测值或观 测手段以进一步提高浮点解精度(Carcanague,2012; Du and Gao, 2012;Kim et al., 2012).不同于此类方法,本文将在仅利用GPS观测数据的情形下进一步挖掘有效约束信息.由于电离层延迟变化是有一定规律的,在一段时间内的变化比较平缓,通过建立电离层变化预报模型,将其预报结果作为约束信息引入到参数估计中,可以改善周跳参数的估计精度,从而提高周跳修复的可能性.3 电离层延迟变化预报模型及定权策略
基于公式(2),L1 频率的电离层延迟 Imnp1 可以表示为
Dai等(2003)对大气延迟的时间相关性开展了大量分析工作,提出了多种预报模型.本文采用线性拟合模型进行短期电离层延迟建模,表示为如下函数: 其中,Imnp1(t)为t 时刻星间差分的电离层延迟,a,b 为线性函数的系数.利用一个时间窗口(本文中取300 s)内的电离层延迟进行拟合计算式(5)的系数,然后预报电离层延迟变化量,并将预报结果作为虚拟观测值引入到参数估计中.虚拟观测方程可以表示为
其中,为电离层延迟变化的预报结果.由于不同情形下的电离层延迟预报精度不同,将上述信息作为虚拟观测值引入到参数估计时,必须合理确定其精度,以保证后续参数估计结果的精度及可靠性.本文根据电离层预报中前后两端预报精度衰减的对称性确定其中误差,具体实施过程如图 1所示.
假定需要预报(T2,T2+Δt2)时段的电离层延迟变化量,并利用(T1,T2)时段的观测值进行拟合;选择 T1 时刻前存在实测观测数据的时刻 T1-Δt1,且使Δt1与Δt2 的差异尽量小,分别预报(T2,T2+Δt2)和(T1,T1-Δt1)时段的电离层延迟变化量,由于(T1,T1-Δt1)时段的电离层延迟变化量可根据实测数据直接计算得到,因此,
预报值 ΔI ^ mnp1(T1,T1-Δt1)的外符合精度 σ0 可以利用预报值与实测值之间的差异近似地描述;在此基础上,考虑电离层预报精度与预报时间长度和卫星高度角之间的关系,则 σI0、σ0 存在如下关系:
其中,p(E)为与高度角相关的比例因子,如式(8)所示: E 为星间差分的两颗卫星在预报时段内的高度角最小值.p(τI0)为与预报时间长度相关的比例因子(Geng et al., 2010),表示为
τ为预报时段长度,以s为单位.由此,根据式(7)可以计算得到预报结果的单位权中误差 σI0. 4 滤波模型及数据处理流程
在实时周跳修复中,采用Kalman滤波方法进行参数估计.联合式(3)和(6),得到线性化后的观测方程:
其中,P,L,I0 分别为改化的伪距、相位观测值和电离层延迟变化预报结果;X,I,CS 为未知参数,包括测站坐标、电离层延迟变化及周跳参数; A,B 为 P 相对于参数 X,I 的系数阵; E 为单位阵; O 表示元素全部为0的矩阵; M 为 L 相对于参数 CS 的系数阵,当所有卫星发生周跳时,该矩阵为单位阵 E ; Γ 为随机误差.状态方程表示为
其中,过程噪声 εX,εI 的中误差设置为无穷大;当相邻历元卫星对未发生周跳时,εCS 的中误差设置为0,否则设置为无穷大.本文中设定周跳修复的最大连续历元数目为3,超过该限差时,认定修复失败,实时PPP重新初始化.针对每个历元,参数估计后提取 CS 的估计结果,采用LAMBDA算法(Teunissen,1995)尝试固定,并采用ratio检验方法(Wang et al., 1998)进行可靠性验证,ratio检验的限差设置为2.尝试固定的流程如图 2所示,首先提取宽巷组合浮点解并尝试固定,若固定失败,则修复失败;否则,利用宽巷固定解分别更新引入或不引入电离层约束时的周跳参数估计结果,并提取窄巷组合的浮点解,尝试固定,若其中一种模式下固定成功,结合宽巷组合固定解,则可以确定两个频率的周跳,进而修复观测数据.
为了分析上述方法在不同区域和不同时期内的具体效果,根据测站分布及数据质量的情况,在中国和欧洲区域分别选择一个跟踪网.其中,中国区域选择了中国地壳运动监测网络(CMONOC,Crustal Movement Observation Network of China)的18个跟踪站,测站分布如图 3所示.根据测站纬度,可将其划分为三个区域,具体如表 1所示;欧洲区域选择了24个IGS站,测站分布如图 4所示.在上述测站中,欧洲区域纬度较高(几乎全部测站的纬度大于45°),中国区域的测站分布范围较广,测站纬度最低约为25°.利用这些观测数据,可以全面分析不同区域内的差异,从而更好的评价本方法的效果.
原始观测数据的采样间隔为30 s,卫星轨道和钟差直接采用IGS最终精密轨道和钟差产品,实际数据处理中可以采用相关研究机构提供的实时轨道和钟差产品(Ge et al., 2012;Zhang et al., 2007;Zhang et al., 2011).在分析过程中,每30 min(00 h 30 min,01 h 00 min,…)模拟一次周跳修复过程,对观测的所有卫星添加周跳,以分析在最严重的情况下的效果;同时,将部分历元的观测数据删除,以模拟不同的中断时间长度(Time Latency).本文分析了30 s、90 s、150 s、210 s及270 s五种情况.为了定量评价修复效果,本文中统计了周跳修复的成功率(Success Rate,修复成功的过程数目与总的周跳修复过程数目的比值)和平均历元数目(Average Epoch,每次周跳修复过程中所利用的观测数据的历元数目的平均值).为了反映修复效果与电离层预报精度的相关性,本文同时统计了每个时刻电离层延迟预报误差的RMS(IONO RMS). 5.2 电离层约束权函数分析
电离层约束权函数是否正确地反映电离层延迟预报信息的真实精度,对周跳修复效果起着关键作用.本节采用中国区域内的18个测站,分别利用2008年(电离层活动低年)第188天以及2012年(电离层活动高年)第89天的实测数据分析电离层权函数所表征的电离层精度信息与电离层延迟预报精度之间的相关性.首先,将预报的电离层延迟与实测的电离层延迟进行对比,对二者作差得到电离层延迟预报值的实测精度;其次,按照第3节所述方法对上述电离层预报值的精度进行预测;最后,通过比较真实精度与预测精度以评价电离层约束权函数的正确性.
图 5给出了代表高、中、低纬地区的BJFS、LUZH以及KMIN站PRN02卫星电离层延迟预报值的预测精度与实测精度,其中,中断时长取30 s.图左和图右分别给出了2008年和2012年的结果.在低纬地区,电离层预报精度优于10 cm,在高纬地区,其精度优于5 cm.总体上看,电离层预报值的预测精度基本能够反映其真实精度,其原因主要是由于利用公式计算电离层预测精度时充分考虑了(T1,T1-Δt1)时段内电离层预报值的实测精度,从而提高了其准确度.
为了进一步分析不同卫星及中断时长情况下电离层延迟预报值的预测精度与实测精度之间的关 系,表 2给出了所有测站两个参数的相关系数统计结果.相关系数越大,预测精度越能够反映电离层延迟预报值的实测精度.从统计结果上看,相同中断时长的情况下,电离层延迟预报值的预测精度与实测精度的相关系数在不同纬度地区基本相当,由此说明第3节中所采用的计算方法能够有效地顾及到不同纬度地区电离层活动差异对预报精度的影响.换言之,该预测精度基本真实地反映了不同地区电离层延迟预报值的实测精度.随着中断时间的延长,二者之间的相关系数逐渐变小,中断30 s的情况下,相关系数大约为0.85,当中断时长变为270 s时,相关系数仅为0.75左右,造成该现象的原因在于随着中断时间逐渐变长,预报窗口两侧电离层延迟变化的对称性逐渐变弱,从而使得基于公式(7)计算得到的电离层延迟预报精度的可靠性逐渐降低.
从上述分析结果来看,采用第3节中建立的电离层延迟精度的预测方法能够基本反映电离层延迟预报值的真实精度,基于该预测精度构造电离层延迟约束权函数,进而辅助周跳修复是基本可行的. 5.3 不同区域周跳修复效果比较
本节比较分析了不同区域的效果差异,选择2008年第188天的观测数据进行模拟分析,欧洲区域和中国区域的统计结果分别如表 3和表 4所示.
由表 3的统计结果可见,当中断时间长度为30 s 时,几乎全部情况下修复成功,成功率达到99.0%,并且除少数时刻外,绝大部分情况下仅利用一个历元的观测数据就可以修复成功,平均历元数目为1.03;进一步增加中断时间长度,成功率仅略有下降,增加至150 s时,成功率仍然在95%以上,平均 历元数目为1.10;由此证明,本方法可以有效提高
坐标参数和周跳参数的估计精度,从而保证了后续宽巷和窄巷组合固定成功的可能性.当中断时间长度达到270 s时,随着电离层预报误差增大,成功率及平均历元数目进一步降低,忽略多路径误差的影响,电离层预报误差是导致修复成功率下降的关键因素;但是由于欧洲区域电离层活动相对比较平静,修复成功率仍然可以达到80%以上,由此验证本方法在欧洲区域可以取得较好的效果.
同时,以BOR1测站为例说明周跳修复算法对实时定位精度的改善效果,其他测站的情形与之类似.当中断时长为30 s时,采用本文的方法进行周跳修复后,BOR1站的实时定位误差时间序列如图 6a所示;作为对比,图 6b绘制了按照传统模式进行实时动态PPP的定位误差时间序列.可以看到,由于在该时段内周跳修复全部成功,且均仅利用一个历元的观测数据,经过起始阶段的初始化以后,在整个时段内一直保持较高的定位精度,而采用传统模式进行实时动态PPP时,需要反复重新初始化,长时间内无法获得连续高精度定位结果,这显然是无法满足实际工程应用需求的.
对于中国区域,当中断时间长度为30 s时,不同纬度的测站的成功率基本相同,而需要的平均历元数目随着纬度的降低略有增加,低纬地区的平均历元数目为1.14,高纬地区仅为1.03,与欧洲区域的平均历元数目相等.随着中断时间长度增大,成功 率不断降低,并且纬度越低,下降幅度越明显;当中断时间长度增加至90 s时,中低纬区域的成功率降低至大约90%,而在高纬区域,成功率仅略微降低;继续增加至270 s时,低纬地区的成功率仅为37.5%,高纬地区可以达到65.1%.由此可见,周跳修复的效果与测站位置存在密切关系.
为了显著反映电离层预报误差的影响,同时计算了不同时刻电离层预报误差精度,结果如图 7所示,并分别统计了不同纬度区域不同时刻修复成功的测站数目,图 8给出了中断时间长度为150、210 s和270 s时的结果.当中断时间较短(30 s)时,不同区域的测站在一天之内的电离层预报误差基本相同,均保持在很小的误差范围内,一方面,可以改善坐标估计精度,提高宽巷组合固定的几率;另一方面,对窄巷组合估计精度的影响较小,从而保证了周跳修复的可行性.随着中断时间长度增大,对于低纬地区的测站,在两个时段(UTC5 ∶ 00~UTC9 ∶ 00,UTC12 ∶ 00~UTC18 ∶ 00)范围内的电离层误差显著增大,中断时间达到270 s时,最大预报误差达到8 cm,即使宽巷组合固定成功,窄巷组合的成功率也会受到极大影响.由图 8的统计结果可见,在这两个时段内,部分时刻全部测站都修复失败;对于中纬地区,仅有一个时段(UTC13 ∶ 00~UTC17 ∶ 00)的电离层预报误差显著增大,导致该时段内的修复效果明显变差;而对于高纬地区,在大约UTC14 ∶ 00时,电离层预报误差略微增大,修复成功数目有一定减少,但是相比于中低纬地区影响较小.
电离层活跃程度与太阳活动存在密切关系.5.3节中分析的观测数据位于太阳活动平静时期,当太阳活动剧烈时,电离层延迟变化规律降低,预报精度将进一步下降.利用中国区域的测站,选择2012年第89天的观测数据按照相同的方式进行模拟分析,统计结果如表 5所示,同时统计了不同时刻的电离层预报误差(图 9)及修复成功的测站数目(图 10).
由图 9及表 5可以看出,当中断时间长度为30 s时,电离层预报误差保持平稳,均在1.5 cm以内,周跳修复成功率均达到97%以上,与太阳活动平静期的结果相当.由此反映出,当中断时间较短时,周跳修复效果受太阳活动影响较小,不同时期的效果基本相同.
随着中断时间的增加,相比于太阳活动平静时 期,电离层变化受到的影响更大,预报精度下降更 快,并且在一天时间之内,电离层预报误差显著偏大 的时段更长,测站纬度越低,上述现象越明显.在低 纬地区,当中断时长超过150 s时,大部分时刻的电离层预报误差大于4 cm.与之相应,中断时长越大,纬度越低,成功率下降越快.当中断时长达到150 s时,高纬地区的修复成功率接近于80%,而低纬地区仅为55%. 6 结论
本文针对实时动态PPP中由于卫星信号中断等原因导致的PPP重新初始化问题,提出了一种新的算法.通过设计电离层延迟预报模型进行电离层延迟变化预报,并给出预报信息的合理定权方法,以实现电离层延迟预报结果的权重的自适应估计,有效顾及不同地区内电离层活动的差异,保证后续参数估计结果的精度及可靠性;在此基础上,利用周跳前后的观测数据,生成历元间/星间差分观测数据,结合电离层延迟变化预报结果进行参数估计,获取周跳浮点解,进而尝试固定,实现PPP快速重新初始化.
基于上述算法,本文利用不同区域、不同时期的 观测数据进行了大量的模拟分析,全面评价了不同 情形下的效果差异.分析结果表明,当中断时间较短时,上述方法在不同区域和太阳活动情形下均可以达到较高的成功率,保证实时PPP的连续高精度定位,由此验证了本方法的效果,可以有效解决实时动态PPP的重新初始化问题,拓展PPP的实际应用范围;中断时间长度进一步增大时,不同区域和时期的周跳修复效果存在显著差异,纬度越低,太阳活动越剧烈,修复成功率越低.
致谢 感谢地壳运动观测网络提供的GPS观测资料,感谢IGS提供的基准站观测数据.[1] | Banville S, Langley R. 2009. Improving real time kinematic PPP with instantaneous cycle slip correction. Proceedings of ION GNSS 2009, Savannah, Georgia, 2470-2478. |
[2] | Bisnath S, Gao Y. 2008. Current state of precise point positioning and future prospects and limitations. Proceedings of IUGG 24th General Assembly, International Association of Geodesy Symposia, 133:615-623. |
[3] | Carcanague S. 2012. Real-time Geometry-based cycle slip resolution technique for single-frequency PPP and RTK. Proceedings of ION GNSS 2012, Nashville, TN, 1136-1148. |
[4] | Dai L, Wang J, Rizos C, et al. 2003. Predicting atmospheric biases for real time ambiguity resolution in GPS/GLONASS reference station networks. J. Geod., 76(11-12): 617-628. |
[5] | Dai Z, Knedlik S, Loffeld O. 2008. Real-time cycle-slip detection and determination for multiple frequency GNSS.//Proceedings of the 5th Workshop on Positioning, Navigation and Communication.Hannover: IEEE, 37-43. |
[6] | De Lacy M C, Reguzzoni M, Sansò F, et al. 2008. The Bayesian detection of discontinuities in a polynomial regression and its application to the cycle-slip problem. J Geod., 82(9): 527-542. |
[7] | Du S, Gao Y. 2012. Inertial aided cycle slip detection and identification for integrated PPP GPS and INS. Sensors, 12(11): 14344-14362. |
[8] | Ge M, Chen J P, Douša J, et al. 2012. A computationally efficient approach for estimating high-rate satellite clock corrections in realtime. GPS Solut., 16(1): 9-17. |
[9] | Geng J H, Meng X L, Dodson A H, et al. 2010. Rapid re-convergences to ambiguity-fixed solutions in precise point positioning. J. Geod., 84(12): 705-714. |
[10] | Geng J, Teferle F N, Meng X, et al. 2011. Towards PPP-RTK: Ambiguity resolution in real-time precise point positioning. Adv. Space Res., 47(10): 1664-1673. |
[11] | Giorgi G, Teunissen P J G, Verhagen S, et al. 2012. Instantaneous ambiguity resolution in global-navigation-satellite-system-based attitude determination applications: A multivariate constrained approach. Journal of Guidance, Control, and Dynamics, 35(1): 51-67. |
[12] | Kim Y, Song J, Yun H, et al. 2012. Cost-effective selection of inertial sensor for cycle-slip detection for land vehicle. Proceedings of ION GNSS 2012, Nashville, TN, 1580-1588. |
[13] | Kouba J, Héroux P. 2001. Precise point positioning using IGS orbit and clock products. GPS Solut., 5(2): 12-28. |
[14] | Li X X, Zhang X H, Ge M R. 2011. Regional reference network augmented precise point positioning for instantaneous ambiguity resolution. J. Geod., 85(3): 151-158. |
[15] | Li X X, Ge M R, Zhang H P, et al. 2013. A method for improving uncalibrated phase delay estimation and ambiguity-fixing in real-time precise point positioning. J. Geod., 87(5): 405-416. |
[16] | Shi C, Gu S F, Lou Y D, et al. 2012. An improved approach to model ionospheric delays for single-frequency precise point positioning. Adv. Space Res., 49(12): 1698-1708. |
[17] | Shi J, Gao Y. 2011. Integer ambiguity resolution to improve accuracy and convergence of PPP-inferred troposphere estimates. Proceedings of ION GNSS 2011, Portland, OR, 588-596. |
[18] | 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. |
[19] | Wang J, Stewart M P, Tsakiri M. 1998. A discrimination test procedure for ambiguity resolution on-the-fly. J Geod., 72(11): 644-653. |
[20] | Xu D Y, Kou Y H. 2011. Instantaneous cycle slip detection and repair for standalone triple-frequency GPS receiver. Proceedings of ION GNSS 2011, 3916-3922. |
[21] | Zhang Q, Moore P, Hanley J, et al. 2007. Auto-BAHN: Software for near real-time GPS orbit and clock computations. Adv. Space Res., 39(10): 1531-1538. |
[22] | Zhang X H, Li X X, Guo F. 2011. Satellite clock estimation at 1Hz for realtime kinematic PPP applications. GPS Solut., 15(4): 315-324. |
[23] | Zhang X H, Li X X. 2012. Instantaneous re-initialization in real-time kinematic PPP with cycle slip fixing. GPS Solut., 16(3): 315-327. |
[24] | Zumberge J F, Heflin M B, Jefferson D C, et al. 1997. Precise point positioning for the efficient and robust analysis of GPSdata from large networks. J Geophys. Res., 102(B3): 5005-5017. |