RTKPPP定位算法流程

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

2、y standard time_t */ double sec; /* fraction of second under 1 s */ gtime_t; 1.2.1 GPST 和 UTC(Universal Time Coordimated) 关系参考【图1】，/仏参考【图21： 图1转换关系公式 精选文档 lAble E，i-1 CPST-LITC(luitil March 2013) Tinii? StitCE? (in TJTC)(s) 1980-01-06 00 00 00 0 1381-07-01 00 00 00 1 1982=07-01 00 00 00 2 1993-07-01

3、 00 00 00 3 1号37-01-01 00 00 00 4 1938-01-01 00 00 00 5 1930-01-01 00 00 00 6 1991-01-01 00 00 00 7 1992-07-01 00 00 00 a 1993-07-01 00 00 00 5 1994-07-01 00 00 00 10 1996-01-01 00 00 00 11 1999-01-01 00 00 00 13 1997-07-01 00 00 00 12 2006-01-01 00 00 00 14 2009-01-01 00 00 00 15 2012-07-01 00 图 00

4、 2 00 IE 通过使用GPS导航信息中的UTC参数，GPST到UTC或者UTC到GPST之前的转换可以用 更准确的表达方式，如【图 3】。 LTC =GPS LS + 心 + 月iUe+ 604800(7/ -fl)j 心用丹are卩。这些参数是由GPS导航消息提供的。 1.2.2 BDT （北斗导航卫星系统时间） BDT （北斗导航卫星系统时间）是一个连续的时间系统，没有润秒。开始历元的时间是 【UTC 2006 年 1 月 1 号 00:00:00 】。 北斗时间计算公式【图4】： UTC和GPST时间转换同上面的GPS样，只不过UTC参数来自与北斗导航信息中。 1.3坐标系统 接收机

5、和卫星的位置在 RTKLIB中表示为在ECEF（地心地固坐标系）坐标系统中的X, Y, Z 组件。 1.3.1 大地坐标到ECEF坐标的转换 转换公式如【图5】。 *=/P- /） /1 -矿sill松广 （节十 A）cos cos + 你 sin 人 ,T（1 一于）510 一 J 第三个公式最后一行有错，应该为：(v(1 -e2)+h)sin此 图5 参数说明：a:地球参考椭球的长半径 f:地球参考椭球的扁平率 h：椭球高度 曲I：纬度 叫经度 当前版本的RTKLIB使用的值为【图6】： (7 = 637$ 13工0 (m) / = 1.0/298.257223563 图7参考椭球体 1.

6、3.2 ECEF坐系到大地坐标的转换 转换公式如【图8】 r = yjx- + j- 如.0 二 丹J+1 二 aKUin 二 札=lim加P 2 = ATANKyx h =. 133 本地坐标到ECEF坐标的转换 在接收机位置的本地坐标，也被称为ENU坐标，通常使用在 GNSS导航处理。ECEF坐 标到本地坐标转换的旋转矩阵表示为【图 9】。 9 0、 cos处 sin 如 J COS 2 -sin 克 gin 如 cos 如 sin 图9 Er旋转矩阵 参数说明：必：接收机位置的纬度 心：接收机位置的经度 通过使用Er和接收机的坐标rr【ECEF】，坐标ecef【ECEF】可以被转换到本地

7、坐标的 坐标riocal，公式如【图10】。 仏】二EK住CEF-fJ 图10 2 RTKPPP定位算法 2.1 单点定位（pntpos ） 1:sat poss 2: est pos 3: estvel 1. 计算计算卫星位置、速度和时钟（sat poss） a）通过广播星历计算卫星钟差（ep hclk） ?卫星编号到卫星系统的转换（satsys） 根据卫星的编号，获取到对应的卫星导航系统。 选择星历（seleph） 1、传入信号传输时间，卫星编号，导航数据等参数。 2、遍历导航数据，遍历导航数据里面的星历数据，判断星历数据的卫星 编号是否和传入的卫星编号相等。 3、如果星历数据的卫星编号和

8、传入的卫星编号相等，就计算星历参考时 间（toe）和信号传输时间的时间差。 如果不相等，继续处理下一条星历数据。 4、判断计算出来的时间差，如果时间差大于了允许的最大时间差，继续 记录星历数据 查找下一个星历数据。否则，判断时间差最小的星历数据, 的位置。 5、返回之前记录出来的星历数据。 使用广播星历计算卫星时钟偏差（ep h2clk） 1、 传入信号发射时刻的时间和星历数据。 2、 计算信号发射时刻的时间和本时段钟差参数参考时间 （星历参数toc） 的时间差。 3、 通过下式计算钟差，这里还没有处理相对论校正项和 tgd :（代码中 有个迭代过程，资料上没看到写） 罕 Atsv 自0 +

9、器 toc) +一Toc) 卫星钟差计算出来之后，信号发射时刻的时间还要加上这个钟差。 b）计算卫星在信号发射时刻的位置、速度和时钟（sat pos） 根据星历表选项来选择不同的处理，如下： ?广播星历（EPHOPT_BRDC : ephpos （广播星历到卫星位置和钟差） 1、根据公式 计算出tk； 2、根据使用的卫星系统，选择使用的地球引力常数（mu）和地球的角速度 (omge) 3、 根据公式 +曲g 计算出平近点角M。 4、 5、 6、 根据以下公式计算卫星在轨道平面内的坐标。 求解开普勒方程 M = E - 小 E，按照以下公式迭代求解。 j siii Kj A/ I- 2、 判断高

