RTKPPP定位算法流程.doc

1 基础知识
1.1 GPS精密单点定位的基本原理
GPS精密单点定位一般采用单台双频GPS接受机，运用IGS提供的精密星历和卫星钟差，基于载波相位观测值进行的高精度定位。观测值中的电离层延迟误差通过双频信号组合消除，对流层延迟误差通过引入未知参数进行估计。
1.2 时间系统
RTKLIB内部使用GPST(GPST时间)用于GNSS的数据解决和定位算法。数据在RTKLIB内部解决之前，需要转换成GPST时间。使用GPST的因素是避免解决润秒。RTKLIB使用以下结构体表达时间：
typedef struct
time_t time; /* time(s) expressed by standard time_t */
double sec; /* fraction of second under 1 s */
} gtime_t;
1.2.1 GPST和UTC(Universal Time Coordinated)
关系参考【图1】，参考【图2】：

图1 转换关系公式

图 2
通过使用GPS导航信息中的UTC参数，GPST到UTC或者UTC到GPST之前的转换可以用更准确的表达方式，如【图3】。

图 3
这些参数是由GPS导航消息提供的。
1.2.2 BDT（北斗导航卫星系统时间）
BDT（北斗导航卫星系统时间）是一个连续的时间系统，没有润秒。开始历元的时间是【UTC 2023年1月1号 00:00:00】。
北斗时间计算公式【图4】：

图 4
UTC和GPST时间转换同上面的GPS同样，只但是UTC参数来自与北斗导航信息中。
1.3 坐标系统
接受机和卫星的位置在RTKLIB中表达为在ECEF(地心地固坐标系)坐标系统中的X, Y, Z组件。
1.3.1 大地坐标到ECEF坐标的转换
转换公式如【图5】。

第三个公式最后一行有错有错，应当为: (v(1 – e2)+h)sin
图 5
参数说明： a : 地球参考椭球的长半径
f : 地球参考椭球的扁平率
h: 椭球高度
: 纬度
: 经度
当前版本的RTKLIB使用的值为【图6】：

图 6

图7 参考椭球体
1.3.2 ECEF坐系到大地坐标的转换
转换公式如【图8】

图 8
1.3.3 本地坐标到ECEF坐标的转换
在接受机位置的本地坐标，也被称为ENU坐标，通常使用在GNSS导航解决。ECEF坐标到本地坐标转换的旋转矩阵表达为【图9】。

图 9 Er旋转矩阵
参数说明：：接受机位置的纬度
：接受机位置的经度
通过使用Er和接受机的坐标rr【ECEF】，坐标recef【ECEF】可以被转换到本地坐标的坐标rlocal，公式如【图10】。

图 10
2 RTKPPP定位算法
2.1 单点定位（pntpos）
1:satposs à 2: estpos à 3: estvel
1. 计算计算卫星位置、速度和时钟（satposs）
a) 通过广播星历计算卫星钟差（ephclk）
Ø 卫星编号到卫星系统的转换（satsys）
根据卫星的编号，获取到相应的卫星导航系统。
Ø 选择星历（seleph）
1、传入信号传输时间，卫星编号，导航数据等参数。
2、遍历导航数据，遍历导航数据里面的星历数据，判断星历数据的卫星编号是否和传入的卫星编号相等。
3、假如星历数据的卫星编号和传入的卫星编号相等，就计算星历参考时间(toe)和信号传输时间的时间差。假如不相等，继续解决下一条星历数据。
4、判断计算出来的时间差，假如时间差大于了允许的最大时间差，继续查找下一个星历数据。否则，判断时间差最小的星历数据，记录星历数据的位置。
5、返回之前记录出来的星历数据。
Ø 使用广播星历计算卫星时钟偏差（eph2clk）
1、 传入信号发射时刻的时间和星历数据。
2、 计算信号发射时刻的时间和本时段钟差参数参考时间(星历参数toc)的时间差。
3、 通过下式计算钟差，这里还没有解决相对论校正项和tgd：（代码中有个迭代过程，资料上没看到写）

卫星钟差计算出来之后，信号发射时刻的时间还要加上这个钟差。
b) 计算卫星在信号发射时刻的位置、速度和时钟（satpos）
根据星历表选项来选择不同的解决，如下：
Ø 广播星历（EPHOPT_BRDC）: ephpos （广播星历到卫星位置和钟差）
1、 根据公式计算出tk；
2、 根据使用的卫星系统，选择使用的地球引力常数(mu)和地球的角速度(omge)
3、 根据公式计算出平近点角M。
4、 求解开普勒方程，按照以下公式迭代求解。
5、 根据以下公式计算出u(改正后的纬度幅角), r(改正后的径向), i(改正后的轨道倾角) 的值。

6、 根据以下公式计算卫星在轨道平面内的坐标。

