RTKLIB的相关模型和算法

主要是RTKLIB ver. 2.4.2 Manual 的翻译.

E.3 GNSS信号测量模型

(1)GNSS信号结构

Figure E.3‐1 展现了典型的GNSS信号架构。GNSS信号一般由载波频率(carrier)、测距码(code)和导航数据(data)的乘积组成。测距码也称为伪随机噪声码(PRN)。GPS、GLONASS、Galileo、QZSS、北斗、SBAS提供的这些GNSS信号的详细规格见附录F。

(2)伪距测量模型

伪距定义为ʺ接收机天线到卫星天线的距离，包括接收机和卫星的时钟偏差(以及其他偏差，如大气延迟)ʺ[9]。$L_i$伪距$P_{r,i}^s$可以用接收时钟测量的信号接收时间${\bar{t}}_r(s)$和卫星时钟测量的信号发送时间${\bar{t}}^s(s)$表示为:
$$P_{r,i}^s=c({\bar{t}}_r-{\bar{t}}^s)------(E.3.1)$$
该方程可以用卫星和接收机天线之间的几何范围${\rho }_r^s$，接收机和卫星时钟偏差$d{t}_{r}d{T}^s$，电离层和对流层延迟$I_{r,i}^s,T_r^s$和测量误差${\epsilon }_P$写成:

(3)载波相位和相位范围测量模型

载波相位为ʺ…实际上是对卫星信号的接收载波与接收机生成的参考频率之间的拍频进行测量"。$L_i$的载波相位${\varphi }_{r,i}^s$可以表示为:

其中$t_0$为初始时间(s)，${\varphi }_{r,i}(t)$为t时刻接收机本振的$L_i$相位(循环)，${\varphi }_i^s$为t时刻发射导航信号$L_i$相位(循环)。${\varphi }_{r,0,i}$为接收机本振在$t_0$时刻$L_i$初始相位(循环)，${\varphi }_{r,0,i}^s$为发射导航信号在$t_0$时刻$L_i$初始相位(循环)。
$L_i$相位范围${\varphi }_{r,i}^s$定义为载波相位乘以载波频率${\lambda }_i$(m)，也可以用载波相位偏置$B_{r,i}^s$和载波相位校正项$d{\Phi }_{(}^{s}r,i)$.包括天线相位中心偏移和变化、固体潮引起的台站位移、相位闭合效应和卫星时钟上的相对论修正为:

$N_{r,i}^s$通常被称为载波相位整数模糊度、载波周期模糊度或简单模糊度。载波相位校正条款的具体制定请参见附录E.9。

(4)接收机与卫星天线的几何距离

几何距离定义为在惯性坐标下卫星天线相位中心位置与接收机天线相位中心位置的物理距离。首先，信号的传输时间$t^s$可以由:
$$t^s=\overline{(}t{)}_r-P_{r,i}^s/c-dT(t^s)---(E.3.7)$$

方程的两边都包含了$t^s$。因此需要多次迭代才能解出方程。几何范围可以通过ECEF(地球中心地球固定)坐标上$t_r时刻和t^S时刻$的接收机和卫星天线相位中心位置$r_r=(x_r,y_r,z_r{)}^T$和$r^s=(x^s,y^s,z^s{)}^T$表示为:
$${\rho }_r^s=\Vert U(t_r)r_r(t_r)-U(t^s)r^s(t^s)\Vert ---(E.3.8)$$
其中$U(t)$为t时刻的ECEF到ECI(地心惯性)坐标变换矩阵。对于ECEF坐标中的表达式，必须考虑地球自转效应，才能得到几何范围。这个方程可以用下列方程中的一个来近似，在1mm水平下有足够的精度。当前版本RTKLIB总是使用方程(F.3.8b)来表示几何范围。(F.3.8b)中的最后一项有时被称为Sagnac效应。

(5)卫星方向的方位角和仰角

从接收机到卫星的单位视线矢量(line - of - sight)可以用ECEF坐标表示为:

方程中忽略了地球自转的影响。从接收站点得到的卫星方向的方位角和仰角$A{z}_r^s$和$E{l}_r^s$由:

其中$E_r$为从ECEF到接收机位置局部坐标的坐标旋转矩阵。有关矩阵的详细构成，请参阅E.2。

E.4 GNSS卫星星历表和时钟

RTKLIB支持GPS、GLONASS、Galileo、QZSS、北斗和SBAS的广播星历和时钟。它还支持精确的星历表和时钟，如SP3‐c[22]和时钟RINEX[15]，包括伽利略，QZSS和北斗的后处理模式。对于实时模式，还支持由SBAS长期快速校正和RTCM 3 SSR(状态空间表示)校正校正的广播星历表和时钟。下面的方程显示了RTKLIB中使用的星历和时钟模型。

(1) GPS、Galileo和QZSS的广播星历表和时钟

GPS、Galileo和QZSS的广播星历和SV时钟参数在导航电文中给出如下:
$$p_{eph}(t_{oe},t_{oc},IDO)=(a,e,i_0,{\Omega }_0,\omega ,M_0,\Delta n,\dot{I},\dot{\Omega },C_{us},C_{uc},C_{rs},C_{rc},C_{is},C_{ic},a{f}_0,a{f}_1,a{f}_2,T_{GD}{)}^T$$
利用这些参数计算出ECEF中的卫星位置(天线相位中心位置)$r^s(t)$，卫星时钟偏差$d{T}^s(t)$和时钟漂移$d{\dot{T}}^s(t)$为:


