RTKLIB Manual之AppendixE Models and Algorithms解读

缩略语：

光速/m/s；

L频段伪距/m；

L频段载波相位/cycle

L频段载波相位/m；

导航信号发送至接收机时间/秒——即：接收机收到信号的时间

卫星发送导航信号时间/s；

卫星和接收机天线之间的几何距离/m；

卫星和接收机天线之间的伪距率/m/s;

在时间t时，ECEF下卫星位置/m；

在时间t时，ECEF下卫星速度/m/s；

在时间t时，ECEF下接收机天线位置/m;

在时间t时，ECEF下接收机天线速度/m/s;

ECEF下，接收机到卫星的单位向量/LOS；

当地坐标系下，接收机到卫星的单位向量；

ECEF坐标系到大地坐标系的旋转变换矩阵；

ECEF坐标系到卫星载体坐标系的旋转变换矩阵；

接收机天线位置的纬度/大地维度；

接收机天线位置的经度/大地经度；

当前位置观测，卫星的方位角/rad；

当前位置观测，卫星的高度角/rad；

时间t，接收机钟差/s；

时间t,接收机钟漂/s/s——频差；

时间t,卫星钟差/s；

时间t,卫星钟漂s/s——频差；

对流层延迟/m;

L频段电离层延迟/m;

L频段载波频率/Hz;

L频段波长/m；

L频段载波偏差/cycle；

L频段载波整周模糊度/cycle；

伪距测量误差/m；

载波相位测量误差/cycle；

载波相位测量误差/m；

地球自转角速度/rad/s；

当地坐标系下，L频段接收机天线相位中心偏差(PCO)/m;

卫星载体坐标系下，卫星天线L波段相位中心偏差(PCO)/m；

L频段接收机天线相位中心变化PCV)/m;

L频段卫星天线相位中心变化PCV)/m;

接收机天线到卫星的天底角/α；

当地坐标系下，接收机天线位置处的地球潮偏移/disp;

相位缠绕效应/cycle;

卫星时钟的相对论改正/Trel;

时间t下，ECI和ECEF的坐标系转换矩阵；

绕X轴旋转θ得到的坐标旋转矩阵；

绕Y轴旋转θ得到的坐标旋转矩阵；

绕Z轴旋转θ得到的坐标旋转矩阵；

卫星端，卫星j和k之间的单差；

地面端，基准站b和移动站/r之间的单差；

对流层天顶总延迟/m；

对流层天顶干延迟/m；

对流层东向梯度分量；

对流层北向梯度分量；

天顶对流层干延迟映射函数；

天顶对流层湿延迟映射函数；

电离层延迟映射函数；

E.1 Time System

       RTKLIB内部采用GPST进行GNSS数据处理和定位；协调世界时/UTC是一个不连续时间，简单说，和GPST有一个“跳秒”的关系。GPST通常采用周数/week和周内秒/TOW来表示；但是对于double变量类型来说，其距离精度只有0.04米，因此，在实际运算时，RTKLIB调用C LIBRARY中的时间结构体：

同时，RTKLIB提供函数能够进行时间的加、减，时间变换等。

• GPST与UTC转换公式：

GPS时间系统起点：00:00:00 UTC on January 6, 1980。

简单计算是跳秒差，严格计算系数，需要在导航电文里面解析相关参数！

• GLONASS与UTC

GLONASS与UTC一样，是一个不连续时间系统，粗略计算，存在三个小时的延迟。

• Galileo与UTC相差 13秒

Galileo时间系统是一个连续时间系统，其与GPST相差1024周，因此，在RINEX文件中，Galileo与GPS提供相同的周数。

Galileo时间系统起点：00:00:00 UTC on August 22, 1999。

• QZSST时间与GPST一致
• BDT与UTC相差 14秒

BDT时间系统是一个连续时间系统，系统起点：00:00:00 UTC on January 1, 2006.
 

E.2 Coordinates System 

利用广播星历计算出的卫星位置是基于ECEF下；

大地坐标转换到ECEF下坐标：