7、 根据不同的卫星系统，做不同的计算。
GPS计算方式：

北斗计算方式：

其中
8、 按照公式计算出时间tc。
9、 按照以下公式计算出钟差和钟漂。

Ø 精密星历（EPHOPT_PREC）: peph2pos
Ø 广播+SBAS（EPHOPT_SBAS）：satpos_sbas
Ø 广播+SSR_APC（EPHOPT_SSRAPC）: satpos_ssr
Ø 广播+SSR_COM（EPHOPT_SSRCOM）：satpos_ssr
Ø QZSS LEX星历（EPHOPT_LEX）：lexeph2pos
2. 使用伪距估算接受机的位置,返回估算状态结果 （estpos）
a) 伪距残差（rescode）
Ø 把ecef坐标系转换成大地坐标系（ecef2pos）
1、 按照以下公式做转换，暂时还没看懂。

Ø 计算几何距离和接受机到卫星的单位矢量（geodist）
1、 用卫星的坐标向量做欧几里德范数，返回值和地球长半轴(WGS84)比较。小于地球长半轴(WGS84)，返回-1；
2、 计算卫星坐标和接受机坐标的差值向量。
3、 用差值向量做欧几里德范数，再用上一步计算出来的差值向量和计算结果做除法，得到视线向量。按照如下公式：

4、 使用以下公式计算几何距离。

Ø 计算卫星方位角/仰角（satazel）
1、 把接受机ecef坐标转换到大地坐标；
2、 判断高度是否大于地球半长轴（WGS84）的负数值；
3、 假如高度小于等于地球半长轴（WGS84）的负数值，方位角为0，仰角为PI/2;
4、 假如高度大于地球半长轴（WGS84）的负数值，把ECEF向量转换到局部坐标。然后对转换出来的坐标做内积。
5、 使用如下公式计算卫星方位角和仰角。

Ø 伪距使用编码残差改正（prange）
暂时没找到相应的文档相应。
Ø 电离层改正（ionocorr）
n 通过广播电离层模型（klobuchar模型）计算出电离层延迟（ionmodel）
1、 校验传入的电离层模型参数，校验失败，使用默认的电离层模型参数；
2、 使用以下公式计算出地球为中心的角度（半圆）；

3、 使用以下公式计算子的电离层的纬度/经度（半圆）；

4、 使用以下公式计算地磁纬度。

5、 计算本地时间，返回值按这个公式【tt-=floor(tt/86400.0)*86400.0;】解决，保证tt的范围（0<=tt<86400）。

6、 计算倾斜因子。

7、 电离层延迟计算，公式如下。

Ø 对流层改正（tropcorr）
n 通过标准大气压和saastamoinen模型计算对流层延迟 （tropmodel）
1、 使用以下公式计算总气压。

2、 使用以下公式计算绝对温度。

3、 使用以下公式计算水蒸汽的分压。代码中没有使用。

4、 使用以下公式计算【Saastamoinen模型】。

