基于基本多项式的GPS精密星历插值方法研究

张 志，廖 瑛，文援兰，杨雅君，李 志,2

(1. 国防科技大学 航天与材料工程学院，湖南 长沙 410073； 2. 中国空间技术研究院，北京 100094)

一、引 言

在GPS精密定位中，尤其是高精度、动态定位中，经常要用到IGS[1](International GNSS Service)提供的精密星历，其中三维坐标的最终精度优于2.5 cm，但GPS精密星历的数据采样间隔为15 min。而GPS接收机的采样率一般为30 s或15 s，甚至更密，因此，要想利用某一时刻的卫星位置，就必须对精密星历进行高精度的插值。在GPS数据后处理中，用户需要选择合理的数学模型对GPS精密星历进行插值或拟合，才能获得任意时刻的卫星坐标，因而研究如何利用精密星历快速地插值出观测瞬间卫星的位置成为实际应用中一项十分必要的工作。

对GPS精密星历的内插法主要有多项式内插法和三角内插法，采用多项式内插法和三角内插法得到的精度基本相当[2]。对于多项式内插法，要想获得理想的精度就要合理地选择多项式的阶数和插值的时间间隔。目前，国内用于精密星历插值较常用的方法有Lagrange插值法、Neville插值法[3-4]及Chebyshev多项式插值法[5-6]等，但对基本多项式插值法的研究很少。Mark Schenewerk[2]和Milan Horemuž[7]均提到了基本多项式插值法用于精密星历的插值计算，前者只给出了概括性的分析说明，后者虽对多项式插值法在精密星历的插值过程中出现的问题进行了分析，但相关文献并未给出具体的插值计算方法。

本文利用Householder变换求解基本多项式拟合系数来避免对病态的法方程系数矩阵直接求逆计算二乘解，并对拟合曲线进行插值计算，为避免高阶基本多项式插值过程中出现的龙格现象的影响，给出了多项式拟合区间长度、有效区间长度和基本插值多项式阶数的选取准则，最后对插值的结果进行比较分析。

二、插值多项式选取原则

对于插值多项式(插值函数)来说，为了使插值函数的数值计算具有可操作性，通常取插值函数类为有限维子空间。为使插值便于实现，插值函数类子空间的一组基函数应当容易计算。一旦确定了插值函数类子空间的一组基函数，插值问题就可以转化为求一组组合系数的问题。

由定理说明组合系数的唯一性：设φi,i=1,2,…，n+1是线性无关的多项式，那么任一多项式p可唯一地由φ1t,φ2t,…,φn+1t的线性组合来表示[8]。

设线性无关的一组基函数1,t,t2,…,tn，则插值多项式(插值函数)pt可唯一地由基函数线性组合表示pt=c1tn+c2tn-1+…+cnt+cn+1，此插值多项式称为基本多项式。同样，Chebyshev插值多项式和Legendre插值多项式等均可以唯一地由一组线性无关的基函数表示。以下将以基本多项式为例进行分析讨论。

三、基本多项式插值方法

给定m个历元时刻t1,t2,…,tm的观测点坐标，假设需要计算n阶基本多项式系数。为了改进数值计算精度，需要将变量t正则化为

式中，t∈t1,t2,…,tm；mean()、std()分别表示向量的均值和标准差。

则GPS卫星坐标(x,y,z)的基本多项式分别为

(1)

式中，n为基本多项式的阶数；cxi、cyi、czi分别为(x,y,z)坐标分量的基本多项式系数，为待求量。

误差方程的矩阵形式表示为

(2)

为最小。此处给出最小二乘问题的法方程组

(3)

式中，R1是为n+1阶上三角矩阵；R中后面m-(n+1)行元素为0。

利用正交矩阵不改变向量长度的性质，又有

这样用QR分解的最小二乘法求解出基本多项式系数cxi，同理，可求得cyi、czi，利用这些系数，根据式(1)便可计算出t∈t1,tm时间区间内任意时刻的卫星坐标。

四、实例分析

为更好地进行试验比对，采用Mark Sche-newerk[2]中的数据，包括时间间隔(步长)为15 min的精密星历ECF_15MI.200，文中称为试验数据；作为插值比较的5 min间隔数据文件ECF_5MIN.200，文中称为真值数据，以上数据均可以从http:∥www.ngs.noaa.gov/gps-toolbox 下载。以试验数据的第一个历元点2002年1月1日00∶00∶00.000作为拟合起始历元，分别采用2 h、3 h、4 h、5 h弧段进行基本多项式拟合，然后对拟合的曲线进行5 min间隔的插值，得到的插值数据与对应历元真值数据之差称为插值误差，以此来说明方法的有效性。插值误差如图1—图4所示。本文给出了PRN01星的比较结果，其他各星结果类似。

