卫星导航定位算法与程序设计-常用参数和公式讲解

1、卫星导航定位算法与程序设计课程常用参数和常用公式一览编制人：刘晖最后更新：2010年11月26日1、常用参考框架的几何和物理参数1.1ITRFyy主要的大地测量常数长半轴a=6.3781366xl06m；地球引力常数（含大气层）GM=3.986004418x10i4m3；地球动力因子J2=1.0826359x10-3；地球自转角速度=7.292115x10-5rad/s。扁率1/f=298.25642；椭球正常重力位U0=6.2636856Ox1O7m2；赤道正常重力Y=9.7803278m/s2；e光速c=2.99792458x108m/s。1.2GTRF主要的大地测量常数长半轴a=6.37

2、813655x106m；地球引力常数GM=3.986004415x1014m3；地球动力因子J2=1.0826267x10-3；扁率1/f=298.25769。WGS84（Gwwww）主要的大地测量常数长半轴a=6.3781370 x106m；地球引力常数（含大气层）GM=3.986004418x1014m3/s2地球自转角速度=7.292115x10-5rad/so扁率1/f=298.257223563；椭球正常重力位U0=62636860.8497m2/s2；赤道正常重力Y=9.7803267714m/s2；e短半轴b=6356752.3142m；引力位二阶谐系数C=-484.16685X

3、10-6；2,0第一偏心率平方e2=0.00669437999013；第二偏心率平方e2=0.006739496742227oPZ90主要的大地测量常数长半轴a=6.378136x106m；地球引力常数GM=3.9860044x1014m3；地球大气引力常数fM=3.5x108m3/s2；a地球自转角速度=7.292115x10-5rad/s。扁率1/f=298.257839303；椭球正常重力位U0=6.2636861074x107m2/s2；赤道正常重力Y=9.780328m/s2；e光速c=2.99792458x108m/s；引力位二阶带谐项系数J0=1.0826257x10-3；2引力

4、位四阶带谐项系数J0=23709x10-6；4海平面上由大气引起的重力改正-0.9m/s2。1.52000国家大地坐标系主要的大地测量常数长半轴a=6378137m；地球引力常数GM=3.986004418x1014m3/s2；地球自转角速度=7.292115x10-5rad/s；扁率f=1/298.257222101。1.61954年北京坐标系主要的大地测量常数长半轴a=6.378245x106m；短半轴b=6.3568630188x106m；扁率1/f=298.3；第一偏心率平方e2=6.693421622966x10-3；第二偏心率平方e2=6.738525414683x10-3。1.7

5、1980西安坐标系主要的大地测量常数长半轴a=6.378140 x106m；地球引力常数（含大气层）GM=3.986005x1014m3；引力位二阶带谐系数J2=1.08263x10-3；地球自转角速度=7.292115x10-5rad/s。扁率1/f=298.257；椭球正常重力位U0=6.2636830 x10-7m2；赤道正常重力Y0=9.78032m/s2。第一偏心率平方e2=6.69438499959x10-3；第二偏心率平方e2=6.73950181947x10-3；/PI/180.0/180/PI/L1FrequencyinHz/L2FrequencyinHz/SpeedofLi

