卫星导航定位算法-常用参数和公式

卫星导航定位常用参数和常用公式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.37813655x106m；地球引力常数GM=3.986004415x1014m3/s2；地球动力因子J2=1.0826267x10-3；扁率1/f=298.25769。WGS84（Gwwww）主要的大地测量常数长半轴a=6.3781370x106m；地球引力常数（含大气层）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.16685X10-6；2,0第一偏心率平方e2=0.00669437999013；第二偏心率平方e'2=0.006739496742227。PZ90主要的大地测量常数长半轴a=6.378136x106m；地球引力常数GM=3.9860044x1014m3®；地球大气引力常数fM=3.5x108m3/s2；a地球自转角速度①=7.292115x10-5rad/so扁率1/f=298.257839303；椭球正常重力位U0=6.2636861074x107m2®；赤道正常重力Y=9.780328m/s2；e光速c=2.99792458x108m/s；引力位二阶带谐项系数J0=1.0826257x10-3；2

引力位四阶带谐项系数J0=23709x10-6；4

海平面上由大气引起的重力改正-0.9m/s2。1.52000国家大地坐标系主要的大地测量常数长半轴a=6378137m；地球引力常数GM=3.986004418x10i4m3/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；第二偏心率平方e'2=6.738525414683x10-3。1.71980西安坐标系主要的大地测量常数长半轴a=6.378140x106m；地球引力常数（含大气层）GM=3.986005x1014m3®；引力位二阶带谐系数J2=1.08263x10-3；地球自转角速度①=7.292115x10-5rad/so扁率1/f=298.257；椭球正常重力位U0=6.2636830x10-7m2®；赤道正常重力Y0=9.78032m/s2o第一偏心率平方e2=6.69438499959x10-3；第二偏心率平方e'2=6.73950181947x10-3；1.8EGM2008计算地球重力场时所使用的参数#define#define#define#define#define#defineGM_Earth398600.4415e+9R_Earth6378136.3GM_Sun1.32712438e+20GM_Moon(GM_Earth/81.300587)P_Sol4.560E-6GM_J2-1.08262617385222255846e-32、课程程序常用常数doublePI=(3.1415926535897932384626433832795);/五十位的圆周率:PI=3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679doubleD2R=(0.017453292519943295769222222222222);//PI/180.0doubleR2D=(57.295779513082320876846364344191);//180/PI//L1FrequencyinHz//L2FrequencyinHz//SpeedofLightm/s//Earthrotation(r/s)doubleFREQ_L1=(1575.42E6);//L1FrequencyinHz//L2FrequencyinHz//SpeedofLightm/s//Earthrotation(r/s)doubleFREQ_L2=(1227.60E6);doubleSPEED_OF_LIGHT=(299792458.0);doubleEARTH_ROTATE=(7.2921151467E-5);doubleGM=3.9860047e14xx3、常用公式扁率:a-bf=a离心率:<a2-b2e=a3.1大地坐标换算到高斯平面坐标高斯投影中，某点的大地坐标到高斯平面坐标的转换公式组如下NNx=X+sinBcosBl2+sinBcos3B(5-12+9n2+4n4)140224N720sinN720sinBcos5B(61—58t2+14)16.（B.1）y=NcosBl+cos3B(1-12+n2)136面cos5B(5-18t2+14+14n2-58n2t2)15式中：X,y——该点的高斯投影坐标，单位为米（m）；子午线弧长，单位为米（m）；X=a(1-e2)AB0TOC\o"1-5"\h\zB'X=a(1-e2)AB0—sin2B+—sin4B-—sin6B+—sin8B2468丿451751102543659A=1+—e2+e4+e6+e8+e1064256163846553672765+e1072765+e1065536—e2+—e4+e6+e8165122048N—B,LC'=D'N—B,LC'=D'=15105220510395+-e4+e6+e8+e10642564096163843531531185e6+e8+e105122048131072315+e816384卯酉圈曲率半径，单位为米（m）3645+e1065536a<1-e2sin2B该点的大地纬度和大地经度，单位为弧度（rad）；L——高斯投影带的中央子午线大地经度，单位为弧度（rad）；01——该点大地经度与投影带中央子午线大地经度的经度差，1=L-L单位为弧度0