Ø 伪距残差
残差值=伪距-(物理距离+dtr-光速*时钟偏差+电离层误差+对流层误差)
Ø 时钟系统和接受器的偏置补偿
Ø 伪距测量误差方差（varerr）
b) 方差权重值 (weight by variance)
c) 最小二乘估计 （lsq）
Ø 最小二乘估计通过求解正规方程（X =（A * A'）^ - 1* A * Y）
1、 计算矩阵A*Y的结果，保存到矩阵Ay；
2、 计算矩阵A*A’的结果，保存到矩阵Q；
3、 求的矩阵Q的逆矩阵,结果保存到矩阵Q；
4、 最后x=Q*Ay
d) 欧几里德范数（norm）
公式：
1、 通过最小二乘法估算出来的参数dx【长度为4】，分别加到位置向量上x【长度为4】。
2、 对这个速度向量做欧几里德范数，返回值同1E-4比较，小于这个值时，就得到估算出的接受机的位置为向量x的值。
e) 验证求解 (valsol)
1、 对伪距残差值做内积，然后同卡方分布的自由度分布值(alpha=0.001)做比较，大于卡方分布的自由度分布值的数据，无效。
2、 计算dops； 暂时还没没找到资料。
3、 用计算出来的dops和配置的最大dops阈值比较，大于配置的最大dops阈值的数据无效。
3. 使用多普勒估算接受机速度（estvel）
a) 多普勒残差（resdop）
Ø 把ecef坐标系转换成大地坐标系（ecef2pos）
Ø 把xyz坐标转换成enu坐标（xyz2enu）
Ø 计算ecef中的瞄准线向量
Ø 计算相对于接受机在ECEF中的卫星速度
Ø 多普勒残差计算
1、速率公式：
b) 最小二乘估计 （lsq）
Ø 最小二乘估计通过求解正规方程（X =（A * A'）^ - 1* A * Y）
1、 计算矩阵A*Y的结果，保存到矩阵Ay；
2、 计算矩阵A*A’的结果，保存到矩阵Q；
3、 求的矩阵Q的逆矩阵,结果保存到矩阵Q；
4、 最后x=Q*Ay
c) 欧几里德范数（norm）
公式：
3、 通过最小二乘法估算出来的参数dx【长度为4】，分别加到速度向量上x【长度为4】。
4、 对这个速度向量做欧几里德范数，返回值同1E-6比较，小于这个值时，就得到估算出的接受机速度为向量x的值。
2.2 精拟定位 （pppos）
1:udstate_ppp à 2: satposs à 3: testeclipse à 4: res_ppp à 5: res_ppp à 6: filter à 7: res_ppp
1. 暂时更新状态（udstate_ppp）
a) 位置更新
b) 时钟更新
c) 对流层参数更新
d) 相位偏差更新
1、 通过LLI检测周跳；
2、 通过Geometry-Free相位跳变检测周跳(假如双频测量值可用)；
2. 计算卫星的位置和速度和时钟（satposs）
c) 通过广播星历计算卫星钟差（ephclk）
Ø 卫星编号到卫星系统的转换（satsys）
根据卫星的编号，获取到相应的卫星导航系统。
Ø 选择星历（seleph）
1、传入信号传输时间，卫星编号，导航数据等参数。
2、遍历导航数据，遍历导航数据里面的星历数据，判断星历数据的卫星编号是否和传入的卫星编号相等。
3、假如星历数据的卫星编号和传入的卫星编号相等，就计算星历参考时间(toe)和信号传输时间的时间差。假如不相等，继续解决下一条星历数据。
4、判断计算出来的时间差，假如时间差大于了允许的最大时间差，继续查找下一个星历数据。否则，判断时间差最小的星历数据，记录星历数据的位置。
5、返回之前记录出来的星历数据。
Ø 使用广播星历计算卫星时钟偏差（eph2clk）
1、 传入信号传输时间和星历数据。
2、 传入信号发射时刻的时间和星历数据。
3、 计算信号发射时刻的时间和本时段钟差参数参考时间(星历参数toc)的时间差。
4、 通过下式计算钟差，这里还没有解决相对论校正项和tgd：（代码中有个迭代过程，资料上没看到写）

d) 计算卫星位置、速度和时钟（satpos）
根据星历表选项来选择不同的解决，如下：
Ø 广播星历（EPHOPT_BRDC）: ephpos （广播星历到卫星位置和钟差）
1、 根据公式计算出tk；
2、 根据使用的卫星系统，选择使用的地球引力常数(mu)和地球的角速度(omge)
3、 根据公式计算出平近点角M。
4、 求解开普勒方程，按照以下公式迭代求解。
5、 根据以下公式计算出u(改正后的纬度幅角), r(改正后的径向), i(改正后的轨道倾角) 的值。

6、 根据以下公式计算卫星在轨道平面内的坐标。

7、 根据不同的卫星系统，做不同的计算。
GPS计算方式：

北斗计算方式：

其中
8、 按照公式计算出时间tc。
9、 按照以下公式计算出钟差和钟漂。

Ø 精密星历（EPHOPT_PREC）: peph2pos
Ø 广播+SBAS（EPHOPT_SBAS）：satpos_sbas
Ø 广播+SSR_APC（EPHOPT_SSRAPC）: satpos_ssr
Ø 广播+SSR_COM（EPHOPT_SSRCOM）：satpos_ssr
Ø QZSS LEX星历（EPHOPT_LEX）：lexeph2pos

3. 排除对被遮蔽卫星的测量（testeclipse）
a) 获得太阳和月亮在ecef中的位置（sunmoonpos）
b) 归一化三维矢量（normv3）
c) 欧几里德范数的向量（norm）
d) 计算太阳-地球-卫星角（dot）
4. 相位和编码残差（res_ppp）
a) 地球潮汐改正（tidedisp）
b) 计算几何距离（geodist）
1、 用卫星的坐标向量做欧几里德范数，返回值和地球长半轴(WGS84)比较。小于地球长半轴(WGS84)，返回-1；
2、 计算卫星坐标和接受机坐标的差值向量。
3、 用差值向量做欧几里德范数，再用上一步计算出来的差值向量和计算结果做除法，得到视线向量。
4、 使用以下公式计算几何距离。

备注：这个公式是地球自转改正公式。
地球自转改正产生地球自转改正的机制在于地固系随地球的自转而旋转, 地固系是非惯性系统, 用地固系中的坐标计算卫星到接受机几何距离时, 由于卫星位置和接受机位置是两个不同时刻的位置矢量, 并且这两个时刻的地固系相对于惯性系是变化的, 因此要考虑地球自转引起的距离变化。其改正公式为:

