GPS 卫星位置的计算 利用 C++编写了一段能计算单一瞬时卫星坐标的程序，在运行程序之前，需做部分准备工作： （1）在 F 盘下建立一名为“单一卫星广播星历”的 txt 文件。 （2）从“广播星历.txt”文件中拷贝从卫星 PRN 号开始的 8 行数据到“单一卫星广播星历.txt” 中 （3）在编辑选项中，将全部的“D”替换为“E”。 下面为我所选取的一个广播星历： 18 06 0.0-2.472363412380E-04-1.023181539495E-12 0.000000000000E+00 8 25 6 0 1.410000000000E+02-1.721875000000E+01 4.502687555010E-09 1.413760604187E+00 -7.990747690201E-07 7.598234573379E-03 1.118145883083E-05 5.153709835052E+03 4.536000000000E+05-1.303851604462E-08-1.095067942661E-01 1.527369022369E-07 9.571235745530E-01 1.640000000000E+02-2.656176299285E+00-8.037477650349E-09 -5.193073455211E-10 1.000000000000E+00 1.389000000000E+03 0.000000000000E+00 2.000000000000E+00 0.000000000000E+00-1.024454832077E-08 1.410000000000E+02 4.464490000000E+05 4.000000000000E+00 程序设计部分： #include #include= 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 m Y = 24133486.931401 m Z = 6032585.919385 m