BDS卫星短时差分码偏差估计与分析

刘冰雨, 王中元, 胡 超, 周圣淇

(1.中国矿业大学 环境与测绘学院,江苏 徐州 221116; 2.安徽理工大学 空间信息与测绘工程学院,安徽 淮南 232001)

2020年6月,最后一颗北斗三号卫星成功组网,标志着北斗卫星导航系统(BeiDou Navigation Satellite System,BDS)星座部署全面完成。目前,BDS由8颗地球静止轨道(Geostationary Earth Orbit,GEO)卫星、11颗倾斜地球同步轨道(Inclined Geosynchronous Orbit,IGSO)卫星和30颗中地球轨道(Medium Earth Orbit,MEO)卫星组成,包含北斗二号全球卫星导航系统(BDS-2)和北斗三号全球卫星导航系统(BDS-3)2种卫星混合星座。相较于BDS-2,BDS-3不仅向下兼容了B1I和B3I信号,还增加了B1C、B2a、B2b和B2(B2a+B2b)信号[1-3]。更多的可用卫星和更加丰富的信号资源使得BDS在导航定位中发挥着更加重要的作用。

差分码偏差(differential code bias,DCB)是指同一时刻不同频率或同一频率上的不同测距信号在发射链路和接收链路中所产生时间延迟的差值[4]。DCB在数量级上可以达到几纳秒甚至几十纳秒,是影响电离层总电子含量(total electron content,TEC)监测和建模的主要误差源。若忽略DCB的影响,不仅会给TEC的计算带来9～30 TECU(TECU为电离层TEC的单位,1 TECU=1×1016个电子/m2)的偏差,也会使得依靠伪距进行定位的测量结果产生数米偏差[5-6]。常用DCB估计方法有2种:① 与电离层参数同步估计[7-8];② 先用电离层产品直接扣除电离层TEC,再解算DCB[9]。目前,多全球卫星导航系统(global navigation satellite system,GNSS)试验跟踪网(Multi-GNSS Experiment,MGEX)向全球用户提供2款北斗DCB产品,分别是德国宇航中心(Deutsches Zentrum für Luft-und Raumfahrt,DLR)提供的DLR DCB产品和中国科学院(Chinese Academy of Sciences,CAS)提供的CAS DCB产品。文献[10]对MGEX发布的北斗DCB产品的精度与稳定性进行质量分析,并初步计算分析了北斗DCB的随机误差特性;文献[11]分析不同太阳活动水平下BDS卫星DCB产品的稳定性变化特性,并建模实现不同太阳活动下BDS卫星DCB短期预报;文献[12-13]首先明确BDS时间群延迟(timing group delay,TGD)和DCB之间的关系,然后推导了适用多种场合的北斗DCB改正模型,并论证在高精度定位中采用合适DCB改正模型的重要性;文献[14]在文献[12]的基础上推导三频无电离层组合DCB改正模型,并用实验分析验证DCB改正对标准单点定位(standard point positioning,SPP)和精密单点定位2种定位模式的影响。

以上关于北斗DCB稳定性变化特性和DCB改正定位模型的分析,都基于DCB为1 d中的常量参数,关于BDS卫星DCB短时变化特性的研究很少。本文分析BDS卫星DCB在1 d中的短时变化,针对单历元数据进行短时DCB估计,即在每个观测历元估计出该时刻对应的DCB估值。首先采用最小二乘和Tikhonov正则化方法,解算BDS卫星在各个历元时刻下的DCB估值,然后对得到的各卫星DCB序列的精度、稳定性进行分析,并采用谱分析方法确定短时DCB周期性,构建模型拟合短时DCB变化,最后通过实验来探讨短时DCB序列改正对SPP的影响。

1 短时DCB估计方法

通过双频无几何组合,可得到伪距和载波相位无几何观测量,表达式为:

(1)

(2)

由于伪距无几何观测量观测噪声较大,常用载波相位对其进行平滑处理,平滑后的伪距观测值可表示为:

