RTKPPP定位算法流程

1.1 1基础知识GPS精密单点定位的基本原理GPS精密单点定位一般采用单台双频仔?$接收机，利用IGS提供的精密星历和卫星钟差，基于载波相位观测值进行的高精度定位。观测值中的电离层延迟误差通过双频信号组合消除，对流层延迟误差通过引入未知参数进行估计。1.2时间系统RTKLIB内部使用GPST(GPST时间)用于GNSS的数据处理和定位算法。数据在RTKLIB内部处理之前，需要转换成GPST时间。使用GPST的原因是避免处理润秒。RTKLIB使用以下结构体表示时间：typedefstructtime_ttime;/*time(s)expressedbystandardtime_t*/1.2.1 doublesec;/*fractionofsecondunder1s*/gtime_t;GPST和UTC(UniversalTimeCoordinated)关系参考【图1】，浇参考【图2】：krrc-!gps-PST=lure+图1转换关系公式TableE.l-1GPST-UTCValues(untilMarch2013)TimeSince(inUTC)1980-01-0600:00:0001981-07-0100:00:0011982-07-0100:00:0021983-07-0100:00:0031987-01-0100:00:0041988-01-0100:00:0051990-01-0100:00:0061991-01-0100:00:0071992-07-0100:00:0081993-07-0100:00:0091994-07-0100:00:00101996-01-0100:00:00111999-01-0100:00:00131997-07-0100:00:00122006-01-01。图(20:0014通过使用GPS导航信息中的UTC参GpST到UtC或者UTC到GpST之前的转换可以用更准确的表达方式，如【图3】。L7TC-Ss-Ms-或)+-S-604800(丽V-所罚)are&上，必Luem这些参数是由GPS导航消息提供的。1.2.2 BDT（北斗导航卫星系统时间）BDT（北斗导航卫星系统时间）是一个连续的时间系统，没有润秒。开始历元的时间是【UTC2006年1月号oo：oo：oo】。北斗时间计算公式【图4】：UTC和GPST时间转换同上面的GPS-样，只不过UTC参数来自与北斗导航信息中。1.3坐标系统接收机和卫星的位置在RTKLIB中表示为在ECEF（地心地固坐标系）坐标系统中的X,Y,Z组件。1.3.1大地坐标到ECEF坐标的转换转换公式如【图5】。*=/(2-/):广S1I1(r+力)co皿cosZ?.牛=(M5Mkv(l-)sm4J第三个公式最后一行有错，应该为：(v(l-e2)+h)sin已图5参数说明：a:地球参考椭球的长半径f:地球参考椭球的扁平率h:椭球高度“：纬度&:经度当前版本的RTKLIB使用的值为【图6】：=6378137.0(m)/=1.0/298.257223563图6图7参考椭球体1.3.2 ECEF坐系到大地坐标的转换转换公式如【图8】r=+.cie2tank二矣血Ilk7、根据不同的卫星系统，做不同的计算。GPS计算方式：r?二口口+(双一叫)一rtoecosncos2-sinwc-os;inQ择(f)=rcosusin1-sinucosicos12sinifsin?北斗计算方式：口=Q-口小-命cos/?-siiifycast或n?K()=R-1如*)L(-V)cositsin/2-sinz/cos?cos/?sin或n?100、cost?sin30|查)=0cos0sin8,5=-sinQcost?0、0sin8853,o011其中8、按照公式二计算出时间tc。9、按照以下公式计算出钟差和钟漂。1. dT*(r)=af+aftc+afc-a/T号inE-bTQDdts(t)=cifx+2af.c精密星历(EPHOPT_PREC):peph2pos广播+SBAS(EPHOPT_SBAS):satpos_sbas广播+SSR_APC(EPHOPT_SSRAPC):satpos_ssr广播+SSR_COM(EPHOPT_SSRCOM):satpos_ssrQZSSLEX星历(EPHOPT_LEX)：lexeph2pos使用伪距估算接收机的位置，返回估算状态结果(estpos)伪距残差(rescode)把ecef坐标系转换成大地坐标系(ecef2pos)暂时还没看懂。，1、按照以下公式做转换再=lim再jA=AT,4N2y.)计算几何距离和接收机到卫星的单位矢量(geodist)1、用卫星的坐标向量做欧几里德范数，返回值和地球长半轴(WGS84)比较。小于地球长半轴(WGS84)，返回-1；2、计算卫星坐标和接收机坐标的差值向量。3、用差值向量做欧几里德范数，再用上一步计算出来的差值向量和计算结果做除法，得到视线向量。按照如下公式：es_k（，）-&（&）4、使用以下公式计算几何距离。用+皿）rW）|+（/好）计算卫星方位角/仰角（satazel）1、把接收机ecef坐标转换到大地坐标；2、判断高度是否大于地球半长轴（WGS84）的负数值；3、如果高度小于等于地球半长轴（WGS84）的负数值，方位角为0，仰角为PI/2;4、如果高度大于地球半长轴（WGS84）的负数值，把ECEF向量转换到局部坐标。然后对转换出来的坐标做内积。5、使用如下公式计算卫星方位角和仰角。AZj.=ATAN2（匕r%）ElI=arcin（%）伪距使用编码残差改正（prange）暂时没找到对应的文档对应。电离层改正（ionocorr）通过广播电离层模型（klobuchar模型）计算出电离层延迟（ionmodel）1、校验传入的电离层模型参数，校验失败，使用默认的电离层模型参数;2、使用以下公式计算出地球为中心的角度（半圆）材=0.013?3+0,11）-0.0223、使用以下公式计算子的电离层的纬度/经度（半圆）；曾=f?-(.7C0S.4Z4=-+.写虬4如COS%4、使用以下公式计算地磁纬度。务=豹-+0.064cos(-L617)5、计算本地时间，返回值按这个公式【tt-=floor(tt/86400.0)*86400.0;】处理，保证tt的范围(0=ttGO5、根据以下公式计算出口（改正后的纬度幅角），（改正后的径向），i（改正后的轨道倾角）的值。_Jl注血EcosE-t?。=arcraila-co(5ii=Cj.jsin2+Ccos2亩=Crssin2。一Crccos203i=qsin2+Ciccos20ii=_Shv=口(1ecosE)。厂i=和-&+i“6、根据以下公式计算卫星在轨道平面内的坐标。Xk=rkcosukk二矣血土7、根据不同的卫星系统，做不同的计算。GPS计算方式：Q=纬十叫)4一叫如cqwcos?-sinifcos?sill尸！f)=rcos?y技-%inhcosi8勺12simfsini北斗计算方式：口=G-云8SMcosx-100Lr-38.45J14、使用以下公式计算LSaastamoinen模型】。r=/2-E：；.e) 卫星天线模型(satantpcv)接收机天线模型，通过天线相位中心参数计算天线偏移(antmodel)相位缠绕校正(windupcorr)电离层和天线相位校正测量(corrmeas)计算卫星时钟和对流层延迟卡尔曼滤波(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)1. 组合残差(res_ppp)解析PPP整周模糊度(pppamb)averageLC(average_LC)固定宽巷模糊(fix_amb_WL)计算LC波长度(m)(lam_LC)宽巷模糊宽巷模糊度的方差fixnarrow-laneambiguity【ARmode:PPP-AR】fixnarrow-laneambiguitybyILS【ARmode:PPP-ARILS】3以下是上面一些算法的参考资料，信息来自于Rtklib的手册3.1接收机和卫星天线之间的几何距离(4)GeometricrangebetweenreceiverandsatelliteantennasThegeometricrangeisdefinedasthephysicaldistancebetweenthesatelliteantennaphasecenterpositionandthereceiverantennaphasecenterpositionintheinertialcoordinates.Atfirst,thesignaltransmissiontimetscanbederivedfrom:尸=E-玲/c-打(L)(E.3.7)geometricrangecanbeexpressedbyusingthereceiverandsatelliteantennaphasecenterpositions弓(弓)=习)atthetimetrandrts)=(xs.ysatthetimetsintheECEF(earthcenterearthfixed)coordinatesas:Pr=|(6kr(6)-跖/)/()(E.3.8)whereU(，)istheECEFtoECI(earthcenterinertial)coordinatestransformationmatrixatthetimet.FortheexpressionintheECEFcoordinates,theearthrotationeffecthastobeincorporatedintoobtainthegeometricrange.Theequationcanbeapproximatedbyoneofthefollowingequationswithadequateprecisionunder1mmlevel.CurrentversionRTKLIBalwaysusestheequation(F.3.8b)forthegeometricrange.Thelasttermin(F.3.8b)issometimescalledasSagnaceffect.(E3.8a)(E3.8b)(E3.8c)Pr*&(弓)-足房*/)prH|r(r；)尸(r)+牛(有-%)+ECEFatsignaltiansmissionECEFatsignalreceptionFigureE.3-1GeometricRangeandEarthRotationCorrection.3.2卫星方向的方位角和仰角(5)AzimuthandelevationanglesofsatellitedirectionTheunitLOS(line-of-sight)vectorfromthereceivertothesatellitecanbeexpressedintheECEFcoordinatesas:W”r)尸(广)一弓0.)(E.3.9)Intheequation,theearthrotationeffectisneglected.TheazimuthandelevationanglesAz.andEl*ofthesatellitedirectionfromthereceiversitecanbederivedfrom:r9enu=ATAN2(ee,en)(E.3.10)(E3.il)El*=arcsin(ew)(E3.12)Pr。(弓)耳0-)扩,珈whereEristhecoordinatesrotationmatrixfromECEFtothelocalcoordinatesatthereceiverpositionReferE.2fordetailedformationofthematrix.x(E)y(N)FigureE.3-2LocalCoordinatesandAzimuthandElevationAngles(E.4.2)(E.4.3)(E.4.4)(E.4.5)(E.4.6)(E.4.7)(E.4.8)(E.4.9)(E.4.10)(F.4.11)(E.4.12)(E.4.13)(E.4.14)(E.4.15)3.3GPS、Galileo和QNSS广播星历和时钟BroadcastephemeridesandclocksforGPS,GalileoandQZSS【5】BroadcastephemerisandSVclockparametersforGPS,GalileoandQZSSaregiveninnavigationmessagesas:“咖（膈膈1）=（，。0，/20，0眼0，小急，。*。”00ThebroadcastephemeridesandclockareappliedincasethattheprocessingoptionSatelliteEphemeris/ClocktoBroadcastaswellasGLONASS,BeiDouandSBAS.3.4北斗广播星历和时钟(3)BroadcastephemeridesandclocksforBeiDou1ForBeiDousatellites,thesimilarephemerisandclockparametersasGPS,GalileoandQZSSareprovidedinthenavigationmessagesas:Pephoe*切*皿0*办”，。,Cus*uc9rs，GrGs，Gc,afq,Tgd)(E.4.28)ForMEOandIGSOsatellitesofBeiDou,thesameformulationsas(1)forGPSephemerisandclock,exceptfor以=3.986004418xlO14,叱=7.2921150xIO-5rad/sandthetimetisexpressedinBDT.Toobtainthesatellitepositionrs(Z)ofBeiDouGEOsatellitesatthetimetinBDT,theequation(E.5.13)and(E.5.14)shouldbereplacedby:Q=ao*htk-饥te(F-4.29)cos/cos/2-sinzzcos/sin/?rs(/)二rR-(coetk)Rx(-5)cos”sin/2+sincoszcos/2(E.4.30)、sinzzsin?)where:qo0、cos。sin。0、)二0coWsinQ,5=sin9cos0、0sin8cos幻、o0b3.5对流层和电流层模型3.5.1对流层模型Thestandardatmospherecanbeexpressedas:p=1013.25x(l-2.2557xl0-5A)525687=15.0-6.5x10-3/7+273.15e=6.108xexp1.?7)-jj(E.5.13)3.6单点定位3.6.1线性最小二乘估计RTKLIBemploysaniteratedweightedLSE(leastsquareestimation)fortheSingle(singlepointpositioning)modewithorwithoutSBAScorrections.(1)LinearLSEAssumeameasurementvectoryaregivenanditcanbemodeledasthefollowinglinearequationsofanunknownparametervector*andarandommeasurementerrorvectorv.y=Hx+v(E.6.1)TheleastsquarecostfunctionJlsisdefinedasthesumofthesquaredmeasurementerrorsas:(E.6.2)Jls=vi+v2+vm2=八(E.6.3)Byusing(E.6.1)and(E.6.2),thecostfunctioncanberewrittenas:=,rry-yJHx-xTHTy+xTHTHxTominimizethecostfunction,thegradientofJsshouldbezero.Then(E.6.4)=0T-yTH(丑弓沪+(HTHx)T+xTHTHdx=-2yTH+2xTHTH=0(E.6.5)HTHx=HTyTosolvethenormalequation,wecangettheestimatedunknownparametervectorxbytheLSEas:x=(HTHylHTy(E.6.6)Iftheweightsofeachmeasurementsaregiven,thecostfunction(E.6.3)canberewrittenbyusingaweightmatrixW.Jwls=Wv(E.6.7)TominimizethecostfunctionJwls,wecanobtaintheestimatedunknownparametervectorbytheweightedLSEbythesimilarwayforthesimpleLSEas:X=(HtWHYxHTWy(E.6.8)TheweightmatrixWfortheweightedLSEisoftengivenas:卬-=g(b2,b2-2,.,b2)wheretisthea-prioristandarddeviationofthei-thmeasurementerror.Tosolvethemeasurementequationtoobtainthefinalestimatedreceiverpositionandthereceiverclockbias,RTKLIBemploystheiteratedweightedLSEas:+1=金+卬丑)TJ，-(&)(E.6.22)FortheinitialparametervectorfortheiteratedweightedLSE,justall0areusedforthefirstepochofthesinglepointpositioning.Onceasolutionobtained,thepositionisusedforthenextepochinitialreceiverposition.FortheweightmatrixIf,RTKLIBusesthefollowingformulasw-枷g(bi-2,b2-2,.,b2)(E.6.23)b=Fb/sinElrj+cepij+弓0初+b.op+bbig(E.6.24)where:Fs:satellitesystemerrorfactor(1:GPS,Galileo,QZSSandBeiDou,1.5:GLONASS,3.0:SBAS)%:code/carrier-phaseerrorratioo:carrier-phaseerrorfactoraandb(m)eph:standarddeviationofephemerisandclockerror(m)Gion:standarddeviationofionospherecorrectionmodelerror(m)Jop:standarddeviationoftropospherecorrectionmodelerror(m)。瞄:standarddeviationofcodebiaserror(m)Forthestandarddeviationofephemerisandclockerror,URA(userrangeaccuracy)orsimilarindicatorsareusedinRTKLIB.Byseveraliterations,thesolutionisconvergedinthenormalcaseandtheestimatedreceiverposition弓andthereceiverclockbiasdtrareobtained.(E.6.25)x二limXj=GF,c瓦VZ00Theestimatedreceiverclockbiasdtrisnotexplicitlyoutputtothesolutionfile.Instead,itisincorporatedinthesolutiontime-tag.Thatmeansthesolutiontime-tagindicatesnotthereceivertime-tagbutthetruesignalreceptiontimemeasuredinGPST.3.6.2估算接收机的位置和钟差EstimationofreceiverpositionandclockbiasForSinglemodeasPositioningMode,thefollowingsingle-pointpositioningprocedureisappliedtoobtainthefinalsolutionbyepoch-by-epochbasis.Foranepochtime,theunknownparametersvectorxisdefinedas:x=(马丁W，(E.6.19)Thepseudorangemeasurementvectorycanbegivenas:)，=(船32,月3,矿产(E.6.20)whereP；isthepseudorangemeasurement.IftheprocessingoptionIonosphereCorrectionsettolono-FreeLC,theionosphere-freeLC(linearcombination)pseudorangedefinedinAppendixE.5(7)isused.Inothercases,justtheApseudorangeisused.FigureE.6-1SatelliteGeometryforSinglePointPositioningThemeasurementequationanditspartialderivativematrixforthesinglepointpositioningareformedas:p).+cdtycdT+】+TyJ官1、+cdtr-cdT2+I+T：一崎1A(x)=pj+cdtycdT，+tH二-&1(F.6.21)uwherethegeometricrangep*andLOSvectoraregivenbyE.3(4)andE.3(5)withthesatelliteandreceiverpositions.TliesatellitepositionsrsandtheclockbiasesdTsarealsoderivedfromtheGNSSsatelliteephemeridesandclocksdescribedinE.4accordingtotheprocessingoptionSatelliteEphemeris/Clock.Tosolvethemeasurementequationtoobtainthefinalestimatedreceiverpositionandthereceiverclockbias,RTKLIBemploystheiteratedweightedLSEas:+1=金+卬丑)TJ，-(&)(E.6.22)FortheinitialparametervectorfortheiteratedweightedLSE,justall0areusedforthefirstepochofthesinglepointpositioning.Onceasolutionobtained,thepositionisusedforthenextepochinitialreceiverposition.FortheweightmatrixIf,RTKLIBusesthefollowingformulasw-枷g(bi-2,b2-2,.,b2)(E.6.23)b=Fb/sinElrj+cepij+弓0初+b.op+bbig(E.6.24)where:Fs:satellitesystemerrorfactor(1:GPS,Galileo,QZSSandBeiDou,1.5:GLONASS,3.0:SBAS)%:code/carrier-phaseerrorratioo:carrier-phaseerrorfactoraandb(m)eph:standarddeviationofephemerisandclockerror(m)Gion:standarddeviationofionospherecorrectionmodelerror(m)Jop:standarddeviationoftropospherecorrectionmodelerror(m)。瞄:standarddeviationofcodebiaserror(m)Forthestandarddeviationofephemerisandclockerror,URA(userrangeaccuracy)orsimilarindicatorsareusedinRTKLIB.Byseveraliterations,thesolutionisconvergedinthenormalcaseandtheestimatedreceiverposition弓andthereceiverclockbiasdtrareobtained.(E.6.25)x二limXj=GF,c瓦VZ00Theestimatedreceiverclockbiasdtrisnotexplicitlyoutputtothesolutionfile.Instead,itisincorporatedinthesolutiontime-tag.Thatmeansthesolutiontime-tagindicatesnotthereceivertime-tagbutthetruesignalreceptiontimemeasuredinGPST.3.6.3估算接收机的速度和钟漂Estimationofreceivervelocityandclock-driftIfDopplerfrequencymeasurementsofGNSSsignalsaregiven,receivervelocitiesandclock-driftscanbeestimatedthefollowingprocedure.Foranepochtime,theunknownparametersvectorxforthevelocityestimationisdefinedas:x=(vf.,cdir)(E.6.26)wherevranddarethereceivervelocityinECEF(m/s)andthereceiverclock-drift(s/s),respectively.Therange-ratemeasurementvectorJcanbegivenas:y=(-&以i-(E.6.27)whereDj.tisthe。Dopplerfrequencymeasurementofthesatellites.RTKLIBalwaysusesDopplerfrequencymeasurements.Tliesemeasurementequationsanditspartialderivativematrixareformedas:+cdty_cdT一#f+cdiy_cdT?昭1+cdiy_cdT，H=-借1+cdir-cdTn=&rrrr(F.6.28)Tlierange-rageVj.betweenthereceiverandthesatelliteintheseequationsisderivedfrom:rr-ei-T()f)+(vrYr+-vbr-/%,)(F.6.29)s,sss、r-Twherev-(vx,vv,vz)andvr二(vrr,vvr).ByusingtheiteratedLSEsimilartotheestimationofthereceiverposition,wecanobtainthereceivervelocityandclock-driftas:a*=limXj=(vrr,cdir)T(E.6.30)wheretheweightmatrixWissettoI(non-weightedLSE).(3) 3.6.4求解验证和RAIMFDESolutionvalidationandRAIMFDETheestimatedreceiverpositionsdescribedin(3)mightincludeinvalidsolutionsduetoumuodeledmeasurementerrors.Totestwhetherthevalidsol