顾及GPS轨道凹凸性的标准化方法改进研究*

赵自强 刘志平

(中国矿业大学国土环境与灾害监测国家测绘地理信息局重点实验室，徐州 221116)

1 引言

IGS 提供的SP3 精密星历产品间隔为15 分钟，远低于导航定位应用的采样率(30 s、15 s 甚至更高)，故需对卫星轨道进行标准化处理［1］。常用的轨道标准化方法有插值算法［2，3］和拟合方法［1，4，5］，在相同阶数的情况下，拟合方法较插值算法的标准化结果更为稳定［5］。常规拟合方法主要有Chebyshev 多项式［6］、Legendre 多项式［1］、三角函数［7］等，当采用多项式拟合方法进行轨道标准化时，过高的阶数易产生“Runger现象”［8］，对此，许多学者［2，9，10］研究滑动式拟合方法，有效地提高了轨道标准化精度及可靠性。但是，滑动式拟合方法采用逐节点移动处理策略，一般只保留滑动区间中心结果，从而降低了数据处理效率。文献［4，11］发现，卫星轨道在天球坐标系下具有较好的周期规律，进而提出了顾及卫星轨道周期特征的滑动式拟合方法。但是，该方法需进行天球坐标与地固坐标的相互转换，轨道标准化精度易受转换矩阵的影响。此外，最佳拟合阶数与轨道弧段特征、节点数及拟合精度紧密相关［12］，但现有研究仅考虑拟合精度并得出拟合阶数为9 阶［5］、10 ～11 阶［1］、11 ～12 阶［9］、15阶［13］等，缺乏合理统一的定阶准则。而且，常规拟合方法中多项式系数的求解多采用等权最小二乘［1，13］，没有利用IGS 事后精密星历产品中的坐标中误差信息。综上分析，本文提出了顾及地固坐标系下轨道凹凸性特征的快速加权拟合方法，并利用实际数据验证了该方法的有效性与正确性。

2 常规拟合方法

由于Chebyshev 多项式拟合、Legendre 多项式拟合方法的轨道标准化精度是等效的［1，5］，故以Chebyshev 多项式为例进行讨论。

Chebyshev 多项式为最佳一致逼近。要在区间［t0，t0+Δt］进行n 阶Chebyshev 多项式拟合，需先对时间变量t∈［t0，t0+Δt］进行变换:

式中，τ∈［－1，1］，t0为初始历元时间，Δt 为拟合区间长度。

卫星X 坐标分量可用q 阶Chebyshev 多项式表示为:

式中，Ci为Chebyshev 多项式系数，Ti(τ)有如下递推规律:

T0(τ)=1，T1(τ)=τ，…，

Tn(τ)=2τTn－1(τ)－Tn－2(τ)

卫星Y、Z 坐标分量的求解方法与X 坐标分量相同。

3 拟合方法的改进

3.1 凹凸区间

拟合区间的选取与轨道弧段特征紧密相关。鉴此，本文根据卫星轨道运动特征，以卫星坐标二次差代数符号的变化(拐点)为特征点进行区间的划分。首先，介绍拐点的确定方法:

卫星的X 坐标分量序列为X(i)，(i=1，2，…，n)，每次取出4 个连续历元的坐标序列X(i－1)，X(i)，X(i+1)，X(i+2)，(i=2，…，n－2)，按式(3)求二次差:

判断式(3)所得两个连续二次差结果k(i－1)与k(i)的代数符号，若二者符号相反，则第i 个历元即为拐点;否则，i=i +1 继续判断。卫星Y、Z 坐标分量拐点的求解方法与X 坐标分量相同。图1 显示了PRN2 卫星轨道拐点求解的结果。

图1 PRN2 卫星轨道拐点Fig.1 Inflection points of PRN2 orbit

事实上，上述二次差结果近似表征了坐标分量加速度，表明卫星坐标分量在相邻两拐点之内的区间具有同向加速度，而该区间在数学图形上表现为上凸、上凹两种特征(图1)。因此，本文将相邻两拐点确定的连续区间称之为凹凸区间。考虑到“Runger 现象”主要发生在拟合区间两端［8］，因此在凹凸区间两端各增加一个节点构成新的区间，本文称之为拟合区间。多项式拟合在上述整个拟合区间内进行，加密星历的内插计算则在凹凸区间内进行。

分析图1 可知，凹凸区间的卫星轨道凹向或凸向同一侧，区间内节点具有一致的轨道运动特征(加速度同向)，因此有利于多项式建模描述。此外，凹凸区间的节点数完全由卫星坐标序列拐点确定(不同凹凸区间节点数不尽相同，X、Y 方向节点数为14 ～23，Z 方向为24 ～27)，表明凹凸区间不同于固定节点数的滑动区间，能够较好地顾及卫星轨道运动特征。需要说明的是，本文所得的凹凸区间一般不能完全覆盖单天所有历元节点，在单天开始与结束历元附近需要搭接相邻的精密星历产品。

3.2 加权拟合