其中, w 为地球自转角速度; C 为光速; x si和x ri分别代表卫星位置矢量和测站位置矢量的分量, i= 1, 2 相应x , y 分量。地球自转对纬度影响甚小, 对经度影响最大, 另一方面是高度。同时地球的自转影响还跟测站的经纬度以及测站与卫星的几何关系有关, 对两极的测站, 影响为零; 对赤道上的测站影响最大; 当卫星在测站子午面内影响为零;卫星在测站东方时, 影响为负, 卫星在测站西方时,影响为正。
c) 卫星方位角/仰角计算（satazel）
1、 把接受机ecef坐标转换到大地坐标；
2、 判断高度是否大于地球半长轴（WGS84）的负数值；
3、 假如高度小于等于地球半长轴（WGS84）的负数值，方位角为0，仰角为PI/2;
4、 假如高度大于地球半长轴（WGS84）的负数值，把ECEF向量转换到局部坐标。然后对转换出来的坐标做内积。
5、 使用如下公式计算卫星方位角和仰角。

d) 对流层延迟修正：
Ø 通过标准大气压和saastamoinen模型计算对流层延迟 （tropmodel）
1、 使用以下公式计算总气压。

2、 使用以下公式计算绝对温度。

3、 使用以下公式计算水蒸汽的分压。代码中没有使用。

4、 使用以下公式计算【Saastamoinen模型】。

e) 卫星天线模型（satantpcv）
f) 接受机天线模型，通过天线相位中心参数计算天线偏移 （antmodel）
g) 相位缠绕校正（windupcorr）
h) 电离层和天线相位校正测量（corrmeas）
i) 计算卫星时钟和对流层延迟
5. 卡尔曼滤波（filter）
卡尔曼滤波按照以下的公式更新状态：
K = P * H * (H' * P * H + R)^-1 ;
xp = x + K * v ;
Pp = (I – K * H') * P ;
参数解释：
x: 状态向量 (n x 1)
p: 状态的协方差矩阵 (n x n)
H: 设计矩阵的转置矩阵 (n x m)
v: 创新(测量—模型) (m x 1)
R: 测量误差的协方差矩阵 (m x m)
n,m: 状态和测量值的数量
xp：状态更新后的向量 (n x 1)
Pp: 状态更新后的协方差矩阵 (n x n)
6. 组合残差 （res_ppp）
7. 解析PPP整周模糊度（pppamb）
a) average LC（average_LC）
b) 固定宽巷模糊（fix_amb_WL）
Ø 计算LC波长度（m）（lam_LC）
Ø 宽巷模糊
Ø 宽巷模糊度的方差
c) fix narrow-lane ambiguity 【AR mode: PPP-AR】
d) fix narrow-lane ambiguity by ILS 【AR mode: PPP-AR ILS】
3 以下是上面一些算法的参考资料，信息来自于Rtklib的手册
3.1 接受机和卫星天线之间的几何距离相应rtklib中的函数
extern double geodist(const double *rs, const double *rr, double *e)

3.2 卫星方向的方位角和仰角相应rtklib的函数
extern double satazel(const double *pos, const double *e, double *azel)

3.3 GPS、Galileo和QNSS广播星历和时钟相应rtklib的函数
extern void eph2pos(gtime_t time, const eph_t *eph, double *rs, double *dts, double *var)

3.4 北斗广播星历和时钟相应rtklib的函数
extern void eph2pos(gtime_t time, const eph_t *eph, double *rs, double *dts, double *var)

3.5 对流层和电流层模型
3.5.1 对流层模型相应rtklib的函数
extern double tropmodel(gtime_t time, const double *pos, const double *azel,double humi)

3.5.2 电离层模型相应rtklib的函数
extern double ionmodel(gtime_t t, const double *ion, const double *pos,const double *azel)

3.6 单点定位
3.6.1 线性最小二乘估计相应rtklib的函数
extern int lsq(const double *A, const double *y, int n, int m, double *x,double *Q)

3.6.2 估算接受机的位置和钟差
static int estpos(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)

3.6.3 估算接受机的速度和钟漂
static void estvel(const obsd_t *obs, int n, const double *rs, const double *dts, const nav_t *nav, const prcopt_t *opt, sol_t *sol, const double *azel, const int *vsat)

3.6.4 求解验证和RAIM FDE
static int valsol(const double *azel, const int *vsat, int n, const prcopt_t *opt, const double *v, int nv, int nx, char *msg)

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)

3.7 精密定位
3.7.1 用于PPP的ZD测量模式

3.7.2 接受机天线相位中心模型

3.7.3 卫星天线相位中心模型

3.7.4 由地球潮汐引起的位移

3.7.5 相位缠绕校正

3.7.6 通过PPP估算接受机的位置