• Geocentric Latitude 地心纬度
• Geodetic Latitude 地理纬度
• Ellipsoidal Height 椭球高

1、从大地坐标:(Latitude、Longitude、Ellipsoidal Height)转到ECEF下的XYZ

2、从ECEF下的坐标XYZ转换到大地坐标系下的LLH：

其中地理纬度的计算，需要经过迭代！

3、从ECEF坐标系转到ENU/当地坐标系：

其中旋转矩阵E中的纬度、经度是接收机在大地坐标系下的位置。

如下公式中，Rr为接收机起始位置在ECEF下的坐标，Recef为任意时刻接收机在ECEF下的坐标：

4、大地水准面模型/Geoid models

由于GPS导航系统，以WGS-84坐标系为基准，因此，求出来的高度为椭球高/ellipsoidal height；但是大地水准面模型更符合真实的地球模型，其高度为 geodetic height，两者之间关系用如下公式表示：

要想使用：必须下载“Geoid Data File ”;同时在设置选项中，做出如下配置：

E.3 GNSS Signal Measurement Models 

1、信号结构：

接收机接收到GNSS卫星的信号，主要为：载波、扩频码、导航数据。信号结构如下图所示：

2、伪距测量模型：

主要包括：卫星和接收机之间的几何距离、卫星端钟差、接收机钟差、电离层延迟误差、对流层延迟误差、噪声等。

两种表示方式：

1、利用时间差求取伪距：

2、误差项展开表示：

3、载波相位/cycle/m测量模型：

Cycle为单位：

m为单位：

其中，ΔΦ包括：卫星端和接收机端的PCO、PCV，地球固体潮极移，相位缠绕，卫星钟相对论效应等

4、卫星和接收机的几何距离：

扣除地球自转的影响，Sagnac Effect，即信号发射和信号接收两个时刻点，不在同一个坐标系，要把这两个信号归到同一个坐标系中。

5、接收机到卫星视距、方位角，高度角：

此处，没有考虑地球自转的影响；

• LOS为：接收机指向卫星

• 方位角和高度角：方位角北偏东为正；

E.4 GNSS Satellite Ephemerides and Clocks

1、GPS、Galileo、Qzss卫星位置、速度、钟差、频移计算：

星历：

options选项配置：

计算过程：

2、利用广播星历计算Glonass

星历信息：

相对论效应已经在星历中（关于卫星钟的星历参数）做了补偿，因此，不需要额外公式计算补偿。

3、利用广播星历计算BeiDou:

星历信息：

对于MEO、IGSO卫星相关信息的计算与GPS导航系统一致，只需要注意三个参数：

对于GEO卫星，在计算参数时有区别！请查阅接口文件。

4、计算 SBAS GEO在ECEF下的位置、钟差:

message type 9

ECEF下SBAS GEO卫星位置、钟差计算：

5、利用SBAS GEO 播发位置、钟差改正项：

在option中，如下配置：

播发参数：

6、精密星历和钟差：

在option中，如下配置：

7、RTCM SSR校正卫星位置、钟差：

星历参数如下：

在option中，如下配置：

E.5 Troposphere and Ionosphere Models 

一、对流层模型：

1、Saastamoinen model/标准模型：

参数：大气压力、绝对温度、海拔高、相对湿度、天顶角（pi/2-高度角）

在RTKLIB中，用椭球高代替海拔高；相对湿度为70%

在option中，如下配置：

2、SBAS 对流层模型：

在option中，如下配置：

3、精密对流层模型：

相对湿度/relative humidity 为0，天顶角为0时，计算天顶对流层干延迟，其映射函数为：NMF，

湿延迟映射函数为:GMF

在option中，如下配置：

Estimate ZTD:

Estimate ZTD+Grad:此种模式下，

二、电离层模型：

1、Klobuchar model

主要针对单频系统，模型参数为：

 计算过程：

在option中，如下配置：

2、SBAS ionosphere model

type18:ionospheric grid point masks

type26: ionospheric delay corrections

在option中，如下配置：

3、Single layer model/单层模型

简化为一个单层模型，进行电离层延迟计算：