rad）；a,b参考椭球的长半轴和短半轴，单位为米（m）；a2-b2e2椭球第一偏心率的平方，e2=，无量纲；a2a2-b2e‘2——椭球第二偏心率的平方，e‘2=，无量纲。b2t该点纬度的正切函数值，t=tanB，无量纲；耳——参变量，耳二eScos2B，无量纲。当AL<3.5时，按式（B.1）进行投影的坐标转换，精度约为0.001m。3.2高斯平面坐标换算到大地坐标某点的高斯平面坐标到大地坐标的转换公式如下y2L=L+y2L=L+0NcosB1-碌(5+n；+3厂91需)+360N4(61+90t2+45t；)y2y41-(1+耳2+2t2)+(5+6q2+28t2+212+24t4)6N2ffl20N4fffffB.2)式（B.2）中：B,L——该点的大地纬度和大地经度，单位为弧度（rad）；L——高斯投影带中央子午线大地经度，AL=L-L，单位为弧度（rad）；00B——横坐标（y）在高斯投影带中央子午线上的垂足点的纬度（底点纬度），单位为弧度（rad）；M——底点纬度B处的子午圈曲率半径，m=.a（1^e2），单位为米（m）；fff（1-e2sin2B）2底点纬度B处的卯酉圈曲率半径，N二，单位为米底点纬度B处的卯酉圈曲率半径，N二，单位为米（m）；ffJl-e2sin2Ba，b参考椭球的长半轴和短半轴，单位为米（m）；a2-b2e2椭球第一偏心率的平方，e2=，无量纲；a2a2-b2e‘2——椭球第二偏心率的平方，e'2=，无量纲。b2t——纬度B的正切函数值，t=tanB，无量纲；fff参变量，1=eScos2B，无量纲。（B.2）中的B可采用以下迭代方式计算:a)b)初值（B）a)b)初值（B）=f0按下式迭代:（B匸flaA(1-e2)+丄Bsin2aA(1-e2)A2(B)-Csin4(B)+Dsin6(B)-Esin8(B)f04f06f08f0上式中：a参考椭球的长半轴和短半轴，单位为米（m）；a2-b2e2——第一偏心率的平方，e2=，无量纲；a2TOC\o"1-5"\h\z3451751102543659A=1+—e2+e4+e6+e8+eio464256163846553672765+eio72765+eio65536e2+e4+e6+e816512204815105220510395+e4+e6+e8+eio642564096163843531531185e6+e8+e1o5122048131072c)E=检查迭代结果c)E=检查迭代结果3153645+e8+e101638465536若|6)-(B)|<O.OOOl”(i=1,2,3,fifi-1）则退出，否则返回b）继续迭代。当AL<3.5时，按式（B.2）进行投影坐标转换，转换精度约为0.0001”。3.3空间直角坐标与大地坐标相互转换的数学模型同一坐标系统的空间直角坐标（X,Y,Z）与大地坐标（B,L,H）的转换关系见下式:*(N+H)cosBcosLY二(N+H)cosBsinL(A.4)Z(N(1-e2)+H)sinBZ+N-e2sinBarctan7(X2+Y2)A.5)YarctanA.5)XMX2+Y2)-NcosB式（A.4）和（A.5）中：[XYZ片——空间直角坐标，单位为米（m）；N——卯酉圈曲率半径，单位为米（m）,N=—1-e2sin2BB——大地纬度，单位为弧度（rad）；L大地经度，单位为弧度（rad）；H大地高，单位为米（m）；a2-b2e椭球第一偏心率，e2=a2a椭球长半轴，单位为米（m）；b——椭球短半轴，单位为米（m）。1X1XX1X1XX3.4地心坐标系与站心坐标系的转换关系设：点A(X,Y,Z),点B(X,Y,Z)AAABBB求：以A为原点的站心地平坐标系下B点的坐标算法：1.AX=X-XAB1.AX=X-XABBAAY=Y-YABBAAZ=Z-Z(X,Y,Z)t〈B,L,H)AAAAA-sinBCOsL

AA-sinLACOsBCOsLAA2.A-sinBsinLAACOsLACOsBsinLAA站心直角坐标系与站心极坐标系之间的转换关系cosBA0sinBAAXABAYABAZABra「fZ=<D丿arCtanra「fZ=<D丿arCtan(y/x)

arCtan((xCOsA+ysinA)/z)