(3)

其中:P4,sm为平滑后的伪距无几何观测量;P4,prd(t)为伪距和载波的组合观测量;ωt为与第t个历元相关的比例系数,ωt=1/t。当t=1时,P4,sm=P4。

电离层延迟通常忽略高阶项影响,采用一阶表达式表示,即

(4)

其中:fj为信号频率,j=1,2;Stec为倾斜TEC。将(4)式代入(3)式,可得:

(5)

由(5)式可得Stec的表达式为:

(6)

为方便电离层建模,常借助一个映射函数FM将倾斜方向Stec转到天顶方向Vtec,表达式为:

(7)

其中:Vtec为天顶方向TEC;z为卫星高度角的余角;R为地球平均半径;H为电离层薄层高度,本文取值为506.7 km;α为比例因子,本文取值为0.978 2。将(7)式代入(6)式,可得:

(P4,sm-cDs-cDr)

(8)

本文采用球谐函数模型对天顶方向Vtec进行建模,表达式为:

[anmcos(ms)+bnmsin(ms)]

(9)

将(8)式、(9)式联立可得:

(10)

(10)式改写成矩阵形式为:

(11)

在解算短时DCB过程中,邻近测站对应穿刺点纬度和日固经度变化较小,系数矩阵中数据存在一定的相关性,在某些历元中直接采用最小二乘法估算会出现秩亏现象。为了减弱或消除病态问题对秩亏历元的影响,本文采用Tikhonov正则化[15]对秩亏历元的参数进行解算。Tikhonov正则化方法是对最小二乘法的改进,正则化最小二乘代价函数为:

(12)

(12)式的解为:

(13)

其中,λ为正则化参数,λ≥0,当λ=0时该解为最小二乘解。为了选取合适的λ,本文采用岭迹法确定秩亏历元中λ参数,对随机选取的50个秩亏历元绘制岭迹图,当λ处于0.001附近时,岭迹趋于平稳,故本文解算过程中λ取0.001。

2 数据处理与实验分析

本文选取2021年年积日第152天至第212天国际GNSS服务(International GNSS Service,IGS)和MGEX测站网BDS的C2I、C6I观测数据,测站分布位置如图1所示。采用4阶球谐函数模型对TEC进行建模,卫星截止高度角设置为10°。为了分析DCB短时变化,同步估计每个观测历元时刻下的DCB参数和球谐函数参数。针对可能出现的极少数历元缺失情况,采用前后2个历元平均值进行替代。

图1 测站位置分布图

以年积日第200天为例,按照本文方法解算出的部分卫星DCB估值序列如图2所示。从图2可以看出:BDS-2中C02、C04、C06、C08卫星在各个历元的估值比较稳定,仅有轻微的波动,1 d中上下波动在0.5 ns左右;BDS-3中C26、C32、C34、C40卫星在某些时段会出现较大的起伏,稳定性与BDS-2中的卫星相比较差,1 d中上下波动在1.2 ns左右。

图2 部分卫星短时DCB序列

2.1 BDS卫星短时DCB精度与稳定性分析

为了验证BDS卫星短时DCB估计精度,将CAS发布的DCB产品作为参考值,以本文方法得到的1 d中各卫星短时DCB序列分别与相应的参考值作差,可获得各卫星DCB估值与参考值的差值序列,对各卫星的差值序列求取平均值,可得到各卫星短时DCB序列在1 d中的精度状况。考虑到仅对比1 d数据具有偶然性,对2021年年积日第152天至第212天共60 d天数据进行绝对值平均。目前,BDS-2卫星有15颗,伪随机噪声(pseudo-random noise,PRN)码编号为C01～C16(C15编号未用);BDS-3卫星PRN码编号为C19～C46,其中C31卫星处于试验阶段,未参与解算。各卫星在60 d中的平均偏差如图3所示。

