BDS是世界上为数不多的使用GEO卫星参与导航的卫星导航系统之一。BDS-1由2颗GEO卫星组成，BDS-2拥有5颗GEO卫星，BDS-3系统拥有3颗GEO卫星。我国GEO卫星基本集中在亚太地区上空，可以进行全天候观测，具有高轨特性，观测效果较为理想。GEO卫星在北斗系统的短报文通信^[1]、实时精密单点定位^[2]、星间链路组建^[3-4]中都起到重要作用，可有效提高北斗系统定位、定轨、授时精度。目前BDS-3卫星已布置完毕。

北斗全星座均可发射多频信号，研究表明，多频信号之间存在系统性偏差，不同观测组合解算的卫星钟差之间也存在差异，称为卫星钟差频间偏差IFCB，IFCB可以对不同频率组合得到的钟差进行改正^[5-7]。目前采用多频无电离层延迟组合的差值计算IFCB，其中相位观测用于变化部分的求解^[8-9]，伪距观测用于常数部分的计算^[10-11]。当未发生周跳时，可以不间断地进行GEO卫星观测，从而为BDS系统的IFCB特性研究提供良好条件。由于BDS-2卫星对应的IFCB具有明显的时变性与周期性^[12]，因此本文分析BDS-3 GEO卫星的IFCB特性，并对BDS-2、BDS-3的IFCB进行对比分析，构建其高精度模型化函数。

1 IFCB解算方法 1.1 相位观测值解算IFCB

将多频BDS观测中的2组无电离层延迟相位组合相减可以得到：

$ \begin{array}{c} \operatorname{DIF}\left(B_1, B_2, B_3\right)=B_1 / B_3-B_1 / B_2= \\ \varphi^{\mathrm{s}}+\text { const }_{1, 3}-\text {const}_{1, 2}+\varepsilon_{1, 3}-\varepsilon_{1, 2} \end{array} $

(1)

式中，φ^s为卫星IFCB的变化部分，const_{1, 3}、const_{1, 2}分别为无电离层相位B₁/B₃、B₁/B₂组合对应的相位模糊度和硬件偏差的常数值，ε_{1, 3}、ε_{1, 2}分别为B₁/B₃、B₁/B₂对应的噪声。对于序号为r的测站，当未发生周跳时，对相邻历元的无电离层延迟相位组合差值进行差分，即对k历元与(k+1)历元数值进行差分，可以得到IFCB相邻历元的变化值:

$ \begin{array}{c} \Delta \varphi^{\mathrm{s}}(k)_r=\operatorname{DIF}\left(B_1, B_2, B_3\right)_r(k+1)- \\ \operatorname{DIF}\left(B_1, B_2, B_3\right)_r k \\ \end{array} $

(2)

若观测站网络中有n个测站，其均值可写为:

$ \Delta \varphi^{\mathrm{s}}(k)=\frac{1}{n} \cdot \sum\limits_{r=1}^n \Delta \varphi^{\mathrm{s}}(k)_r $

(3)

通过该方法可以计算得到每个历元中IFCB变化部分φ^s的变化量。由于IFCB剩余的常数部分为模糊度、硬件偏差等参数构成的常数值，因此该变化量也应考虑在内。

假设初始历元的IFCB^s为IFCB₀^s，m历元的IFCB变化量为Δφ^s(k)_m，则i历元的IFCB可表示为：

$ \mathrm{IFCB}_i^\mathrm{s}=\mathrm{IFCB}_0^\mathrm{s}+\sum\limits_{m=1}^i \Delta \varphi^{\mathrm{s}}(k)_m $

(4)

通过变化量可以推算出每个历元的IFCB值，从而生成时间序列。

1.2 伪距观测值解算IFCB

将伪距无电离层延迟观测组合相减可得:

$ \begin{array}{c} \operatorname{DIF}\left(C_1, C_2, C_3\right)=C_1 / C_3-C_1 / C_2= \\ \operatorname{IFCB}^{\mathrm{r}}-\left(\varphi^{\mathrm{s}}+\operatorname{IFCB}^{\mathrm{s}}\right)+\omega_{1, 3}-\omega_{1, 2} \end{array} $

(5)

式(5)也可改写为：

$ \begin{array}{c} \operatorname{DIF}\left(C_1, C_2, C_3\right)=C_1 / C_3-C_1 / C_2= \\ \left(a_3 b_1^{\mathrm{r}}-a_4 b_3^{\mathrm{r}}\right)-\left(a_3 b_1^{\mathrm{s}}-a_4 b_3^{\mathrm{s}}\right)-\left(a_1 b_1^{\mathrm{r}}-\right. \\ \left.a_2 b_2^{\mathrm{r}}\right)+\left(a_1 b_1^{\mathrm{s}}-a_2 b_2^{\mathrm{s}}\right)+\omega_{1, 3}-\omega_{1, 2}= \\ a_2 \mathrm{DCB}^{\mathrm{s}}\left(C_1-C_2\right)-a_4 \mathrm{DCB}^{\mathrm{s}}\left(C_1-C_3\right)+ \\ a_4 \mathrm{DCB}_{\mathrm{r}}\left(C_1-C_3\right)-a_2 \mathrm{DCB}_{\mathrm{r}}\left(C_1-C_2\right)+ \\ \omega_{1, 3}-\omega_{1, 2} \end{array} $

