GPS卫星定位坐标讲义计算与程序设计
卫星导航定位算法与程序设计_常用参数和公式讲解
《卫星导航定位算法与程序设计》课程常用参数和常用公式一览编制人：刘晖最后更新：2010年11月26日1、常用参考框架的几何和物理参数1.1 ITRFyy 主要的大地测量常数长半轴a=6.3781366×106m；地球引力常数（含大气层）GM=3.986004418×1014 m3/s2；地球动力因子J2=1.0826359×10-3；地球自转角速度ω=7.292115×10-5 rad/s。
扁率1/f =298.25642；椭球正常重力位U0=6.26368560×107 m2/s2；γ=9.7803278 m/s2；赤道正常重力e光速c=2.99792458×108 m/s。
1.2 GTRF主要的大地测量常数长半轴a=6.37813655×106 m；地球引力常数GM=3.986004415×1014 m3/s2；地球动力因子J2=1.0826267×10-3；扁率1/f =298.25769。
1.3 WGS84（Gwwww）主要的大地测量常数长半轴a=6.3781370×106 m；地球引力常数（含大气层）GM=3.986004418×1014 m3/s2；地球自转角速度ω=7.292115×10-5 rad/s。
扁率1/f =298.257223563；椭球正常重力位U0=62636860.8497 m2/s2；γ=9.7803267714m/s2；赤道正常重力e短半轴b=6356752.3142m；引力位二阶谐系数2,0C=-484.16685×10-6；第一偏心率平方2e=0.00669437999013；e'=0.006739496742227。
第二偏心率平方21.4 PZ90 主要的大地测量常数长半轴a=6.378136×106m；地球引力常数GM=3.9860044×1014 m3/s2；fM=3.5×108 m3/s2；地球大气引力常数a地球自转角速度ω=7.292115×10-5 rad/s。
GPS卫星坐标计算ppt课件
第二章 GPS卫星位置的计算
2）轨道摄动九参数 卫星平均运动角速度与计算值之差⊿n，或称⊿n为平近地点角速度的改正数；升交点赤经的变化率 ；轨道倾角的变化率 ；升交角距的正余弦调和改正项之振幅Cus、Cuc；轨道倾角的正余弦调和改正项之振幅Cis、Cic；轨道向径的正余弦调和改正项之振幅Crs、Crc。
第二章 GPS卫星位置的计算
*
最后得到各空间坐标系之间的转换关系：
第二章 GPS卫星位置的计算
*
第二章 GPS卫星位置的计算
在卫星运行轨道平面上，以地心为原点，至升交点方向为x轴方向， y轴垂直于x轴指向地极北方向，z轴向上垂直于轨道平面，构成的右手坐标系。由任意观测时刻t的卫星的开交角距为 t时刻卫星的三维坐标为
*
第二章 GPS卫星位置的计算
3) 2个时间参数 从星期日子夜零点开始度量的星历参考时刻toe ；星历表的数据龄期AODE, 即预报星历的外报时间间隔。
*
2.31卫星位置计算推导过程
在用GPS信号进行导航定位以及制定观测计划时，都必须已知GPS卫星在空间的瞬时位置。卫星位置的计算是根据卫星电文所提供的轨道参数按一定的公式计算的。 1、计算卫星运行的平均速度n 根据开普勒第三定律，卫星运行的 平均角速度n0可以用下式计算
*
3、地心空间直角坐标系
以地球自转轴为Z轴，而以地球赤道面为其基准面。由于地球坐标系是随同地球自转的旋转坐标系，显然X轴的方向不能指向与地球自转无关的春分点，而是采用地球赤道面与格林尼威治子午面交线的方向。
第二章 GPS卫星位置的计算
平均经度零点X`T(t)
GAST
平均经度零点X`T(t)
第二章 GPS卫星坐标的计算
2.3 卫星在轨瞬时位置计算
GPS卫星定位坐标计算及程序设计_97
2 2 2 ( X X ) ( YY ) ( Z Z ) c s p s p s p
3.2绝对定位基本方程
X X Y Y Z Z s p 0 s p 0 s p 0 d = d X d Y d Z d t R R R 0 0 0
a
s
2
式中：n和as 为广播星历参数。
用广播星历参数计算卫星位置
• 2、计算归化时间
t t toe
• 3、计算t时刻卫星的平近点角 M
s
M0为广播星历中给出的参考时刻toe的平近点角度 。 Es • 4、计算t时刻卫星的偏近点角 E Me s i n E
s
M Mnt ( t ) s 0 o e
2.3GPS卫星的信号
• 导航电文 导航电文是包含有关卫星的星历、卫星工作状态 时间系统、卫星钟运行状态、轨道摄动改正、大 气折射改正和C/A码捕获P码等导航信息的数据码 （或D码），是利用GPS进行定位的数据基础。 导航电文的内容包括遥测码（TLW）、转换码（ HOW）、第一数据块、第二数据块和第三数据块5 部分。
• 以2009年5月7日南京工业大学江浦校区控 制网20号控制点观测数据为例，来说明如 何利用该程序计算卫星坐标和地面点的近 似坐标。该数据利用华测GPS接收机观测， 观测时间为2小时。
5.实例计算和精度分析
• 利用本程序计算得到的20号控制点的坐标如表5-1 所示：
表5-1 本程序计算的20号点坐标
• 11、计算观测时刻升交点的经度 ＝ － G A S T
( t t ) 0 e 0 e
则：
G A S T = G A S T ( t ) ( tt ) 0 e 0
GPS卫星位置的计算(C++程序计算)
GPS卫星位置的计算利用C++编写了一段能计算单一瞬时卫星坐标的程序，在运行程序之前，需做部分准备工作：（1）在F盘下建立一名为“单一卫星广播星历”的txt文件。
（2）从“广播星历.txt”文件中拷贝从卫星PRN号开始的8行数据到“单一卫星广播星历.txt”中（3）在编辑选项中，将全部的“D”替换为“E”。
下面为我所选取的一个广播星历：18 06 8 25 6 0 0.0-2.472363412380E-04-1.023*********E-12 0.000000000000E+001.410000000000E+02-1.721875000000E+01 4.502687555010E-09 1.413760604187E+00-7.990747690201E-07 7.598234573379E-03 1.118145883083E-05 5.153709835052E+034.536000000000E+05-1.303851604462E-08-1.095067942661E-01 1.527369022369E-079.571235745530E-01 1.640000000000E+02-2.656176299285E+00-8.0374********E-09-5.193073455211E-10 1.000000000000E+00 1.389000000000E+03 0.000000000000E+002.000000000000E+00 0.000000000000E+00-1.024*********E-08 1.410000000000E+024.464490000000E+05 4.000000000000E+00程序设计部分：#include<stdio.h>#include<math.h>int main(){int i = 0;double n[50], n0, nn, t, tk, Mk, Ek, Vk, Yk, Gu, Gr, Gi, uk, rk, ik, xk, yk, zk, X, Y, Z, Lk, UT, yy, mm, JD, gpsz;FILE *fp;fp = fopen("F:\\单一卫星广播星历.txt", "r");if (fp == NULL){printf ("文件打开失败!\n");return 0;}while (! feof (fp)){fscanf(fp, "%lf", &n[i]);i++;}n0 = (sqrt(3986005E+8))/pow(n[17], 3);nn = n0 + n[12];/*计算卫星运行的平均角速度*/UT = n[4] + (n[5] / 60) + (n[7] / 3600);/*民用日的时分秒化为实数时*/if (n[1] >= 80)/*广播星历中年只有后两位，化为4位，参考1980年1月6日0点*/ {if (n[1] == 80 && n[2] == 1 && n[3] < 6){n[1] = n[1] + 2000;}n[1] = n[1] + 1900;}else{n[1] = n[1] + 2000;}if (n[2] <= 2){yy = n[1] - 1;mm = n[2] + 12;}if (n[2] > 2){yy = n[1];mm = n[2];}JD = (int)(365.25 * yy) + (int)(30.6001 * (mm + 1)) + n[3] + (UT / 24) + 1720981.5;/*化为儒略日*/gpsz = (int)((JD - 2444244.5) / 7);/*计算GPS周*/t = (JD - 2444244.5 - 7 * gpsz) * 24 * 3600;/*得出GPS秒*/tk = t - n[18];/*tk1为中间值，用以判断tk与正负302400的关系，然后返回到tk上*/ while (tk > 302400 || tk < -302400){if (tk > 302400){tk = tk - 604800;}else{tk = tk + 604800;}}/*计算归化观测时间*/Mk = n[13] + nn * tk;/*观测时刻的卫星平近点角*/Ek = Mk;Ek = Mk + n[15] * sin(Ek);Ek = Mk + n[15] * sin(Ek);/*迭代两次计算观测时刻的偏近点角*/Vk = atan(sqrt(1 - n[15] * n[15]) * sin(Ek)) / (cos(Ek) - n[15]);/*真近点角*/Yk = Vk + n[24];/*升交距角*/Gu = n[14] * cos(2 * Yk) + n[16] * sin(2 * Yk);Gr = n[23] * cos(2 * Yk) + n[11] * sin(2 * Yk);Gi = n[19] * cos(2 * Yk) + n[21] * sin(2 * Yk);/*摄动改正项*/uk = Yk + Gu;rk = n[17] * n[17] * (1 - n[15] * cos(Ek)) + Gr;ik = n[22] + Gi + n[26] * tk;/*经摄动改正后的升交距角、卫星矢径、轨道倾角*/xk = rk * cos(uk);yk = rk * sin(uk);zk = 0;/*卫星在轨道坐标系的坐标*/Lk = n[20] + (n[25] - 7.29211515E-5) * tk - 7.29211515E-5 * n[18];/*观测时刻t的升交点经度*/X = xk * cos(Lk) - yk * cos(ik) * sin(Lk);Y = xk * sin(Lk) + yk * cos(ik) * cos(Lk);Z = yk * sin(ik);/*卫星在WGS-84坐标系的坐标*/printf("该卫星在WGS-84坐标系中的坐标为：\nX = %lf m\nY = %lf m\nZ = %lf m\n", X, Y, Z);fclose(fp);return 0;}计算结果：该卫星在WGS-84坐标系中的坐标为：X = 9223153.692525 mY = 24133486.931401 mZ = 6032585.919385 m。
卫星导航定位算法与程序设计讲义
第1章
1.1 GPS 算法与程序设计 .......................................................................................... 1 1.2 MATLAB ................................................................................................................ 1 1.3 课程目的及内容 ................................................................................................. 3 第2章 程序设计基础........................................................................................ 4 2.1 程序设计基本思想 ............................................................................................. 4
GPS卫星坐标计算分解
主要内容 2.1卫星坐标系简介 2.2偏近点角E与真近点角f的关系公式推导 2.3卫星位置计算推导过程 2.4二体问题的星位置计算解析
第二章 GPS卫星位置的计算
2.1卫星坐标系简介
一、WGS-84大地坐标系 1、WGS-84大地坐标系定义
WGS-84(World Geodetic System，1984年)是美国国防 部研制确定的大地坐标系，其 坐标系的几何定义是：
顶为正），以子午线方向为x轴（向北为正 ），y轴与x、z轴垂直（向东为正）。
站心赤道直角坐标系与站心地平直角坐 标系之间的关系
X sin B cos L
Y
sin
B
sin
L
Z
站赤
cos B
sin L cos L
0
cos B cos L x
cos
B
sin
L
y
sin B z 地平
站心地平直角坐标系与球心空间直角坐
椭球第一偏心率： e2=0.00669437999013 地球引力常数： GM=（39860050.6）108（m3/s2）
正常化二阶带谐系数：J2=（–484.166851.30）10–9（rad/s） 地球自转角速度： ω=（72921150.1500）10–11（rad/s）
国际大地测量与地球物理联合会(IUGG)——International Union of Geodesy and Geophysics
1 e2 sin E
cos f
,sin f
1 e cos E
1 e cos E
第二章 GPS卫星坐标的计算 2.3 卫星在轨瞬时位置计算
2.3.1 广播星历
卫星导航定位算法与程序设计_单点定位程序流程
GPS 单点定位程序流程亠、计算流程1、读取RlNEX N 文件，将所有星历放到一个列表（数组）ephlst 中2、读取RINEX C 文件,读取一个历元观测值epoch3、数据预处理根据epoch 中的卫星号和历元时刻T R 在ephlst 查找相应的卫星星历 准则 TR-ToEl < 3600.0 s 。
4、程序初始化，置测站概略位置为X r ，接收机钟差初值dt r5、选择epoch 中一颗卫星S i 观测值，设其伪距为T S& 计算卫星S i 的信号发射的概略时刻T S方法如下:a ） 卫星S i 的信号传播时间：∙°S= S i /c -dt r ∙ dt S i；dt S i为卫星钟差，需要进行相对论改正；b ） 卫星S i 的信号发射时刻：T Si=T R - O S i；C ）卫星 S i 在T S i时刻的位置 X S （T Si）=（X S iY S iZS i I T S；d ） 对卫星位置X S iT S i进行地球自转改正，得到X W ST S；e ） 根据X W i（T S i 闲测站概略位置X r 计算卫星和测站的几何距离 R Sf ） 根据几何距离R S i求信号传播时间∙ιS = R S i/ c 。
I X 〕∣Xr1Y-Cdt r 0Zcdt r 一X 0 =第一次迭代，取X 。
=■01 0∣g）如果”I S70S l <10-,则退出迭代。
T l Si=T R―补即为卫星信号发射时刻h）否则-S = .i Si,回带到b）进行迭代。
7、求卫星S i方向余弦SSS(X-X i)bi S i (Y-Y i)-S i (Z-Z i)R S i R S R S8、求卫星S i在观测方程式中的余数项:C dt -∙d trop -'d iono 'D RTCM其中:J S ------ 卫星S i的伪距观测值；R S i――卫星S i到测站的几何距离；C dt S以米表示的卫星S i的钟差；d tro——对流层延迟改正量，单位米，用简化的hopfi模型计算；d iono――电离层延迟改正量，单位米，采用无电离层伪距组合观测值时，此项为0；D RTC M——对伪距的差分改正值，此处为°；9、选择epoch中下一颗卫星S j观测值，设其伪距为10、重复第6—9步，计算每颗卫星的系数和余数项11、将所有卫星的系数组成误差方程，以x,y, Z I Cdt r为未知参数进行求解，形式应该是：AX=Lb S0b S0b2S° 1A= b0 bl: b2 1X = [χy z Cdtl T L=[卢IS …I S I.b S b S i b S 1ji = 0,1, ,svnu mT12、求解法方程— A T PA J A T PL，求出定位结果X 。
GPS卫星位置的计算(C++程序计算)
GPS卫星位置的计算利用C++编写了一段能计算单一瞬时卫星坐标的程序，在运行程序之前，需做部分准备工作：（1）在F盘下建立一名为“单一卫星广播星历”的txt文件。
（2）从“广播星历.txt”文件中拷贝从卫星PRN号开始的8行数据到“单一卫星广播星历.txt”中（3）在编辑选项中，将全部的“D”替换为“E”。
下面为我所选取的一个广播星历：18 06 8 25 6 0 0.0-2.472363412380E-04-1.023*********E-12 0.000000000000E+001.410000000000E+02-1.721875000000E+01 4.502687555010E-09 1.413760604187E+00-7.990747690201E-07 7.598234573379E-03 1.118145883083E-05 5.153709835052E+034.536000000000E+05-1.303851604462E-08-1.095067942661E-01 1.527369022369E-079.571235745530E-01 1.640000000000E+02-2.656176299285E+00-8.0374********E-09-5.193073455211E-10 1.000000000000E+00 1.389000000000E+03 0.000000000000E+002.000000000000E+00 0.000000000000E+00-1.024*********E-08 1.410000000000E+024.464490000000E+05 4.000000000000E+00程序设计部分：#include<stdio.h>#include<math.h>int main(){int i = 0;double n[50], n0, nn, t, tk, Mk, Ek, Vk, Yk, Gu, Gr, Gi, uk, rk, ik, xk, yk, zk, X, Y, Z, Lk, UT, yy, mm, JD, gpsz;FILE *fp;fp = fopen("F:\\单一卫星广播星历.txt", "r");if (fp == NULL){printf ("文件打开失败!\n");return 0;}while (! feof (fp)){fscanf(fp, "%lf", &n[i]);i++;}n0 = (sqrt(3986005E+8))/pow(n[17], 3);nn = n0 + n[12];/*计算卫星运行的平均角速度*/UT = n[4] + (n[5] / 60) + (n[7] / 3600);/*民用日的时分秒化为实数时*/if (n[1] >= 80)/*广播星历中年只有后两位，化为4位，参考1980年1月6日0点*/ {if (n[1] == 80 && n[2] == 1 && n[3] < 6){n[1] = n[1] + 2000;}n[1] = n[1] + 1900;}else{n[1] = n[1] + 2000;}if (n[2] <= 2){yy = n[1] - 1;mm = n[2] + 12;}if (n[2] > 2){yy = n[1];mm = n[2];}JD = (int)(365.25 * yy) + (int)(30.6001 * (mm + 1)) + n[3] + (UT / 24) + 1720981.5;/*化为儒略日*/gpsz = (int)((JD - 2444244.5) / 7);/*计算GPS周*/t = (JD - 2444244.5 - 7 * gpsz) * 24 * 3600;/*得出GPS秒*/tk = t - n[18];/*tk1为中间值，用以判断tk与正负302400的关系，然后返回到tk上*/ while (tk > 302400 || tk < -302400){if (tk > 302400){tk = tk - 604800;}else{tk = tk + 604800;}}/*计算归化观测时间*/Mk = n[13] + nn * tk;/*观测时刻的卫星平近点角*/Ek = Mk;Ek = Mk + n[15] * sin(Ek);Ek = Mk + n[15] * sin(Ek);/*迭代两次计算观测时刻的偏近点角*/Vk = atan(sqrt(1 - n[15] * n[15]) * sin(Ek)) / (cos(Ek) - n[15]);/*真近点角*/Yk = Vk + n[24];/*升交距角*/Gu = n[14] * cos(2 * Yk) + n[16] * sin(2 * Yk);Gr = n[23] * cos(2 * Yk) + n[11] * sin(2 * Yk);Gi = n[19] * cos(2 * Yk) + n[21] * sin(2 * Yk);/*摄动改正项*/uk = Yk + Gu;rk = n[17] * n[17] * (1 - n[15] * cos(Ek)) + Gr;ik = n[22] + Gi + n[26] * tk;/*经摄动改正后的升交距角、卫星矢径、轨道倾角*/xk = rk * cos(uk);yk = rk * sin(uk);zk = 0;/*卫星在轨道坐标系的坐标*/Lk = n[20] + (n[25] - 7.29211515E-5) * tk - 7.29211515E-5 * n[18];/*观测时刻t的升交点经度*/X = xk * cos(Lk) - yk * cos(ik) * sin(Lk);Y = xk * sin(Lk) + yk * cos(ik) * cos(Lk);Z = yk * sin(ik);/*卫星在WGS-84坐标系的坐标*/printf("该卫星在WGS-84坐标系中的坐标为：\nX = %lf m\nY = %lf m\nZ = %lf m\n", X, Y, Z);fclose(fp);return 0;}计算结果：该卫星在WGS-84坐标系中的坐标为：X = 9223153.692525 mY = 24133486.931401 mZ = 6032585.919385 m。