本文内容主要包括卫星位置计算原理、所需文件及其格式内容、数据结构、计算步骤、编程流程、算例及结果分析的阐述，以及笔者对编程过程中一些常见问题和注意事项的总结。
（注：该程序单纯用于卫星位置计算的学习，在GPS伪距单点定位程序中用到的卫星位置计算模块需要在此基础上做一些调整和改进）
代码详见：
https://download.csdn.net/download/m0_58307078/88622825?spm=1001.2014.3001.55031 实验任务及目的
了解星历文件的基本格式及内容，编写程序，实现基于广播星历的GNSS卫星位置计算，使用PPP精密单点定位原理实验平台或orbit.txt文件检验结果的正确性，完成实验报告，内容包括：1）数据来源；2）处理过程及结果；3）精度评价分析；4）存在问题及解决方法等。
2 数据来源及编程测试环境
（1）数据来源 ：GNSS广播星历文件
（2）编程环境：MATLAB 2016a
（3）测试环境：MATLAB 2016a及PPP精密单点定位原理实验平台
3 实验内容
3.1 实验准备
（1）了解星历文件格式及内容
表3-1 RINEX格式的广播星历数据
其中的数据含义见表3-2：
表3-2主要符号说明
（2）计算原理
利用卫星星历计算卫星位置的原理及计算流程如下：
图3-1 计算流程示意
Step1计算卫星运行的平均角速度
1）卫星的平均角速度：
其中，地球引力常数GM取398600.4km3/s2 （WGS-84坐标系中）
2）卫星运行的平均角速度：
其中，Δn为平均角速度改正数，由星历文件给出。
Step2计算观测瞬间卫星的平近点角：
Ms=M0+n(t-toe)
Step3解开普勒方程求偏近点角：
Es=Ms+esinEs
Step4计算真近点角fs
Step5计算升交距角：
u0=ω+fs
Step6计算轨道摄动改正项：
Step7计算摄动改正值：
Step8计算卫星在轨道坐标系中的坐标：
Step9计算观测瞬间t升交点的经度：
其中，地球自转角速度ω=7.2921151467×10-5rad/s
Step10计算卫星在WG-S84坐标系中的空间直角坐标:
3.2 编程流程
程序编写思路及步骤如下图：
图3-2 编程流程示意
（1）确认星历中所需数据的变量符号
程序所需的星历数据变量定义如表3-3所示：
表3-3 星历中所需变量符号
| id | / | / | / | |
| / | Crs | dlta_n | M0 | |
| Cuc | e | Cus | sqrta | |
| toe | Cic | Omega0 | Cis | |
| i0 | Crc | omega | dOmega | |
| di | / | / | / | |
| / | / | / | / | |
| / | ||||
说明：
- 符号在表中位置与对应参数表3-1 RINEX格式的广播星历数据中的位置一致；‘/’表示该程序中不需要此位置对应的这些参数；
- id：卫星的PRN号+星历发布时刻（年月日时分秒），如：'G05 2021 11 12 23 59 44'。
（2）读取星历文件、导入所需参数
思路：考虑到封装函数、算法复杂度要求及避免浪费资源的需要，不进行一次性读取并导入整个星历文件；而是通过识别某些标识，一次仅从文件中读取并存储所需的8行星历数据（即：指定卫星在指定星历的数据）。
算法支持：MATLAB自带的textscan()函数适合读入大文件，不仅可以从文件的任何位置开始读入，也可以从上一次textscan()调用结束的位置开始读入，且其返回值为元胞数组，适合存放既有文本又有数字的星历数据。
具体实现：
图3-3 对应代码截图（部分）
如图3-3代码所示，打开文件后，通过使用textscan()函数跳过第1-第id_r-1行，按指定格式读取第id_r行并转换为数值型矩阵R1；由于再次调用了具有相同文件标识符（fid）的textscan()恢复文件的文本扫描（间图3-4），textscan()会在上次终止读取的位置（第id_r行，即所需数据开始的第一行）自动恢复读取，接着读取所需数据的第二行，以此类推，读取并存储完所需的8行数据。
图3-4 对应代码截图（续 图3-3）
所需数据导入结束，将表3-3中所需星历参数按对应位置分别赋值给相应变量：
图3-5 变量赋值
（3）编写计算程序
按照3.1(1)中计算步骤及公式编写程序，其中需要特别注意的几点如下：
- 当前时刻t与参考时刻toe时间差的计算：
- 星历文件中的toe本身为GPS时；
- 需要输入的当前时刻t可直接输入GPS时，若输入为世界协调时，则需要进行时间转换；
- 计算时间差tk=t-toe时需要进行归化，以保证tk在区间[-302400，302400]内。
- 计算偏近点角Es求解开普勒方程Es=Ms+esinEs可用迭代法或直接使用MATLAB的fsolve()函数。
（4）初步测试（测试计算程序）
封装函数前首先要确保程序计算结果正确，调试过程中发现的主要问题为：计算升交点经度时，地球自转角速度精度影响卫星坐标精度。
取ω=7.2921151467×10-5rad/s时，卫星XYZ坐标误差在2米甚至1米以内，而取ω=7.292115×10-5rad/s（而许多资料只精确到这个程度）时，误差有时高达10米！
参考位置（视作真实位置）：（-18295994.4110,-8270444.0440,17998338.1510）
图3-6 ω不同取位G07卫星的计算结果与真实位置的对比
原因分析：ω用于计算升交点经度L，而时间差t-toe与toe的数量级达105，在公式 L=Ω0+Ω-ωt-toe-ωtoe的作用下，由ω取位误差造成的L的误差数量级约为0.1（rad），而具有0.1（rad）误差的L在cosL的作用下继续被放大，导致旋转矩阵存在误差，进而使得最终计算得到的卫星坐标误差高达10米。
（5）封装函数
思路：实现任意时刻t卫星位置的计算至少需要输入卫星的PRN号和时刻t两个参数，另外还需要选择一个星历，由于星历的参考时刻toe与其PRN号后的年月日时分秒一一对应，因此可将卫星的PRN号+年月日时分秒作为识别所需的8行星历数据位置的唯一标识（id,格式：'G05 2021 11 12 23 59 44'）。
算法支持：将MATLAB中的find()函数与strcmp()函数联用可以简洁地实现在元胞数组中查找特定元素位置（索引）的功能。
具体实现：id_r=find(strcmp(index,id)) 其中，strcmp(index,id)返回一个与元胞数组同维度大小的逻辑数组，再用find()找到其中的非零元素索引即为唯一标识id所在行数。其中，index为星历文件的第一列，如图3-7：
图3-7 index（部分）
封装后函数形参：id,t；返回值：卫星坐标XYZ列向量
图3-8 函数头
（6）测试函数
思路：1）输入id时格式要求较为严格，易出错，因此不仅需要编写输入格式错误时触发的操作，最好还能提供可以输入的选择；2）函数计算结果输出的显示形式应尽量清晰明了，需要多次调试；
具体实现：
图3-9 对应代码（部分）
交互效果：
图3-10函数交互效果
4 实验成果
由于星历文件obs_binex与参考位置文件orbit不完全对应（orbit文件中只有某些特定时刻的GPS卫星位置），此处选择在两个文件中都数据齐全的卫星作为代表，分析其在t=579600时刻的位置计算结果及精度。
4.1 算例展示
所选卫星在t=579600时刻的算例如表4-1：
表4-1 t=579600.0时刻的算例结果
| PRN_ID | X/m | Y/m | Z/m |
| G01 2021 11 13 16 00 00 | -13126065.600 | 10724714.960 | 20077603.740 |
| G04 2021 11 13 16 00 00 | -18593526.340 | 4111319.775 | -18545523.930 |
| G07 2021 11 13 16 00 00 | -6216287.487 | 23647611.720 | -9365135.821 |
| G08 2021 11 13 16 00 00 | -24265269.020 | -4672456.072 | 10150426.500 |
| G09 2021 11 13 16 00 00 | -7865809.128 | 13269809.260 | -21696100.150 |
| G14 2021 11 13 16 00 00 | -2992350.580 | 17589754.360 | 19655268.810 |
| G16 2021 11 13 16 00 00 | -20580243.250 | -2818188.912 | -17067771.050 |
| G17 2021 11 13 16 00 00 | 5639176.624 | 18361592.970 | 18774852.980 |
| G22 2021 11 13 16 00 00 | -23027408.010 | 8017556.105 | 10447590.330 |
所选卫星在t=579600时刻的位置如下图所示：
图4-1 t=579600.0时刻的算例卫星位置
4.2 分析与评价
编程计算出的卫星位置与orbit文件提供的参考位置之间的偏差用vX、vY、vZ表示。由表4-2，Y、Z坐标误差均值较小，但方差和峰值较大且分布分散（范围大），综合考虑均值和方差，X坐标计算值精度要优于Y、Z坐标。
坐标误差分布如图4-2所示：显然，除个别点外，计算误差都在2米以内，大多数在1.5米以内，结果表明计算结果精度正常，符合误差在2米左右的要求。
| 表4-2 坐标误差的描述性统计 | |||||||
| 个案数 | 范围 | 最小值 | 最大值 | 平均值 | 标准差 | 方差 | |
| vX | 9 | 2.41 | -0.42 | 1.99 | 0.9113 | 0.72601 | 0.527 |
| vY | 9 | 3.16 | -0.95 | 2.21 | 0.2302 | 1.11536 | 1.244 |
| vZ | 9 | 4.33 | -1.94 | 2.39 | 0.0786 | 1.40444 | 1.972 |
| v | 27 | 4.33 | -1.94 | 2.39 | 0.4067 | 1.13493 | 1.288 |
图4-3 t=579600.0时刻算例卫星位置误差
表4-2、表4-3、表4-5分别为卫星X、Y、Z坐标值及对应误差：
表4-4 t=579600.0时刻算例的X坐标误差
| PRN_ID | 参考值X0/m | 计算值X/m | 误差vX/m |
| G01 | -13126066.390 | -13126065.600 | 0.790 |
| G04 | -18593527.350 | -18593526.340 | 1.010 |
| G07 | -6216289.476 | -6216287.487 | 1.989 |
| G08 | -24265270.410 | -24265269.020 | 1.390 |
| G09 | -7865810.160 | -7865809.128 | 1.032 |
| G14 | -2992351.088 | -2992350.580 | 0.508 |
| G16 | -20580242.830 | -20580243.250 | -0.420 |
| G17 | 5639175.021 | 5639176.624 | 1.603 |
| G22 | -23027408.310 | -23027408.010 | 0.300 |
表4-5 t=579600.0时刻算例的Y坐标误差
| PRN_ID | 参考值Y0/m | 计算值Y/m | 误差vY/m |
| G01 | 10724715.910 | 10724714.960 | -0.950 |
| G04 | 4111319.709 | 4111319.775 | 0.066 |
| G07 | 23647611.240 | 23647611.720 | 0.477 |
| G08 | -4672457.422 | -4672456.072 | 1.350 |
| G09 | 13269810.070 | 13269809.260 | -0.810 |
| G14 | 17589755.250 | 17589754.360 | -0.890 |
| G16 | -2818191.121 | -2818188.912 | 2.209 |
| G17 | 18361593.350 | 18361592.970 | -0.380 |
| G22 | 8017555.105 | 8017556.105 | 1.000 |
表4-6 t=579600.0时刻算例的Z坐标误差
| PRN_ID | 参考值Z0 | 计算值Z | 误差vZ/m |
| G01 | 20077604.620 | 20077603.740 | -0.880 |
| G04 | -18545524.350 | -18545523.930 | 0.420 |
| G07 | -9365137.438 | -9365135.821 | 1.617 |
| G08 | 10150426.790 | 10150426.500 | -0.290 |
| G09 | -21696101.180 | -21696100.150 | 1.030 |
| G14 | 19655269.240 | 19655268.810 | -0.430 |
| G16 | -17067773.440 | -17067771.050 | 2.390 |
| G17 | 18774854.190 | 18774852.980 | -1.210 |
| G22 | 10447592.270 | 10447590.330 | -1.940 |
代码详见：
注：该程序单纯用于卫星位置计算的学习，在GPS伪距单点定位程序中用到的卫星位置计算模块在此基础上做了一些调整和改进：
1705
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