图1 2 h拟合区间，8阶基本多项式插值误差

图2 3 h拟合区间，12阶基本多项式插值误差

图3 4 h拟合区间，16阶基本多项式插值误差

图4 5 h拟合区间，20阶基本多项式插值误差

由图1—图4可以看出，高阶基本多项式插值结果并不稳定，插值结果存在明显的龙格现象[11]，即插值的两边缘部分的精度较低，且存在一定的震荡性。但图2—图4清晰地表明，基本多项式拟合区间的中间2 h区域插值精度为毫米级，因此可以称拟合区间中间2 h区域为有效区间，这样，只要每次插值计算都进入有效区间，则能有效避免龙格现象的影响，从而使拟合达到毫米级的精度。对2 h的有效区间，最大插值误差不超过5 mm，这对于厘米级的精密星历来说，是完全可以接受的。

综合考虑计算量和拟合阶数，本文选取阶数为16，拟合区间为4 h的多项式，在2 h有效区间进行插值计算，如图5所示，对试验数据进行4 h拟合，在一天内有效区间插值数据与对应历元真值数据的差值如图6所示。

由图6可以看出，试验数据在时间段00∶00∶00.000—23∶45∶00.000内的插值计算结果表明，有效区间插值数据只能得到1～21 h的数据结果，若希望得到试验数据中24 h的有效插值结果，则需要提供前一天和后一天的精密星历数据。试验结果显示，插值数据精度可达到毫米级，误差小于3 mm。由于精密单点定位对卫星星历误差的要求为毫米级，因此这并不影响多项式拟合法插值出的数据在实际中的应用。用两个统计量Max和RMS来反映试验数据所有卫星基本多项式插值结果，表1列出了PRN01、PRN15和PRN31号卫星的分析结果。其中，Max指插值结果三维位置误差绝对值的最大值；RMS指内插结果三维位置误差的均方根误差。

由表1可知，采用2 h的有效区间基本多项式插值方法，试验数据3颗卫星x、y、z3个方向的插值误差最大值均不超过3 mm，3个方向的均方根小于0.6 mm。对于厘米级的精密星历来说，试验结果是比较理想的。

图5 拟合区间和有效区间选取准则

图6 一天有效区间插值结果与真值数据差值

MaxRMSxyzxyzPRN012.55102.15321.98610.48390.48880.4776PRN151.61571.84242.36110.43260.44210.5290PRN311.73111.73941.86420.41360.44300.4534

五、结束语

本文对插值多项式的选取原则进行了探讨，根据选取原则设计了一种简单的基本插值多项式，然后利用基本多项式拟合得到插值多项式系数，并对精密星历进行插值测试。试验结果表明，在插值边缘部分出现明显的龙格现象， 用基本多项式插值精密星历时，4 h拟合区间、2 h有效区间、多项式拟合阶数选取16能够有效消除龙格现象的影响。根据以上选取准则，对试验数据进行一天的插值，给出有效区间之内的插值结果，插值误差小于3 mm，能够满足用户定位精度要求。

参考文献：

[1] GERHARD B, MOORE A W, MUELLER I I. The International Global Navigation Satellite Systems Service (IGS):Development and Achievements[J]. Journal of Geodesy，2009, 83(3-4): 297-307.

[2] SCHENEWERK M. A Brief Review of Basic GPS Orbit Interpolation Strategies[J]. GPS Solutions，2003, 6(4): 265-267.

[3] 刘迎春，林宝军，张晓坤. 一种卫星精密星历的插值方法[J]. 飞行器测控学报，2004, 23(4): 43-46.

[4] 王晓明，成英燕，刘立. 基于阶次组合的GPS精密星历插值研究[J]. 大地测量与地球动力学，2011, 31(4): 103-106.

[5] 彭泽泉. GPS精密星历拟合方法的研究[J]. 测绘科学，2010, 35(S1): 63-65.

[6] 孔巧丽. 用切贝雪夫多项式拟合GPS卫星精密坐标[J]. 测绘通报，2006(8): 1-3.

[7] HOREMUŽ M,ANDERSSON J V. Polynomial Interpolation of GPS Satellite Coordinates[J]. GPS Solutions，2006, 10(1): 67-72.

[8] 关治，陆金甫. 数值方法[M]. 北京: 清华大学出版社, 2006： 192-193.

[9] CRASSIDIS J L,JUNKINS J L. Optimal Estimation of Dynamic Systems[M].Boca Raton,Florida: CRC Press, 2004：7-13.

[10] 金凯德，切尼，王国荣. 数值分析[M].3版.北京: 机械工业出版社, 2005：221-224.

[11] 吉长东，徐爱功，冯磊. GPS精密星历拟合与插值中龙格现象的处理方法[J]. 测绘科学，2011, 36(6): 169-171.