RTKPPP定位算法流程

1、1 基础知识1.1GPS 精密单点定位的基本原理GPS精密单点定位一般采用单台双频GPS接收机，利用 IGS 提供的精密星历和卫星钟差，基于载波相位观测值进行的高精度定位。观测值中的电离层延迟误差通过双频信号组合消除，对流层延迟误差通过引入未知参数进行估计。1.2时间系统RTKLIB 内部使用 GPST(GPST时间 )用于 GNSS的数据处理和定位算法。数据在RTKLIB内部处理之前，需要转换成GPST时间。使用 GPST的原因是避免处理润秒。RTKLIB 使用以下结构体表示时间：typedef structtime_t time;/* time(s) expressed by standa

2、rd time_t */double sec;/* fraction of second under 1 s */ gtime_t;1.2.1GPST 和 UTC(Universal Time Coordinated)关系参考【图1】，参考【图 2】：图 1 转换关系公式图 2通过使用 GPS导航信息中的UTC 参数， GPST到 UTC 或者 UTC 到 GPST之前的转换可以用更准确的表达方式，如【图3】。图 3这些参数是由 GPS导航消息提供的。1.2.2BDT （北斗导航卫星系统时间）BDT （北斗导航卫星系统时间）是一个连续的时间系统，没有润秒。开始历元的时间是【UTC 2006 年

3、 1月 1号 00:00:00 】。北斗时间计算公式【图4】：图 4UTC 和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.2ECEF

4、坐系到大地坐标的转换转换公式如【图8】图 81.3.3本地坐标到 ECEF 坐标的转换在接收机位置的本地坐标，也被称为ENU 坐标，通常使用在GNSS导航处理。 ECEF坐标到本地坐标转换的旋转矩阵表示为【图9】。图 9 Er 旋转矩阵参数说明：接收机位置的纬度：接收机位置的经度通过使用 Er 和接收机的坐标r r【 ECEF】，坐标 recef【 ECEF】可以被转换到本地坐标的坐标 rlocal ，公式如【图 10】。图 102RTKPPP定位算法2.1单点定位（ pntpos ）1:satposs2: estpos3: estvel1. 计算计算卫星位置、速度和时钟（satposs）a)

5、 通过广播星历计算卫星钟差（ ephclk ）? 卫星编号到卫星系统的转换（ satsys）根据卫星的编号，获取到对应的卫星导航系统。? 选择星历（ seleph）1、传入信号传输时间，卫星编号，导航数据等参数。2、遍历导航数据，遍历导航数据里面的星历数据，判断星历数据的卫星编号是否和传入的卫星编号相等。3、如果星历数据的卫星编号和传入的卫星编号相等，就计算星历参考时间(toe)和信号传输时间的时间差。如果不相等， 继续处理下一条星历数据。4、判断计算出来的时间差，如果时间差大于了允许的最大时间差，继续查找下一个星历数据。否则，判断时间差最小的星历数据，记录星历数据的位置。5、返回之前记录出来

6、的星历数据。?使用广播星历计算卫星时钟偏差（eph2clk）1、 传入信号发射时刻的时间和星历数据。2、 计算信号发射时刻的时间和本时段钟差参数参考时间的时间差。3、 通过下式计算钟差，这里还没有处理相对论校正项和有个迭代过程，资料上没看到写）(星历参数 toc)tgd：（代码中卫星钟差计算出来之后，信号发射时刻的时间还要加上这个钟差。b) 计算卫星在信号发射时刻的位置、速度和时钟（ satpos）根据星历表选项来选择不同的处理，如下：?广播星历（ EPHOPT_BRDC） : ephpos （广播星历到卫星位置和钟差）1、 根据公式计算出 tk；2、 根据使用的卫星系统，选择使用的地球引力常

7、数(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_SS

8、RAPC）: satpos_ssr? 广播 +SSR_COM（ EPHOPT_SSRCOM）： satpos_ssr? QZSS LEX星历（ EPHOPT_LEX）： lexeph2pos2. 使用伪距估算接收机的位置 ,返回估算状态结果 （estpos）a) 伪距残差（ rescode）?把 ecef 坐标系转换成大地坐标系（ecef2pos）1、 按照以下公式做转换，暂时还没看懂。? 计算几何距离和接收机到卫星的单位矢量（geodist ）1、 用卫星的坐标向量做欧几里德范数，返回值和地球长半轴(WGS84) 比较。小于地球长半轴(WGS84) ，返回 -1；2、 计算卫星坐标和接收机坐

