近几年来GNSS系统发展迅速，多系统多频率是将来的发展趋势^[1-2]。BDS系统于2012年正式投入运行，服务范围覆盖亚太地区，在2015年中期开始了第三阶段的建设，预计将于2020年完成全球组网，最终形成5颗同步卫星和30颗非同步卫星的星座构成^[3-4]。GPS系统正在实施现代化计划，主要内容包括使用新的民用信号、发射新型号的卫星以及对地面控制部分进行更新，目前已有12颗BLOCK-IIF卫星支持L5频率信号。

随着多个全球和区域卫星导航系统的建成，多系统PPP能够有效提高定位的精确性、可用性和可靠性^[5-8]。然而不同系统的信号在同一台接收机中的硬件延迟并不相同，导致了硬件延迟偏差的出现^[9-13]。系统间偏差ISB包括了系统时间偏差和硬件延迟偏差^[14]。在进行多系统PPP数据处理时必须要考虑ISB的影响^[15]。ISB建模有助于深入了解ISB的误差特性，对多系统PPP研究具有一定的参考价值。通过对ISB进行精确建模和预报，采用ISB预报值进行先验约束能够提高多系统PPP的定位精度。因此有必要对ISB建模和预报方法进行研究。

文献[16]指出系统间偏差与接收机内部的信号时延有关。文献[17]发现在多系统PPP中Galileo、BDS与GPS的ISB同接收机类型相关，不同类型接收机之间GPS/BDS的ISB差异最大超过了100 ns。文献[18]对GPS/Galileo紧组合当中的差分ISB进行了研究，利用双差组合分析了差分ISB的稳定性及其对定位结果的影响，试验结果表明差分ISB的长期稳定性较好，大小与接收机类型有关。文献[19]研究了IRNSS相对GPS、Galileo、QZSS在L5/E5a频率上的差分ISB及其应用。文献[20]采用谱分析和LS拟合的方法对ISB进行了短期预报，取得了良好的成果。需要注意的是，参考文献[1, 18, 19]中提到的系统间偏差为差分ISB，与PPP中的ISB在定义上存在本质的差别。

文献[20]中LS方法对拟合数据是等权处理的，而在实际情况中随着时间的推移早期数据与预报时刻的时空相关性会逐渐减弱，近期数据对ISB预报更有价值。因此本文提出一种改进的ISB建模方法，采用Kalman滤波估计模型参数，并根据时空相关性对ISB拟合数据的方差进行调整。试验结果表明该方法可以有效提高ISB的预报精度。

1 ISB估计和建模 1.1 ISB估计方法

采用消电离层组合消去电离层延迟误差一阶项，PPP的观测方程为^[21-22]

式中，上标G和C分别表示GPS和BDS系统；P和L分别表示伪距观测量和载波相位观测量，单位为m；ρ^S为在信号发射和接收时刻卫星与接收机天线相位中心的几何距离；c为光在真空中的速度；dt_r为接收机钟差；dt^S为卫星钟差；τ^C为BDS系统时相对GPS系统时的时间偏差；d_rS、d^S和b_rS、b^S分别为接收机、卫星的码硬件延迟和相位硬件延迟；T^S为对流层延迟误差；N^S为相位模糊度，单位为m；ε_P^S和ε_L^S分别为伪距和载波相位观测误差包括观测噪声和多路径误差。

采用精密轨道和钟差产品后，消去了卫星钟差和卫星码硬件延迟，观测方程为

对式(2) 进行整理得到

式中，和分别为吸收了偏差后的接收机钟差和相位模糊度；ISB为BDS/GPS系统间偏差。对流层干延迟部分采用经验模型进行改正，湿延迟部分作为未知参数进行估计。相位缠绕、天线相位中心和固体潮等误差根据相应模型进行改正。具体处理策略见表 1。

1.2 ISB建模方法

分析中心在计算精密钟差产品时，天与天之间基准钟的调整会导致卫星钟差发生跳变，并引起ISB序列的跳变。该跳变会影响ISB模型的拟合程度，因此在进行建模之前需要对拟合数据进行预处理。如果一段ISB序列超出了正常的阈值范围，则判断此处发生了跳变，并将该段ISB序列减去跳变的幅度。跳变修复的公式为