图3 BDS卫星短时DCB(C2I、C6I)平均偏差

从图3可以看出,BDS-2卫星DCB短时估计的精度小于BDS-3卫星。BDS-2卫星短时DCB估值与参考值在60 d中的平均偏差为0.34 ns,且各卫星平均偏差之间相差较大,C04卫星与参考值的平均偏差(0.79 ns)最大,C16卫星与参考值的平均偏差(0.11 ns)最小。BDS-3卫星短时DCB在60 d 中的平均偏差为0.14 ns,各卫星短时DCB平均偏差相差较小,C38卫星平均偏差(0.35 ns)最大, C34卫星平均偏差(0.08 ns)最小,其他卫星平均偏差大多在0.10 ns左右。

为了分析卫星短时DCB在日内变化的稳定性,本文统计各观测日1 d内卫星短时DCB的标准差(standard deviation,STD)σ,并将其作为反映卫星短时DCB稳定性的指标[9]。σ表达式为:

(14)

图4 BDS卫星平均日内变化稳定性

BDS-2中GEO卫星包括C01～C05卫星,BDS-2中非GEO卫星包括C06～C16卫星(C15编号未用),BDS-3中非GEO卫星为C19～C46卫星。从图4可以看出:BDS-2和BDS-3中的非GEO卫星(MEO卫星和IGSO卫星)的平均日内稳定性非常接近,其STD均值分别为0.42、0.43 ns；BDS-2中的GEO卫星日内变化稳定性更高,其STD均值为0.24 ns。

2.2 短时DCB周期性分析与拟合建模

通过对年积日第152天至第212天各卫星DCB短时序列进行谱分析,可发现DCB短时估值的周期性变化。本文使用快速傅里叶变换对DCB序列进行分析,得到部分卫星在各观测日中的频谱图,如图5所示。

从图5a可以看出,BDS的GEO卫星DCB估值在频率为2.035×10-6、1.170×10-5、2.289×10-5Hz会达到一个极大值,此时三者对应的周期分别为5.6 d、23.7 h、12.1 h,即GEO卫星会表现出以5.6 d、23.7 h、12.1 h为周期的周期性变化。

从图5b可以看出:IGSO卫星在频率为1.157×10-5、2.314×10-5、3.471×10-5Hz或4.641×10-5Hz会达到一个极大值；在频率为2.314×10-5Hz时振幅都会达到最大,在1.157×10-5Hz时振幅次之,上述2个频率对应的周期分别为12.0、24.0 h;个别卫星在频率为3.471×10-5Hz 或4.641×10-5Hz时也会达到一个相对较小的峰值,上述2个频率对应的周期分别为8.0、6.0 h。这说明IGSO卫星除了会有以12.0、24.0 h为周期的强周期性变化,也会在某卫星上表现出以6.0、8.0 h为周期的弱周期性变化。

从图5c可以看出:MEO卫星在频率为9.918×10-6、2.149×10-5、3.319×10-5Hz或4.311×10-5Hz会达到一个极大值；在频率为9.918×10-6、4.311×10-5Hz时振幅较为明显,上述2个频率对应的周期分别为28.0、6.4 h;个别卫星在频率为2.149×10-5Hz或3.319×10-5Hz时也会达到一个相对较小的峰值,此时对应的周期分别为12.9、8.4 h。因此,MEO卫星有以6.4、28.0 h为周期的强周期性变化,也会出现以12.9、8.4 h为周期的弱周期性变化。

图5 BDS部分卫星短时DCB序列频谱

由上述周期性分析可知,卫星短时DCB在1 d内存在周期性变化,考虑到此特性,可将短时DCB的拟合模型表示为一个二次函数与几个周期函数的总和。

表达式为:

y=A0T2+B0T+C0+

(15)