(6)

式中，DCB_r和DCB^s分别为接收机端和卫星端的差分码偏差DCB，b_i^r与b_i^s分别为接收机和卫星在i(i=1, 2, 3)频段内的硬件误差。设f₁、f₂、f₃分别为B₁、B₂、B₃频点对应的频率，则有：

$ \left\{\begin{array}{l} a_1=\frac{f_1^2}{f_1^2-f_2^2}, a_2=\frac{f_2^2}{f_1^2-f_2^2} \\ a_3=\frac{f_1^2}{f_1^2-f_3^2}, a_4=\frac{f_3^2}{f_1^2-f_3^2} \end{array}\right. $

(7)

当采用式(6)进行偏差估计时，接收机、卫星端对应的偏差无法分离，因此引入约束:

$ \sum\limits_{j=1}^d\left(b_{1, 2, 3}^{\mathrm{r}}\right)_j=0 $

(8)

式中，b_{1, 2, 3}^r)_j为测站j对应的偏差，d为测站总数。线性化后可得:

$\begin{array}{c} \left[\begin{array}{cc} \boldsymbol{e}_{1 \times d} & \bf{0}_{1 \times t} \\ \boldsymbol{I}_{d \times d} \bigotimes \boldsymbol{e}_{t \times 1} & \boldsymbol{e}_{d \times 1} \bigotimes \boldsymbol{I}_{t \times t} \end{array}\right]\left[\begin{array}{c} \left(b_{1, 2, 3}^{\mathrm{r}}\right)_1 \\ \cdots \\ \left(b_{1, 2, 3}^{\mathrm{r}}\right)_d \\ \left(b_{1, 2, 3}^{\mathrm{s}}\right)_1 \\ \cdots \\ \left(b_{1, 2, 3}^{\mathrm{s}}\right)_t \end{array}\right]=\\ \left[\begin{array}{c} 0 \\ \left(\operatorname{DIF}\left(C_1, C_2, C_3\right)\right)_{1, 1} \\ \cdots \\ \left(\operatorname{DIF}\left(C_1, C_2, C_3\right)\right)_{1, t} \\ \cdots \\ \left(\operatorname{DIF}\left(C_1, C_2, C_3\right)\right)_{d, t} \end{array}\right] \end{array} $

(9)

式中，(b_{1, 2, 3}^s)_t为卫星t对应的偏差，I_d×d与I_t×t分别为单位矩阵，e_1×d为每个元素均为1的向量，0_1×t为每个元素均为0的向量，$\otimes $为克罗内克积算子，t为卫星个数。在式(9)的基础上采用最小二乘方法进行卫星、接收机端伪距IFCB估计。由于添加了约束，解算的IFCB会包含基准偏差，由此需采用星间差分方法进行消除。选择卫星l作为卫星k的参考星，则其IFCB为：

$ \Delta\left(b_{1, 2, 3}^{\mathrm{s}}\right)_{k, l}=\left(b_{1, 2, 3}^{\mathrm{s}}\right)_k-\left(b_{1, 2, 3}^{\mathrm{s}}\right)_l $

(10)

对每个历元执行上述计算，得到Δ(b_{1, 2, 3}^s)_{k, l}值作为该历元卫星k的IFCB结果。

2 数据处理与分析

BDS-2共包含5颗GEO卫星，PRN分别为C01、C02、C03、C04、C05；BDS-3共包含3颗GEO卫星，PRN分别为C59、C60、C61。截至2022-01-15，C61卫星仍处于在轨测试阶段，IGS测站无观测数据，因此分析对象仅为C01、C02、C03、C04、C05(BDS-2)与C59、C60(BDS-3)7颗GEO卫星(图 1)。由于卫星数量较多，仅展示部分卫星的分析结果。利用BDS-2和BDS-3分别计算B1、B2、B3与B1、B2b、B3频点组合的IFCB。

为分析BDS-2、BDS-3系统GEO卫星的IFCB结果，采用186个IGS站对2022年doy1~15的观测数据进行处理。根据解算的IFCB结果分析对比BDS-2、BDS-3 GEO卫星IFCB的异同，并研究其变化特性。在伪距观测结果分析中，BDS-2系统采用C01卫星作为参考星；BDS-3系统采用C59卫星作为参考星，以此消除基准对结果的影响。

IFCB采用基于谐波函数的混合函数和多项式进行模型化处理。基于谐波函数的混合函数为：

$ \begin{array}{l} \operatorname{IFCB}_i^{\mathrm{s}}=a+e \cdot t+ \\ \sum\limits_{x=1}^n y_x \sin \left(\frac{2 \pi}{T_x} \cdot t+\theta_x\right) \end{array} $

(11)

式中，a为常数值，e为线性项系数，x为谐波函数序号，T_x为周期值，θ_x为初始相位, y_x为振幅。n次多项式为：