式中，ISB_(t)为t时刻的ISB值；t_i和t_j分别为ISB序列超出阈值的起始和结束时刻；t_i－1为ISB序列超出阈值起始时刻的前一历元；mean和std分别为ISB拟合数据的平均值和标准差；d为跳变的幅度。

ISB拟合数据可以看作是由趋势项和周期项两部分组成的。首先利用FFT(fast fourier transform)分析ISB序列并确定拟合数据的周期，然后分别使用二次曲线和三角函数对拟合数据的趋势项和周期项进行拟合，建模的数学模型为

式中，t=t_i－t₀，t₀为起始历元，t_i为第i个历元；n为周期项的个数；P_i为已知的周期；A、B、C、D_i和E_i为待估参数。

在实际情况中，随着时间的推移早期数据与预报时刻的时空相关性会逐渐减弱，近期数据对ISB预报更有价值。因此提出一种改进的ISB建模方法，采用Kalman滤波方法估计模型参数，并根据时空相关性对ISB拟合数据的方差进行调整。Kalman滤波的增益矩阵为^[22-23]

式中，下标t为t时刻；m为待估参数的个数；K_t为增益矩阵；为模型参数预报向量的协方差；A_t为设计矩阵；R_t为ISB拟合值的方差。

为了增加靠近预测时刻的拟合数据对模型参数的贡献，用ISB序列末端多个拟合数据的平均值作为确定拟合数据方差的标准，考虑到在ISB波动较大的情况下采用多天数据的平均值能够增强预测的稳定性，因此选择使用拟合数据最后3 d的平均值为标准。确定ISB拟合数据方差的具体方法为

式中，为7 d ISB拟合数据中最后3 d的平均值。

2 试验分析 2.1 ISB建模和预报试验

考虑到BDS和GPS系统均要有足够的可用卫星，选择MGEX(multi-GNSS experiment)网中观测条件较好的CUT0、MAYG、MRO1和NNOR 4个测站2015年年积日87—93天共7 d的数据进行建模。首先利用PPP估计ISB序列，在修正ISB跳变后分别采用Kalman滤波和等权LS方法进行ISB建模，根据所建模型预报年积日94天的ISB值，最后利用当天的ISB估计值来验证预报结果的精度，并对Kalman滤波和等权LS方法进行了对比分析。为了验证ISB序列跳变修正的效果，ISB跳变修正前和修正后均采用等权LS方法进行建模，并对ISB模型的拟合程度进行了比较。表 2给出了测站在E、N、U方向的定位误差。4个测站7 d静态PPP的定位精度在东向和高程方向均为厘米级，北向定位精度为毫米级，较高的定位精度保证了ISB的估计质量。图 1给出了跳变修正前后的ISB估计值和模型值。图 2给出了年积日88—89天精密卫星钟差的历元差。

由图 1可以看出，ISB跳变发生在年积日第88—89天和第91—92天，跳变是由于分析中心在计算精密钟差产品时天与天之间基准钟的调整导致的。以年积日第88—89天为例，将年积日第88天23:50至第89天00:05共4个历元的精密卫星钟差依次做差。由图 2可知，第88天23:50至23:55之差和第89天00:00至00:05之差接近于0 ns，而第88天23:55和第89天00:00之差约为-20 ns，发生了很大的跳变。

根据图 1可知，在跳变修正之前，跳变的幅度干扰了LS拟合模型对ISB序列振幅的拟合，导致拟合模型的振幅过大，不利于对ISB的预报。在修正跳变之后，效果最为明显的是NNOR站，LS拟合模型的振幅变小，周期性也更加明显，模型的拟合精度得到了提高。CUT0、MRO1和NNOR 3个站，跳变修正前模型与原始ISB序列的拟合精度分别为0.62 ns、0.63 ns、0.60 ns，跳变修正后模型与修正ISB序列的拟合精度分别为0.28 ns、0.58 ns、044 ns，拟合精度均有较大提高。

为了对Kalman滤波和等权LS方法进行比较，试验均采用跳变修正后的ISB数据，其中MAYG站ISB序列没有超过跳变阈值，未作跳变修正。使用LS方法估计模型参数时，对ISB拟合数据采用等权处理的方式。采用Kalman滤波方法估计模型参数时，以拟合数据最后3 d的ISB平均值为标准确定拟合数据的方差，如式(8) 所示。在确定模型参数后利用ISB模型预报年积日94天的ISB值，并用当天的ISB估计值对预报结果进行验证。Kalman滤波方法的ISB模型参数和周期项见表 3。此外，对ISB模型参数估值进行了显著性检验。给定显著水平α=0.05，自由度为326，对应的检验临界值为t_α/2=1.967。周期P_i采用FFT方法确定为已知量，对A、B、C、D_i、E_i等参数构造t统计量。参数估值的显著性检验结果见表 4。

