RTKPPP定位算法流程.docx
《RTKPPP定位算法流程.docx》由会员分享，可在线阅读，更多相关《RTKPPP定位算法流程.docx（42页珍藏版）》请在冰点文库上搜索。
RTKPPP定位算法流程
1基础知识
2
2.1GPＳ周密单点定位的基本原理
2.2
ＧPS周密单点定位一般采纳单台双频GPS接收机,利用ＩGS提供的周密星历与卫星钟差,基于载波相位观测值进行的高精度定位。
观测值中的电离层延迟误差通过双频信号组合消除,对流层延迟误差通过引入未知参数进行估计。
2.3时间系统
2.4
ＲTＫＬＩＢ内部使用GPSＴ（GＰＳＴ时间）用于GNSS的数据处理与定位算法、数据在RTＫＬIＢ内部处理之前,需要转换成GＰSＴ时间。
使用GPＳT的缘故是幸免处理润秒、RTKLＩB使用以下结构体表示时间:
ｔyｐｅdeｆｓtrｕct
time_ｔ ｔime; /＊ time（s） expｒessedｂyｓtanｄard time_t*/
ｄｏｕbｌe sec; /*fraｃtｉonofｓecoｎduｎdeｒ１s*／
｝gtｉｍe_ｔ;
2.4.1GPＳT与ＵTC（ＵｎiｖｅｒsalTｉmｅ Ｃoordinａted）
2.4.2
关系参考【图1】,参考【图2】:
图1转换关系公式
图２
ﻩ通过使用GＰS导航信息中的UTＣ参数,ＧPST到ＵＴＣ或者UTC到GPＳＴ之前的转换能够用更准确的表达方式,如【图3】、
图3
这些参数是由GPS导航消息提供的、
2.4.3BＤＴ（北斗导航卫星系统时间）
2.4.4
ＢＤＴ（北斗导航卫星系统时间）是一个连续的时间系统,没有润秒、开始历元的时间是【ＵTＣ2０06年1月１号0０:
０0:
00】、
北斗时间计算公式【图4】:
图 4
ﻩUTＣ与GPＳT时间转换同上面的GＰS一样,只只是UTC参数来自与北斗导航信息中、
2.5坐标系统
2.6
接收机与卫星的位置在ＲTＫＬＩB中表示为在EＣEF（地心地固坐标系）坐标系统中的X,Y, Z组件、
2.6.1大地坐标到ＥＣＥF坐标的转换
2.6.2
转换公式如【图5】。
第三个公式最后一行有错,应该为:
（v（１ –e2）+h）sin
图５
参数讲明:
a:
地球参考椭球的长半径
f:
地球参考椭球的扁平率
ｈ:
椭球高度
:
纬度
:
经度
当前版本的RTKLIＢ使用的值为【图６】:
图6
图7参考椭球体
2.6.3ECEF坐系到大地坐标的转换
2.6.4
转换公式如【图8】
图 8
2.6.5本地坐标到EＣＥＦ坐标的转换
2.6.6
在接收机位置的本地坐标,也被称为ENＵ坐标,通常使用在GＮSＳ导航处理。
ECEF坐标到本地坐标转换的旋转矩阵表示为【图9】。
图9Er旋转矩阵
参数讲明:
:
接收机位置的纬度
:
接收机位置的经度
ﻩ通过使用Ｅr与接收机的坐标rr【EＣＥF】,坐标recef【ECＥF】能够被转换到本地坐标的坐标rlocal,公式如【图１0】。
图
10
3RTKＰPP定位算法
4
4.1单点定位（pntpos）
4.2
1:
ｓａtpoｓs → 2:
eｓｔpos → 3:
estvｅl
1.计算计算卫星位置、速度与时钟（ｓaｔｐｏss）
2.
a）通过广播星历计算卫星钟差（epｈｃｌk）
b）
Ø卫星编号到卫星系统的转换（ｓaｔsys）
Ø
依照卫星的编号,获取到对应的卫星导航系统、
Ø选择星历（ｓｅｌｅph）
Ø
1、传入信号传输时间,卫星编号,导航数据等参数、
2、遍历导航数据,遍历导航数据里面的星历数据,判断星历数据的卫星编号是否与传入的卫星编号相等。
3、假如星历数据的卫星编号与传入的卫星编号相等,就计算星历参考时间（ｔoｅ）与信号传输时间的时间差。
假如不相等,接着处理下一条星历数据、
４、判断计算出来的时间差,假如时间差大于了允许的最大时间差,接着查找下一个星历数据、否则,判断时间差最小的星历数据,记录星历数据的位置、
5、返回之前记录出来的星历数据、
Ø使用广播星历计算卫星时钟偏差（ｅph2clk）
Ø
1、传入信号发射时刻的时间与星历数据、
2、
3、计算信号发射时刻的时间与本时段钟差参数参考时间（星历参数ｔoc）的时间差、
4、
5、通过下式计算钟差,这个地方还没有处理相对论校正项与tgd:
（代码中有个迭代过程,资料上没看到写）
6、
卫星钟差计算出来之后,信号发射时刻的时间还要加上这个钟差。
c）计算卫星在信号发射时刻的位置、速度与时钟（sａtpoｓ）
d）
依照星历表选项来选择不同的处理,如下:
Ø广播星历（EＰHOPT_ＢＲDC）:
ｅｐhpoｓ（广播星历到卫星位置与钟差）
Ø
1、依照公式计算出tk;
2、
3、依照使用的卫星系统,选择使用的地球引力常数（ｍu）与地球的角速度（omｇｅ）
4、
5、依照公式计算出平近点角M。
6、
7、求解开普勒方程,依照以下公式迭代求解、
8、
9、依照以下公式计算出u（改正后的纬度幅角）,ｒ（改正后的径向）, i（改正后的轨道倾角）的值、
10、
11、依照以下公式计算卫星在轨道平面内的坐标。
12、
13、依照不同的卫星系统,做不同的计算、
14、
ＧPＳ计算方式:
北斗计算方式:
其中
15、依照公式计算出时间tc。
16、
17、依照以下公式计算出钟差与钟漂。
18、
Ø周密星历（EPHＯPT_PRＥＣ）:
pｅph2pｏｓ
Ø
Ø广播+SＢＡS（EＰHＯPT＿ＳＢAS）:
sａｔpos_sｂａｓ
Ø
Ø广播＋ＳSＲ_ＡPC（EＰHOＰT＿SSＲAPC）:
ｓatpos_ｓｓｒ
Ø
Ø广播+SＳR_（EＰHOＰＴ＿SSＲ）:
sａｔpｏs_ssr
Ø
ØQＺSSLＥX星历（EＰHOPT_LEX）:
lexｅpｈ2ｐos
Ø
3.使用伪距估算接收机的位置,返回估算状态结果 （eｓtpos）
4.
a）伪距残差（rｅscｏde）
b）
Ø把ecｅf坐标系转换成大地坐标系（eceｆ2poｓ）
Ø
1、依照以下公式做转换,暂时还没看明白、
2、
Ø计算几何距离与接收机到卫星的单位矢量（geodisｔ）
Ø
1、用卫星的坐标向量做欧几里德范数,返回值与地球长半轴（WGS8４）比较。
小于地球长半轴（WＧＳ84）,返回—1;
2、
3、计算卫星坐标与接收机坐标的差值向量。
4、
5、用差值向量做欧几里德范数,再用上一步计算出来的差值向量与计算结果做除法,得到视线向量。
依照如下公式:
6、
7、使用以下公式计算几何距离。
8、
Ø计算卫星方位角／仰角（sataｚel）
Ø
1、把接收机ｅcｅｆ坐标转换到大地坐标;
2、
3、判断高度是否大于地球半长轴（ＷGS８４）的负数值;
4、
5、假如高度小于等于地球半长轴（WGＳ８4）的负数值,方位角为0,仰角为ＰＩ/２;
6、
7、假如高度大于地球半长轴（WGS84）的负数值,把EＣＥF向量转换到局部坐标。
然后对转换出来的坐标做内积、
8、
9、使用如下公式计算卫星方位角与仰角。
10、
Ø伪距使用编码残差改正（prange）
Ø
暂时没找到对应的文档对应、
Ø电离层改正（ioｎocoｒｒ）
Ø
⏹通过广播电离层模型（klｏbｕｃｈａr模型）计算出电离层延迟（ｉonmodel）
⏹
1、校验传入的电离层模型参数,校验失败,使用默认的电离层模型参数;
2、
3、使用以下公式计算出地球为中心的角度（半圆）;
4、
5、使用以下公式计算子的电离层的纬度/经度（半圆）;
6、
7、使用以下公式计算地磁纬度。
8、
9、计算本地时间,返回值按这个公式【ｔt-=flｏoｒ（ｔt/86４0０、0）*8６400、０;】处理,保证tt的范围（0＜=ｔt<86４00）。
10、
11、计算倾斜因子、
12、
13、电离层延迟计算,公式如下、
14、
Ø对流层改正（tｒopcorr）
Ø
⏹通过标准大气压与ｓaastaｍｏｉｎｅｎ模型计算对流层延迟 （ｔｒｏｐmoｄel）
⏹
1、使用以下公式计算总气压。
2、
3、使用以下公式计算绝对温度、
4、
5、使用以下公式计算水蒸汽的分压。
代码中没有使用、
6、
7、使用以下公式计算【Ｓaastａｍoｉｎen模型】。
8、
Ø伪距残差
Ø
残差值=伪距—（物理距离＋ｄtr—光速＊时钟偏差+电离层误差+对流层误差）
Ø时钟系统与接收器的偏置补偿
Ø
Ø伪距测量误差方差（varｅrr）
Ø
c）方差权重值（weｉghtbｙvariaｎｃe）
d）
e）最小二乘估计（ｌｓq）
f）
Ø最小二乘估计通过求解正规方程（X=（Ａ*Ａ’）＾ -1*A*Ｙ）
Ø
1、计算矩阵Ａ*Ｙ的结果,保存到矩阵Aｙ;
2、
3、计算矩阵A＊A’的结果,保存到矩阵Ｑ;
4、
5、求的矩阵Q的逆矩阵,结果保存到矩阵Q;
6、
7、最后x=Ｑ＊Ａｙ
8、
g）欧几里德范数（ｎorm）
h）
公式:
1、通过最小二乘法估算出来的参数dｘ【长度为4】,分不加到位置向量上x【长度为４】。
2、
3、对这个速度向量做欧几里德范数,返回值同１E—4比较,小于这个值时,就得到估算出的接收机的位置为向量x的值。
4、
i）验证求解 （ｖalsｏl）
j）
1、对伪距残差值做内积,然后同卡方分布的自由度分布值（ａｌpha=０、0０1）做比较,大于卡方分布的自由度分布值的数据,无效、
2、
3、计算dops;暂时还没没找到资料。
4、
5、用计算出来的dops与配置的最大dｏps阈值比较,大于配置的最大ｄoｐs阈值的数据无效、
6、
5.使用多普勒估算接收机速度（eｓｔｖｅl）
6.
a）多普勒残差（resｄoｐ）
b）
Ø把ｅｃef坐标系转换成大地坐标系（ecｅｆ2ｐos）
Ø
Ø把xyz坐标转换成enu坐标（xyｚ2enu）
Ø
Ø计算ecｅｆ中的瞄准线向量
Ø
Ø计算相关于接收机在ECＥF中的卫星速度
Ø
Ø多普勒残差计算
Ø
１、速率公式:
c）最小二乘估计（lsq）
d）
Ø最小二乘估计通过求解正规方程（X =（A ＊ A’）^ -１* Ａ＊Y）
Ø
1、计算矩阵A＊Y的结果,保存到矩阵Aｙ;
2、
3、计算矩阵A＊A'的结果,保存到矩阵Ｑ;
4、
5、求的矩阵Ｑ的逆矩阵,结果保存到矩阵Q;
6、
7、最后x=Ｑ*Ａy
8、
e）欧几里德范数（ｎorｍ）
f）
公式:
5、通过最小二乘法估算出来的参数dｘ【长度为4】,分不加到速度向量上x【长度为4】。
6、
7、对这个速度向量做欧几里德范数,返回值同1Ｅ—6比较,小于这个值时,就得到估算出的接收机速度为向量x的值。
8、
4.3精确定位（ｐｐpos）
4.4
1:
udstａtｅ＿ppp →2:
ｓａｔposs → 3:
ｔeｓｔeclｉｐse→４:
res_pｐp→５:
ｒｅｓ_ｐpp→ 6:
ｆiｌｔer → 7:
res_pｐp
1.暂时更新状态（ｕｄsｔａte_ｐpp）
2.
a）位置更新
b）
c）时钟更新
d）
e）对流层参数更新
f）
g）相位偏差更新
h）
1、通过LＬI检测周跳;
2、
3、通过Geomeｔry—Frｅe相位跳变检测周跳（假如双频测量值可用）;
4、
3.计算卫星的位置与速度与时钟（satpｏss）
4.
e）通过广播星历计算卫星钟差（eｐｈｃｌｋ）
f）
Ø卫星编号到卫星系统的转换（satsｙｓ）
Ø
依照卫星的编号,获取到对应的卫星导航系统。
Ø选择星历（sｅlｅpｈ）
Ø
1、传入信号传输时间,卫星编号,导航数据等参数、
2、遍历导航数据,遍历导航数据里面的星历数据,判断星历数据的卫星编号是否与传入的卫星编号相等。
3、假如星历数据的卫星编号与传入的卫星编号相等,就计算星历参考时间（tｏe）与信号传输时间的时间差、假如不相等,接着处理下一条星历数据。
4、判断计算出来的时间差,假如时间差大于了允许的最大时间差,接着查找下一个星历数据。
否则,判断时间差最小的星历数据,记录星历数据的位置、
５、返回之前记录出来的星历数据、
Ø使用广播星历计算卫星时钟偏差（ｅpｈ２clｋ）
Ø
7、传入信号传输时间与星历数据、
8、
9、传入信号发射时刻的时间与星历数据。
10、
11、计算信号发射时刻的时间与本时段钟差参数参考时间（星历参数toc）的时间差、
12、
13、通过下式计算钟差,这个地方还没有处理相对论校正项与ｔｇd:
（代码中有个迭代过程,资料上没看到写）
14、
g）计算卫星位置、速度与时钟（ｓatｐos）
h）
依照星历表选项来选择不同的处理,如下:
Ø广播星历（EPHＯPT_BＲDＣ）:
ephpoｓ（广播星历到卫星位置与钟差）
Ø
1、依照公式计算出ｔk;
2、
3、依照使用的卫星系统,选择使用的地球引力常数（mu）与地球的角速度（omｇe）
4、
5、依照公式计算出平近点角Ｍ、
6、
7、求解开普勒方程,依照以下公式迭代求解。
8、
9、依照以下公式计算出ｕ（改正后的纬度幅角）,ｒ（改正后的径向）,i（改正后的轨道倾角）的值。
10、
11、依照以下公式计算卫星在轨道平面内的坐标、
12、
13、依照不同的卫星系统,做不同的计算。
14、
ＧPS计算方式:
北斗计算方式:
其中
15、依照公式计算出时间tｃ。
16、
17、依照以下公式计算出钟差与钟漂、
18、
Ø周密星历（EPHＯＰT_PＲEＣ）:
peｐh2pos
Ø
Ø广播＋SＢAS（EPHOＰT＿SＢAS）:
ｓａｔpｏｓ_sbas
Ø
Ø广播+ＳSR_APＣ（EＰHOＰＴ_SＳＲAPＣ）:
ｓatpｏs_ｓｓr
Ø
Ø广播+SSR＿（EPHOＰT_SＳＲ）:
saｔpos_ｓｓr
Ø
ØQＺＳSＬEＸ星历（ＥPＨOＰT_LEX）:
lexepｈ2pｏs
Ø
5.排除对被遮蔽卫星的测量（ｔｅsteclｉｐsｅ）
6.
a）获得太阳与月亮在ｅcｅf中的位置（ｓunmoonpos）
b）
c）归一化三维矢量（normv3）
d）
e）欧几里德范数的向量（noｒm）
f）
g）计算太阳—地球－卫星角（ｄｏｔ）
h）
7.相位与编码残差（ｒes＿ppｐ）
8.
a）地球潮汐改正（tｉdediｓｐ）
b）
c）计算几何距离（gｅodiｓt）
d）
1、用卫星的坐标向量做欧几里德范数,返回值与地球长半轴（ＷGS8４）比较、小于地球长半轴（WＧＳ８4）,返回—1;
2、
3、计算卫星坐标与接收机坐标的差值向量、
4、
5、用差值向量做欧几里德范数,再用上一步计算出来的差值向量与计算结果做除法,得到视线向量。
6、
7、使用以下公式计算几何距离。
8、
ﻩ备注:
这个公式是地球自转改正公式。
地球自转改正产生地球自转改正的机制在于地固系随地球的自转而旋转, 地固系是非惯性系统,用地固系中的坐标计算卫星到接收机几何距离时,由于卫星位置与接收机位置是两个不同时刻的位置矢量, 而且这两个时刻的地固系相关于惯性系是变化的,因此要考虑地球自转引起的距离变化。
其改正公式为:
其中, w 为地球自转角速度;C为光速;xsi与xri分不代表卫星位置矢量与测站位置矢量的分量,i=1,2对应ｘ,y分量。
地球自转对纬度影响甚小, 对经度影响最大,其次是高度、同时地球的自转影响还跟测站的经纬度以及测站与卫星的几何关系有关, 对两极的测站, 影响为零; 对赤道上的测站影响最大;当卫星在测站子午面内影响为零;卫星在测站东方时,影响为负,卫星在测站西方时,影响为正。
e）卫星方位角/仰角计算（satａｚeｌ）
f）
1、把接收机eｃef坐标转换到大地坐标;
2、
3、判断高度是否大于地球半长轴（ＷＧＳ84）的负数值;
4、
5、假如高度小于等于地球半长轴（WGS84）的负数值,方位角为０,仰角为ＰI／2;
6、
7、假如高度大于地球半长轴（ＷGS８4）的负数值,把ＥＣＥF向量转换到局部坐标。
然后对转换出来的坐标做内积、
8、
9、使用如下公式计算卫星方位角与仰角。
10、
g）对流层延迟修正:
h）
Ø通过标准大气压与ｓaastａｍｏｉnen模型计算对流层延迟 （ｔroｐmodel）
Ø
1、使用以下公式计算总气压。
2、
3、使用以下公式计算绝对温度。
4、
5、使用以下公式计算水蒸汽的分压。
代码中没有使用。
6、
7、使用以下公式计算【Ｓａasｔamoinｅn模型】、
8、
i）卫星天线模型（ｓatantpcv）
j）
k）接收机天线模型,通过天线相位中心参数计算天线偏移 （antｍｏdeｌ）
l）
m）相位缠绕校正（wｉｎdｕpcｏrr）
n）
o）电离层与天线相位校正测量（corrｍeaｓ）
p）
q）计算卫星时钟与对流层延迟
r）
9.卡尔曼滤波（ｆiｌtｅｒ）
10.
卡尔曼滤波依照以下的公式更新状态:
K ＝P＊H＊（Ｈ'*Ｐ ＊Ｈ+R）^-1;
xp =x +Ｋ＊v;
Ｐｐ＝（I– K＊ H'）* P ;
参数解释:
x:
状态向量 （nｘ1）
ｐ:
状态的协方差矩阵（ｎ ｘｎ）
H:
设计矩阵的转置矩阵（ｎ x m）
v:
创新（测量—模型）（mx １）
Ｒ:
测量误差的协方差矩阵（m xm）
n,m:
状态与测量值的数量
xp:
状态更新后的向量 （ｎｘ 1）
Ｐp:
状态更新后的协方差矩阵（nxn）
11.组合残差（res＿ppｐ）
12.
13.解析PPP整周模糊度（pppａmｂ）
14.
a）ａvｅｒagｅ ＬC（aｖeｒage_LＣ）
b）
c）固定宽巷模糊（fix＿ａmb_WL）
d）
Ø计算LC波长度（m）（lam＿LC）
Ø
Ø宽巷模糊
Ø
Ø宽巷模糊度的方差
Ø
e）fiｘ narrｏw－lane ａｍbｉｇｕｉty 【ＡRmoｄe:
ＰPP－ＡR】
f）
g）ｆixnａｒｒow－lａnｅ ambｉgｕity ｂyIＬS【ARmode:
PＰP—ＡRILS】
h）
5ﻬ以下是上面一些算法的参考资料,信息来自于Rtkｌib的手册
6
6.1接收机与卫星天线之间的几何距离
6.2
6.3卫星方向的方位角与仰角
6.4
6.5GPS、Galｉｌeｏ与QＮSS广播星历与时钟
6.6
6.7北斗广播星历与时钟
6.8
ﻩ
6.9对流层与电流层模型
6.10
6.10.1对流层模型
6.10.2
6.10.3电离层模型
6.10.4
6.11单点定位
6.12
6.12.1线性最小二乘估计
6.12.2
6.12.3估算接收机的位置与钟差
6.12.4
6.12.5估算接收机的速度与钟漂
6.12.6
6.12.7求解验证与RAIM FDE
6.12.8
6.13周密定位
6.14
6.14.1用于PPP的ＺD测量模式
6.14.2
6.14.3接收机天线相位中心模型
6.14.4
6.14.5卫星天线相位中心模型
6.14.6
6.14.7由地球潮汐引起的位移
6.14.8
6.14.9相位缠绕校正
6.14.10
6.14.11通过ＰＰP估算接收机的位置
6.14.12
升级会员 