9、标的差值向量。3、 用差值向量做欧几里德范数， 再用上一步计算出来的差值向量和计算结果做除法，得到视线向量。按照如下公式：4、 使用以下公式计算几何距离。? 计算卫星方位角 / 仰角（ satazel）1、 把接收机 ecef坐标转换到大地坐标；2、 判断高度是否大于地球半长轴（WGS84 ）的负数值；3、 如果高度小于等于地球半长轴（WGS84 ）的负数值，方位角为0，仰角为 PI/2;4、 如果高度大于地球半长轴（WGS84 ）的负数值，把ECEF向量转换到局部坐标。然后对转换出来的坐标做内积。5、 使用如下公式计算卫星方位角和仰角。? 伪距使用编码残差改正（ prange）暂时没找到对应

10、的文档对应。? 电离层改正（ ionocorr ）通过广播电离层模型 （ klobuchar 模型）计算出电离层延迟 （ ionmodel ）1、 校验传入的电离层模型参数，校验失败， 使用默认的电离层模型参数；2、 使用以下公式计算出地球为中心的角度（半圆）；3、 使用以下公式计算子的电离层的纬度/经度（半圆）；4、 使用以下公式计算地磁纬度。5、 计算本地时间， 返回值按这个公式【 tt-=floor(tt/86400.0)*86400.0; 】处理，保证 tt的范围（ 0=tt86400 ）。6、 计算倾斜因子。7、 电离层延迟计算，公式如下。? 对流层改正（ tropcorr ）通过标

11、准大气压和saastamoinen 模型计算对流层延迟（ tropmodel ）1、 使用以下公式计算总气压。2、 使用以下公式计算绝对温度。3、 使用以下公式计算水蒸汽的分压。代码中没有使用。4、使用以下公式计算【Saastamoinen模型】。? 伪距残差残差值 =伪距 -(物理距离 +dtr- 光速 *时钟偏差 +电离层误差 +对流层误差 )?时钟系统和接收器的偏置补偿?伪距测量误差方差（ varerr ）b)方差权重值(weight by variance)c)最小二乘估计（ lsq）?最小二乘估计通过求解正规方程（X =（A * A ） - 1* A * Y ）1、 计算矩阵 A*Y

12、 的结果，保存到矩阵Ay ；2、 计算矩阵 A*A 的结果，保存到矩阵Q；3、 求的矩阵 Q的逆矩阵 ,结果保存到矩阵Q；4、 最后 x=Q*Ayd) 欧几里德范数（ norm ）公式 ：1、 通过最小二乘法估算出来的参数dx【长度为 4】，分别加到位置向量上x【长度为 4】。2、 对这个速度向量做欧几里德范数，返回值同1E-4比较，小于这个值时，就得到估算出的接收机的位置为向量x的值。e) 验证求解 (valsol)1、 对伪距残差值做内积，然后同卡方分布的自由度分布值(alpha=0.001) 做比较，大于卡方分布的自由度分布值的数据，无效。2、 计算 dops； 暂时还没没找到资料。3、

13、 用计算出来的dops和配置的最大dops阈值比较，大于配置的最大dops阈值的数据无效。3. 使用多普勒估算接收机速度（ estvel ）a) 多普勒残差（ resdop）?把 ecef 坐标系转换成大地坐标系（ecef2pos）?把 xyz 坐标转换成enu 坐标（ xyz2enu）?计算 ecef 中的瞄准线向量?计算相对于接收机在ECEF中的卫星速度?多普勒残差计算1、速率公式：b)最小二乘估计（ lsq）?最小二乘估计通过求解正规方程（1、 计算矩阵 A*Y 的结果，保存到矩阵X =（A * A ） - 1* A * Y ）Ay；2、 计算矩阵 A*A的结果，保存到矩阵Q；3、 求的

14、矩阵 Q 的逆矩阵 ,结果保存到矩阵Q；4、 最后 x=Q*Ayc) 欧几里德范数（ norm ）公式 ：3、 通过最小二乘法估算出来的参数dx【长度为 4】，分别加到速度向量上x【长度为 4】。4、 对这个速度向量做欧几里德范数，返回值同 1E-6比较，小于这个值时，就得到估算出的接收机速度为向量 x的值。2.2精确定位（ pppos ）1:udstate_ppp2: satposs3: testeclipse4: res_ppp5: res_ppp6: filter7:res_ppp1. 暂时更新状态（ udstate_ppp ）a) 位置更新b) 时钟更新c) 对流层参数更新d) 相位偏