首先计算电离层刺穿点的经纬度：其中Re=6378137，H=350000。

如果电子总量知道，则得到电离层延迟：

在option中，如下配置：目前仅用于后处理模式中！

4、电离层无关LC/线性组合 ：无电离层延迟观测值LC

满足下列方程的线性组合皆为无电离层延迟的观测值：

nf1+mf2=0

在option中，如下配置：目前仅用于单点定位和PPP定位模式！

E.6 Single Point Positioning

       RTKLIB采用迭代加权最小二乘进行单点定位！注意，在RTKLIB中所有变量的存储是以“列”优先，因此，在一些公式的描述中和教材会出现“转置”的情况！

加权最小二乘解为：

1、单点定位：接收机位置和钟差

对于卫星定位来说，其方程为非线性，因此需要在某一个点进行线性化，然后进行加权最小二乘定位！

流程如下：

初始参数x0

观测模型求偏导数矩阵/设计矩阵：

经过泰勒展开后：

线性化后，方程为：

经过加权最小二乘，最终得到方程解：

循环，迭代：

RTKLIB中，设置迭代10次，判断结束条件为：状态变化量小于1e-4

单点定位，在option中，如下配置：

其中，在单点定位中，系统默认采用的是L1频点伪距；如果设置了无电离层模型，L1&L2用于GPS、GLONASS、QZSS；L1&L5用于Galileo.

实现流程：

对于第一个历元系统状态可以设为0；后面的历元初始状态直接用前一个历元的状态。

权矩阵如何设置：

最终求出接收机位置和钟差：

注意：求解出接收机位置的时刻是经过接收机钟差补偿后的GPST时；

2、接收机速度和频差：

注意点：

1. 在接收机速度计算过程中，同样存在地球自转现象，导致卫星信号发射时刻和接收信号时刻不是同一个坐标系，因此需要转换。
2. 加权矩阵为单位阵。
3. RTKLIB中用L1频段的多普勒作为测量值进行测速。

3、Solution Validation and Raim FDE

首先进行Solution Validation 判断：要同时保证几个条件成立：

1、卡方分布；

2、定位方程个数（包括了人为添加与系统间时间相关的观测方程）大于求解状态参数个数；

3、保证参与定位卫星构成的Gdop小于gdop阈值；其中在求Gdop时，已经将高度角小于预设最小高度角除去。

在option中，如下配置：

如果判断失败，则进行Raim FDE:

判断条件：计算残差平方和最小，为最终结果：

程序如下：

