北斗/GNSS码观测值在卫星端和接收机端都受到与信号类型相关的硬件延迟影响，而不同码类型的硬件延迟并不相同。这种不同码类型之间的硬件延迟差异即差分码偏差(differential code bias, DCB)，是电离层建模和高精度定位中不可忽略的误差，影响可达数米^[1]。多GNSS星座播发的信号可依照频率、调制方式和接收机跟踪方式等分类，其中MGEX(multi-GNSS experiment) 网多系统观测数据所涉及的信号类型近40种。理论上不同信号类型间的硬件延迟均存在差异，因而为实现多系统多频组合高精度定位，需要精确处理各类码偏差^[2]。

DCB可按不同码类型的信号频率是否相同分为频内DCB和频间DCB。频内DCB一般可由同频伪距观测值作差估计，频间DCB一般基于电离层建模理论估计获得。部分学者提出采用精密单点定位的方法估计频间DCB^[3-4]，所获得的码偏差与IGS码偏差产品具有良好的一致性。但此方法需要依赖精密轨道和钟差产品，且处理站网的运算效率较低。CODE、DLR和CAS等IGS机构基于电离层建模理论估计出各系统码偏差，并将其作为产品供用户使用^[5-7]。其中，由于DCB格式在应对各种码类型组合时修正方法较为复杂，CODE和CAS也推出独立定义每种原始观测量偏差的OSB(observable-specific signal bias)格式产品^[8]。

本文拟采用相位平滑伪距法，联合多系统双频观测值实现电离层建模，并进一步将电离层模型回代解算各系统多频DCB，最终转换为OSB格式。本文算法利用MGEX网测站数据与广播星历，不依赖于外部事后产品，能以较快的时效性和较高的精度解算码偏差。

1 码偏差估计模型 1.1 DCB估计

将接收机r从卫星s观测的不同类型伪距观测值$P_{\mathrm{r}, i}^{\mathrm{s}} $和$ P_{\mathrm{r}, j}^{\mathrm{s}}$互差，可以获得卫星端DCB和接收机端DCB之和SPRDCB(satellite-plus-receiver differential code bias)。忽略噪声、电离层延迟高阶项和多路径效应影响，其表达式为：