$ \operatorname{IFCB}_i^\mathrm{s}=\sum\limits_{z=1}^n a_z \cdot t^z $

(12)

式中，a_z为z次项对应的系数。

2.1 IFCB结果

采用相位观测数据进行IFCB解算，结果如图 2所示。由图可见，BDS-2系统GEO卫星的相位IFCB变化部分具有明显的时变性与周期性，BDS-3系统GEO卫星的相位IFCB变化幅度约为mm级，周期性与时变性不明显。对每颗卫星1 d内的IFCB进行快速傅里叶变换(fast Fourier transform, FFT)，得到每颗卫星的8个周期值，结果如图 3所示。

由图 3可知，出现频率最高的8个周期值依次为：12 h、8 h、6 h、4.8 h、4 h、3.428 6 h、3 h、2.666 7 h，说明BDS-3 GEO卫星的相位IFCB变化存在周期性。选取BDS-2、BDS-3参考星后，其对应的伪距IFCB结果如图 4所示(BDS-2结果中doy10观测数据质量较差，不作展示)。

由图 4可见，北斗GEO卫星伪距IFCB围绕某一固定值变化，不同卫星的变化范围存在差异。若忽略噪声影响，则BDS-2 GEO卫星伪距IFCB具有明显的时变性与周期性，BDS-3 GEO卫星伪距IFCB的时变性与周期性不明显。使用FFT对结果进行周期分析，部分周期结果如图 5所示。

由图 5可知，出现频率最高的8个周期值依次为：12 h、8 h、6 h、4.8 h、4 h、3.428 6 h、3 h、2.666 7 h，与相位IFCB的结果一致，说明BDS-2与BDS-3 GEO卫星的IFCB变化具有显著的周期性，且BDS-3的周期值稳定度有所下降，不利于伪距IFCB的长期拟合与预报。

2.2 模型化

针对BDS-2 GEO卫星，分别采用多项式、谐波函数混合函数对§2.1中FFT变换得到的8个周期值进行相位、伪距IFCB结果模型化分析。其中，2022年doy1相位、伪距对应的IFCB以及部分模型化对应结果分别如图 6、7所示。

从图 6、图 7可以看出，使用3次及以上的多项式和4阶、8阶谐波函数混合函数均可以较好地对相位IFCB进行模型化处理；使用4阶、8阶谐波函数混合函数对伪距IFCB进行模型化的效果优于多项式函数。上述拟合效果的差异与IFCB曲线的极值点数量有关。

为定量比较模型化效果，以d为单位计算模型化函数拟合的均方根误差(RMS)和平均值，结果见表 1(单位mm)和表 2(单位m)。

由表 1、表 2可以看出，随着阶数或次数的增加，谐波函数混合函数和多项式函数拟合精度均有所提升，但多项式函数拟合精度始终低于谐波函数混合函数。使用8阶谐波函数混合函数的拟合效果最好：相位IFCB拟合RMS平均值最低为1.239 4 mm，最高为3.437 0 mm；伪距IFCB拟合RMS平均值最低为0.069 3 m，最高为0.123 8 m。

对于BDS-3 GEO卫星，虽然FFT得到的周期结果证明其相位和伪距IFCB具有周期性，但变化数值较小，因此在模型化中可以忽略其时变性。采用均值作为常数值对其进行模型化处理，结果如图 8所示。计算每日拟合RMS的平均值，结果见表 3。

由图 8和表 3可知，使用均值作为常数值对IFCB进行模型化处理的精度较高：相位IFCB拟合RMS平均值分别为0.372 8 mm和0.533 5 mm；伪距IFCB对应的平均值为0.035 2 m。

3 结语

本文围绕BDS-2系统和部分BDS-3系统GEO卫星的三频观测值偏差展开讨论。使用2022年doy1~15的186个IGS观测站数据，对北斗GEO卫星IFCB变化特性进行分析。对相位和伪距IFCB进行解算发现，BDS-2系统GEO卫星IFCB的变化具有明显的周期性与时变性，而BDS-3系统GEO卫星无明显时变性。采用快速傅里叶变换对BDS-2与BDS-3 GEO卫星的周期进行分析并对周期结果进行统计，结果表明，相位和伪距IFCB出现频率最高的8个周期值均为12 h、8 h、6 h、4.8 h、4 h、3.428 6 h、3 h、2.666 7 h。

基于上述周期，使用4阶、8阶谐波函数混合函数和3~6次多项式对BDS-2 GEO卫星的IFCB时间序列进行模型化处理。结果表明，在相位IFCB模型化处理中，谐波函数混合函数和多项式拟合效果较为相近；在伪距IFCB模型化处理中，谐波函数混合函数效果明显优于多项式函数。对模型化拟合RMS进行分析可知，随着多项式次数的增加，拟合效果有所提升，但多项式拟合精度始终低于谐波函数混合函数，8阶谐波函数混合函数的拟合效果最佳。由于BDS-3 GEO卫星的时变性不明显，因此采用均值作为常数值对其IFCB进行模型化处理，结果表明，该模型化方法精度较高。