6、ghtm/s/Earthrotation(r/s)/地球自转角速度2、课程程序常用常数doublePI=(3.1415926535897932384626433832795);doubleD2R=(0.017453292519943295769222222222222)doubleR2D=(57.295779513082320876846364344191);doubleFREQ_L1=(1575.42E6);doubleFREQ_L2=(1227.60E6);doubleSPEED_OF_LIGHT=(299792458.0);doubleEARTH_ROTATE=(7.2921151467

7、E-5);doubleGM=3.9860047e14;doubleomega=7.2921151467E-5;3、坐标系转换公式大地坐标换算到高斯平面坐标高斯投影中，某点的大地坐标到高斯平面坐标的转换公式组如下NNx=X+sinBcosBl2+sinBcos3B(5一12+9n2+4n4)l40224N+sinBcos5B(61-58t2+14)16.（.（B.1）Ny=NcosBl+cos3B(1-12+n2”36N+cos5B(5-18t2+14+14n2-58n2t25120式中：x，y该点的高斯投影坐标，单位为米（m）；X子午线弧长，单位为米（m）；0X=a(1-e2)AB0X=a(1

8、-e2)AB0sin2B+sin4B-sin6B+sin8B2468丿45175110254365964256163846553615525220572765e4+e6+e8+e101651220486553615105220510395一e4+e6+e8+e10642564096163843531531185e6+e8+e1051220481310723153645+e8+e1016384+65536aA=1+e2+e4+e6+e8+eio16384+65536a4BB=e2+D=D=E=N卯酉圈曲率半径，单位为米（m）,N=1一e2sin2BB,L该点的大地纬度和大地经度，单位为弧度（ra

9、d）；L高斯投影带的中央子午线大地经度，单位为弧度（rad）；心L心L-Lo单位为弧度l该点大地经度与投影带中央子午线大地经度的经度差，rad）；a,b参考椭球的长半轴和短半轴，单位为米（m）；a2-b2e2椭球第一偏心率的平方，e2=，无量纲；a2a2-b2b2e2椭球第二偏心率的平方，e2=-b2xxt该点纬度的正切函数值，t二tanB，无量纲；耳参变量，耳二eScos2B，无量纲。当AL3.5时，按式（B.1）进行投影的坐标转换，精度约为0.001m。高斯平面坐标换算到大地坐标某点的高斯平面坐标到大地坐标的转换公式如下：y2L=L+y2L=L+0NcosB1-济(5+n2+3t2一9n2

10、t2)+360N4(61+90t廿45t2)y2y41一(1+耳2+2t2)+(5+6耳2+28t2+212+24t4)6N2ffl20N4fffffB.2)式（B.2）中：B,L该点的大地纬度和大地经度，单位为弧度（rad）；L高斯投影带中央子午线大地经度，AL=L-L，单位为弧度（rad）；00Bf横坐标（y）在高斯投影带中央子午线上的垂足点的纬度（底点纬度），单位为弧度（rad）；M底点纬度B处的子午圈曲率半径，m=.，单位为米（m）；fff（1-e2sin2B）2N底点纬度B处的卯酉圈曲率半径，N=N底点纬度B处的卯酉圈曲率半径，N=，单位为米（m）；fff1-e2sin2Ba,b参考

11、椭球的长半轴和短半轴，单位为米（m）；a2-b2e2椭球第一偏心率的平方，e2=，无量纲；a2e2椭球第二偏心率的平方，e2=ab，无量纲。b2t纬度B的正切函数值，t二tanB，无量纲；ffff耳参变量，耳=e2cos2B，无量纲。a)b)（B.2）中的B可采用以下迭代方式计算:f初值（a)b)（B.2）中的B可采用以下迭代方式计算:f初值（B）=f0按下式迭代:（B）=flaA(l-e2)+sin2aA(1-e2)A2Br(B)Csin4(B)+Dsin6(B)-sin8(B)04f06f08f0上式中：a参考椭球的长半轴和短半轴，单位为米（m）；a2-b2e2第一偏心率的平方，e2=a2

12、，无量纲；TOC o 1-5 h z451751102543659A=1+e2+e4+e6+e8+eio64256163846553672765+eio72765+eio65536e2+e4+e6+e816512204815105220510395+e4+e6+e8+eio642564096163843531531185e6+e8+e1o5122048131072c)E=检查迭代结果c)E=检查迭代结果3153645+e8+e101638465536若)-(B)|0.0001”(i=1,2,3,fi/i-1）则退出，否则返回b）继续迭代。当AL3.5时，按式（B.2）进行投影坐标转换，转换精度

13、约为0.0001”。O3.3空间直角坐标与大地坐标相互转换的数学模型同一坐标系统的空间直角坐标（X,Y,Z）与大地坐标（B,L,H）的转换关系见下式:X(N+H)cosBcosLY=(N+H)cosBsinL(A.4)Z(N(1-e2)+H)sinBZ+N-e2sinBarctanX2+Y2)A.5)YarctanA.5)X/X2+Y2)-NcosB式（A.4）和（A.5）中：XYZT空间直角坐标，单位为米（m）；N卯酉圈曲率半径，单位为米（m）,N=1一e2sin2BB大地纬度，单位为弧度（rad）；L大地经度，单位为弧度（rad）；H大地高，单位为米（m）；a2一b2e椭球第一偏心率，e2

