RTKPPP定位算法流程.docx
《RTKPPP定位算法流程.docx》由会员分享，可在线阅读，更多相关《RTKPPP定位算法流程.docx（34页珍藏版）》请在冰豆网上搜索。
RTKPPP定位算法流程
1基础知识
1.1GPS精密单点定位得基本原理
GPS精密单点定位一般采用单台双频GPS接收机,利用IGS提供得精密星历与卫星钟差,基于载波相位观测值进行得高精度定位。
观测值中得电离层延迟误差通过双频信号组合消除,对流层延迟误差通过引入未知参数进行估计。
1.2时间系统
RTKLIB内部使用GPST（GPST时间）用于GNSS得数据处理与定位算法。
数据在RTKLIB内部处理之前,需要转换成GPST时间。
使用GPST得原因就是避免处理润秒。
RTKLIB使用以下结构体表示时间:
typedefstruct
time_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】:
图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.2ECEF坐系到大地坐标得转换
转换公式如【图8】
图8
1.3.3本地坐标到ECEF坐标得转换
在接收机位置得本地坐标,也被称为ENU坐标,通常使用在GNSS导航处理。
ECEF坐标到本地坐标转换得旋转矩阵表示为【图9】。
图9Er旋转矩阵
参数说明:
:
接收机位置得纬度
:
接收机位置得经度
通过使用Er与接收机得坐标rr【ECEF】,坐标recef【ECEF】可以被转换到本地坐标得坐标rlocal,公式如【图10】。
图10
2RTKPPP定位算法
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_（EPHOPT_SSR）:
satpos_ssr
ØQZSSLEX星历（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）
⏹通过广播电离层模型（klobuchar模型）计算出电离层延迟（ionmodel）
1、校验传入得电离层模型参数,校验失败,使用默认得电离层模型参数;
2、使用以下公式计算出地球为中心得角度（半圆）;
3、使用以下公式计算子得电离层得纬度/经度（半圆）;
4、使用以下公式计算地磁纬度。
5、计算本地时间,返回值按这个公式【tt-=floor（tt/86400、0）*86400、0;】处理,保证tt得范围（0<=tt<86400）。
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*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）
4、传入信号传输时间与星历数据。
5、传入信号发射时刻得时间与星历数据。
6、计算信号发射时刻得时间与本时段钟差参数参考时间（星历参数toc）得时间差。
7、通过下式计算钟差,这里还没有处理相对论校正项与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_（EPHOPT_SSR）:
satpos_ssr
ØQZSSLEX星历（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为光速;xsi与xri分别代表卫星位置矢量与测站位置矢量得分量,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:
状态向量（nx1）
p:
状态得协方差矩阵（nxn）
H:
设计矩阵得转置矩阵（nxm）
v:
创新（测量—模型）（mx1）
R:
测量误差得协方差矩阵（mxm）
n,m:
状态与测量值得数量
xp:
状态更新后得向量（nx1）
Pp:
状态更新后得协方差矩阵（nxn）
6.组合残差（res_ppp）
7.解析PPP整周模糊度（pppamb）
a）averageLC（average_LC）
b）固定宽巷模糊（fix_amb_WL）
Ø计算LC波长度（m）（lam_LC）
Ø宽巷模糊
Ø宽巷模糊度得方差
c）fixnarrow-laneambiguity【ARmode:
PPP-AR】
d）fixnarrow-laneambiguitybyILS【ARmode:
PPP-ARILS】
3
以下就是上面一些算法得参考资料,信息来自于Rtklib得手册
3.1接收机与卫星天线之间得几何距离
3.2卫星方向得方位角与仰角
3.3GPS、Galileo与QNSS广播星历与时钟
3.4北斗广播星历与时钟
3.5对流层与电流层模型
3.5.1对流层模型
3.5.2电离层模型
3.6单点定位
3.6.1线性最小二乘估计
3.6.2估算接收机得位置与钟差
3.6.3估算接收机得速度与钟漂
3.6.4求解验证与RAIMFDE
3.7精密定位
3.7.1用于PPP得ZD测量模式
3.7.2接收机天线相位中心模型
3.7.3卫星天线相位中心模型
3.7.4由地球潮汐引起得位移
3.7.5相位缠绕校正
3.7.6通过PPP估算接收机得位置
升级会员 