$ \left\{\begin{array}{l} P_{\mathrm{r}, i}^{\mathrm{s}}-P_{\mathrm{r}, j}^{\mathrm{s}}=B_{i j}^{\mathrm{s}}+B_{\mathrm{r}, i j}, f_i=f_j \\ P_{\mathrm{r}, i}^{\mathrm{s}}-P_{\mathrm{r}, j}^{\mathrm{s}}=B_{i j}^{\mathrm{s}}+B_{\mathrm{r}, i j}+\left(\frac{40.3}{f_i^2}-\frac{40.3}{f_j^2}\right) \cdot \\ \quad \mathrm{STEC}_{\mathrm{r}}^{\mathrm{s}}, f_i \neq f_j \end{array}\right. $

(1)

式中，f_i和f_j为不同码类型i和j的观测值频率；B_ij^s和B_{r, ij}分别为卫星端和接收机端的DCB，一般在单天内视为常数；STEC_r^s为卫星和接收机间信号路径的总电子含量。

为估计和消除电离层延迟项，可通过高精度的STEC观测值将STEC模型化。由于伪距无几何组合无法精确反映电离层延迟的变化，且单系统卫星观测值弧段不够饱满，故采用多系统相位平滑伪距作为STEC的观测值^[9]：

$ \left\{\begin{array}{l} \widetilde{P}_{\mathrm{r,GF}}^{\mathrm{G,s},}=A^{\mathrm{G}, \mathrm{s}} \mathrm{STEC}_{\mathrm{r}}^{\mathrm{G}, \mathrm{s}}+B^{\mathrm{G}, \mathrm{s}}+B_{\mathrm{r}}^{\mathrm{G}} \\ \widetilde{P}_{\mathrm{r}, \mathrm{GF}}^{\mathrm{R,s}}=A^{\mathrm{R}_k, \mathrm{s}} \mathrm{STEC}_{\mathrm{r}}^{\mathrm{R}, \mathrm{s}}+B^{\mathrm{R}, \mathrm{s}}+B_{\mathrm{r}}^{\mathrm{R}} \\ \widetilde{P}_{\mathrm{r}, \mathrm{GF}}^{\mathrm{E,s}}=A^{\mathrm{E}, \mathrm{s}} \mathrm{STEC}_{\mathrm{r}}^{\mathrm{E}, \mathrm{s}}+B^{\mathrm{E}, \mathrm{s}}+B_{\mathrm{r}}^{\mathrm{E}} \\ \widetilde{P}_{\mathrm{r}, \mathrm{GF}}^{\mathrm{C,s} }=A^{\mathrm{C}, \mathrm{s}} \mathrm{STEC}_{\mathrm{r}}^{\mathrm{C}, \mathrm{s}}+B^{\mathrm{C}, \mathrm{s}}+B_{\mathrm{r}}^{\mathrm{C}} \end{array}\right. $

(2)

式中，G、R、E和C分别表示GPS、GLONASS、Galileo和BDS；$ \widetilde{P}_{\mathrm{r}, \mathrm{GF}}^{*, \mathrm{~s}}$为各系统的无几何组合的相位平滑伪距观测值，其中*为卫星系统的通配符；B_r^*和B^{*, s}为各系统的接收机端和卫星端的DCB，作为天内常数估计$ A^{*, s}=40.3 / f_{*, i}^2-40.3 / f_{*, j}^2$，依据各系统卫星双频信号的频率确定，其中A^{R_k, s}需要根据每颗GLONASS卫星的频道号k确定；STEC_r^{*, s}可根据CODE的电离层建模方法^[10]，基于单层假设，采用球谐函数将其模型化：

$ \left\{\begin{array}{l} \operatorname{STEC}(\varphi, \lambda, z)=M(z) \cdot \operatorname{VTEC}(\varphi, \lambda) \\ \operatorname{VTEC}(\varphi, \lambda)=\sum\limits_{n=0}^{n_{\max }} \sum\limits_{m=0}^n \\ {\left[\widetilde{P}_{n m}(\sin \varphi)\left(\widetilde{C}_{n m} \cos (m \lambda)+\widetilde{S}_{n m} \sin (m \lambda)\right)\right]} \end{array}\right. $

(3)

式中，φ和λ为日固地磁坐标系下的纬度和经度；z为卫星位置相对于接收机位置的天顶距，M(z)为薄层穿刺点处的垂直方向电子含量VTEC到斜路径电子含量STEC(φ, λ, z)的投影函数，具体计算方法参考CODE的MSLM(modified single-layer model)投影函数；$ \widetilde{P}_{n m}$为n次m阶正交化缔合勒让德函数；n_max为球谐函数最大阶数，本文取15；$ \widetilde{C}_{n n}$和$ \widetilde{S}_{n m}$分别为各阶次的球谐系数。

联立式(1)和式(3)，可以获得各历元频间SPRDCB估计值；频内SPRDCB由同频伪距观测值直接作差获得各历元估计值。依据各历元卫星高度角可计算出SPRDCB的加权平均值：

$ \left\{\begin{array}{l} B_{\mathrm{r}, i j}^{\mathrm{s}}(t)=P_{\mathrm{r}, i}^{\mathrm{s}}(t)-P_{\mathrm{r}, j}^{\mathrm{s}}(t)- \\ \left(\frac{40.3}{f_i^2}-\frac{40.3}{f_j^2}\right) \mathrm{STEC}_{\mathrm{r}}^{\mathrm{s}}(t), f_i \neq f_j \\ B_{\mathrm{r}, i j}^{\mathrm{s}}(t)=P_{\mathrm{r}, i}^{\mathrm{s}}(t)-P_{\mathrm{r}, j}^{\mathrm{s}}(t), f_i=f_j \\ \hat{B}_{\mathrm{r}, i j}^{\mathrm{s}}=\frac{\sum\limits_{t=1}^n w(t) B_{\mathrm{r}, i j}^{\mathrm{s}}(t)}{\sum\limits_{t=1}^n w(t)} \\ \sigma_{\hat{B}}=\sqrt{\frac{\sum\limits_{t=1}^n\left(B_{\mathrm{r}, i j}^{\mathrm{s}}(t)-\hat{B}_{\mathrm{r}, i j}^{\mathrm{s}}\right)^2}{n-1}} \\ w(t)=\frac{1}{1+\cos ^2 e(t)} \end{array}\right. $

(4)

式中，w(t)为卫星高度角e(t)依据历元t所确定的权重$ \hat{B}_{\mathrm{r}, i j}^{\mathrm{s}}$为多历元SPRDCB的加权平均值，$ \sigma_{\hat{B}}$为其标准差；STEC_r^s(t)可在电离层建模后由式(3)计算获得。根据式(4)可以获得各系统任意两种码类型之间的SPRDCB；再通过施加卫星端DCB零均值约束，可以分离卫星端DCB。

1.2 多频OSB估计

各系统多频伪距观测值间的卫星端DCB可转换为多频伪距OSB。根据定义，可以认为DCB即两种码类型OSB之差，另外对双频基准类型OSB附加无电离层组合零值约束，可实现OSB的分离。综上，同一GNSS的非钟差基准类型以及钟差基准类型的OSB求解方式为：

$ \left[\begin{array}{c} B_{m n} \\ B_{m k} \\ \vdots \\ B_{m l} \\ 0 \end{array}\right]=\left[\begin{array}{ccccc} 1 & -1 & 0 & \cdots & 0 \\ 1 & 0 & -1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ 1 & 0 & 0 & \cdots & -1 \\ \alpha_{m n} & \beta_{m n} & 0 & \cdots & 0 \end{array}\right]\left[\begin{array}{c} \mathrm{OSB}_m \\ \mathrm{OSB}_n \\ \mathrm{OSB}_k \\ \vdots \\ \mathrm{OSB}_l \end{array}\right] $

(5)

式中，OSB_m和OSB_n为钟差基准类型的伪距OSB；$ \alpha_{m n}=f_m^2 /\left(f_m^2-f_n^2\right)$和$ \beta_{m m}=-f_n^2 /\left(f_m^2- f_n^2\right)$为无电离层组合系数；OSB_k和OSB_l为非钟差基准类型的伪距OSB。

2 多频码偏差处理

MGEX网的测站所观测到的多系统多频数据可以按照RINEX3.04格式标准区分，区分依据主要是信号频率、跟踪方式和通道。本文参考DLR和CAS码偏差产品所涉及的码类型，统计处理的码偏差类型如表 1。

数据处理步骤为：1)电离层建模。首先对多系统观测数据进行周跳探测和弧段划分，并形成无几何组合的相位平滑伪距观测值；然后利用广播星历计算卫星位置和信号路径在近似薄层的穿刺点位置，并根据球谐函数和DCB参数模型进行最小二乘平差；最后对各系统卫星端DCB引入零均值基准，解算出球谐系数。2)多频码偏差解算。首先利用同频但不同码类型的伪距观测值互差，求得频内SPRDCB；然后利用双频伪距无几何组合扣除由球谐函数模型计算的电离层延迟，获得频间SPRDCB；之后通过卫星端DCB零均值基准约束从SPRDCB中分离卫星端DCB；最后对钟差基准类型OSB施加无电离层组合约束，由多种卫星端DCB计算多种码类型的OSB。

上述步骤的具体策略见表 2。多频码偏差解算时，为削弱卫星低高度角时的电离层延迟模型化误差和多路径效应影响，设置截止高度角为20°。

3 结果与分析 3.1 SPRDCB稳定性

频间SPRDCB可由各历元伪距无几何观测值扣除电离层项提取，而其多历元稳定性可在一定程度上反映电离层模型精度。本文首先验证多系统电离层建模后提取的频间SPRDCB估值的稳定性。取2022-01-01 MGEX网全部测站数据，分别在单G系统、双系统和多系统电离层建模后，获得各观测历元的电离层延迟，再按照式(4)计算卫星观测弧段内的多历元SPRDCB估值及标准差。以北斗C26(C211)卫星为例，统计全球可观测此卫星的测站与其形成的SPRDCB标准差(图 1)。

由图 1可见，采用G、GR、GEC和GREC等系统组合电离层建模后，所得的C26与测站间的SPRDCB稳定性的多站平均值分别是1.482 ns、1.215 ns、1.197 ns和1.185 ns。使用单G系统电离层建模时，北美、澳洲和欧洲区域的大部分测站相应的SPRDCB日估值精度优于1.5 ns。另外，赤道区域电离层活跃程度较高，而极地和海洋等测站稀疏区域的观测弧段较少，因而这些区域的电离层模型精度较差，相关测站SPRDCB的估值标准差较大。双系统或多系统组合电离层建模相比于单G系统，可以提升大陆边缘、海洋等测站稀疏区域的SPRDCB估值精度。相比于众多电离层分析中心常使用的单G或GR组合电离层建模，GEC组合和GREC组合对SPRDCB的估值精度更高。进一步地，SPRDCB多历元标准差也将作为其分离站星DCB的定权依据。

另外，本文分别通过GREC四系统组合电离层建模和CODE电离层产品(采用GR双系统)解算SPRDCB，将两种电离层模型所得的SPRDCB均值差异绘制成图 2。G074、R802、E210和C201与各测站之间的SPRDCB均值差异平均值分别为0.35 ns、0.49 ns、0.45 ns和0.35 ns，并且大部分测站与卫星间的SPRDCB均值差异优于1 ns。这种参考CODE电离层产品解算结果的SPRDCB均值差异即为SPRDCB的外符合精度，与各区域的电离层模型精度相关。

各测站的SPRDCB估值精度主要与其地理分布、观测条件等因素有关。北美、澳洲和欧洲等区域的测站质量较好，且密度较高，所以这些区域电离层模型精度也相对较高，相关测站的频间SPRDCB的外符合精度也表现较好。南极区域电离层穿刺点投影误差较大^[11]，且测站分布较为稀疏，因而此区域电离层模型精度较低，相关测站的SPRDCB外符合精度表现较差。同时，由于E210卫星和C201卫星在南美洲测站的观测时段处于单天的始末时间，且时长较短，而始末时间中分段线性球谐函数对电离层模型的拟合效果较差，因而这些区域测站的SPRDCB外符合精度较差。另外，CODE电离层产品的VTEC精度在赤道及测站稀疏区域较差。综上，本文的SPRDCB估值外符合精度总体可靠。

3.2 码偏差稳定性分析

选取部分C2I类OSB值长期稳定的BDS-3卫星作为基准卫星，根据它们每天的日解OSB值形成拟稳基准，其时间序列如图 3(a)。在此基准下，绘制另一些C2I类型OSB不稳定卫星的OSB时间序列，如图 3(b)。图 3(b)中，OSB在2022年下半年发生跳变，这可能与卫星升级维护有关。在卫星升级维护过程中，其伪距偏差作为天内常数估计不一定准确，且可能发生显著变动，故应使用当天或者日期较近的OSB产品。

通过GREC四系统组合电离层建模，取2019-01~2023-07 MGEX网测站数据，按天解解算多系统多频伪距OSB，并统计各系统码偏差平均周稳定性，如图 4所示。图中，横轴为PRN和对应的SVN。总体上，各系统的多种OSB平均周稳定性可达0.07~0.26 ns；同一系统内，信号频率较高的码偏差较稳定；信号频率相同或相近的码偏差稳定性较为相近。

GPS和GLONASS的各类OSB较为稳定，平均周稳定性分别可达0.06~0.12 ns和0.10~0.17 ns。而Galileo和BDS-3的相同频点的OSB稳定性相近，为0.13~0.21 ns。其中BDS-3中PRN较大的部分卫星由于测站数较少，其OSB稳定性相对较低。BDS-2的各类OSB总体稳定性为0.17~0.26 ns，其中C06~C10、C13和C16等IGSO卫星的多频OSB稳定性较好，其余卫星的多频OSB稳定性较差，这与BDS-2的星上多路径效应及测站分布和观测条件有关。QZSS的J07属于GEO卫星，稳定性相对较差，QZSS其余卫星稳定性良好。

3.3 码偏差一致性分析

对于码偏差一致性分析，选取2019-01~2023-07 DLR和CAS发布的未经天线改正的系列码偏差产品作为参考，将本文解算结果与二者按天解统一基准后作比较，互差结果如表 3所示，其中本文解算的码偏差产品以WUM标记，“-”表示DLR或CAS无数据。

总体上，WUM-DLR和WUM-CAS在各类DCB上的一致性均可达0.3 ns。其中对于同频DCB类型，如GPS的C1C-C1W、C2W-C2L和C2W-C2X，QZSS的C1C-C1X和GLONASS的C1C-C1P，WUM与DLR具有良好的一致性，互差RMS为0.02~0.07 ns，且其与CAS的互差RMS也可保持在0.04~0.11 ns。对于GPS、GLONASS、Galileo和QZSS系统的各类不同频码类型间DCB，WUM与DLR、CAS均有较好的一致性，互差RMS基本在0.08~0.16 ns。因BDS-2系统频间DCB估值稳定性相对于其他系统较差，WUM与DLR、CAS一致性水平也相对较差。BDS的C2I-C6I类型DCB一致性较差，这可能与本文将BDS-2和BDS-3分为两个系统进行处理的策略有关。其中，随着BDS-3的MGEX网测站数量逐渐增多，WUM-DLR的BDS-3系统的C2I-C6I类型DCB一致性逐年趋近，2019~2022年每年的总体互差RMS依次为0.36 ns、0.33 ns、0.17 ns和0.13 ns。BDS-3中C2I与C1P、C5P、C1X、C5X、C7Z、C8X等新信号码类型形成的频间DCB不受是否合并BDS-2和BDS-3系统估计的策略差异影响，与DLR一致性较好，其互差RMS为0.15~0.18 ns。

4 结语

1) 基于GREC四系统组合电离层建模可实现最高稳定性的SPRDCB估计；利用GREC电离层建模与CODE电离层产品所得的SPRDCB均值一致性可达0.35~0.49 ns。

2) BDS-3卫星C2I类型OSB一般可长期保持稳定，但可能会因卫星升级维护发生较大变化。

3) 本文利用2019-01~2023-07的MGEX网数据逐日解算多系统多频码偏差，各类OSB平均周稳定性为0.07~0.26 ns，而各类DCB与DLR、CAS的一致性优于0.3 ns。其中，BDS-3与Galileo的相同频率OSB的稳定性相近，平均周稳定性为0.13~0.21 ns。

在上述工作的基础上，仍存在一些问题：1)对于多种码类型间的DCB，应如何选取DCB组合估计或考虑多种DCB的闭合差，并生成相应的OSB，从而改善码偏差对定位精度的影响；2)如何对接收机码偏差的短时变化精确建模。