RTKPPP定位算法流程

1基础知识1.1GPS精密单点定位的基本原理GPS精密单点定位一般采用单台双频GPS接收机，利用IGS提供的精密星历和卫星钟差，基于载波相位观测值进行的高精度定位。观测值中的电离层延迟误差通过双频信号组合消除，对流层延迟误差通过引入未知参数进行估计。1.2时间系统RTKLIB内部使用GPST(GPST时间)用于GNSS的数据处理和定位算法。数据在RTKLIB内部处理之前，需要转换成GPST时间。使用GPST的原因是避免处理润秒。RTKLIB使用以下结构体表示时间：typedefstructtime_ttime;/*time(s)expressedbystandardtime_t*/doublesec;/*fractionofsecondunder1s*/}gtime_t;1.2.1GPST和UTC(UniversalTimeCoordinated)关系参考【图1】，参考【图2】：图1转换关系公式图2通过使用GPS导航信息中的UTC参数，GPST到UTC或者UTC到GPST之前的转换可以用更准确的表达方式，如【图3】。图3这些参数是由GPS导航消息提供的。1.2.2BDT（北斗导航卫星系统时间）BDT（北斗导航卫星系统时间）是一个连续的时间系统，没有润秒。开始历元的时间是【UTC2006年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坐系到大地坐标的转换转换公式如【图8】图81.3.3本地坐标到ECEF坐标的转换在接收机位置的本地坐标，也被称为ENU坐标，通常使用在GNSS导航处理。ECEF坐标到本地坐标转换的旋转矩阵表示为【图9】。图9Er旋转矩阵参数说明：：接收机位置的纬度：接收机位置的经度通过使用Er和接收机的坐标rr【ECEF】，坐标recef【ECEF】可以被转换到本地坐标的坐标rlocal，公式如【图10】。图102RTKPPP定位算法2.1单点定位（pntpos）1:satposs2:estpos3:estvel1.计算计算卫星位置、速度和时钟（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QZSSLEX星历（EPHOPT_LEX）：lexeph2pos2.使用伪距估算接收机的位置,返回估算状态结果（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）通过广播电离层模型（klobuchar模型）计算出电离层延迟（ionmodel）1、校验传入的电离层模型参数，校验失败，使用默认的电离层模型参数；2、使用以下公式计算出地球为中心的角度（半圆）；3、使用以下公式计算子的电离层的纬度/经度（半圆）；4、使用以下公式计算地磁纬度。5、计算本地时间，返回值按这个公式【tt-=floor(tt/86400.0)*86400.0;】处理，保证tt的范围（0=tt86400）。6、计算倾斜因子。7、电离层延迟计算，公式如下。对流层改正（tropcorr）通过标准大气压和saastamoinen模型计算对流层延迟（tropmodel）1、使用以下公式计算总气压。2、使用以下公式计算绝对温度。3、使用以下公式计算水蒸汽的分压。代码中没有使用。4、使用以下公式计算【Saastamoinen模型】。伪距残差残差值=伪距-(物理距离+dtr-光速*时钟偏差+电离层误差+对流层误差)时钟系统和接收器的偏置补偿伪距测量误差方差（varerr）b)方差权重值(weightbyvariance)c)最小二乘估计（lsq）最小二乘估计通过求解正规方程（X=（A*A'）^-1*A*Y）1、计算矩阵A*Y的结果，保存到矩阵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、用计算出来的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*Ayc)欧几里德范数（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_ppp1.暂时更新状态（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QZSSLEX星历（EPHOPT_LEX）：lexeph2pos3.排除对被遮蔽卫星的测量（testeclipse）a)获得太阳和月亮在ecef中的位置（sunmoonpos）b)归一化三维矢量（normv3）c)欧几里德范数的向量（norm）d)计算太阳-地球-卫星角（dot）4.相位和编码残差（res_ppp）a)地