其中:y为DCB值;T为相对起始历元的时间间隔;A0、B0、C0分别为二次函数的二次项系数、一次项系数和常数项;Dk、Ek为周期函数在第k个周期的系数;Qk为周期函数的周期。为了得到每日短时DCB的拟合函数,首先对每日卫星DCB短时估值进行谱分析,确定周期函数的周期后,采用最小二乘法确定周期函数和二次函数的系数。

以年积日第200天为例,部分卫星DCB拟合函数中系数与周期函数周期分别见表1、表2所列。

表1 拟合模型中二次函数对应系数

表2 拟合模型中周期函数对应系数与周期

为了评估模型的拟合精度,本文将拟合模型得到的拟合值与解算出的DCB估值进行比较。部分卫星短时DCB采用(15)式的拟合值与解算值对比如图6所示。拟合值与解算值差值的相关统计特性见表3所列。表3中,RMS表示差值的均方根(root mean square)值。

从图6可以看出,各卫星DCB拟合值与DCB解算值具有很好的一致性。

从表3可以看出,各卫星DCB拟合值与解算值的平均差值几乎为0 ns,所选6颗卫星的RMS值分别为0.006、0.017、0.053、0.006、0.018、0.047 ns,表明本文提出的拟合模型具有较高的精度,可以很好地拟合卫星DCB在日内的短时变化,也从侧面反映了采用谱分析方法分析周期性的可行性。

图6 短时卫星DCB拟合结果

表3 解算值与拟合值之间差值统计特性单位:ns

3 DCB改正对SPP的影响

为了分析短时DCB估计对SPP的影响,选用IGS测站网下ABPO、ULAB、LEIJ和WUH2测站在2021年年积日第200天的数据进行SPP实验。对BDS的B1频点观测量分别采用3种不同的处理策略,即不改正DCB、CAS产品改正、本文方法改正,分析单频SPP结果。DCB改正的单频SPP模型见文献[12]。4个测站在3种策略下单频SPP三维(3D)残差结果如图7所示。

从图7可以看出:

(1) 经过DCB改正后,单频SPP的定位精度提升明显,基本在5 m以内。

(2) 使用CAS产品与本文方法对单频SPP改正,3D残差的量级和历元间变化趋势大致相当。在不同测站上,两者定位精度互有优劣,互差在cm级,改正效果大致相当。

4个测站单频SPP在不同方向上的定位结果见表4所列。

由表4可知:

(1) CAS产品改正和本文方法改正对SPP的影响为m级,CAS产品改正的平面方向、U方向和3D方向平均RMS分别为1.118、3.162、3.407 m,本文方法改正的平面方向、U方向和3D方向平均RMS分别为1.112、3.143、3.387 m。可以看出,两者改正后SPP定位精度互差在cm级,在平面方向、U方向和3D方向,本文方法改正与CAS产品改正相比,略微提升0.6、1.9、2.0 cm,其原因可能是短时DCB参数更新周期短,能更好反映DCB在各个历元的真实状况。

(2) CAS产品改正的SPP在E、N、U和3D方向定位精度平均提升55.5%、58.6%、42.5%、45.0%；本文方法改正的SPP在E、N、U和3D方向定位精度平均提升56.9%、58.2%、43.0%、45.3%。

图7 4个测站DCB改正前、后SPP定位3D残差对比

表4 4个测站3种DCB处理策略下SPP定位平均RMS值

4 结 论

通过在单历元间估算DCB,得到卫星DCB短时序列,可以对DCB短时序列变化特性进行分析,也能充分考虑DCB短时变化改正对导航定位的影响。

本文在各个历元间采用最小二乘和Tikhonov正则化方法同步估算电离层参数和BDS的DCB短时序列,分析卫星DCB短时精度、稳定性和周期性,最后通过SPP实验分析DCB短时序列对定位的影响。结果表明：施加DCB短时改正和施加CAS DCB产品都能显著提高SPP定位精度,定位精度提升都在40%以上;施加短时DCB改正和CAS产品改正,在不同测站上SPP定位精度各有优劣,两者差异在cm级,改正效果大致相当。