IGS 事后精密星历产品不仅给出了轨道坐标，也给出了相应坐标的中误差［14］。若以加权最小二乘融合上述中误差信息σi，则所得结果将更加合理。因此，利用IGS 产品中的坐标中误差构造权阵P，进而获得加权Chebyshev 多项式拟合系数C 如下:

式中，

其中，σXi表示X 坐标分量中误差。卫星Y、Z 坐标分量的求解方法与X 坐标分量相同。

3.3 最佳阶数

按照上述方法对15 分钟间隔的卫星轨道数据进行凹凸区间划分，拟合区间节点数记为m，各区间以8 ～m－1 阶进行逐阶加权拟合，将拟合中误差小于0.001 m 时的最小阶数作为相应区间的最佳拟合阶数。图2 显示了32 颗卫星轨道区间节点数与最佳拟合阶数的统计关系，表1 为PRN2 卫星一个区间的计算结果。

为保证各区间拟合结果都能满足上述精度要求，故取图2各处最大值，得出最佳拟合阶数与区间节点总数的经验关系式(阶数曲线如图2 所示):

图2 拟合阶数与节点数关系Fig.2 Relation between fitting order and node number

表1 PRN2 拟合误差统计结果Tab.1 Statistical results of PRN2 fitting errors

式(5)中，q 为拟合阶数，m 为拟合区间的节点总数，int表示按四舍五入原则取整。

4 应用结果与分析

4.1 外符合精度分析

为验证本文方法的正确性和可靠性，以文献［11］提供的2002年1月1日的15 分钟精密星历ECF_15MIN.200(其中PRN12、16、19、24、32 无数据)内插为5 分钟间隔的低频数据，并与5 分钟间隔的精密星历ECF_5MIN.200 比较，以内插误差和中误差为指标进行外符合精度分析。其中，PRN2 卫星各分量误差结果如图3(a)所示，所有卫星的点位中误差如图3(b)所示。

图3(a)结果显示，PRN2 号卫星X、Y、Z 方向内插结果的最大误差均不超过3 mm，且基本服从正态分布，各方向的中误差分别为0.6、0.5 和0.7 mm。其他卫星的内插误差分布序列与图3(a)类似，不再展示。图3(b)统计了所有卫星内插结果的点位中误差，其最大值不超过1.2 mm，平均值为0.9 mm。由此，说明本文方法所得卫星位置精度可达到1 mm水平，验证了外符合精度的可靠性。

图3 外符合精度统计结果Fig.3 Statistical results of external accuracy

4.2 内符合精度分析

为进一步验证所提方法的优越性和可靠性，采用不同的Chebyshev 拟合法将ECF_15MIN.200 和ECF_5MIN.200 分别内插为1 s 采样率的数据，并以两者对应历元坐标的最大互差和中误差为指标进行差异分析。具体设计方案如下:

方案一:采用常规固定长度的Chebyshev 拟合法［1］，每次取3 小时长度的数据段，以10 阶内插得到1s 采样率的数据序列。

方案二:采用常规8 阶滑动式Chebyshev 拟合法［2］，其余处理策略同方案一。

方案三:采用本文方法，其余处理策略同方案一。

上述三个方案中，PRN2 卫星坐标分量误差的统计结果见表2，所有卫星坐标点位误差的统计结果如图4 所示。

表2 PRN2 坐标分量误差统计结果比较Tab.2 Results comparison of PRN2 coordinate component errors

图4 卫星坐标点位误差统计结果比较Fig.4 Results comparison of satellite position errors

由表2 可知，以X 坐标分量为例，方案一、二、三的最大互差依次为0.031 2、0.004 7、0.000 7 m，呈现一种明显的递减规律，表明所提出方法的最大互差均较其他方法的小;以中误差为指标，方案三比方案一的精度提高7 倍，比方案二的精度提高3 倍。从而表明，以X 坐标分量进行高频内插时，本文方法较其他方法更具稳定性。其他坐标分量也呈现类似的统计特性。

图4(a)中径向代表最大互差，切向代表卫星的PRN 序列;图4(b)中径向代表中误差，切向代表卫星的PRN 序列。图4 统计结果显示，三个方案的最大点位互差及其平均中误差依次为:方案一分别为5.1 cm 和2.3 mm，方案二分别为8 mm 和1.5 mm，方案三的分别为3 mm 和0.5 mm。对比结果说明，对于不同星历产品的轨道标准化误差，固定长度拟合方法存在很大差异，且不同卫星序列间差异显著，而滑动式拟合算法有所改善，但它们均不及本文所提的方法。由此表明，利用本文方法进行轨道标准化时，所得拟合轨道形态具有较好的一致性，较其他方法具有更高的精度和稳定性。

为说明本文方法标准化效率的优越性，比较分析上述三个方案程序计算所消耗的时间。程序基于MATLAB2009a 平台，在主频为2.93 GHz，内存为3.5 G 的计算机中运行，所得CPU 时间如表3 所示。

表3 三个方案程序执行的CPU 时间Tab.3 Execution CPU time of program with three methods