15、差更新1、 通过 LLI 检测周跳；2、 通过 Geometry-Free 相位跳变检测周跳(如果双频测量值可用) ；2.计算卫星的位置和速度和时钟（satposs）c) 通过广播星历计算卫星钟差（ ephclk ）? 卫星编号到卫星系统的转换（ satsys）根据卫星的编号，获取到对应的卫星导航系统。? 选择星历（ seleph）1、传入信号传输时间，卫星编号，导航数据等参数。2、遍历导航数据，遍历导航数据里面的星历数据，判断星历数据的卫星编号是否和传入的卫星编号相等。3、如果星历数据的卫星编号和传入的卫星编号相等，就计算星历参考时间(toe)和信号传输时间的时间差。如果不相等，继续处理下一

16、条星历数据。4、判断计算出来的时间差，如果时间差大于了允许的最大时间差，继续查找下一个星历数据。否则，判断时间差最小的星历数据，记录星历数据的位置。5、返回之前记录出来的星历数据。?使用广播星历计算卫星时钟偏差（eph2clk）1、 传入信号传输时间和星历数据。2、 传入信号发射时刻的时间和星历数据。3、 计算信号发射时刻的时间和本时段钟差参数参考时间的时间差。4、 通过下式计算钟差，这里还没有处理相对论校正项和有个迭代过程，资料上没看到写）(星历参数 toc)tgd：（代码中d) 计算卫星位置、速度和时钟（ satpos）根据星历表选项来选择不同的处理，如下：?广播星历（ EPHOPT_BR

17、DC） : ephpos （广播星历到卫星位置和钟差）1、 根据公式计算出 tk；2、 根据使用的卫星系统，选择使用的地球引力常数(mu)和地球的角速度(omge)3、 根据公式计算出平近点角M 。4、 求解开普勒方程，按照以下公式迭代求解。5、 根据以下公式计算出 u( 改正后的纬度幅角 ), r( 改正后的径向 ), i(改正后的轨道倾角 ) 的值。6、 根据以下公式计算卫星在轨道平面内的坐标。7、 根据不同的卫星系统，做不同的计算。GPS计算方式：北斗计算方式：其中8、 按照公式计算出时间 tc。9、 按照以下公式计算出钟差和钟漂。? 精密星历（ EPHOPT_PREC）: peph2p

18、os? 广播 +SBAS（ EPHOPT_SBAS）： satpos_sbas? 广播 +SSR_APC（ EPHOPT_SSRAPC）: satpos_ssr? 广播 +SSR_COM（ EPHOPT_SSRCOM）： satpos_ssr? QZSS LEX星历（ EPHOPT_LEX）： lexeph2pos3. 排除对被遮蔽卫星的测量（ testeclipse）a) 获得太阳和月亮在 ecef 中的位置（ sunmoonpos ）b) 归一化三维矢量（ normv3 ）c) 欧几里德范数的向量（ norm ）d) 计算太阳 -地球 -卫星角（ dot ）4. 相位和编码残差（ res_

19、ppp）a) 地球潮汐改正（ tidedisp ）b) 计算几何距离（ geodist ）1、 用卫星的坐标向量做欧几里德范数，返回值和地球长半轴(WGS84) 比较。小于地球长半轴(WGS84) ，返回 -1；2、 计算卫星坐标和接收机坐标的差值向量。3、 用差值向量做欧几里德范数，再用上一步计算出来的差值向量和计算结果做除法，得到视线向量。4、 使用以下公式计算几何距离。备注：这个公式是地球自转改正公式。地球自转改正产生地球自转改正的机制在于地固系随地球的自转而旋转, 地固系是非惯性系统, 用地固系中的坐标计算卫星到接收机几何距离时, 由于卫星位置和接收机位置是两个不同时刻的位置矢量, 而

20、且这两个时刻的地固系相对于惯性系是变化的, 因此要考虑地球自转引起的距离变化。其改正公式为 :其中 , w 为地球自转角速度; C 为光速 ; x si和 x ri 分别代表卫星位置矢量和测站位置矢量的分量, i= 1, 2 对应 x , y 分量。地球自转对纬度影响甚小, 对经度影响最大 , 其次是高度。 同时地球的自转影响还跟测站的经纬度以及测站与卫星的几何关系有关 , 对两极的测站, 影响为零 ; 对赤道上的测站影响最大; 当卫星在测站子午面内影响为零;卫星在测站东方时, 影响为负 , 卫星在测站西方时,影响为正。c)卫星方位角 / 仰角计算（ satazel）1、 把接收机 ecef坐

21、标转换到大地坐标；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;x