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 expr 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 与 UTCUn i versa I Tim e Coord i n at e d 关系参考【图1】，参考【图2】图1转换关系公式Table E.l-1 GPST-UTC Values until Marcli 2013Time Since in LTTCrLS s1980-01-0600 000001981-07-0100000011982-07-0100000021983-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之前得转换可以用 更准确得表达方式，如【图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号000000】。北斗时间计算公式【图4】图4UTC与GPST时间转换同上而得GPS-样,只不过UTC参数来自与北斗导航信息中。1.3坐标系统接收机与卫星得位苣在RTKLIB中表示为在ECEF（地心地固坐标系）坐标系统中得X.Y.Z组件、1.3.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,xf rah 1cos Jl_/snF0.图 8133本地坐标到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 ntpos1 s atposs -2 e s tp o s - 3 estvel1.计算计算卫星位置、速度与时钟sa t p ossa通过广播星历计算卫星钟差ephclk 卫星编号到卫星系统得转换sa t sys根据卫星得编号，获取到对应得卫星导航系统。
选择星历seleph1、传入信号传输时间，卫星编号，导航数据等参数、2、遍历导航数据，遍历导航数据里而得星历数据，判断星历数据得卫星编号 就是否与传入得卫星编号相等。3、如果星历数据得卫星编号与传入得卫星编号相等，就计算星历参考时间 t oc与信号传输时间得时间差。如果不相等,继续处理下一条星历数据、4、判断计算出来得时间差，如果时间差大于了允许得最大时间差，继续査找下一个星历数据、否则，判断时间差最小得星历数摇，记录星历数拯得位宜、5、返回之前记录出来得星历数据乙 使用广播星历计算卫星时钟偏差e ph2c 1 k1、传入信号发射时刻得时间与星历数拯。2、讣算信号发射时刻得时间与本时段钟差参数参考时间星历参数sc 得时间差。3、通过下式讣算钟差，这里还没有处理相对论校正项与tgd代码中有个 迭代过程，资料上没瞧到写。卫星賈差汁算出來之信号发射时时间还要加上这个钟.b计算卫星在信号发射时刻得位宜、速度与时钟satp o s 根据星历表选项来选择不同得处理，如下 广播星历E PHO PT_BRDCephpos 广播星历到卫星位置与钟差1、根据公式计算出tk;2、根据使用得卫星系统，选择使用得地球引力常数muj地球得角速度omge3、根据公式计算出平近点角M、4、求解开普勒方程，按照以下公式迭代求解。5、根据以下公式汁算出11改正后得纬度幅角,I改正后得径向,i改正后 得轨逍倾角得值、J1-以 sinEacosE - o arctail dll CU5 sin 20 Cuc cos 1/dr Crs sill 20 Crc cos 2pSi C15 sin 2 Clc cos 2u p 6u7 dl-ccosE 丹i z0 Jz6、根据以下公式计算卫星在轨道平面内得坐标、7、根据不冋得卫星系统，做不同得讣算、GPS计算方式C Q _叫0 Cte Lcosu co/2 -siinz cos sin I尸5r r cos/ sin/ sin// cosf cosZ2 sin u sini北斗计算方式Q 人一叫 tQecos cos一sill cos sin2厂心（-5。）com sin/ sin;/ cosj cosQ.sinw sin/00cos 9sinZ0COS2、汁算卫星坐标与接收机坐标得差值向量。3、用差值向量做欧几里徳范数，再用上一步计算出来得差值向量与计算结果做除法，得到视线向量。按照如下公式计-也rsts-rrt.4、使用以下公式计算几何距离。Pr H b - r W | -yC-Vr -ySXr计算卫星方位角/仰角satazel1、把接收机ccef坐标转换到大地坐标；2、判断髙度就是否大于地球半长轴WGS8 4得负数值；3、如果髙度小于等于地球半长轴WGS84得负数值，方位角为0.仰角为P 1/24、如果高度大于地球半长轴WGS84得负数值，把ECEF向量转换到局 部坐标。然后对转换出来得坐标做内枳。5、使用如下公式计算卫星方位角与仰角、er,emi 四 ，勺勺Az j.ATAN2ee,euElj.arcsin处理，保证tt得范用0 t t 8 6400o6、计算倾斜因子。7、电离层延迟计算，公式如下。3兀2力-50400/工0网/；0Fx5x10 一 Fx5xl0924r x x 1 224H1.57悴 1.57对流层改正tropcorr 通过标准大气压与saas t amoi n en模型计算对流层延迟t r o pm o del 1、使用以下公式计算总气压。2、使用以下公式汁算绝对温度、3、使用以下公式汁算水蒸汽得分压、代码中没有使用、4、使用以下公式计算Saa s ta mo i nen模型】丁S 0.002277 f 1255 门“、、耳p 0.05 e-tan-z cosz l T丿 伪距残差残差值伪距一物理距离dtr-光速冬时钟偏差电离层误差对流层误差时钟系统与接收器得偏置补偿 伪距测呈误差方差v arerrb 方差权重值weigh t b y va r i a neec 最小二乘估计1 sq 最小二乘估计通过求解正规方程X A * A* a-1*A* Y1、计算矩阵A*Y得结果.保存到矩阵Ay;2、讣算矩阵A* A得结果，保存到矩阵Q;3、求得矩阵Q得逆矩阵.结果保存到矩阵Q4、最后xQ* Ayd 欧几里徳范数norm公式1、通过最小二乘法估算出来得参数dx【长度为4】，分别加到位置向量上x【长 度为4】。2、对这个速度向量做欧几里徳范数，返回值同1E-4比较，小于这个值时,就得到 估算出得接收机得位置为向量X得值。e 验证求解valsol1、对伪距残差值做内积，然后同卡方分布得自由度分布值alph a0 .00 1 做比较，大于卡方分布得自由度分布值得数据,无效。2、计算d ops ；暂时还没没找到资料。3、用计算出来得dop s与配置得最大d o ps阈值比较，大于配置得最大d o p s 阈值得数据无效。3.使用多普勒估算接收机速度es t vela 多普勒残差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 *Y1、计算矩阵A*Y得结果，保存到矩阵A y ;2、计算矩阵A * “得结果，保存到矩阵Q;3、求得矩阵Q得逆矩阵，结果保存到矩阵Q;4、最后 xQ * Ayc 欧几里徳范数norm公式3、通过最小二乘法估算岀来得参数dx【长度为4】，分别加到速度向量上x【长度为4】。4、对这个速度向屋做欧几里徳范数,返回值同1E 6比较,小于这个值时，就得到估算岀得接收机速度为向量x得值。2.2精确定位pppos1uds 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 pa 位置更新b 时钟更新c 对流层参数更新d 相位偏差更新1、通过LLI检测周跳；2、通过Geomet ryF r X相位跳变检测周跳如果双频测量值可用；2.计算卫星得位程与速度与时钟satp o ssc通过广播星历计算卫星钟差ephclk 卫星编号到卫星系统得转换sat s ys根摇卫星得编号，获取到对应得卫星导航系统、选择星历（s eleph）1、传入信号传输时间，卫星编号，导航数据等参数。2、遍历导航数据.遍历导航数据里而得星历数据，判断星历数据得卫星编号 就是否与传入得卫星编号相等。3、如果星历数据得卫星编号与传入得卫星编号相等，就计算星历参考时间 （toe）与信号传输时间得时间差。如果不相等,继续处理下一条星历数据。4、判断计算出来得时间差，如果时间差大于了允许得最大时间差,继续查 找下一个星历数据。否则,判断时间差最小得星历数据，记录星历数拯得位 置。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、根据使用得卫星系统，选择使用得地球引力常数（mu）与地球得角速度 （omge）3、根据公式计算出平近点角4、求解开普勒方程，按照以下公式迭代求解。5、根据以下公式讣算出u（改正后得纬度幅角）,1（改正后得径向）,i（改正后 得轨逍倾角）得值、cosE-o0 arctana e 创CU5 sin 20 Cuc cos lfdr Cs sill 20 Crc cos 205i Cs 血 2 Clc cos 2 u7 dl-ccos 77 Z0 十刃 6、根据以下公式计算卫星在轨道平而内得坐标匚7、根据不同得卫星系统，做不同得讣算。GPSil-算方式x /q /2-de tk -oe roecos” cos 2 - sitn/ cos sin/2 r r r cosi/ sinjQ sin// cosz cos2、sin m sin ZJ北斗计算方式G 0 Gg 一叫 loecos cos/-siii cosz sin/2rJ f rRzder0Ev-5 cosm sin/ sinz/ cosz cos/2.sinw sin其中100、cos9sine0尽0COS2、计算卫星坐标与接收机坐标得差值向S。3、用差值向量做欧几里徳范数，再用上一步汁算出来得差值向量与计算 结果做除法，得到视线向量、4、使用以下公式计算几何距离。Pr -6 -心产| - 5备注这个公式就是地球自转改正公式。地球自转改正产生地球自转改正得机制在于地固系随地球得自转而旋转.地 固系就是非惯性系统，用地固系中得坐标计算卫星到接收机几何距离时.由于 卫星位置与接收机位置就是两个不同时刻得位宜矢疑，而且这两个时刻得地 固系相对于惯性系就是变化得，因此要考虑地球自转引起得距离变化、其改 正公式为心，二彳-“兀 1 一 兀 r XS2- Xsl- Xr2XSJ其中，w为地球自转角速度;C为光速Xsi与Xri分别代表卫星位置矢绘与测 站位巻矢屋得分i 1,2对应X,y分量。地球自转对纬度影响甚小，对经度影 响最大，英次就是髙度、同时地球得自转影响还跟测站得经纬度以及测站与卫星 得几何关系有关，对两极得测站，影响为零；对赤道上得测站影响最大；当卫星在 测站子午而内影响为零;卫星在测站东方时，影响为负，卫星在测站西方时,影响为正。c 卫星方位角/仰角计算satazel1、把接收机ccef坐标转换到大地坐标；2、判断高度就是否大于地球半长轴WGS84得负数值；3、如果髙度小于等于地球半长轴WGS 8 4f负数值，方位角为0，仰角 为 PI/24、如果高度大于地球半长轴W G S84得负数值，把ECE F向量转换到局 部坐标。然后对转换出来得坐标做内积。5、使用如下公式计算卫星方位角与仰角。r.emi 四 ee en euAZj.ATAN2ee,eJJEll arcsiiiQ“d 对流层延迟修正 通过标准大气压与saastamo i n en模型计算对流层延迟t r o pm odel1、使用以下公式计算总气压、2、使用以下公式计算绝对温度、3、使用以下公式计算水蒸汽得分压、代码中没有使用。z 1255 J、P 0.05e 一 tair z1T,4、使用以下公式计算Saas tamo i ncn模型】。0.002277COS 2e 卫星天线模型satantpcvf 接收机天线模型，通过天线相位中心参数计算天线偏移a n tmo d elg 相位缠绕校正windupc o r r h 电离层与天线相位校正测量corrmea s i 计算卫星时钟与对流层延迟5.卡尔曼滤波f i Iter卡尔曼滤波按照以下得公式更新状态K P * H *H;* P * H R1;xp x K * v;Pp I- K *Hr *P;参数解释X状态向虽n X 1p状态得协方差矩阵n x nII设讣矩阵得转1矩阵n x mv 创新测量模型m x 1R测虽渓差得协方差矩阵m x mn.m状态与测量值得数量xp状态更新后得向星nx 1 Pp状态更新后得协方差矩阵n x n6.组合残差r e s_pp p 7.解析PPP整周模糊度ppp a mba aver a g e LC average_LCb 固左宽巷模糊f ix_a mb_WL 计算LC波长度mlam_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接收机与卫星天线之间得几何距离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 derived fromE.3.7The 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 56v； F at the time lr and r/s .v.y5.ziT at the time ts in the ECEF earth center earth fixed coordinates asE-3.8P；二帆曲“一附where rz is the ECEF to ECI earth center inertial coordinates transation 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 can 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 effectE.3.8aPr 一 Rz /c” r ||pr E lim EE418E.4.19E4.20Tlie broadcast ephemerides and clock are applied in case that the processing option Satellite Ephemeris/aockM to Broadcast** as well as GLONASS,BeiDou and SB AS.3.4北斗广播星历与时钟3 Broadcast ephemerides and docks for BiDou 门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/.7pE.4.28For MEO and IGSO satellites of BeiDou,the same ulations as 1 for GPS ephemeris and clock,except for 3.98600441Sxl014,7.2921150 x 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 equationE.5.13 and E.5.14 should be replaced by2 2 Cg -咲 Qcos/ cosP-sun/ cof sm/2 rst cos/ sin2 sin” cos/ cos/2sin// sin/F.4.29E.4.30whereu00COSsm00cos 6sin,R二一 sincos 800X一 sin0COSZ0X0-3.5对流层与电流层模型3.5.1对流层模型The standard atmosphere can be expressed asE.5.1E5.2E53p 1013.25x1 - 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