(xcosA+ysinA)sinZ+zcosZ丿3.5三参数三维转换数学模型某点从坐标系统1到坐标系统2的三参数数学模型如下T2T1ZZY0AYY：03XXb_LXXAX0Y=Y+AY0ZZAZ21——该点在坐标系统2下的坐标，单位为米(m)；——该点在坐标系统1下的坐标，单位为米(m)；AZ卜——坐标系统2原点相对于坐标系统1原点的平移量,0单位为米3.6七参数三维转换数学模型某点从坐标系统1至坐标系统2的七参数模型如下+(1+m+(1+m)-R(8)-R(8)-R(8)1X2Y3ZXYZA.2)sinsin8XCOS810R(8)=0COS8XXX0-sin80—sin£Y100—sin£Y100cos£Y式中[X[X-T2-TzY0A——该点在坐标系统2下的坐标，单位为米（m）；——该点在坐标系统1下的坐标，单位为米（m）；AZ卜——坐标系统2原点相对于坐标系统1原点的平移量，单位为米0COS£YR(£)=0YYsin£cos£sin£ZZR(£)=一sin£cos£ZZZZ00Y£,£,£——坐标系统2和坐标系统1间的三个欧拉角，单位为弧度（rad）；XYZm——两个坐标系统间的尺度参数，10-9（ppb）。一般情况下，£X，£y，£Z土匀为微小角，因此式（A.2）可简化为下式:-XY-AX01£z-XY-AX01£z—£YX=AY+(1+m)—£1£Y0zXAZ£—£1Z20YX—11A.3)式中[X[XJ7k-zAY——该点在坐标系统2下的坐标，单位为米（m）；——该点在坐标系统1下的坐标，单位为米（m）；AZ卜——坐标系统2原点相对于坐标系统1原点的平移量，单位为米0£,£,£——坐标系统2和坐标系统1间的三个欧拉角，单位为弧度（rad）；XYZm——两个坐标系统间的尺度参数，10-9（ppb）。3.7二维转换数学模型某点从平面坐标系统1到平面坐标系统2的二维相似变换数学模型如下=(1+m)cosasin=(1+m)cosasina—sinacosa1A.3)式中：[xy卜——该点在坐标系统1下的坐标，单位为米（m）；11[xy卜——该点在坐标系统2下的坐标，单位为米（m）；22[axAyJT坐标系统2原点相对于坐标系统1原点的平移量，单位为米（m）；00a——坐标系统1与坐标系统2间的夹角，单位为弧度（rad）；m——两个坐标系统间的尺度参数，10-9（ppb）。3.8通用时转换到儒略日JD二INT【365.25y]+INT[30.6001(m+1)]+D+UT/24+1720981.5MJD二INT【365.25y]+INT[30.6001(m+1)]+D+UT/24-679019.0其中：如果M<2,贝Uy=Y-1,m=M+12；如果M>2,贝Uy=Y,m=M；JD为儒略日，实数；MJD为新儒略日，实数；Y为年，M为月，D为日，UT为日内小时数；INT[]表示取实数的整数部分。3.9儒略日转换到通用时a=INT[jD+0.5]或者a=INT[MJD+2400000.5+0.5]b=a+1537c=INT[(b-122.1)/365.25]d=INT[365.25xc]e=INT[(b-d),30.6001]D=b-d-INT[30.6001xe]+FRAC[JD+0.5](日)M=e-1-12-INT[e/14](月)Y=c-4715-INT[(7+M)；10](年)N=modIiNT[JD+0.5],7}(星期几。N=0,星期一；N=1,星期二;…)其中：JD为儒略日，实数；MJD为新儒略日，实数；Y为年，M为月，D为日，UT为日内小时数(时)；INT[]表示取实数的整数部分；FRAC[]表示取实数部分的小数部分。3.10儒略日转换到GPS时间GPSWEEK=INT[(JD-2444244.5)/7]GPSSECOND=(JD-2444244.5)*3600*24-GPSWEEK*3600*24*7其中：JD为儒略日，实数；如果采用新儒略日计算，贝\MJD=JD-2400000.53.11通用时转换到GPS时间第一步：通用时转换到儒略日（Q）或新儒略日（MJD）第二步：儒略日（JD）或新儒略日（MD）转换到GPS时间3.12新儒略日（MJD）与儒略日（JD）的关系MJD=JD-2400000.53.13天球坐标系转换关系（XYZ-©,九,r）Xcos©ECos九y=rcos©ESin九zsin©X2+y2+z2

arctanyxarctan一Zx2+y24.123456789、4、根据广播星历计算卫星位置和速度1卫星位置计算流程计算轨道长半轴AA丄AI,其中\：'A来自广播星历参数n0计算平均运动角速度n0n0，其中GM为地球引力常数，本课程中设定参数值为3・9860047e14。计算相对于星历参考历元的时间tkt—604800,当t>302400kkt=t—t，koet=\t+604800,当tt=t—t，koekkkt,其它情况kt为所计算卫星位置的时刻，t为星历中的星历参考时刻。oe对平均运动角速度进行改正nn二n+An，其中An来自广播星历。0计算平近点角MkM二M+nUt，其中M0来自广播星历。k0k0计算偏近点角EkM二E—eDsinE，其中e来自广播星历kkk计算真近点角frJ1—e2sinEf=arctancosE—e\丿计算升交角距u'u'=/,其中°来自广播星历（omega）计算升交角距改正数5u，向径改正数5r，轨道倾角改正数5iTOC\o"1-5"\h\zkkku=Csin2u'+Ccos2ukusuc<5r=Csin2u'+Ccos2u'krsrci=Csin2u'+Ccos2u'kisic

10、向径%和轨道倾角ik计算改正后的升交角距10、向径%和轨道倾角iku=u'+6ukkr二A(1—e-cosE)+6r,kkki二i+6i+(IDOT)•t0kk其中e,i0,IODT来自广播星历11、计算卫星在轨道平面上的位置（兀3丿x=rcosuTOC\o"1-5"\h