《RTKPPP定位算法流程.doc》由会员分享，可在线阅读，更多相关《RTKPPP定位算法流程.doc（49页珍藏版）》请在三一文库上搜索。
1、1基础知识1.1 GPS精密单点定位得基本原理GPS精密单点定位一般采用单台双频GPS接收机，利用IGS提供得精密星历与卫星钟差， 基于载波相位观测值进行得高精度宦位。观测值中得电离层延迟误差通过双频信号组合消除, 对流层延迟误差通过引入未知参数进行估计。1.2时间系统RTK L IB内部使用GPST( G P S T时间)用于G N S S得数据处理与泄位算法。数据在RT K LIB内部处理之前，需要转换成GPST时间。使用GPST得原因就是避免处理润秒、RTKLIB 使用以下结构体表示时间：t y pedef s tru c tt ime_ t time;/ * tim e ( s ) e
2、xpr e ss e d by s ta n d ar d time_t */double sec;/* fr a ctio n of second under 1 s */ gtime_t;1.2.1 GPST 与 UTC(Un i versa I Tim e Coord i n at e d )关系参考【图1】，参考【图2】:图1转换关系公式Table E.l-1 GPST-UTC Values (until Marcli 2013)Time Since (in LTTC)rLS (s)1980-01-0600 : 000001981-07-0100000011982-07-0100000
3、021983-07-0100000031987-01-0100000041988-01-0100000051990-01-0100000061991-01-0100000071992-07-0100000081993-07-0100000091994-07-01000000101996-01-01000000111999-01-01000000131997-07-01000000122006-01-01000000142009-01-01000000152012-07-0100000016图3。通过使用GPS导航信息中得UTC参数,GPST到UTC或者UTC到GPST之前得转换可以用 更准确得
4、表达方式，如【图3】、fUTC =fGPS LS + -% + 4血 Iot + 604800（W7vt-仍：）图3月0, A, fE, fot,are开K这些参数就是由GPS导航消息提供得。1.2.2 BDT（北斗导航卫星系统时间）BDT（北斗导航卫星系统时间）就是一个连续得时间系统，没有润秒。开始历元得时间就 是【UTC 2006年1月1号00:00:00】。北斗时间计算公式【图4】:图4UTC与GPST时间转换同上而得GPS-样,只不过UTC参数来自与北斗导航信息中。1.3坐标系统接收机与卫星得位苣在RTKLIB中表示为在ECEF（地心地固坐标系）坐标系统中得X. Y. Z组件、1.3.
5、1大地坐标到ECEF坐标得转换转换公式如【图5】。= 7（2- /）_ aJl _ g亠修广f、（卩 +力）COS给 COS /i”rr （v + 力）cos 0 sin 第三个公式最后一行有错，应该为：（v（l 一 c2）+h）s i n参数说明:“：地球参考椭球得长半径图6/:地球参考椭球得扁平率 h:椭球高度.纬度:经度当前版本得RTKLIB使用得值为【图6】:图7参考椭球体1.3.2 ECEF坐系到大地坐标得转换转换公式如【图8】加0 = 0= arctan/-ae1 tail z舛=lim 0VA = ATAN2Qv,x)f rah =1=cos J(l_/)snF0.图 8133本
6、地坐标到ECEF坐标得转换在接收机位置得本地坐标，也被称为ENU坐标，通常使用在GNS S导航处理。ECEF坐标 到本地坐标转换得旋转矩阵表示为【图9】。-sinA,.cos 切0、Er = -sinr cos/ir -sin sin2r cosg 、cos0cos/ cossin sin丿图9 7旋转矩阵参数说明:接收机位置得纬度:接收机位置得经度通过使用与接收机得坐标ECEF，坐标咒理f【ECEF】可以被转换到本地坐标 得坐标rlocal.公式如【图10】、图102 RTKPPP定位算法2.1单点定位(p ntpos)1: s atposs -2: e s tp o s - 3: estv
7、el1.计算计算卫星位置、速度与时钟(sa t p oss)a)通过广播星历计算卫星钟差(ephclk) 卫星编号到卫星系统得转换(sa t sys)根据卫星得编号，获取到对应得卫星导航系统。 选择星历(seleph)1、传入信号传输时间，卫星编号，导航数据等参数、2、遍历导航数据，遍历导航数据里而得星历数据，判断星历数据得卫星编号 就是否与传入得卫星编号相等。3、如果星历数据得卫星编号与传入得卫星编号相等，就计算星历参考时间 (t oc)与信号传输时间得时间差。如果不相等,继续处理下一条星历数据、4、判断计算出来得时间差，如果时间差大于了允许得最大时间差，继续査找下一个星历数据、否则，判断时
8、间差最小得星历数摇，记录星历数拯得位宜、5、返回之前记录出来得星历数据乙 使用广播星历计算卫星时钟偏差(e ph2c 1 k)1、传入信号发射时刻得时间与星历数拯。2、讣算信号发射时刻得时间与本时段钟差参数参考时间(星历参数sc) 得时间差。3、通过下式讣算钟差，这里还没有处理相对论校正项与tgd:(代码中有个 迭代过程，资料上没瞧到写)。卫星賈差汁算出來之信号发射时时间还要加上这个钟.b)计算卫星在信号发射时刻得位宜、速度与时钟(satp o s) 根据星历表选项来选择不同得处理，如下： 广播星历(E PHO PT_BRDC):ephpos (广播星历到卫星位置与钟差)1、根据公式计算出tk
9、;2、根据使用得卫星系统，选择使用得地球引力常数(mu)j地球得角速度(omge)3、根据公式计算出平近点角M、4、求解开普勒方程，按照以下公式迭代求解。5、根据以下公式汁算出11(改正后得纬度幅角),I(改正后得径向),i(改正后 得轨逍倾角)得值、J1-以 sinEa=cosE - o= arctail +dll = CU5 sin 20 + Cuc cos 1(/)dr = Crs sill 20 + Crc cos 2(pSi = C15 sin 2+ Clc cos 2u =(p + 6u7 = d(l-ccosE) +丹i =z0 + Jz6、根据以下公式计算卫星在轨道平面内得坐标
10、、7、根据不冋得卫星系统，做不同得讣算、GPS计算方式:C = Q +(_叫)0 Cte Lcosu co$/2 -siinz cos? sin I?尸5(r) = r cos/ sin/? + sin/ cosf cosZ2 sin u sini北斗计算方式:Q =人一叫 tQecos cos?一sill cos? sin2厂心（-5。）com sin/? +sin;/ cosj cosQ.sinw sin?/00cos fA sin E - bTGDdTs(t) = afi2af2tc 精密星历(EPH 0 PT_P R EC): p eph 2 po s 广播+S B AS ( E P
11、HOPT_S B A S):s a tpos_ s bas 广播+SSR_ A PC( E P 110 PT_SS RAPC): satpos_ s s r 广播 + S SR_COM(EPH O P T_SSR):sa t pos_ s sr Q ZSS LEXM5(EPHO PT_LEX): 1 ex e ph 2 p os2.使用伪距估算接收机得位置，返回估翕状态结果(est pos)a)伪距残差(re s c o de) 把ecef坐标系转换成大地坐标系(e c ef2p o s )1、按照以下公式做转换，暂时还没瞧懂。(pT 弁1 = arctanae2 tail 如彳1+(1 -外
12、丿,0 = lini 如A = ATAN.x)ah=cos 族 J(l_F)sii2g计算几何距离与接收机到卫星得单位矢量(geodi s t)1、用卫星得坐标向量做欧几里徳范数，返回值与地球长半轴(WGS 8 4 )比较、小于地球长半轴(WGS84),返回一 1;2、汁算卫星坐标与接收机坐标得差值向量。3、用差值向量做欧几里徳范数，再用上一步计算出来得差值向量与计算结果做除法，得到视线向量。按照如下公式:计)-也)rs(ts)-rr(t.)4、使用以下公式计算几何距离。Pr H b ( ) - r W )| + -yC-Vr -ySXr)计算卫星方位角/仰角(satazel)1、把接收机cc
13、ef坐标转换到大地坐标；2、判断髙度就是否大于地球半长轴(WGS8 4)得负数值；3、如果髙度小于等于地球半长轴(WGS84)得负数值，方位角为0.仰角为P 1/2:4、如果高度大于地球半长轴(WGS84)得负数值，把ECEF向量转换到局 部坐标。然后对转换出来得坐标做内枳。5、使用如下公式计算卫星方位角与仰角、er,emi =四 =(%，勺勺)Az j. = ATAN2(ee,eu)Elj. = arcsin(&“)伪距使用编码残差改( P ran g e)暂时没找到对应得文档对应。电离层改正(iono corr) 通过广播电离层模型(klo b uchar模型)计算岀电离层延迟(i o n
14、model)1、校验传入得电离层模型参数，校验失败.使用默认得电离层模型参数；2、使用以下公式计算岀地球为中心得角度(半圆)；3、使用以下公式汁算子得电离层得纬度/经度(半圆)；4、使用以下公式计算地磁纬度。5、计算本地时间，返回值按这个公式tt-=floor(t t / 86400o 0)*86400、0;处理，保证tt得范用(0= t t1.57)悴 1.57)对流层改正(tropcorr) 通过标准大气压与saas t amoi n en模型计算对流层延迟(t r o pm o del) 1、使用以下公式计算总气压。2、使用以下公式汁算绝对温度、3、使用以下公式汁算水蒸汽得分压、代码中没
15、有使用、4、使用以下公式计算Saa s ta mo i nen模型】丁S 0.002277 f 1255 门“、 、耳=cosz l T丿 伪距残差残差值=伪距一(物理距离+dtr-光速冬时钟偏差+电离层误差+对流层误差)时钟系统与接收器得偏置补偿 伪距测呈:误差方差(v arerr)b) 方差权重值(weigh t b y va r i a nee)c) 最小二乘估计(1 sq) 最小二乘估计通过求解正规方程(X =(A * A* )a-1*A* Y)1、计算矩阵A*Y得结果.保存到矩阵Ay;2、讣算矩阵A* A得结果，保存到矩阵Q;3、求得矩阵Q得逆矩阵.结果保存到矩阵Q：4、最后x=Q*
16、 Ayd) 欧几里徳范数(norm)公式：1、通过最小二乘法估算出来得参数dx【长度为4】，分别加到位置向量上x【长 度为4】。2、对这个速度向量做欧几里徳范数，返回值同1E-4比较，小于这个值时,就得到 估算出得接收机得位置为向量X得值。e) 验证求解(valsol)1、对伪距残差值做内积，然后同卡方分布得自由度分布值(alph a=0 . 00 1) 做比较，大于卡方分布得自由度分布值得数据,无效。2、计算d ops ；暂时还没没找到资料。3、用计算出来得dop s与配置得最大d o ps阈值比较，大于配置得最大d o p s 阈值得数据无效。3.使用多普勒估算接收机速度(es t vel
17、)a) 多普勒残差(r e sd o p) 把ecef坐标系转换成大地坐标系(ecef2p o s ) 把xyz坐标转换成enu坐标(xyz2enu) 计算ecef中得瞄准线向量 计算相对于接收机在ECE F中得卫星速度 多普勒残差计算1、速率公式：b) 最小二乘估计(I S q) 最小二乘估计通过求解正规方程(X =(A * A) -1* A *Y)1、计算矩阵A*Y得结果，保存到矩阵A y ;2、计算矩阵A * “得结果，保存到矩阵Q;3、求得矩阵Q得逆矩阵，结果保存到矩阵Q;4、最后 x=Q * Ayc) 欧几里徳范数(norm)公式：3、通过最小二乘法估算岀来得参数dx【长度为4】，分
18、别加到速度向量上x【长度为4】。4、对这个速度向屋做欧几里徳范数,返回值同1E 6比较,小于这个值时，就得到估算岀得接收机速度为向量x得值。2.2精确定位(pppos)1:uds t a t e _ppp - 2: sat p oss - 3 : t e s tec 1 i p se - 4: res_ppp - 5 : r e s_p p p 6: filt e r 7: r e s_ p pp1. 暂时更新状态(u dst a te_p p p)a) 位置更新b) 时钟更新c) 对流层参数更新d) 相位偏差更新1、通过LLI检测周跳；2、通过Geomet ryF r X相位跳变检测周跳(如
19、果双频测量值可用)；2. 计算卫星得位程与速度与时钟(satp o ss)c)通过广播星历计算卫星钟差(ephclk) 卫星编号到卫星系统得转换(sat s ys)根摇卫星得编号，获取到对应得卫星导航系统、 选择星历（s eleph）1、传入信号传输时间，卫星编号，导航数据等参数。2、遍历导航数据.遍历导航数据里而得星历数据，判断星历数据得卫星编号 就是否与传入得卫星编号相等。3、如果星历数据得卫星编号与传入得卫星编号相等，就计算星历参考时间 （toe）与信号传输时间得时间差。如果不相等,继续处理下一条星历数据。4、判断计算出来得时间差，如果时间差大于了允许得最大时间差,继续查 找下一个星历数
20、据。否则,判断时间差最小得星历数据，记录星历数拯得位 置。5、返回之前记录出来得星历数据。 使用广播星历计算卫星时钟偏差（eph 2 c Ik）1、传入信号传输时间与星历数据。2、传入信号发射时刻得时间与星历数据。3、il算信号发射时刻得时间与本时段钟差参数参考时间（星历参数sc） 得时间差、4、通过下式汁算钟差，这里还没有处理相对论校正项与tgd:（代码中有个 迭代过程，资料上没瞧到写）d）计算卫星位置、速度与时钟（satp os）根据星历表选项来选择不同得处理，如下： 广播星历（E PHOPT_BR 0 C）: e p h p os （广播星历到卫星位置与钟差）1、根据公式计算出tk；2、
21、根据使用得卫星系统，选择使用得地球引力常数（mu）与地球得角速度 （omge）3、根据公式计算出平近点角4、求解开普勒方程，按照以下公式迭代求解。5、根据以下公式讣算出u（改正后得纬度幅角）,1（改正后得径向）,i（改正后 得轨逍倾角）得值、cosE-o0 = arctana + e 创=CU5 sin 20 + Cuc cos l(f)dr = Cs sill 20 + Crc cos 205i = Cs 血 2+ Clc cos 2 u7 = d(l-ccos) + 77 =Z0 十刃 +6、根据以下公式计算卫星在轨道平而内得坐标匚7、根据不同得卫星系统，做不同得讣算。GPSil-算方式：
22、x? = /?q + (/2-d)e) tk -oe roecos” cos (r) = r cosi/ sinjQ +sin/ cosz cos 精密星历(EPHOPT_P R E C): peph 2pos 广播+ SB A S(EPHOPT_SBAS):sa t p o s_ s bas 广播+ SS R_A P C( E P H O P T_ S SRA P C): satpo s _ s s r 广播+SS R_(EPH OP T_SS R C OM):satpos_ s s r QZSS LEX 星历(EPHOPT_LEX):lexeph 2 pos3. 排除对被遮蔽卫星得测量(t
23、esteci ipse)a) 获得太阳与月亮在ecef中得位K(s u nm o onpos)b) 归一化三维矢量(n ormv3)c) 欧几里徳范数得向量(norm)d) 计算太阳一地球-卫星角(dot)4. 相位与编码残差(r es_ p pp)a) 地球潮汐改正(tide d i s p)b) 计算几何距离(ge o dist)1、用卫星得坐标向量做欧几里徳范数，返回值与地球长半轴(WGS 84)比 较。小于地球长半轴(WG S 8 4 ),返回- 1 ;2、计算卫星坐标与接收机坐标得差值向:S：。3、用差值向量做欧几里徳范数，再用上一步汁算出来得差值向量与计算 结果做除法，得到视线向量
24、、4、使用以下公式计算几何距离。Pr -(6 )-心产)| +- 5)备注:这个公式就是地球自转改正公式。地球自转改正产生地球自转改正得机制在于地固系随地球得自转而旋转.地 固系就是非惯性系统，用地固系中得坐标计算卫星到接收机几何距离时.由于 卫星位置与接收机位置就是两个不同时刻得位宜矢疑，而且这两个时刻得地 固系相对于惯性系就是变化得，因此要考虑地球自转引起得距离变化、其改 正公式为：心，二彳-“兀 $1 一 兀 r) XS2- ( Xsl- Xr2)XSJ其中，w为地球自转角速度;C为光速:Xsi与Xri分别代表卫星位置矢绘与测 站位巻矢屋得分i= 1,2对应X,y分量。地球自转对纬度影响
25、甚小，对经度影 响最大，英次就是髙度、同时地球得自转影响还跟测站得经纬度以及测站与卫星 得几何关系有关，对两极得测站，影响为零；对赤道上得测站影响最大；当卫星在 测站子午而内影响为零;卫星在测站东方时，影响为负，卫星在测站西方时,影响为正。c) 卫星方位角/仰角计算(satazel)1、把接收机ccef坐标转换到大地坐标；2、判断高度就是否大于地球半长轴(WGS84)得负数值；3、如果髙度小于等于地球半长轴(WGS 8 4)f#负数值，方位角为0，仰角 为 PI/2:4、如果高度大于地球半长轴(W G S84)得负数值，把ECE F向量转换到局 部坐标。然后对转换出来得坐标做内积。5、使用如下
26、公式计算卫星方位角与仰角。r.emi =四 (ee en eu)AZj. = ATAN2ee,eJJ)Ell = arcsiii(Q“)d) 对流层延迟修正： 通过标准大气压与saastamo i n en模型计算对流层延迟(t r o pm odel)1、使用以下公式计算总气压、2、使用以下公式计算绝对温度、3、使用以下公式计算水蒸汽得分压、代码中没有使用。z 1255 J、P+ 0.05e 一 tair z1T,4、使用以下公式计算(Saas tamo i ncn模型】。0.002277COS 2e) 卫星天线模型(satantpcv)f) 接收机天线模型，通过天线相位中心参数计算天线偏移
27、(a n tmo d el)g) 相位缠绕校正(windupc o r r )h) 电离层与天线相位校正测量(corrmea s )i) 计算卫星时钟与对流层延迟5. 卡尔曼滤波(f i Iter)卡尔曼滤波按照以下得公式更新状态：K = P * H *(H;* P * H+ R)1;xp =x + K * v;Pp= (I- K *Hr )*P;参数解释：X：状态向虽(n X 1)p：状态得协方差矩阵(n x n)II：设讣矩阵得转1矩阵(n x m)v :创新(测量模型)(m x 1)R：测虽:渓差得协方差矩阵(m x m)n.m:状态与测量值得数量xp:状态更新后得向星：(nx 1 )P
28、p：状态更新后得协方差矩阵(n x n)6. 组合残差(r e s_pp p )7. 解析PPP整周模糊度(ppp a mb)a) aver a g e LC (average_LC)b) 固左宽巷模糊(f ix_a mb_WL) 计算LC波长度(m)(lam_LC) 宽巷模糊 宽巷模糊度得方差c) f i x narrow-lane ambigui t y AR mode: PPPARd) f i x narrow Ian e ambig u i t y by ILS 【AR mode: PPP AR ILS 以下就是上面一些算法得参考资料，信息来自于Rtklib得手册3.1接收机与卫星天线
29、之间得几何距离(4) Geometric range between receiver and satellite antennasTh。geometric range is defined as th。physical distance botveen the satellite anteima phasp center position and the receiver antenna phase center position in the inertial coordinates. At first, the signal transmission time P can be deri
30、ved from：(E.3.7)The both sides in the equation includes . So several iterations are needed to solve the equation. The geometric range can be expressed by using the receiver and satellite antenna phase center positions 5(6)=(v； )F at the time lr and r(/s) = (.v$.y5.zi)T at the time ts in the ECEF (ea
31、rth center earth fixed) coordinates as：(E-3.8)P；二帆曲(“)一附)where ?r(z) is the ECEF to ECI (earth center inertial) coordinates transformation matrix at the time / For the expression in the ECEF coordinate the earth rotation effect has to be incorporated in to obtain the geometric range. The equation ca
32、n be approximated by one of the following equations with adequate precision under 1 nun level. Current versioiA RTKLIB always uses the equation (F.3.8B) for the geometric range The last term in (F.3.8b) is sometimes called as Sagiac effect(E.3.8a)Pr 一 Rz(%# /c” (r )|pr &|rr(rr)-(训 + 手(丹r-ysxr)(E.3.$
33、b)Pr 彳-冬(E.3.SC)+ +ECEF at signal transinission ECEF at signal receptionFiguie E.3-1 Geoxnetiic Range and. Earth Rotation Correction3.2卫星方向得方位角与仰角Th。unit LOS (lin-of-sight) vector from the receiver to th。satellite can bo expressed in th。ECEFcoordinates as：(E.3.9)In tlie equation, the earth rotation
34、effect is neglected. The azimuth and elevation angles and 7； of the satellite direction from the receiver site can be derived from：(E.3.10)(E.3.11)(E.3.L2)= E 启=(s 竝)丁AZy = ATAN2(ee,etl)Elf. = arcsint)where Er is the coordinates rotation matrix from ECEF to the local coordinates at the receiver posi
35、tionRefer E.2 for detailed fonnation of the matrix.Figure E.3-2 Local Coordinates and Azimuth and Elevation An吕les3.3 GPS、Gali I eo与QNSS广播星历与时钟(1) Broadcast ephemerides and clocks for GPS, Galileo and QZS5 问Broadcast ephemeris and SV dock parameters for GPS, Galileo and QZSS are given in navigationm
36、essages as：By using tliese parameteTS, the satellite position anteixna phase center position) rs (/) in ECEFf the satellite clock bias d厂(/) and dock drift dT5 (?) are computed as:(E42)M=% + | + d fk(E.4.3)M = E - e sin E(E44)siiia =(E45)cos-e = arctaiia +(y(E.4.6)6“ = CtiS sin 20 + Citc cos 2$(E.4.
37、7)6r- Crs sin 20 + Cfc cos 20(E.4.8)5i = Cis sin 辿 + Cic cos 20(E.4.9)u - 0十加(E.4.10)r = a(-ecosE) + 6r(F.4.11)/ = /0 4- $/ 4- hk(E.4.12)Q =山-畋尬(E.4.13)cos” cos J? - sill cos 7 sin I?r5(r) = r cos” sin2 +sill co cos = /i2 / fi2 for 厶 pseudorangeTgd : group delay parameters for GPS and QZSS, BOD for
38、Galileo (s)Tlie Kepler equation (E.4.4) can be solved the following iteration by Newtons methodQ - esinQ - Ml-ecos;E = lim E(E418)(E.4.19)(E4.20)Tlie broadcast ephemerides and clock are applied in case that the processing option Satellite Ephemeris/aockM to Broadcast* as well as GLONASS, BeiDou and
39、SB AS.3.4北斗广播星历与时钟(3) Broadcast ephemerides and docks for B?iDou 门For BeiDou satellites, tlie siuiilar ephemeris and clock parameters as GPS, Galileo and QZSS areprovided in the navigation messages as:.CWJ.CwcrCrj.Cr.Cj5.Crc，r7/.n/.q/.7p)(E.4.28)For MEO and IGSO satellites of BeiDou, the same formul
40、ations as (1) for GPS ephemeris and clock, except for =3.98600441Sxl014,= 7.2921150x 10-5 rad/s and the tinw is expressed in BDT.1b obtain the satellite position rs (/) of BeiDou GEO satellites at the time / in BDT, the equation(E.5.13) and (E.5.14) should be replaced by:2 = 2+ Cg -咲 Qcos?/ cosP-sun
41、/ co$f sm/2 rs(t) =cos?/ sin(2 + sin” cos/ cos/2sin/ sin/(F.4.29)(E.4.30)where:u00COS&sm00cos 6?sin&,R：二一 sincos 800X一 sin0COS&Z0X0-3.5对流层与电流层模型3.5.1对流层模型The standard atmosphere can be expressed as：(E.5.1)(E5.2)(E53)p = 1013.25x(1 - 2.2557 x 10-5/)52568T = 15.0-6.5x10-3/j + 273.15r ino J7.157-4684.0： hrel e = 6.108x expvr-38.45；1000.002277cos 二+ 0.05 ；(J-tanwhere P is the total pressure (liPa), T is the absolute temperature (K) of the air, h is t