10、度是否大于地球半长轴（WGS84 ）的负数值; 3、 如果高度小于等于地球半长轴（ WGS84 ）的负数值，方位角为 0,仰 角为PI/2; 4、 如果高度大于地球半长轴（ WGS84 ）的负数值，把ECEF向量转换到 局部坐标。然后对转换出来的坐标做内积。 5、 使用如下公式计算卫星方位角和仰角。 4二：-ATAN2e.ej) El：二 arcsiii(cJ 伪距使用编码残差改正（pran ge） 暂时没找到对应的文档对应。 电离层改正（ionocorr） 1、 通过广播电离层模型（klobuchar模型）计算出电离层延迟（ionmodel） 校验传入的电离层模型参数， 校验失败，使用默认的

11、电离层模型参数; 2、 使用以下公式计算出地球为中心的角度（半圆） = 0.0137/(/ + OJl)- 0;022 3、 使用以下公式计算子的电离层的纬度/经度（半圆）； 丹=P 十Az 4 二A + 护will卫二丿cos殍 4、 使用以下公式计算地磁纬度。 殉= . -1.617) 5、 计算本地时间，返回值按这个公式【tt-=floor(tt/86400.0)*86400.0;】处 理，保证tt的范围(0=tt 1.57) r = C sin20 + C枕逊20 tf = G sill 20 + C祀 cos 20 r = n(l-ecos) + tir 6、根据以下公式计算卫星在轨