其中:
$\mu$:地球引力常数( $3.9860050\ast 10^{14}{m}^3/s^2$for GPS and QZSS, $3.986004418\ast 10^{14}{m}^3/s^2$for Galileo
${\omega }_e$:地球角速度(7.2921151467*10^-5rad/s)
$b=f_1^2/f_i^2$for $L_i$的伪距
$T_{GD}$:GPS and QZSS的群延迟参数,对$B_{GD}$ for galileo
开普勒方程(E.4.4)可以用牛顿ʹ的方法进行以下迭代求解。

如果处理选项为ʺSatellite Ephemeris/ clockʺ到ʺbroadcastʺ，包括GLONASS、北斗和SBAS，则应用广播星历表和时钟。

（2）GLONASS广播星历表和时钟参数

在导航消息中给出如下:
$$p_{eph}(t_b)=(x,y,z,v_x,v_y,v_z,a_x,a_y,a_z,{\tau }_n,{\gamma }_n)$$
ECEF (PZ90.02)中卫星位置$r^s(t)=(x,y,z{)}^T$和速度$v^s(t)=(v_x,v_y,v_z{)}^T$的差分方程可以表示为：

式中:
$a_e$:地球半长轴(6378136.0 m)
$\mu$:地球引力常数(9 398600.4410⁹m³/s²)
${\omega }_e$:地球角速度(7.29211510^-5rad/s)
$J_2$:地球位势的二次纬向谐波(1082625.7*10^-9)
$$r=\sqrt{x^2+y^2+z^2}$$
注意，GLONASS ICD 5.1[4]的A.3.1.2中的两个勘误表已在上述模型中进行了修正。
在时间t时卫星的位置$r^s(t)$和速度$v^s(t)$可以通过RK4(龙格库塔第四阶和阶段)推导来解决这些微分方程数值积分与最初的卫星位置$r^s(t_b)$和速度$v^s(t_b)$的参考时间$t_b$。在历元时刻t的卫星时钟偏差$d{T}^s(t)$ 和漂移$d{\dot{T}}^s(t)$也可推导为:

GLONASS时钟参数中包含了卫星时钟中的相对论效应。所以在这种情况下，相对论修正是不适用的。

（3）北斗的广播星历表和时钟

北斗卫星导航信息中提供了与GPS、Galileo、QZSS类似的星历表和时钟参数如下:

（4）SBAS卫星的广播星历表和时钟

SBAS地理卫星的导航消息参数在SBAS消息(消息类型9)中给出如下:

(5) SBAS轨道和时钟校正

SBAS轨道和时钟校正定义为以下参数。

(6)精确星历表和时钟

GPS、GLONASS、伽利略、QZSS和北斗的精确星历表通常以SP3‐c文件的形式每15分钟或5分钟提供一次，包含卫星位置和速度(可选)。为了得到t时刻的卫星位置，需要进行适当的插值。RTKLIB使用固定度(10 n)多项式插值，通过Newton‐Nevilleʹ算法为:

(7) RTCM SSR轨道和时钟校正

RTCM SSR轨道和时钟校正定义为以下参数。

E.5 对流层和电离层模型

(1)对流层模型

标准大气可以表示为:

(2) SBAS对流层模型

如果将处理选项ʺ对流层校正ʺ设置为ʺSBASʺ，则应用SBAS接收机规范中定义的SBAS对流层模型。该模型通常称为ʺMOPS模型ʺ。详情请参考[8]A.4.2.4。

(3)精确对流层模型

如果将处理选项ʺ对流层校正ʺ设置为ʺEstimate ZTDʺ或ʺEstimate ZTD+Gradʺ，则应用更精确的对流层模型，严格映射函数为:

(4)广播电离层模型

对于单频GNSS用户，GPS和QZSS导航数据包含以下广播电离层参数。

如果处理选项ʺ电离层校正ʺ设置为ʺ广播ʺ或ʺQZSS广播ʺ，则应用广播电离层模型进行校正。

(5) SBAS电离层模型

对电离层延迟的SBAS修正由18型电离层格点掩码和26型电离层延迟修正提供。如果处理选项ʺ电离层校正ʺ设置为ʺSBASʺ，RTKLIB将使用SBAS电离层校正，并且这些SBAS消息将在输入文件中提供。模型算法及IGPs(电离层网格点)定义参见SBAS接收机规范[8]的A.4.4.9和A.4.4.10。

(6)单层模型

电离层通常被建模为简单的单层模型，如图E.5‐1所示。单层模型也称为薄壳模型。

(7)无电离层LC(线性组合)

为了消除GNSS信号测量中的电离层效应，GNSS数据处理中经常使用双频测量的线性组合。iL和jl伪橙和相位范围的无电离层LC表示为:

E.6单点定位

RTKLIB采用迭代加权LSE(最小二乘估计)用于ʺSingleʺ(单点定位)模式，有或没有SBAS校正。

(1)线性LSE

假设有一个测量向量y，它可以被建模为一个未知参数向量x和一个随机测量误差向量v的线性方程。
$$y=Hx+v---(E.6.1)$$
最小二乘损失函数$J_{JS}$定义为测量误差平方和:
$$J_{LS}=v_1^2+v_2^2+...+v_m^2=v^{T}v---(E.6.2)$$
利用(E.6.1)和(E.6.2)，损失函数可以改写为:
$$\begin{array}{rl}{J}_{LS}& =(y-Hx{)}^T(y-Hx)\\ & =y^{T}y-y^{T}Hx-x^T{H}^{T}y+x^T{H}^{T}Hx---(E.6.3)\end{array}$$
为了使代价函数最小化，$J_{LS}$的梯度应为零。然后
$$\begin{array}{rl}\frac{\partial J_{LS}}{\partial x}& =0^T-y^{T}H-(H^{T}y{)}^T+(H^{T}Hx{)}^T+x^T{H}^{T}H\\ & =-2y^{T}H+2x^T{H}^{T}H=0---(E.6.4)\end{array}$$
它给出了所谓的ʺ正规方程ʺ为:
$$H^{T}Hx=H^{T}y---(E.6.5)$$
对于法向方程，通过LSE估计得到的未知参数向量$\hat{x}$为:
$$\hat{x}=(H^{T}H{)}^{-1}{H}^{T}y---(E.6.6)$$
如果给出了每个测量值的权重，则成本函数(E.6.3)可以用一个权重矩阵W。
$$J_{WLS}=v^{T}Wv---(E.6.7)$$
为了最小化代价函数$J_{WLS}$，我们可以通过加权LSE得到估计的未知参数向量，方法与简单LSE类似:
$$\hat{x}=(H^{T}WH{)}^{-1}{H}^{T}Wy---(E.6.8)$$
加权LSE的权矩阵W通常为:
$$W=diag({\sigma }_1^{-2},{\sigma }_2^{-2},...,{\sigma }_m^{-2})$$
其中${\sigma }_i$是第i个测量误差的先验标准差。

(2)非线性LSE的高斯-牛顿迭代

如果测量结果不是线性模型，则测量方程可以用一般的非线性向量函数表示为:
$$y=h(x)+v---(E.6.9)$$
式中$h(x)$是参数向量x的测量向量函数，对初始参数向量$x_0$采用泰勒级数可将方程推广为:
$$h(x)=h(x_0)+H(x-x_0)+...---(E.6.10)$$
式中，H为$h(x)$对x在x=x_0处的偏导数矩阵:
$$H=\frac{\partial h(x)}{\partial x}{|}_{x=x_0}---(E.6.11)$$
假设初始参数充分接近真实值，泰勒级数的第二项和进一步项可以忽略。我们可以将(E.6.9)近似为:
$$y\approx h(x_0)+H(x-x_0)+v---(E.6.12)$$
我们可以得到以下线性方程:
$$y-h(x_0)=H(x-x_0)+v---(E.6.13)$$
对(E.6.13)应用线性加权LSE (E.6.8)，可以得到非线性加权LSE的法向方程:
$$H^{T}WH(\hat{x}-x_0)=H^{T}W(y-h(x_0))---(E.6.14)$$
因此，我们可以得到估计的未知参数向量$\hat{x}$:
$$\hat{x}=x_0+(H^{T}WH{)}^{-1}{H}^{T}W(y-h(x_0))---(E.6.15)$$
如果初始参数x₀在真实值附近不够，我们可以迭代改进估计的参数，如:
$$\begin{gathered} \widehat{x_0}=x_0---(E.6.16) \\ {\hat{x}}_{i+1}={\hat{x}}_i+(H^{T}WH{)}^{-1}{H}^{T}W(y-h({\hat{x}}_i))---(E.6.17) \end{gathered}$$
如果迭代收敛，我们可以得到最终的估计参数为:
$$\hat{x}=\underset{i\to \infty }{\lim }{\hat{x}}_i---(E.6.18)$$
迭代LSE通常被称为高斯-牛顿法。需要注意的是，这种迭代并不总是通过简单的高斯牛顿法收敛，特别是对于具有大非线性的病态测量方程。在这种情况下，我们应该对这种非线性LSE采用另一种策略。最常用的非线性LSE方法是LM (Levenberg‐Marquardt)方法。

(3)接收机位置和时钟偏差的估计

对于ʺ单点定位ʺ模式为ʺ定位模式ʺ，采用以下单点定位程序逐历获取最终解。对于一个历元时间，定义未知参数向量x为:

$$x=(r_r^T,cd{t}_r{)}^T---(E.6.19)$$
伪距测量向量y可表示为:
$$y=(P_r^1,P_r^2,P_r^3,...,P_r^m{)}^T---(E.6.20)$$
其中$P_r^s$为伪距测量值。如果处理选项ʺ电离层校正ʺ设置为ʺ无电离层LCʺ，则使用附录E.5(7)中定义的无电离层LC(线性组合)伪距。在其他情况下，只使用L1伪距。

单点定位测量方程及其偏导数矩阵为:
$$h(x)=\left(\begin{array}{c}{\rho }_r^1+cd{t}_r-cd{T}^1+I_r^1+T_r^1\\ {\rho }_r^2+cd{t}_r-cd{T}^2+I_r^2+T_r^2\\ {\rho }_r^3+cd{t}_r-cd{T}^3+I_r^3+T_r^3\\ ⋮\\ {\rho }_r^m+cd{t}_r-cd{T}^m+I_r^m+T_r^m\end{array}\right)H=\left(\begin{array}{cc}-e_r^{1T}& 1\\ -e_r^{2T}& 1\\ -e_r^{3T}& 1\\ ⋮& ⋮\\ -e_r^{mT}& 1\end{array}\right)---(F.6.21)$$
其中几何范围${\rho }_t^s$和LOS矢量$e_r^s$由E.3(4)和E.3(5)与卫星和接收机位置给出。卫星位置$e^s$和时钟偏差$d{T}^s$也来自于E.4中描述的GNSS卫星星历和时钟，根据处理选项ʺ卫星星历/时钟ʺ。
为了求解测量方程，得到最终估计的接收机位置和接收机时钟偏差，RTKLIB采用迭代加权LSE为:
$${\hat{x}}_{i+1}={\hat{x}}_i+(H^{T}WH{)}^{-1}{H}^{T}W(y-h({\hat{x}}_i))---(E.6.22)$$
对于迭代加权LSE的初始参数向量x₀，所有0都用于单点定位的第一阶。一旦得到一个解，该位置用于下一个epoch初始接收机位置。对于权值矩阵W, RTKLIB使用以下公式
$$\begin{gathered} W=diag({\sigma }_1^{-2},{\sigma }_2^{-2},...,{\sigma }_m^{-2})---(E.6.23) \\ {\sigma }^2=F^s{R}_r(a_{\sigma }^2+b_{\sigma }^2/sinE{l}_r^s)+{\sigma }_{eph}^2+{\sigma }_{bias}^2+{\sigma }_{ion}^2+{\sigma }_{trop}^2---(E.6.24) \end{gathered}$$
其中：
$$\begin{array}{rl}{F}^s& :卫星系统误差因子(1:GPS、Galileo、QZSS、北斗，1.5:GLONASS,3.0:SBAS)\\ R_r& :码/载波相位误差比率\\ a_{\sigma },b_{\sigma }& :载波相位误差因子a和b(m)\\ {\sigma }_{eph}& :星历表标准偏差和时钟误差(m)\\ {\sigma }_{ion}& :电离层校正模型误差标准偏差(m)\\ {\sigma }_{trop}& :对流层修正模型误差标准偏差(m)\\ {\sigma }_{bias}& :码偏误差标准差(m)\end{array}$$
RTKLIB使用市建局(用户距离精度)或类似的指标来衡量星历表的标准偏差和时钟误差。经过多次迭代，求解在正常情况下收敛，得到估计的接收机位置${\hat{r}}_r$和接收机时钟偏置$d{\hat{t}}_r$。
$$\hat{x}=\underset{i\to \infty }{\lim }{\hat{x}}_i=({\hat{r}}_r^T,{\widehat{cdt}}_r{)}^T---(E.6.25)$$
估计的接收机时钟偏差$d{\hat{t}}_r$不会显式输出到解决方案文件。相反，它被包含在解决方案的时间标签中。这意味着解决时标不是指接收器时标，而是指在GPST中测量的真实信号接收时间。

（4）估计接收器速度和时钟漂移

如果给出了GNSS信号的多普勒频率测量，则可以按照以下步骤估算接收器的速度和时钟漂移。 在某个时间段内，将用于速度估计的未知参数向量x定义为：
$$x=(v_r^T,cd{\dot{t}}_r{)}^T---(E.6.26)$$
其中$v_r^T$和$cd{\dot{t}}_r$分别为接收机速度，单位为ECEF (m/s)和接收机时钟漂移(s/s)。量程速率测量矢量y可以表示为:
$$y=(-{\lambda }_i{D}_{r,i}^1,-{\lambda }_i{D}_{r,i}^2,-{\lambda }_i{D}_{r,i}^3,...,-{\lambda }_i{D}_{r,i}^m{)}^T---(E.6.27)$$
其中，$D_{r,i}^s$为卫星s的L_i多普勒频率测量值。RTKLIB总是使用L₁多普勒频率测量值。这些测量方程及其偏导数矩阵构成如下:
$$h(x)=\left(\begin{array}{c}{r}_r^1+cd{\dot{t}}_r-cd{\dot{T}}^1\\ r_r^2+cd{\dot{t}}_r-cd{\dot{T}}^2\\ r_r^3+cd{\dot{t}}_r-cd{\dot{T}}^3\\ ⋮\\ r_r^m+cd{\dot{t}}_r-cd{\dot{T}}^m\end{array}\right)H=\left(\begin{array}{cc}-e_r^{1T}& 1\\ -e_r^{2T}& 1\\ -e_r^{3T}& 1\\ ⋮& ⋮\\ -e_r^{mT}& 1\end{array}\right)---(E.6.28)$$
在这些方程中，接收器和卫星之间的距离‐范围$r_r^s$来自于:
$$r_r^s=e_r^{sT}(v^{s}s(r^s)-v_r)+\frac{{\omega }_e}{c}(v_y^s{x}_r+y^s{v}_{x,r}-v_x^s{y}_r-x^s{v}_{y,r})---(E.6.29)$$
式中，$v^s=(v_x^s,v_y^s,v_z^s,{)}^T,v_r=(v_{x,r},v_{y,r},v_{z,r}{)}^T$。利用类似于接收机位置估计的迭代LSE，我们可以得到接收机速度和时钟漂移如下:
$$\hat{x}=\underset{i\to \infty }{\lim }{\hat{x}}_i=({\dot{v}}_r^T,c\hat{d}{\dot{t}}_r{)}^T---(E.6.30)$$
其中权矩阵W设为I(非加权LSE)。

(5)方案验证和RAIM FDE

(3)中所描述的估计接收机位置可能包含由于未建模测量误差而无效的解。RTKLIB在得到估计的接收机位置后，为了测试是否有效并剔除无效的解，采用了以下验证。如果验证失败，该解决方案会被警告消息拒绝。
$$\begin{gathered} v_s=\frac{(P_r^s-({\hat{\rho }}_r^s+c\hat{d}{t}_r-cd{T}^s+I_r^s+T_r^s))}{{\sigma }_s}---(E.6.31) \\ v=(v_1,v_2,v_3,...v_m{)}^T---(E.6.32) \\ \frac{v^{T}v}{m-n-1}<{\chi }_{\alpha }^2(m-n-1)---(E.6.33) \\ GDOP<GDO{P}_{thres}---(E.6.34) \end{gathered}$$
其中n为估计参数的个数，m为测量值的个数。${\chi }_{\alpha }^2(n)$是自由度n的卡方分布和$\alpha =0.001(0.1$。GDOP是几何DOP(精度稀释)。GDOP_thres可以设置为处理选项ʺGDOP的拒绝阈值ʺ。
除上述方案验证外，2.4.2版本中还增加了RAIM(接收机自主完整性监测)FDE(故障检测和排除)功能。如果处理选项ʺRAIM FDEʺ被启用，(E.6.33)的卡方测试失败，RTKLIB通过排除一个又一个可见卫星来重新估计。在所有重试后，选择归一化平方残差最小的估计接收机位置tv v作为最终解。在这种方案中，由于卫星故障、接收机故障或大多径导致的无效测量被排除为离群值。请注意，对于两个或多个无效度量，此特性是无效的。它还需要两颗冗余的可见卫星，这意味着至少需要6颗可见卫星才能得到最终的解决方案。

E.7 动态，静态和移动基线

RTKLIB利用EKF(扩展卡尔曼滤波)，结合附录E.3中的GNSS信号测量模型和附录E.5中的对流层和电离层模型，获得了DGPS/DGNSS、静态、运动和移动基线模式下的最终解。

(1) EKF公式

利用EKF，可以用测量向量y_k在历元t_k处估计未知模型参数的状态向量x及其协方差矩阵P:
$$\begin{array}{rl}{\hat{x}}_k(+)& ={\hat{x}}_k(-)+K_k(y_k-h({\hat{x}}_k(-)))---(E.7.1)\\ P_k(+)& =(I-K_{k}H({\hat{x}}_k(-)))P_k(-)---(E.7.2)\\ K_k& =P_k(-)H({\hat{x}}_k(-))(H({\hat{x}}_k(-))P_k(-)H({\hat{x}}_k(-){)}^T+R_k{)}^{-1}---(E.7.3)\end{array}$$
其中${\hat{x}}_k$和P_k为在历元时间t_k时的估计状态向量及其协方差矩阵。(-)和(+)表示EKF测量更新前后。$h(x),H(x)and{R}_k$分别为测量模型向量、偏导数矩阵和测量误差协方差矩阵。假设系统模型是线性的，EKF的状态向量及其协方差矩阵的时间更新表示为:
$$\begin{array}{rl}{\hat{x}}_{k+1}(-)& =F_k^{k+1}{\hat{x}}_k(+)---(E.7.4)\\ P_{k+1}(-)& =F_k^{k+1}{P}_k(+)F_k^{k+1T}+Q_k^{k+1}---(E.7.5)\end{array}$$
其中，$F_k^{k+1}$和$Q_k^{k+1}$是系统噪声从历时t_k到t_k+1的转移矩阵和协方差矩阵。

(2)DD短基线双差测量模型

对于探测器r和基站b之间的短基线(< 10 km)载波相对定位，通常使用以下DD测量方程来计算iL相位范围和伪距。在这些方程中，利用DD技术几乎消除了卫星和接收机的时钟偏差、电离层和对流层效应和其他小的修正项。
$$\begin{array}{rl}{\Phi }_{rb,i}^{jk}& ={\rho }_{rb}^{jk}+{\lambda }_i(B_{rb,i}^j{-}_{rb,i}^k)+d{\Phi }_{r,i}^s+{\epsilon }_{\Phi }\\ P_{rb,i}^{jk}& ={\rho }_{rb}^{jk}+{\epsilon }_P---(E.7.6)\end{array}$$

其中，$d{\Phi }_{r,i}^s$是载波相位校正项，在短基线情况下可以忽略，除了带有不同天线的接收机PCV项。为了得到方程中的几何范围${\rho }_r^s$，除移动基线情况外，将基站位置$r_b$固定为预定值。
请注意，接收机之间的SD最好是在同一历元时间的测量之间进行。然而，由于接收机的时钟偏差不同，接收机并不是完全同步的。在某些典型情况下，漫游者的采样间隔与基站不同，如10hz和1hz。为了控制SD, RTKLIB采用一个简单的标准来选择一个测量对。RTKLIB只是选择在探测器测量历元时间之前或等于这个时间的最后一次测量。漫游者和基站之间的历时差有时被称为ʺ差异年龄ʺ。随着时间差的增大，由于卫星时钟漂移和电离层时延的变化，求解精度逐渐降低。为了仅补偿卫星时钟漂移，RTKLIB利用广播SV时钟参数对SD测量进行校正。最大ʺAge of Differentialʺ设置为处理选项ʺMAX Age of Diffʺ。
在卫星侧SD生成方面，RTKLIB逐次选择一颗具有最大俯仰角的参考卫星。请注意，不同导航系统的卫星之间不会产生卫星侧SD，比如GPS和GLONASS之间。这是因为对于不同导航系统的信号，即使它们具有相同的载波频率，接收机通常也有不同的群时延。接收机的组延迟差称为接收机ISB(系统间偏置)。
假设移动站和基站均采用三频GPS/GNSS接收机，则待估计的未知状态向量x可定义为:
$$x=(r_r^T,v_r^T,B_1^T,B_2^T,B_5^T{)}^T---(E.7.7)$$
其中$B_i=(B_{rb,i}^1,B_{rb,i}^2,...,B_{rb,i}^m)$为L_iSD(单差)载波相位偏差(循环)。作为RTKLIB实现，它内部使用SD载波相位偏差而不是DD，以避免参考卫星的麻烦移交处理。SD偏差还有助于解决GLONASS FDMA信号中的整数歧义。
测量矢量y也被DD相位范围和伪距测量定义为:
$$y=({\Phi }_1^T,{\Phi }_2^T,{\Phi }_5^T,P_1^T,P_2^T,P_5^T{)}^T---(E.7.8)$$
其中：
$$\begin{gathered} {\Phi }_i=({\Phi }_{rb,i}^{12},{\Phi }_{rb,i}^{13},{\Phi }_{rb,i}^{14},...,{\Phi }_{rb,i}^{1m}{)}^T \\ {P}_i=(P_{rb,i}^{12},P_{rb,i}^{13},P_{rb,i}^{14},...,P_{rb,i}^{1m}{)}^T \end{gathered}$$

(3)短基线EKF的测量更新

由式(E.7.6)，测量模型向量$h(x)$，偏导数矩阵$H(x)$，测量误差R的协方差矩阵可记为:
$$h(x)=({h_{\Phi ,1}}^T,{h_{\Phi ,2}}^T,{h_{\Phi ,5}}^T,{h_{P,1}}^T,{h_{P,2}}^T,{h_{P,5}}^T{)}^T---(E.7.9)$$

$$H(x)=\frac{\partial h(x)}{\partial x}{|}_{x=\hat{x}}=\left(\begin{array}{ccccc}-DE& 0& {\lambda }_{1}D& 0& 0\\ -DE& 0& 0& {\lambda }_{2}D& 0\\ -DE& 0& 0& 0& {\lambda }_{5}D\\ -DE& 0& 0& 0& 0\\ -DE& 0& 0& 0& 0\\ -DE& 0& 0& 0& 0\end{array}\right)---(E.7.10)$$

$$R=\left(\begin{array}{cccccc}D{R}_{\Phi ,1}{D}^T& & & & & \\ & D{R}_{\Phi ,2}{D}^T& & & & \\ & & D{R}_{\Phi ,5}{D}^T& & & \\ & & & D{R}_{P,1}{D}^T& & \\ & & & & D{R}_{P,2}{D}^T& \\ & & & & & D{R}_{P,5}{D}^T\end{array}\right)---(E.7.11)$$
其中：
$$h(x)=\left(\begin{array}{c}{\rho }_{rb}^{12}+{\lambda }_i(B_{rb}^1-B_{rb}^2)\\ {\rho }_{rb}^{13}+{\lambda }_i(B_{rb}^1-B_{rb}^3)\\ ⋮\\ {\rho }_{rb}^{1m}+{\lambda }_i(B_{rb}^1-B_{rb}^m)\end{array}\right)H=\left(\begin{array}{c}{\rho }_{rb}^{12}\\ {\rho }_{rb}^{13}\\ ⋮\\ {\rho }_{rb}^{1m}\end{array}\right)$$
$$\left(\begin{array}{ccccc}1& -1& 0& \cdots & 0\\ 1& 0& -1& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& 0& 0& \cdots & -1\end{array}\right):SD(单差)矩阵$$
$$\begin{array}{rl}E& :(e_r^1,e_r^2,...,e_r^m{)}^T\\ R_{\Phi ,i}& :diag({2{\sigma }_{\Phi ,i}^1}^2,{2{\sigma }_{\Phi ,i}^2}^2,...,{2{\sigma }_{\Phi ,i}^m}^2)\\ R_{P,i}& :diag({2{\sigma }_{P,i}^1}^2,{2{\sigma }_{P,i}^2}^2,...,{2{\sigma }_{P,i}^m}^2)\\ {\sigma }_{\Phi ,i}^s& :L_i相位范围测量误差的标准偏差(m)\\ {\sigma }_{P,i}^s& :L_i伪距测量误差的标准偏差(m)\end{array}$$
利用这些方程求解EKF公式(E.7.1)，得到了漫游者天线的位置、速度和浮动SD载波相位偏差的历元时间t_k。

(4)EKF时间更新

对于RTKLIB中带有接收机动力学的运动学定位模式(positioning mode = kinematic and REC dynamics = ON)， EKF (E.7.2)的时间更新表示为:
$$F_k^{k+1}=\left(\begin{array}{ccc}{I}_{3\times 3}& I_{3\times 3}{\tau }_r& \\ & I_{3\times 3}& \\ & & I_{(3m-3)\times (3m-3)}\end{array}\right),Q_k^{k+1}=\left(\begin{array}{ccc}{0}_{3\times 3}& & \\ & Q_v& \\ & & 0_{(3m-3)\times (3m-3)}\end{array}\right)---(E.7.12)$$
其中：
$$Q_v={E_r}^{T}diag({{\sigma }_{ve}}^2{\tau }_r,{{\sigma }_{vn}}^2{\tau }_r,{{\sigma }_{vu}}^2{\tau }_r)E_r$$
而且${\tau }_r=t_{k-1}-t_k$为GPS/GNSS接收机采样间隔(s)，$({\sigma }_{ve},{\sigma }_{vn},{\sigma }_{vu})$为探测器速度系统噪声东、北、上分量的标准差($m/s/\sqrt{s}$)。
对于没有接收机动力学的纯动态模式(Positioning mode = kinematic and REC dynamics = OFF)，公式(E.7.9)代入:
$$F_k^{k+1}=\left(\begin{array}{ccc}{I}_{3\times 3}& & \\ & I_{3\times 3}& \\ & & I_{(3m-3)\times (3m-3)}\end{array}\right),Q_k^{k+1}=\left(\begin{array}{ccc}{\infty }_{3\times 3}& & \\ & 0_{3\times 3}& \\ & & 0_{(3m-3)\times (3m-3)}\end{array}\right)---(E.7.13)$$
为了避免数字不稳定，请在接收器位置的方差中添加无限的处理噪声，接收器位置状态将重置为每个时期的初始猜测值，并且将足够大的过程噪声($10^4{m}^2$)添加到RTKLIB的方差中。 初始位置来自单点定位过程，该过程用于避免非线性信号测量模型的迭代。
在静态模式下(Positioning mode = static)，公式(F.7.10)可简单替换为:
$$F_k^{k+1}=\left(\begin{array}{ccc}{I}_{3\times 3}& & \\ & I_{3\times 3}& \\ & & I_{(3m-3)\times (3m-3)}\end{array}\right),Q_k^{k+1}=\left(\begin{array}{ccc}{0}_{3\times 3}& & \\ & 0_{3\times 3}& \\ & & 0_{(3m-3)\times (3m-3)}\end{array}\right)---(E.7.14)$$
在瞬时模糊度分辨模式(整数模糊度分辨=瞬时)中，SD载波相位偏差B_i的时间更新处理方式与上述略有不同。在这种模式下，载波相位偏置状态的值不会被EKF时间更新后续到下一个epoch。在每个时点，将偏差重置为初始的猜测值，并在方差中加入足够大的过程噪声(10⁴ m²)。如果在测量数据中检测到一个周跳，相应的SD载波相位偏置状态也被重置为初始值。RTKLIB通过输入测量数据中的LLI(丢失锁定指示器)检测周期滑移，如果双频测量可用，则通过几何自由LC(线性组合)检测相位跳变。周期滑移阈值可以通过处理选项ʺslip Thresʺ来更改。

(5)整数模糊度求解

一旦EKF测量中得到的估计状态更新，浮动载波相位模糊可以被分解成整数值，以提高精度和收敛时间(整模糊度Res =连续、瞬时或固定保持)。首先，将估计状态及其协方差矩阵转化为DD形式:
$$\begin{gathered} {\hat{x}}_k^{\prime }=G{\hat{x}}_k(+)=({{\hat{r}}_r}^T,{{\hat{v}}_r}^T,{\hat{N}}^T{)}^T---(E.7.15) \\ {P}_k^{\prime }=G{P}_k(+)G^T=\left(\begin{array}{cc}{Q}_R& Q_{NR}\\ Q_{RN}& Q_N\end{array}\right)---(E.7.16) \end{gathered}$$
其中:
$$G=\left(\begin{array}{cccc}{I}_{6\times 6}& & & \\ & D& & \\ & & D& \\ & & & D\end{array}\right):SD到DD的变换矩阵$$
在这种转换中，SD载波相位偏差被转移到DD载波相位形式，以消除接收机初始相位项，获得整数模糊$\hat{N}$及其协方差Q_N。在这些公式中，通过求解整数最小二乘问题得到最适合整数二义性的整数向量$\breve{N}$，表达式为:

$$\breve{N}=argmin((N-\hat{N}{)}^T{Q}_N^{-1}(N-\hat{N}))(N\in Z)---(E.7.17)$$
为了解决盲降问题，RTKLIB采用了一种众所周知的高效搜索策略LAMBDA[66]及其扩展MLAMBDA[67]。LAMBDA和MLAMBDA提供了一个线性变换的组合，以缩小整数向量搜索空间和一个转换空间中熟练的树搜索程序。通过以下简单的ʺRatio‐Testʺ来验证这些程序的整数向量解。在ʺRatio‐Testʺ中，ratio‐factor R，定义为第二优解${\breve{N}}_2$的残差平方根的加权和与第一优解$\breve{N}$的加权和之比，用于检验解的可靠性。验证阈值R_thres可以通过处理选项ʺ最小修正歧义比率ʺ来设置。当前版本RTKLIB只支持固定的阈值。
$$R=\frac{({\breve{N}}_2-\hat{N}{)}^T{Q_N}^{-1}({\breve{N}}_2-\hat{N})}{(\breve{N}-\hat{N}{)}^T{Q_N}^{-1}(\breve{N}-\hat{N})}>R_{thres}---(E.7.18)$$
验证后，通过求解下式得到rover天线位置和速度${\breve{r}}_r$和${\breve{v}}_r$的ʺFIXEDʺ解。如果验证失败，RTKLIB输出ʺFLOATʺ解决方案${\hat{r}}_r$和${\hat{v}}_r$。
$$\left(\begin{array}{c}{\breve{r}}_r\\ {\breve{v}}_r\end{array}\right)=\left(\begin{array}{c}{\hat{r}}_r\\ {\hat{v}}_r\end{array}\right)-Q_{RN}{Q_N}^{-1}(\hat{N}-\breve{N})---(E.7.19)$$
如果处理选项设置为ʺFix and Holdʺ模式(整数模糊度解决方案=固定和保持)，并且固定的解决方案经之前的测试正确验证，DD载波相位偏置参数被严格约束为固定的整数值。为此，RTKLIB向EKF输入以下ʺ伪ʺ测量值，并通过(F.7.1)更新EKF。
$$\begin{gathered} y=\breve{N} \\ h(x)=Gx \\ H(x)=G \\ R=diag({{\sigma }_c}^2,{{\sigma }_c}^2,{{\sigma }_c}^2,...) \end{gathered}$$
其中：
$$\begin{gathered} G=\left(\begin{array}{cccc}0& D& & \\ 0& & D& \\ 0& & & D\end{array}\right):SD到DD的变换矩阵 \\ {\sigma }_c:对固定整数二义性的约束(=0.001周期) \end{gathered}$$
ʺFix and Holdʺ模式在RTKLIB版本中首次引入。2.4.0，以提高固定率，特别是在运动模式下跟踪移动接收机。

(6) Long baseline DD measurement model

对于移动站r与基站b之间的长基线处理，可以形成类似于短基线DD模型的DD测量方程:
$$\begin{array}{rl}{\Phi }_{rb,i}^{jk}& ={\rho }_{rb}^{jk}-I_{rb,k}^{jk}+T_{rb}^{jk}+{\lambda }_i(B_{rb,i}^j{-}_{rb,i}^k)+d{\Phi }_{r,i}^s+{\epsilon }_{\Phi }\\ P_{rb,i}^{jk}& ={\rho }_{rb}^{jk}+I_{rb,i}^{jk}+T_{rb}^{jk}+{\epsilon }_P---(E.7.24)\end{array}$$

其中，$I_{ri}^s$和$T_r^s$分别作为电离层延迟(m)和对流层延迟(m)添加到短基线DD模型中。应使用卫星位置的精确星历表来减轻100公里以上基线的广播星历表误差。在载波相位校正项$d{\Phi }_{r,i}^s$中，对于基线大于500公里的情况，应考虑固体潮效应。为了消除电离层项，有时会形成无电离层LC(线性组合)。然而，RTKLIB没有使用这种明确的LC，而是使用EKF的双频或三频测量来直接估计电离层项，用于基线处理。

长基线情况下的未知状态向量x也可求解为:
假设移动站和基站均采用三频GPS/GNSS接收机，则待估计的未知状态向量x可定义为:
$$x=(r_r^T,v_r^T,Z_r,G_{N,r},G_{E,r},Z_b,G_{N,b},G_{E,b},I^T,B_1^T,B_2^T,B_5^T{)}^T---(E.7.25)$$
其中Z_rand Z_b为移动站和基站站点的ZTD(天顶总延迟)，$G_{N,r},G_{E,r},G_{N,b},G_{E,b}$是对流层梯度的北部和东部分量。$I=(I_{rb}^1,I_{rb}^2,...,I_{rb}^m)$也是在频率为L₁ (m)的SD垂直电离层延迟。
测量模型向量h(x)和偏导数矩阵H(x)可表示为:
由式(E.7.6)，测量模型向量$h(x)$，偏导数矩阵$H(x)$，测量误差R的协方差矩阵可记为:
$$h(x)=({h_{\Phi ,1}}^T,{h_{\Phi ,2}}^T,{h_{\Phi ,5}}^T,{h_{P,1}}^T,{h_{P,2}}^T,{h_{P,5}}^T{)}^T---(E.7.26)$$

$$h_{\Phi ,i}=\left(\begin{array}{c}{\rho }_{rb}^{12}+T_{rb}^{12}-{\gamma }_k(m_I^1{I}_{rb}^1-m_I^2{I}_{rb}^2)+{\lambda }_i(B_{rb}^1-B_{rb}^2)+d{\Phi }_{rb,i}^{12}\\ {\rho }_{rb}^{13}+T_{rb}^{13}-{\gamma }_k(m_I^1{I}_{rb}^1-m_I^3{I}_{rb}^3)+{\lambda }_i(B_{rb}^1-B_{rb}^3)+d{\Phi }_{rb,i}^{13}\\ ⋮\\ {\rho }_{rb}^{1m}+T_{rb}^{1m}-{\gamma }_k(m_I^1{I}_{rb}^1-m_I^m{I}_{rb}^m)+{\lambda }_i(B_{rb}^1-B_{rb}^m)+d{\Phi }_{rb,i}^{1m}\end{array}\right)---(E.7.27)$$

$$h_{P,i}=\left(\begin{array}{c}{\rho }_{rb}^{12}+T_{rb}^{12}+{\gamma }_k(m_I^1{I}_{rb}^1-m_I^2{I}_{rb}^2)\\ {\rho }_{rb}^{13}+T_{rb}^{13}+{\gamma }_k(m_I^1{I}_{rb}^1-m_I^3{I}_{rb}^3)\\ ⋮\\ {\rho }_{rb}^{1m}+T_{rb}^{1m}+{\gamma }_k(m_I^1{I}_{rb}^1-m_I^m{I}_{rb}^m)\end{array}\right)---(E.7.27)$$