/* raim fde (failure detection and exclution) -------------------------------*/
static int raim_fde(const obsd_t *obs, int n, const double *rs,const double *dts, const double *vare, const int *svh,const nav_t *nav, const prcopt_t *opt, sol_t *sol,double *azel, int *vsat, double *resp, char *msg)
{obsd_t *obs_e;sol_t sol_e={{0}};char tstr[32],name[16],msg_e[128];double *rs_e,*dts_e,*vare_e,*azel_e,*resp_e,rms_e,rms=100.0;int i,j,k,nvsat,stat=0,*svh_e,*vsat_e,sat=0;trace(3,"raim_fde: %s n=%2d\n",time_str(obs[0].time,0),n);if (!(obs_e=(obsd_t *)malloc(sizeof(obsd_t)*n))) return 0;rs_e = mat(6,n); dts_e = mat(2,n); vare_e=mat(1,n); azel_e=zeros(2,n);svh_e=imat(1,n); vsat_e=imat(1,n); resp_e=mat(1,n); for (i=0;i<n;i++) {/* satellite exclution */for (j=k=0;j<n;j++) {if (j==i) continue; //每次去除一颗卫星obs_e[k]=obs[j];matcpy(rs_e +6*k,rs +6*j,6,1);matcpy(dts_e+2*k,dts+2*j,2,1);vare_e[k]=vare[j];svh_e[k++]=svh[j];}/* estimate receiver position without a satellite */if (!estpos(obs_e,n-1,rs_e,dts_e,vare_e,svh_e,nav,opt,&sol_e,azel_e,vsat_e,resp_e,msg_e)) {trace(3,"raim_fde: exsat=%2d (%s)\n",obs[i].sat,msg);continue;}for (j=nvsat=0,rms_e=0.0;j<n-1;j++) {if (!vsat_e[j]) continue;rms_e+=SQR(resp_e[j]);nvsat++;}if (nvsat<5) {trace(3,"raim_fde: exsat=%2d lack of satellites nvsat=%2d\n",obs[i].sat,nvsat);continue;}rms_e=sqrt(rms_e/nvsat);trace(3,"raim_fde: exsat=%2d rms=%8.3f\n",obs[i].sat,rms_e);if (rms_e>rms) continue;/* save result */for (j=k=0;j<n;j++) {if (j==i) continue;matcpy(azel+2*j,azel_e+2*k,2,1);vsat[j]=vsat_e[k];resp[j]=resp_e[k++];}stat=1;*sol=sol_e;sat=obs[i].sat;rms=rms_e;vsat[i]=0;strcpy(msg,msg_e);}if (stat) {time2str(obs[0].time,tstr,2); satno2id(sat,name);trace(2,"%s: %s excluded by raim\n",tstr+11,name);}free(obs_e);free(rs_e ); free(dts_e ); free(vare_e); free(azel_e);free(svh_e); free(vsat_e); free(resp_e);return stat;
}

在option中，如下配置：

E.7 Kinematic, Static and Moving‐Baseline

在DGPS/DGNSS、Static、Kinematic and Moving-base 定位模式中，RTKLIB利用EKF求解其估计参数结果。

一、EKF的五个公式：

首先利用某一历元的观测量求解待估的状态量、状态协方差矩阵：

其中：

h(x)为观测模型观测向量；

H(x)为偏导数矩阵/设计矩阵；

R(x)为观测误差协方差矩阵；

（-）与（+）代表测量更新前与后

状态更新方程为：状态和状态协方差矩阵更新

二、短基线下的双差观测模型：

对于基线距离小于10km，移动站和基准站的L频段的载波相位和伪距双差方程构成观测量：

其中，卫星和接收机的钟差、电离层、对流层和其他小误差量认为已经扣除。

       在第一个方程中，第三项为载波相位观测值改正项，短基线可以忽略，但使用不同 PCV的接收机天线的基线需要考虑。除了Moving-base定位模式，其他模式下，需要提前预知基准站的位置。

1、站间单差：

       在基准站和移动站间的单差中，由于两者的采样频率不同，一般基准站为1HZ，移动站为10HZ，为进行站间单差 SD， RTKLIB 采用一种简单的标准来选择同步观测数据，只选择比流动站观测历元更早或相等的最后一次基准站数据，这种流动站与基准站间的时间偏差即所谓的“差分龄期”。 随着龄期的增大，由于卫星时钟漂移和电离层延迟的变化，解的精度会随之降低。

2、星间单差：

       至于星间单差 SD，RTKLIB 逐历元选取高度角最高的卫星作为参考卫星。但GPS和GLONASS之间不能进行星间差分，因为不同类型导航系统即便相同频率的信号， 接收机的群延迟也是不同的。这种群延迟差异被称为接收机系统内偏差（ISB，inter system bias）。

在option中，如下配置：

1、设置基准站位置/Base Station

2、设置差分数据龄期：

假设流动站、 基准站均为三频 GPS/GNSS 接收机， 未知参数状态向量 x 可定义为：

其中，B1/B2/B5是三频站间载波相位单差模糊度/cycle，RTKLIB 之所以不使用“双差模糊度”，是为避免频繁更换参考卫星影响。 同时， 站间“单差模糊度”也有助于解决 GLONASS 频分多址信号模糊度解算的问题。

双差载波、伪距观测量为：

三、短基线EKF的测量更新：

在如下方程的基础上：

忽略载波相位改正项，其观测模型向量h(x)，设计矩阵H(x)，协方差矩阵R为：

根据上述的方程，通过EKF能够计算出任一历元的接收机位置、速度、单差模糊度浮点解。

四、EKF的时间更新/状态更新

1、Positioning Mode = Kinematic and RECDynamics = ON

在option中，如下配置：

在以上设置下，状态转移矩阵和状态噪声矩阵分别如下：

其中τ为时间间隔/采样时间。

2、Positioning Mode = Kinematic and REC Dynamics = OFF

在option中，如下配置：

在以上设置下，状态转移矩阵和状态噪声矩阵分别如下：

为避免数值不稳定性对定位的影响， RTKLIB 对各历元接收机位置方差的初始值重置为较大值（10e4）。 为避免非线性信号观测模型相互影响， RTKLIB 使用单点定位计算初始位置。

3、Positioning Mode = Static

在option中，如下配置：

在以上设置下，状态转移矩阵和状态噪声矩阵分别如下：

在模糊度求解方式为：Integer Ambiguity Resolution = Instantaneous 时，其站间单差模糊度值并不随着时间更新，更新到下一个历元，而是，在mode=Instantaneous下，单差模糊度在每个历元初始化初始值，并且其方差噪声为10e4。

如果 RTKLIB 探测到观测值发生周跳， 站间单差模糊度同样会被初始化。 目前， RTKLIB 探测周跳的方法包括： LLI（失锁标志）、电离层残差法（几何无关 LC 线性组合）。 用户可通过ʺSlip Thresʺ设置周跳探测阈值。

五、整周模糊度解算：

一旦EKF测量更新得到状态变量的估计值，载波相位模糊度浮点解可进一步用于求解模糊度整数解从而提高精度和收敛时间(Integer Ambiguity Res = Continuous,Instantaneous or Fix and Hold)。 首先， 将状态估值、方差矩阵由单差 SD 形式转为双差 DD 形式， 具体如下：

 式中对模糊度参数由单差模糊度SD转双差模糊度DD（G 对位置、速度使用的单位矩阵，未改变它们的形式；仅模糊度由单差转为双差形式），消除了接收机初始相位，得整周模糊度N和方差Q。通过式 E.7.17 解决整数最小二乘方式问题，得整周模糊度的最优估计向量N。

RTKLIB 采用知名的 LAMBDA/扩展 LAMBDA(MLAMBDA)方法是解决整数最小二乘问题。 LAMBDA 和 MLAMBDA 结合了线性变换来缩小整数向量搜索空间， 是变换空间中一个巧妙的树搜索过程。 该方法通常使用ʺRatio‐Testʺ来检测整数解的可靠性， ʺRatio‐Testʺ的 Ratio 值 R 即次优解的方差残差之和 N2 与最优解的方差残差之和 N 的比值。 用户通过ʺMin Ratio to Fix Ambiguityʺ可以设置Ratio 值，但目前 RTKLIB 版本仅支持 ratio 定值。

当 ratio 检验通过后， 由式 E.7.19 得接收机天线位置、 速度的“固定解”， 若检验失败， 则 RTKLIB 输出“浮点解”。

当模糊度解算选项为ʺFix and Holdʺ (Integer Ambiguity Resolution = Fix and Hold)， 经上述模糊度检验已获得了“固定解”时， 双差模糊度 DD 值则直接被固定为整数值。 因此，RTKLIB 根据 F.7.1 引入下列“伪”观测值到 EKF 并更新EKF。

RTKLIB V2.4.0 初次采用ʺFix and Holdʺ， 尤其对动态模式动态接收机跟踪，提升 ratio 值发挥了作用。

六、长基线双差测量模型

       其中长基线与短基线相比，L波段中双差载波、伪距的电离层、对流层延迟需要考虑；对于baseline>100km，需要考虑精密星历以消除广播星历存在的误差；对于baseline>500km，需要考虑载波相位改正项（ISB）、潮汐效应。为确定电离层延迟值可以采用无电离层 LC 组合的方法（ linear combination），但RTKLIB基线解算并不使用这种 LC 方法，而是使用 EKF， 双频或三频观测值直接确定电离层延迟。

待估计状态参数为：

Zr、Zb表示流动站、基准站的天顶对流层延迟；

GN,r、GE,r、GN,b和GE,b表示流动站/基准站的对流层梯度在北方向、东方向的分量；

I=(I1rb, I2rb,…, Imrb,)T表示 L1 频点垂直方向的单差电离层延迟值（单位:m）；

观测模型向量和设计矩阵如下：

 

时间更新如上所示：

       式中：QT、QI表示对流层、电离层项的解算噪声矩阵。方程认为流动站和基准站的天顶对流层延迟和梯度值、各卫星的单差电离层延迟值为随机游走（random‐walk）。 为确定电离层、对流层延迟，RTKLIB 对 2.4.1版本的长基线解算引入ʺPartial fixingʺ特性， 即部分模糊度固定的方法。 该方法对所有模糊度中的一部分进行模糊度固定，而其他的模糊度仍使用浮点值。 RTKLIB 采用卫星高度角作为某个模糊度是否固定的标准，如果卫星低于高度角阈值，则不固定该卫星的模糊度，只有当卫星满足高度角限制时才解算整周模糊度。 ʺFix and Holdʺ特性中，通过数据处理选项ʺMin Elevation to Fix Ambʺ和ʺMin Elevation to Hold Ambʺ设置整周模糊度的高度角阈值。

七、动态基线模型：

       动态基线通常用于流动站和基准站接收机都处于运动状态、 只要求获取流动站相对基准站相对位置的情况。 当一个运动平台上安装有两台天线后， 可用动态基线模式计算姿态角。 用户通过ʺMoving‐Baseʺ选项选择动态基线解算模式。

       在动态模型中，基准站位置由每个历元的单点定位得到。当获得基准站位置后，移动站位置通过短基线模式下的Kinematic模式，计算EKF的五个方程得到。因此， 动态基线模式下流动站和基准站的位置精度只相当于单点定位的精度， 只有相对位置关系才有意义。

       在动态基线模式中， RTKLIB 改正了基准站、 流动站间的时间偏差。 由于流动站和基准站时间不同步， 接收机时钟偏差最大可达 2ms， 会导致快速移动平台的精度下降，因此利用基准站坐标进行基线解算前，要先对基准站坐标进行偏差改正，具体方法如下：

式中， tr 和 tb 是流动站、基准站通过单点定位后得到的信号接收时刻。 Vb 表示通过多普勒观测值获得的基准站速度。 另外， RTKLIB 可使用ʺBaseline Length Constraintʺ数据处理选项打开基线长约束，解算动态基线方位角， 该约束应用于EKF观测值更新的表达式如下：

式中， 基线长 rbaseline 表示基线长设定值（单位： m）， 表示基线长方差值（ 单位：m）。 为处理基线长度很短情况下的非线性问题， 可通过ʺNumber of Filter Iterationʺ选项设置为大于 1 的值， 实现 EKF 观测量的迭代更新。 

E.8 PPP (Precise Point Positioning)

RTKLIB中，采用L波段的载波、伪距非差观测数据，通过EKF方式，得到参数估计。

1、PPP非差测量模型：

对载波相位和伪距采用无电离层LC组合，

2、接收机天线相位中心模型

PCO：相对于天线参考点ARP (antenna reference point).
PCV：与卫星高度角和方位角有关；

计算：

3、卫星天线相位中心模型：

PCO相对卫星质量中心，并且参数是基于卫星载体坐标系：因此，需要进行坐标转换，变到ECEF下面

卫星载体坐标系：

4、Site displacement by earth tides/固体潮汐
 

option中，配置如下：

如果选择Solid 则调用相关的函数求出固体潮变化；

如果选择Solid/OTL,则系统同时进行Solid和OTL改正，其中第一个为调用函数，第二项改正还需要添加 ʺOTL BLQ Fileʺ.文件：

5、相位缠绕效应/Phase windup correction

6、接收机参数估计 by PPP

待估参数为：位置、速度、钟差、天顶对流层延迟、东向、北向梯度、非差无电离层观测值模糊度

非差LC组合载波和伪距观测量如下：

测量模型向量、设计矩阵：

option配置如下：

参考：

RTKLIB MANUAL