12、道平面内的坐标。 耳kCOSllk 7、根据不同的卫星系统，做不同的计算。 GPS计算方式: Q = Qq +（d -叫山-匸気 (? COSH cos/? sin/ cos/ sin/? 北斗计算方式: 丿（F吳血唧鸟（-聾） 其中 8、 9、 COWH 弓in/? Yin” cOA? cos/2 sinw 5111/ cosi/ cos/? -sinJ/ cos? sin/? j COS2/ sin/? siiij/ cos / cos/? I sillsin; 1 0 0 siiit? O 比（叭= 0 cos日 山1旧 -siiit? costf 0 ：0 Sill + 笳r -二-何

13、SQ 精密星历(EPHOPT_PREC : peph2pos 广播 +SBAS ( EPHOPT_SBAS: satpos_sbas 广播 +SSR_APC( EPHOPT_SSRAP)C: satpos_ssr 广播 +SSR_COM(EPHOPT_SSRCOM satpos_ssr QZSS LE)星历(EPHOPT_LEX： lexeph2pos 3. 排除对被遮蔽卫星的测量（testecli pse） a） b） c） d） 获得太阳和月亮在 ecef中的位置（sunmoonpos） 归一化三维矢量（normv3） 欧几里德范数的向量（norm） 计算太阳-地球-卫星角（dot） 4.

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

15、位置是两个不同时刻的位置矢量，而且这两个时刻的地固系 相对于惯性系是变化的，因此要考虑地球自转引起的距离变化。其改正公式为： 二十 f 2 - XfI 丿心 2 一 ( Xs2- ”Hj.Bl 丿 其中,w为地球自转角速度；C为光速；X si和X ri分别代表卫星位置矢量和测 站位置矢量的分量，i= 1,2对应x , y分量。地球自转对纬度影响甚小，对经度影响 最大，其次是高度。同时地球的自转影响还跟测站的经纬度以及测站与卫星的几何 关系有关，对两极的测站，影响为零；对赤道上的测站影响最大；当卫星在测站子 午面内影响为零；卫星在测站东方时，影响为负，卫星在测站西方时，影响为正。 C）卫星方位角

16、/仰角计算（satazel） 1、 把接收机ecef坐标转换到大地坐标; 2、 判断高度是否大于地球半长轴（WGS84 ）的负数值; 3、 如果高度小于等于地球半长轴（ WGS84 ）的负数值，方位角为 0,仰 角为PI/2; 4、 如果高度大于地球半长轴（ WGS84 ）的负数值，把ECEF向量转换到 局部坐标。然后对转换出来的坐标做内积。 5、 使用如下公式计算卫星方位角和仰角。 AZj.二 ATAN2(e,ej) d） 对流层延迟修正： ?通过标准大气压和 saastamoinen模型计算对流层延迟（tropmodel） 1、 使用以下公式计算总气压。 P = 1013.25 X (1-

17、2.2557 y 1 护可 2、 使用以下公式计算绝对温度。 3、使用以下公式计算水蒸汽的分压。代码中 100 没有使用。 0 = 6 J OS X exp* 4、 使用以下公式计算【Saastamoi nen模型】。 cosr 5. 6. 7. e) f) g) h) i) 卫星天线模型( 接收机天线模型， 相位缠绕校正( 电离层和天线相位校正测量( 计算卫星时钟和对流层延迟 sata ntpcv) 通过天线相位中心参数计算天线偏移 windupcorr) corrmeas) 卡尔曼滤波(filter) 卡尔曼滤波按照以下的公式更新状态: xp (antmodel ) P * H * (H

18、* P * H + R)A-1 Pp = (I -K * H) * P 参数解释: X：状态向量 (n x 1) p:状态的协方差矩阵 (n x n) H:设计矩阵的转置矩阵(n X m) V：创新(测量一模型)(m X 1) R:测量误差的协方差矩阵(m x m) n,m:状态和测量值的数量 xp :状态更新后的向量(n x 1) Pp:状态更新后的协方差矩阵 (n x n) 组合残差 (res_ppp) 解析PPP整周模糊度(PPPamb) a) b) c) d) average LC( average_LC) 固定宽巷模糊(fix_amb_WL) ? 计算LC波长度(m)( lam_LC

19、) ?宽巷模糊 ? 宽巷模糊度的方差 fix narrow-la ne ambiguity【AR mode: PPP-AR fix narrow-la ne ambiguity by ILS 【AR mode: PPP-AR ILS 以下是上面一些算法的参考资料，信息来自于 Rtkiib 的 手册 3.1 接收机和卫星天线之间的几何距离 Geometric range betivRen receiver and sateilite anieiuias The geometric tfinge is defined as the physicfd diand ?： of lhe satellit

20、e direction front I he receiver site tan be derived from： (E.3.L0) tE.3.11) (E.3.L2) 匕；,咖=E严;=代灼-弓F -上;=卫TLlN工（P甘片） E = arcsiu（fjJ wliBre p is 一he coordinates 1,G罗iion InarliM from ECEF fo 一he lac巴 Lroordil-ElrFS 算 fhe r巾ceNer posi-iQn Refer E 2 fOT d阵比云d fornwfinll of Hw s 畀rix. PgLire E,32 Local c

