1. 引言
在GPS定位中，为了确定地面点的坐标，必须先确定卫星的坐标。卫星坐标的计算主要依靠卫星的星历。GPS星历分为精密星历和广播星历两类。它们各有优缺点，前者精度较高，但是存在延时，主要用于事后处理，例如IGS精密星历，通常采用SP3格式，以卫星在地固坐标系中离散点的坐标(和速度)提供给用户使用；后者的经度较低，但是能够实时转发给用户，用于实时的导航和定位。广播星历由GPS的地面运行控制中心站计算得到，以广播星历参数的形式提供给用户使用 [1] 。
目前，许多国家都在建设自己的卫星导航系统，导航卫星的广播星历拟合是卫星导航系统中一个重要的技术环节 [2] 。它的选择和设计不仅决定了导航卫星星历所能达到的精度，而且还决定了导航用户算法的复杂程度 [3] 。本文用精密星历进行广播星历参数的拟合，并对计算过程的方程病态性加以分析、并利用奇异值截断法进行计算改进，得到一个可靠的结果。
2. 广播星历表达式
GPS广播星历具有16个参数，包括星历参考时刻
；6个开普勒根数：卫星轨道半长轴平方根
、卫星轨道偏心率
、参考时刻
的轨道倾角
、当前时间周开始时刻的升交点赤经
、近地点角距
；参考时刻
的平近点角
；9个摄动参数：平近地点角速度的改正数
、升交点赤经的变化率
、轨道倾角的变化率
、升交角距的正、余弦调和改正项振幅
、
、轨道倾角的正、余弦调和改正项振幅
、
、地心距的正、余弦调和改正项振幅
、
。
GPS星历表中的时间和坐标分别属于GPS时间系统和WGS-84坐标系。在GPS数据处理中，常利用卫星的星历表参数来计算卫星天线相位中心在WGS-84坐标系中的位置，其具体计算过程及相应的表达式如下 [1] [4]
1) 半长轴：
2) 卫星平均角速度：
3) 瞬时历元到参考历元的时间差：
4) 平角速度：
5) 平近点角：
6) 偏近点角
：
7) 真近点角
：
8) 维度参数：
9) 周期改正项：
10) 改正后的纬度参数：
11) 改正后的向径：
12) 改正后的倾角：
13) 升交点的经度：
14) 卫星在轨道平面内的坐标：
15) 卫星在WGS-84坐标系中的坐标：
其中
和
分别为地球引力常数和地球自转角速度。
3. 广播星历拟合方法
3.1. 计算模型
将上节中卫星的坐标表示成位置向量：
(1)
因此广播星历计算卫星坐标可以表示为如下形式 [5] ：
(2)
其中
是要计算卫星坐标的时刻，已知。参数
也是作为已知量给出的，因此只需对余下15个参数求偏导数便可得到线性化的观测方程：
(3)
由(3)式可得，单历元的误差方程为：
(4)
其中
令
其中
为观测历元总数。由此可得总的观测误差方程为：
(5)
式中
含有15个参数，因此只要
即可根据最小二乘原理 [1] [6] ，利用间接平差解得广播星历参数的改正数：
(6)
(7)
根据(5)、(6)、(7)式进行间接平差迭代求解。迭代的收敛条件为：
(8)
其中，
是预先给定的阈值，取1.0e−05，
是第
次迭代的单位权方差。
3.2. 数值导数计算系数矩阵B
B阵为(3)中观测方程系数的矩阵形式，可以通过解析法求偏导计算，但解析法较为复杂不利于计算机运算，本文采用数值导数的方式计算。
由偏导数的定义：
(9)
当
足够小的时候，这一导数可以近似写成数值导数的形式
(10)
结合本文系数阵B阵的实际，则有
(11)
因此误差方程系数阵
第一列为：
(12)
其他偏导系数同理可依次求出，此处略。
4. 算例分析
本文采用GPS星期第1765周第一天即2013年11月03日的IGS精密星历作为观测值，拟合当日格林尼治时间12 h参考时刻的广播星历参数。IGS精密星历提供的卫星坐标时间间隔为15分钟。由计算模型可知，只要历元数大于5便能根据最小二乘原理进行间接平差。本文用不同的历元数进行拟合计算，并比较不同拟合历元数所得广播星历的精度。将拟合的广播星历参数计算参考历元时刻的卫星坐标与精密星历对应时刻的坐标差异作为广播星历精度的评价指标。最后将最优的拟合参数与当日卫星所发送的广播星历进行比较，说明实际的应用意义。
15个广播星历参数中9个摄动参数的初值设为0，6个轨道根数的初值如表1所示。
本文数值导数的增量为1.0e−08。
表1. 六个轨道根数的初值
4.1. 不同拟合历元个参数变化
不同拟合历元，各参数的变化情况略有不同，但在某一特定历元大部分参数出现了曲线的突变，这样的变化有助于我们找出最佳的拟合历元。下面以卫星11号，2013年11月03日格林尼治12:00的广播参数拟合为例，说明各参数不同拟合历元的参数变化，各参数的变化曲线见图1。
从各参数的变化情况可以看出，参数变化的出现在两个历元出。其一为9个拟合历元处，另一个为17个拟合历元处。其中在9个拟合历元处出现转折的参数为：
，而在17个拟合历元处出现转折的参数分别为：
。从参数的变化情况来看，可以初步断定最佳的拟合历元数应该在9到17个历元之间。参数在特定的拟合历元数发生突变的原因有待于进一步的研究。
4.2. 拟合参数与实际广播星历参数的比较
为了比较拟合广播星历的实际应用价值，下面将拟合计算所得的参数与相同参考时刻相同卫星发布的广播星历进行比较，比较结果见表2。
表2中给出了卫星11号2013年11月03日格林尼治12：00为参考时刻的，15个拟合历元数过拟合的结果与当日卫星广播的结果的比较。从表中可以看出，拟合参数均达到了较高的精度，证明算法的正确性以及使用性。同时结果中存在不足的是，对于极小量级的参数的拟合效果不佳。这个主要原因为计算机浮点数的精度有限，因此为了提高参数的拟合参数因该加大浮点数的位数，以获得高精度的拟合结果。
5. 结论
本文推导了拟合广播星历的方法，并且在求取误差方程系数阵时引入数值导数的方法。最后通过算例证明了方法是正确的。
通过算例的分析可知，不同的拟合历元数所得参数的结果具有差异，并且参数值在特定的拟合历元数出现明显的突变与转折。该转折的具体原因有待进一步的研究。不过笔者断定转折有助于我们决定最
Table 2. Comparison of fitting parameters and parameters of satellite broadcasting
表2. 拟合参数与卫星广播参数的比较
佳的拟合历元数。最后本文将拟合的结果与实际的广播参数进行了比较。从比较的结果中可以看出，应用本文所述的方法得到了较高精度的广播星历参数，具有很高的实际运用价值。但是极小量级的参数收到了计算机浮点数精度的影响，未能达到很高的精度。所以应该增加浮点数的位数，以获得高精度小量级参数。
基金项目
国家自然科学基金，编号：41201425。