表3 结果显示，三个方案中，方案三计算用时最短，方案一次之，方案二最长。其主要原因是本方法大大延拓了内插窗口长度，减少了多项式拟合系数的求解次数。尽管本文方法需预先判断拐点，但总体上并未使结果计算时间延长。对比结果表明，本文所提方法较现有标准化方法效率更高。

5 结论

顾及地固坐标系下GPS 卫星轨道运动特征，并利用精密星历产品的坐标中误差信息，构建了一种基于轨道凹凸性的快速加权标准化方法。所提的加权准则与最佳阶数经验公式使得整个凹凸区间内的标准化结果更准确可靠，有效地提高了轨道标准化效率。大量计算结果表明，该方法的内外符合精度分别达到0.5 mm、1 mm，较Chebyshev 拟合方法更加稳定可靠。应指出的是，本文的凹凸区间划分等改进研究可直接应用于其他拟合类型的轨道标准化方法。

1 冯炜，等.两种常用GPS 星历拟合方法的精度分析［J］.大地测量与地球动力学，2010，(1):145－149.(Feng Wei，et al.Accuracy assessment of two fitting methods for GPS precise ephemeris［J］.Journal of Geodesy and Geodynamics，2010，30(1):145－149)

2 洪樱，欧吉坤，彭碧波.GPS 卫星精密星历和钟差三种内插方法的比较［J］.武汉大学学报(信息科学版)，2006，31(6):516－518，556.(Hong Ying，Ou Jikun and Peng Bibo.Three interpolation methods for precise ephemeris and clock offset of GPS satellite［J］.Geomatics and Information Science of Wuhan University，2006，31(6):516－518，556)

3 张守建，等.两种IGS 精密星历插值方法的比较分析［J］.大地测量与地球动力学，2007，(2):80－83.(Zhang Shoujian，et al.Comparative analysis on two methods for IGS precise ephemeris interpolation［J］.Journal of Geodesy and Geodynamics，2007，(2):80－83)

4 刘伟平，郝金明.一种新的IGS 精密星历插值算法［J］.武汉大学学报(信息科学版)，2011，(11):1320－1323.(Liu Weiping and Hao Jinming.A new interpolation method for IGS precise ephemeris［J］.Geomatics and Information Science of Wuhan University，2011，(11):1320－1323)

5 李明峰，江国焰，张凯.IGS 精密星历内插与拟合法精度的比较［J］.大地测量与地球动力学，2008，(2):77－80.(Li Mingfeng，Jiang Guoyan and Zhang Kai.Comparison between the accuracies with interpolating and fitting precise ephemeris［J］.Journal of Geodesy and Geodynamics，2008，(2):77－80)

6 杨学锋，等.利用切比雪夫多项式拟合卫星轨道坐标的研究［J］.测绘通报，2008，(12):1－3.(Yang Xuefeng，et al.On fitting satellite orbit coordinate by using Chebyshev polynomial［J］.Bulletin of Surveying and Mapping，2008，(12):1－3)

7 Feng Y and Zheng Y.Efficient interpolations to GPS orbits for precise wide area applications［J］.GPS Solutions，2005，9(4):273－282.

8 吉长东，徐爱功，冯磊.GPS 精密星历拟合与插值中龙格现象的处理方法［J］.测绘科学，2011，36(6):169－171.(Ji Changdong，Xu Aigong and Feng Lei.Treatment of Runge’s phenomenon in fitting ＆ interpolating GPS precise ephemeris［J］.Science of Surveying and Mapping，2011，36(6):169－171)

9 万亚豪，张书毕，侯东阳.GPS 精密星历插值法与拟合法的精度分析［J］.全球定位系统，2011，36(2):40－44.(Wan Yahao，Zhang Shubi and Hou Dongyang.Accuracy analysis on the interpolation and fitting of GPS precise ephemeris［J］.GNSS World of China，2011，36(2):40－44)

10 常亮，何秀凤.基于移动区间的GPS 轨道标准化方法［J］.大地测量与地球动力学，2009，(1):110－113.(Chang Liang and He Xiufeng.Standardization of GPS satellite orbits based on moving intervals［J］.Journal of Geodesy and Geodynamics，2009，(1):110－113)

11 Schenewerk M.A brief review of basic GPS orbit interpolation strategies［J］.GPS Solutions，2003，6(4):265－267.

12 高伟，姜水生.分段曲线拟合与离散度加权的数据误差处理方法［J］.中国测试技术，2005，31(6):55－56.(Gao Wei and Jiang Shuisheng.Method of deal with data error by subsection curve fitting and discrete degree weight［J］.China Measurement Technology，2005，31(6):55－56)

13 孔巧丽.用切贝雪夫多项式拟合GPS 卫星精密坐标［J］.测绘通报，2006(8):1－3.(Kong Qiaoli.Using Chebyshev polynomial to fit the precise satellite ephemeris［J］.Bulletin of Surveying and Mapping，2006，(8):1－3)

14 The extended standard product 3 orbit format(SP3-c)［EB/OL］.http://www.igscb.jpl.nasa.gov/igscb/data/format/sp3c.txt.