21、grdi二Htem and Azimu+h EleTaHOIl Allglew 精选文档 29 3.3 GPS、Galileo和QNSS广播星历和时钟 Broadcast eplienieiides and flocks for GPS Galileo find QZSS 习同 Broadcast pphfineris and SV clock paranieteTS for GPSj. Galilio and QZSS are 耳iven in Tiavi耳xtico messages as： %忒3、汕2巧=(2心岛网丄心如匚也5工时5沁話工耐0耐卫上5卩 (E.4.1) By using

22、 these paranieters, die satellite pMition (antennd pbase center posiHon) h(f iii ECEE the satellite clock bias dF/) arid dcxk drift dtr are computed as： (E-4-3) (E.4.3) Af = E - 劭iu E Jl -/ siuf a = cos -e. (E.4.4) (E.4.5) tSfj = Cum sill 20 + 佃 cos 辿 (E.4.6) (E47) (E.4.e) (E.屯刃 r = f7(l -fcos) + I

23、f = /g + 血 + J =+(/?-性)4 -饨 Q uc cosC-shiH cosi sin J? cosM sill /? + sill ti COST cos/? 11LM silt/ (E.4.10) (F.4.11) (E.4.12) (E.414 (E.4J5 (E.4.17) (F.4.I6) H丁=吠 +血*亠-Sill E -tJcD where； 馱-.parth gMTit日Kona】constaiit ( 3.9S60050 x10* ni/ffoT GPS and QZS勺 5.9SfiOCM41 沁】。口 ni for Galileo) : earth an

24、guUr vEofity ( 7.2921151467 k 10 rad/sj 5= fC 音 for A pseuilorai吧e 花Q group delay parameters for GFS aiid QZSS,比d for Galileo The Keplpr equation (F.4.4 can bp plveni thp folloivin itprahon by Xjewlons method. E严订 T R Ef -ehinEi-J /+i (E.4.18) -pens 為 (E.4.19) F = liiii Ef 豆二斑+f”-氓,Q cos” cos?-siiif

25、t cos; sill？ (zcisM siiiJ? + SIU ff COST cos I? sin/ sill), (E4.29) rp) (E.4.30) (E.4.2G) Ihe broadcast epliemerides and clock are Applied in case that the processing optioii Satellite FphpiiLfTis/Clock to Broadcast re wpII as GLON ASS, BeiDon nd SBAS. 3.4北斗广播星历和时钟 (3) Broadcast pplieTiieridps and c

26、locks for Be-iDcii PI For BciDoii atellites the similar t?pbcnieTi5 卍id clock parameters as GPS, Galileo and QZSS are proded ill the iidvigatioii iiiessdges as; fF.4.28) Peph (阳 aJ 二仗占-Q 坯-似 Q, JJI, 7./?. f. Gr 匚 I 心T，5-jf 蟻卫GD F For MEO and IGSO satellites of BeiDott, the same forniulaticns as (1)

27、for GFS ephenwris and clock, except for P = .986004118x10= 7.2921150 x IC rad/s and the time f is expressed in UDI lb obtAuL the satellite position 严 of KeiDou GEO satellites at tJw time 卞 in BDX tlw equation (E.5.13) and (E.5.H) shotild be replaced by： where; 1 0 0、 cos(9 sill 日 0 %(= 0 cosH !iin00

28、2277 F eosz where - is th? zenith angle (rad) as ： = x I - E!；- The standard atmosphere and the Sflastainoinen model are applied in case that the processing option Tropcisphere Correction is set to Saaslamoiiien, where Hie geode He lieiglit is approxiniated by the ellipsoitlal heiglit and the relati

29、ve huiniility is fixed to 70 %. 3.5.2电离层模型 Bioadceist ionosphere model For ioaosphErg corretioi for single frqutncy GNSS iisrs, GP5 nnil QZSS navigation tlata iiichi的 tl follcuving broadcast ionospheric paianieters. (E.5.5) By lining these ionospheric parainrs, tlie A ionospheric deky I： (in) can 坯

30、derived the folloivii procedure 一 llie model is often called as Klobiicbar model. ;/= O.OH7 (7 + 01) -D.022 的二clor -V nd a riindom nipflsiirpiupnl erroT vocter F . y = I/x+ r The least square cost function丄 is defuietl as tlw sluii of the squared nieasureutent errors as： (E.6.2) By using (E.C.l) and

31、 (E-丘2)the cost function Cdii be rewritten rts； -yy-yHx-+ yW% (E.63) To minimize (he cost function, the giadient of丄j should be zero.门wn = 一己甘(H丁、屮 + xHH dx =-打+ xd lQdfkifr a旳 H -i-s 少 H二JAJr rs Ee亘 wh伦re khe w伦ighf mMIrrix H iw wet 3 F monWEgh需d L5E)， 3 精选文档 364求解验证和RAIM FDE (5) Solution va lidati

32、on aiid R A IM FDE The estiinrttfftl receiver positicniE described in (3) might iiichiLle invalid soJiitioiis due to iiimiodeleLl nipasurPTlKiit prroi-s. Ib tpsi whptbpT flip valid woliiticm or not and rpjpvt the inwUi 規 solnhonS RTKIJR applies the fcilowiiig valitLatioii after obtaiiwig tlie receie

33、r position estiiiiateLl. If the Vdlidation failed, the sohition is rejected with warning messages. (E6.31) (用 -3； + 血.-cdP +1： +7；?) 叫- (E.6.32) (E,633) (E.6.34) 亠一伽-D 川-JJ - 1 GDOF GDOP仙範 where is hbe niniihcr of pstiiiifitrd parflinctors diid 刖 is the iinniber of iiicasiircnicnts.is clii-square di

34、stribtition of the degree of freedom and a =0.001 (0.1%). GDOP is geometric DCF (dilution of precisiciTi), GDD耳人怦 can be set as the processing option Reject Threslield of GDOP In Rd dition to the solution validation dcscrihcd dbove, RAIM (recouPT antoncunoiis integrity monitoring) FDE (fault detecti

35、on aiid exclusion) function is added in ver.2-4.2- Jf the processing option RAIM FDE is eEiabled and the chi-sqviared test by (E-6.33) is faikd, KTKLIB retries the 附im吊tion by excluding one by one of the visible satellites. After all of retries the estimated receiver position with the inii-iiiimiii

36、normalized stjLiflTed residudls vv is selected as the soliition. In such scheme, an invalid measurement which might be due to satellite malRuiction, receiver fault or largeis excluded as 曰n ouHier. Note that tliis ftaUire dots not eJfectivt Avith two or more invalid measurenents It also needs two re

37、cAiuiciant visible satellifes, tlmt nienns at least which are also provided as tlw ANTEX file typically by IGS international GNSS service) like the rewiver aiilenna KTKLIB also tan import such the antenna model for SAtelliies. The satellite dnk-iuLa PCO is usually given in satellite fixed coordiiial

38、es with respect to the satellite CoM (cent4?r of nia (E.S.3) (E-3.6) (E.8J) r； = X W E* =，白 (E.0.S) (E朋) where fsjfj? IS the siin positicn in ECEF CQardiiiatcs The nadir angle E can also be dcrivrd as: ill flip PPP modp. Hie pffect should be iiicoTporaJeci iu thm mode becHuae the aiiiplitiide someti

39、mes neaclics several 10 cm as the vertical coinpoiwiit. Figure E.8-2 shows the example of the site displdceiiwnt by earth tides- e u- 20 -2U i8 20 -4n 则I MNI Tin 3410Alta IHI -20 Figiirp E.乩2 Di which is deterniiiied as to avoiding cycle-jumps. 3.7.6通过PPP估算接收机的位置 (6) Fstinifltion of rpcpiver position by PPP By u