随着时间的推移，早期数据与预报时刻的时空相关性会逐渐减弱，近期数据对ISB预报更有价值。Kalman滤波能够根据拟合数据的方差调整增益矩阵，而等权LS模型无法利用该特点。图 3给出了Kalman滤波和等权LS方法估计的ISB模型及其ISB预报值。由图 3可知，使用Kalman滤波方法估计模型参数时，由于最后3 d拟合数据的方差较小，第5—7天的数据对模型参数的贡献更大，避免了模型参数特别是振幅参数受到早期数据大幅波动的影响，使得拟合曲线在时序末端与ISB序列更吻合。

由于Kalman滤波方法估计的ISB模型更侧重于ISB时序末端(第5—7天)的拟合程度，因此比较两种ISB模型在最后3天的拟合精度比较客观。表 5给出了Kalman滤波和等权LS方法ISB模型的预报精度和最后3 d的建模精度。由表 5可知，CUT0、MAYG、MRO1、NNOR 4个站使用LS方法第5—7天的拟合精度分别为0.22 ns、2.25 ns、0.31 ns、0.24 ns，使用Kalman滤波方法的拟合精度分别为0.22 ns、2.34 ns、0.24 ns、0.22 ns。只有MAYG站Kalman滤波方法的拟合精度有所降低，从图 3中可以看出，这是MAYG站第5天ISB波动较大所导致的。因为本文拟合数据方差的确定方法能够降低ISB模型对于大幅波动的敏感性，因此MAYG站在第5天的拟合精度较低，而在第6天和第7天Kalman方法的拟合程度依然较好。

从ISB预报精度来看，由表 5可知CUT0、MAYG、MRO1、NNOR 4个站等权LS模型的预报精度分别为0.37 ns、0.87 ns、0.23 ns、1.00 ns，Kalman滤波模型的预报精度分别为0.26 ns、0.77 ns、0.13 ns、0.68 ns，预报精度分别提高了29.7%、11.5%、43.5%、32.0%。图 4给出了Kalman滤波和等权LS方法ISB模型的预报精度。

2.2 PPP定位精度分析

为了验证ISB预报值的可用性，采用第8天的数据进行静态PPP试验。分别将Kalman滤波和等权LS两种模型的ISB预报值作为先验约束，当ISB模型第5—7天的拟合精度小于1 ns时，先验精度取1 ns，当拟合精度大于1 ns时，先验精度取拟合精度。图 5给出了不加ISB先验约束和使用ISB先验约束PPP的定位误差，表 6给出了相应的定位误差统计结果。

根据表 6可知CUT0、MAYG、MRO1、NNOR 4个站静态PPP在E、N、U方向的定位误差RMS值，不加ISB先验约束时为1.51 cm、1.14 cm、2.38 cm，使用等权LS模型ISB先验约束时为1.39 cm、0.84 cm、2.32 cm，分别提高了7.9%、26.3%、2.5%，使用Kalman滤波模型ISB先验约束时为1.35 cm、0.84 cm、2.30 cm，分别提高了10.6%、26.3%、3.4%。与等权LS方法相比，Kalman滤波方法在E和U方向上分别多提高了2.7%和0.9%，在N方向上表现相同。

3 结论

针对等权LS方法估计ISB模型参数时忽略了拟合数据时空相关性的问题，提出采用Kalman滤波方法估计模型参数，并根据时空相关性调整拟合数据的方差。分别使用上述两种方法进行ISB建模和预报，对ISB预报精度和附加ISB先验约束的PPP定位结果进行了验证分析。试验结果表明，4个测站Kalman拟合模型的ISB预报精度比等权LS模型分别提高了29.7%、11.5%、43.5%、32.0%；采用Kalman拟合模型ISB先验约束的定位精度与不加先验约束相比在E、N、U方向分别提高了10.6%、26.3%、3.4%，比使用等权LS模型ISB先验约束在E和U方向上分别多提高了2.7%和0.9%，在N方向上表现相同。