14、=a2a椭球长半轴，单位为米（m）；b椭球短半轴，单位为米（m）。地心坐标系与站心坐标系的转换关系设：点A（X,Y,Z）,点B（X,Y,Z）AAABBB求：以A为原点的站心地平坐标系下B点的坐标算法：AX=X-XABBA1.AY=Y-YABBAAZ=Z-Z2.ABBA(X,Y,2.ABBA(X,Y,Z)T(B-sinBcosLAA-sinLAcosBcosLAAA-sinBsinLAAcosLAcosBsinLAAcosBA0sinBAAXABAYABAZAB站心直角坐标系与站心极坐标系的转换关系SABAABEABJN2+E2+U2SABAABEABABABAB(E=arctanABINAB.

15、(U=arcsmabISAB其中：S为A到B的距离；ABA为A到B的方位角;ABE为A到B的高度角。AB三参数三维转换数学模型某点从坐标系统1到坐标系统2的三参数数学模型如下XXAX0Y=Y+AY0ZZAZ210式中：XYZT该点在坐标系统2下的坐标，单位为米（m）XYZ该点在坐标系统1下的坐标，单位为米（m）；axAYAZt坐标系统2原点相对于坐标系统1原点的平移量，单位为米000m）；七参数三维转换数学模型某点从坐标系统1至坐标系统2的七参数模型如下11XAXnXY=0AY+（1+m）-R（）-R（）-R（）YZ20AZ01X2Y3ZZm）YY：0式中Xm）YY：0式中XXA式FT2FT1

16、ZZARX（X）=Ry（y）=RZ（Z）=0cosX0一sinsinXcosXXcos0-sinYY010，sin0cosYYcossin0ZZ-sincos0ZZ该点在坐标系统2下的坐标，单位为米（m）；该点在坐标系统1下的坐标，单位为米（m）；AZ片一一坐标系统2原点相对于坐标系统1原点的平移量，单位为米0,坐标系统2和坐标系统1间的三个欧拉角，单位为弧度（rad）；XYZm两个坐标系统间的尺度参数，10-9（ppb）。一般情况下，,均为微小角，因此式（A.2）可简化为下式：XYZ+（1+m）+（1+m）1-zYz1-X式中：XYxY式中：XYxY二（1+m）cosasina-sinaco

17、saZF该点在坐标系统2下的坐标，单位为米（m）；ZF该点在坐标系统1下的坐标，单位为米（m）；AXAYAZ卜一一坐标系统2原点相对于坐标系统1原点的平移量，单位为米000m）；,坐标系统2和坐标系统1间的三个欧拉角，单位为弧度（rad）；XYZm两个坐标系统间的尺度参数，10-9（ppb）。3.8二维转换数学模型某点从平面坐标系统1到平面坐标系统2的二维相似变换数学模型如下式中：Yb该点在坐标系统1下的坐标，单位为米（m）；1Y卜一一该点在坐标系统2下的坐标，单位为米（m）；2Ayr坐标系统2原点相对于坐标系统1原点的平移量，单位为米（m）；0-坐标系统1与坐标系统2间的夹角，单位为弧度（r

18、ad）；m两个坐标系统间的尺度参数，10-9（ppb）。4、不同时间标识法之间的转换公式4.1通用时转换到儒略日JD二INT【365.25y+INT30.6001(m+1)+D+UT/24+1720981.5MJD二INT【365.25y+INT30.6001(m+1)+D+UT/24-679019.0其中：如果M2,贝Uy=Y,m=M；JD为儒略日，实数；MJD为新儒略日，实数；Y为年，M为月，D为日，UT为日内小时数；INT表示取实数的整数部分。4.2儒略日转换到通用时a=INTjD+0.5或者a=INTMJD+2400000.5+0.5b=a+1537c=INT(b-122.1)/365

19、.25d=INT【365.25xce=INT(b-d)30.6001D=b-d-INT【30.6001xe+FRAC【JD+0.5(日)M=e-1-12-INTe/14(月)Y=c-4715-INT(7+M)10(年)N=modIiNTJD+0.5,7(星期几。N=0,星期一；N=1,星期二;)其中：JD为儒略日，实数；MJD为新儒略日，实数；Y为年，M为月，D为日，UT为日内小时数(时)；INT表示取实数的整数部分；FRAC表示取实数部分的小数部分。4.3儒略日转换到GPS时间GPSWEEK=INT（JD-2444244.5”7GPSSECOND=（JD-2444244.5）*3600*24

20、-GPSWEEK*3600*24*7其中：JD为儒略日，实数；如果采用新儒略日计算，贝则MJD=JD-2400000.54.4通用时转换到GPS时间第一步：通用时转换到儒略日（Q）或新儒略日（MJD）第二步：儒略日（JD）或新儒略日（MD）转换到GPS时间4.5新儒略日（MJD）与儒略日（JD）的关系MJD=JD-2400000.55.11、2、3、4、5、6、7、8、9、5、根据广播星历计算卫星位置和速度卫星位置计算流程计算轨道长半轴AA=（；A，其中龙A来自广播星历参数计算平均运动角速度n0n再，其中GM为地球引力常数，本课程中设定参数值为3-9860047e14计算相对于星历参考历元的时

21、间tktt=tt，koet604800,kt=302400k当t302400k其它情况t为所计算卫星位置的时刻，t为星历中的星历参考时刻。oe对平均运动角速度进行改正nn=n+An，其中An来自广播星历。0计算平近点角MkMk=M0+电七，其中Mo来自广播星历。计算偏近点角EkM=E一sinE，其中e来自广播星历TOC o 1-5 h zkkk计算真近点角才Q1-e2sinEf=arctanI4，其中e来自广播星历，Ek来自第6步。 HYPERLINK l bookmark103 o Current Document （cosE-e丿k计算升交角距uu=f+3，其中来自广播星历（omega）计

22、算升交角距改正数5u，向径改正数5r，轨道倾角改正数5ikkku=Csin2u+Ccos2ukusuc5r=Csin2u+Ccos2ukrsrci=Csin2u+Ccos2ukisic10、计算改正后的升交角距叫，向径rk和轨道倾角i10、u-u+5uTOC o 1-5 h zk/k、r-A（1一ecosE）+5rkkki-i+5i+（IDOT）tk0kk其中e,i0,IODT来自广播星历11、计算卫星在轨道平面上的位置（叫Jx-rcosuvkkky-rsinukkk12、计算改正后的升交点经度LkLO+（Q）t-Qt，k0earthoe其中:Q来自星历（OMEGA）,Q来自星历（OMEGAD

23、ot）;0为地球自转角速度，数值为7.292115e10-5,eartht为计算的历元时刻（周秒），t为星历参考时刻（周秒）。oe13、计算卫星在地心坐标系下的坐标（&Z）X-xcosL一ycosisinLkkkkkvY-xsinL+ycosicosLkkkkkZ-ysinikk根据广播星历计算卫星速度计算卫星速度时，需要用到计算卫星位置的一些量，下述公式中符号的含义与4.1相同。】、计算偏近点角的时间导数Ekn1一ecosEk2、计算真近点角的时间导数f/1/1+e）1/21e丿cos2-2_Ecos2-2r3、计算升交角距的时间导数u4、计算升交角距向径和轨道倾角的时间导数u=u(1+2C

24、cos2u-2Csin2u)kusucr=EAesinE+2u(Ccos2u一Csin2u)kkkrsrcI=IODT+2u(Ccos2u-Csin2u)kisic其中，其中IDOT,a，CcCus，Cuc,Crs，Crc来均自广播星历cosLksinLk0cosLksinLk0-sinLcosikkcosLcosikksinik-(xsinL+ycosLcosi)kkkkcosL-ysinLcosiysinLsinikkk-ycosLsinikkkycosikk5、计算卫星轨道平面位置的时间导数xcosu一rsinurk=kkkkysinurcosuukLkkkJLk6、计算卫星在地心坐标系

25、下的速度_X_xkY=RykLZk1丄Ik其中：Lk=Q-earth，其中来自广播星历，%为地球自转角速度。计算卫星高度角和方位角设：原点位于测站P,Z轴与P点的椭球法线重合，X轴垂直于Z轴并指向椭球短轴,Y轴垂直于XPZ平面向东，构成左手系。1、卫星S在以P点为原点的站心平面直角坐标的坐标（N,E,U）S表示为H=-sinBH=-sinBcosLpp-sinLpcosBcosLpp-sinBsinLppcosLpcosBsinLppcosBp0sinBp一N_AX_(XX、E=HAY=HY-YUcAZcnZcZQ丿其中：XZ卜为卫星的地心坐标系坐标;XZ卜为卫星的地心坐标系坐标;sXZ卜为观

26、测站的地心坐标系坐标;pB,L为测站的大地纬度和大地经度。pp2,卫星S相对于测站P的方位角和高度角、距离为:A=arctan色E=arctanE=arctanU：N2+E2=arcsin其中：r为卫星到测站P的距离，A为方位角，E为高度角6、地球自转改正由于地球自转，卫星位置将发生变化，设(Xs，Ys，Zs)为卫星在某时刻的位置，则由于地球自转而对卫星坐标带来的变化量(AX,AY,AZ):SSS_AX_S-0WT0_一XSAY=-WT00YSSAZ000ZS11S对星站距离带来的变化叽为:Ap=(X-X)Y-(Y-Y)X其中:WcSRSSRS其中:为地球自转角速度；(X,Y,Z)为卫星在信号

27、发射时刻的地心坐标；SSS1(X,Y,Z)为接收机的地心坐标；RRR(AX,AY,AZ)为信号到到用户接收机的时刻(t2),卫星位置的改正数;SSS2T-t-1，为信号传播时间。21其中:XSY其中:XSYSCOS(WT)-sin(WT)0sin(WT)Cos(WT)0XSYSZSw为地球自转角速度;(X,Y,Z)为卫星在信号发射时刻的地心坐标；SSS1(X,YZ)为卫星在信号接收时刻匚的改正后的地心坐标;SSS2T-t-t，为信号传播时间。217、对流层改正模型7.1Hopfield模型直接得到视线方向路径上的对流层改正量Trop=A+AwK=155.2-10-7P(h-H)dTdssK=1

28、55.2.10-748e(h-H)wT2swssh=40136+148.72(T-273.16)md0h=11000mw其中：A测站到卫星视线方向对流层总延迟量，单位：米；TropH测站高度，单位：米；sT测站上的绝对温度，单位：开氏温度；sP测站上的大气压，单位：毫巴；sRH测站的相对湿度，无单位；se测站上的水汽压，单位：毫巴；sT参考面的绝对温度，单位：开氏温度；0E卫星相对于测站的高度角，单位：度。投影函数在编程中应按照下式计算（式中E的单位为度）：E2+6.25E2+6.25丿=sinE2+6.25*兀/180.0sinsinG2+2.25)=sinC;E2+2.25)*k/180.0)若给出某参考面上的温度、气压和湿度，则测站处的气象元素由以下公式计算得到：T=T-0.0065.(H-H)TOC o 1-5 h zs0s0P=P-0.000022