由卫星定位基本原理可知，确定卫星的轨道位置是一项必要的基础工作。卫星位置一般通过卫星星历获取。精密星历提供的卫星三维坐标精度较高，被广泛应用于大地测量、地球动力学等精密领域^[1]。但是精密星历是离散数据，它一般提供15 min或5 min时间间隔的卫星轨道数据和精密钟差信息^[2]，而在卫星定位数据处理中常需要时间间隔为秒级的卫星位置数据。因此需要利用插值拟合数学模型，内插得到任一时刻的卫星三维坐标，同时要保证较高的插值精度和效率。

目前，对于全球定位系统(Global Positioning System, GPS)精密星历的内插研究较多，所用方法有拉格朗日(Lagrange)插值法、牛顿插值法、切比雪夫多项式拟合、最小二乘曲线拟合等等。文[3]利用滑动式拉格朗日多项式对15 min的全球定位系统精密星历进行了内插研究；文[4]利用拉格朗日插值法、切比雪夫多项式拟合和内维尔算法对全球定位系统精密星历内插，并比较了3种算法的优劣；文[5]利用切比雪夫多项式拟合法对全球定位系统精密星历的最佳拟合阶数进行了研究；文[6]利用广义延拓插值法对全球定位系统精密星历进行内插与外推研究；文[7]利用滑动式切比雪夫多项式对5 min间隔的全球定位系统精密星历进行了拟合研究；文[8]利用傅里叶级数对全球定位系统精密星历进行了内插研究；文[9]提出了一种顾及卫星运动特点的全球定位系统精密星历插值新方法等等。但是以上研究只是针对全球定位系统精密星历，对我国的北斗导航系统精密星历内插研究相对很少。本文分别利用非滑动式和滑动式拉格朗日多项式插值法，对5 min时间间隔的北斗导航系统精密星历进行内插精度对比分析。

1 非滑动式与滑动式拉格朗日插值法 1.1 拉格朗日插值法

对于已知节点x_i(i=0, 1, …, n)和对应函数值y_i(i=0, 1, …, n)，构造一个n次插值多项式：

$ {L_n}\left( x \right) = \sum\limits_{k = 0}^n {{y_k}{l_k}\left( x \right)}, $

(1)

其中，l_k(x)(k=0, 1, …, n)为基本插值多项式，使它满足在该点上取值为1，其它点上取值为0，即

$ {l_k}({x_i}) = \left\{ \begin{array}{l} 1, \;\;\;i = k\\ 0, \;\;\;i \ne k \end{array} \right., $

(2)

由(2)式可见，n个节点x_i(i≠k)都是n次多项式l_k(x)的零点，因此：

$ {l_k}({x_i}) = {A_k}\prod\limits_{\begin{array}{*{20}{c}} {i = 0}\\ {i \ne k} \end{array}}^n {(x-{x_i})}, $

(3)

由条件l_k(x_k)=1可得

$ {l_k}\left( x \right) = \prod\limits_{\begin{array}{*{20}{c}} {i = 0}\\ {i \ne k} \end{array}}^n {\frac{{(x-{x_i})}}{{({x_k}-{x_i})}}}, \;k = 0, \;\;1, \ldots {\rm{ }}n $

(4)

将(4)式代入(1)式，得它的n阶拉格朗日插值多项式表达形式如下^[10]：

$ {L_n}\left( x \right) = \sum\limits_{k = 0}^n {{y_k}} \frac{{(x-{x_0}) \ldots (x-{x_{k-1}})(x - {x_{k + 1}}) \ldots (x - {x_n})}}{{({x_k} - {x_0}) \ldots ({x_k} - {x_{k - 1}})({x_k} - {x_{k + 1}}) \ldots ({x_k} - {x_n})}}. $

(5)

当利用拉格朗日插值法对北斗导航系统卫星精密星历进行插值时，历元就成为自变量，卫星位置为因变量，具体表达形式如下：

$ \left\{ \begin{array}{l} {L_X}\left( t \right) = \sum\limits_{k = 0}^n {{X_k}} \frac{{(t-{t_0}) \ldots (t-{t_{k-1}})(t - {t_{k + 1}}) \ldots (t - {t_n})}}{{({t_k} - {t_0}) \ldots ({t_k} - {t_{k - 1}})({t_k} - {t_{k + 1}}) \ldots ({t_k} - {t_n})}}\\ {L_Y}\left( t \right) = \sum\limits_{k = 0}^n {{Y_k}} \frac{{(t - {t_0}) \ldots (t - {t_{k - 1}})(t - {t_{k + 1}}) \ldots (t - {t_n})}}{{({t_k} - {t_0}) \ldots ({t_k} - {t_{k - }}_1)({t_k} - {t_{k + 1}}) \ldots ({t_k} - {t_n})}}\\ {L_Z}\left( t \right) = \sum\limits_{k = 0}^n {{Z_k}} \frac{{(t - {t_0}) \ldots (t - {t_{k - 1}})(t - {t_{k + 1}}) \ldots (t - {t_n})}}{{({t_k} - {t_0}) \ldots ({t_k} - {t_{k - 1}})({t_k} - {t_{k + 1}}) \ldots ({t_k} - {t_n})}} \end{array} \right., $

(6)

其中，t_i(i=0, 1, 2, … n)为已知历元；t为待求历元；X_k, Y_k, Z_k分别为历元k时的已知卫星X, Y, Z坐标分量；L_X(t), L_Y(t), L_Z(t)分别为内插历元t时的卫星坐标分量值。

1.2 非滑动式与滑动式插值算法

非滑动式插值算法不限定待插值点的所处区间位置，待内插点只要处于插值区间内即可，故可以内插出较多数据，插值效率较高。为了提高插值精确度，常利用高阶插值多项式，但是此时处于插值区间两端的待插值点不稳定，易发生震荡或跳跃现象，这在数值计算上被称为“龙格”现象。

所谓滑动式插值法，就是使待内插点始终处于插值区间的中间，将插值区间作为一个固定“窗口”，“窗口”每向后移动一个单位，内插点也随之向后移动一个单位。例如当插值节点数为10，那么第1次插值时，插值节点为第1个到第10个，待插值点位于第5个和第6个节点之间；第2次插值时，插值节点为第2个到第11个，待插值点位于第6个和第7个节点之间，以此类推。该法能保持较高的插值精度，避免“龙格”现象的发生，但是相对非滑动式插值法，它单次的内插数据相对较少，插值效率相对较低。

2 北斗导航系统精密星历内插算例分析

从网址ftp://cddis.gsfc.nasa.gov/pub/gps/products/mgex/下载了2018年1月27日的混合精密星历文件“gbm19856”，文件包含了14颗健康的北斗导航系统卫星位置信息，这14颗北斗卫星分别为5颗地球同步轨道(Geostationary Earth Orbit, GEO)卫星、5颗倾斜地球同步轨道(Inclined Geosynchronous Orbit, IGSO)卫星和4颗中地球轨道(Medium Earth Orbit, MEO)卫星，采样时间间隔为5 min，即采样时间为0 min, 5 min, 10 min, ……, 1 435 min，故每颗卫星共含有288个历元下的位置数据。

把从10 min历元间隔(即0 min, 10 min, 20 min, ……, 1 430 min，共144个)下的北斗导航系统卫星坐标构成的数组里，每次插值选取部分作为已知节点(如3阶插值，则每次取4个数据点)，分别利用非滑动式和滑动式拉格朗日多项式插值法，对北斗导航系统精密星历中的这14颗北斗卫星进行内插。对于非滑动式插值，其内插点位于插值区段内；对于滑动式插值，其内插点位于插值区段中间，比如取3阶插值，则第1次插值节点为0 min，10 min，20 min，30 min时的卫星坐标，那么非滑动式法的内插点为5 min，15 min，25 min时的卫星坐标，而滑动式插值的内插点为15 min时的卫星坐标。把精密星历提供的10 min历元间隔(即5 min，15 min，25 min，……，1 425 min，共143个)下的已知卫星坐标作为“真值”，将内插出来的卫星坐标与同一历元和卫星下的“真值”对应求差，对残差统计，求出均方根差，以分析其插值精度，求出其绝对值的最大值，分析插值准确度。图 1为插值设计流程图。

本文主要从3方面对北斗导航系统精密星历内插精度进行分析，分别是非滑动式拉格朗日多项式插值法在不同插值阶数下的插值精确度，滑动式拉格朗日插值法在不同插值阶数下的插值精确度，以及非滑动式和滑动式拉格朗日插值法在同一插值阶数下的插值精确度。

2.1 不同插值阶数下的非滑动式拉格朗日插值精度分析

表 1为非滑动式拉格朗日插值法在不同插值阶数下的误差统计结果。从表 1可以看出，插值阶数与插值精度存在一定的规律性：(1)当插值阶数取值较小时，其插值误差较大，如取4阶插值时，其三维坐标误差最大值达到了米级，5阶插值误差最大值处于分米级，这远不能满足插值精度的要求。这是由于利用较少已知数据建立插值模型，难以完全表达研究对象的特性；(2)在一定范围内，插值精确度随插值阶数的增加而提高，即它们为正相关关系，如插值阶数从4阶增加到8阶时，其插值误差最大值从米级持续降低到毫米级，即插值精度提高了3个数量级；(3)当插值阶数取值较高时，插值精度不升反降，即出现了“龙格”现象，如插值阶数从8阶增加到18阶时，其插值误差最大值从毫米级一直增加到分米级，即插值精度降低了2个数量级，这是由于插值多项式经过太多含有观测误差的插值节点，导致误差累积过多，插值模型精度降低。

2.2 不同插值阶数下的滑动式拉格朗日插值精度分析

表 2为滑动式拉格朗日插值法在不同插值阶数下的误差统计结果。从表 2可以看出，类似于非滑动式拉格朗日插值法，滑动式拉格朗日插值法的插值阶数较小时，同样其插值精度较差，如4阶插值下的X，Y，Z坐标分量误差最大值分别达到米、米、分米级。当插值阶数从5阶增加到6阶时，插值精度迅速提高，三维坐标误差最大值进入毫米级。插值阶数从7阶开始，插值精度处于一个稳定状态，不再随插值阶数的增加而变化，X，Y，Z坐标分量误差最大值稳定在1.1 mm左右，误差的均方根差稳定在0.4 mm。这是由于处于插值区间中间位置的待插值点稳定，即使插值阶数取值较高时，也不会出现震荡或跳跃问题。这也验证了滑动式插值算法的高度稳定性和准确性。

2.3 同一插值阶数下的非滑动式与滑动式拉格朗日插值精度对比分析

本小节主要从两个角度对同一插值阶数下的非滑动式与滑动式拉格朗日插值精度进行对比分析：(1)对3类不同轨道卫星的独立分析，即把3类不同轨道类型的地球同步轨道卫星、倾斜地球同步轨道卫星、中地球轨道卫星作为3个独立的分析对象，求出它们的插值误差最大值；(2)对3类不同轨道卫星的整体分析，即把3类不同轨道卫星看为一个整体，求出它们的插值误差最大值和均方差。通过这两个角度的精度分析，对比得出非滑动式与滑动式拉格朗日插值法的优劣。

2.3.1 对3类不同轨道卫星的独立分析

表 3、表 4、表 5分别给出了地球同步轨道、倾斜地球同步轨道、中地球轨道卫星在4~18阶下的非滑动式与滑动式拉格朗日插值误差最大值。从表 3~5可以看出，对于3类不同轨道类型的北斗卫星来说，当取相同的插值阶数时，滑动式拉格朗日插值法的插值误差最大值总是小于非滑动式拉格朗日插值法的插值误差最大值，即滑动式拉格朗日插值法的插值精度总是更高些。而且随着插值阶数的增大，非滑动式拉格朗日插值法的“龙格”问题越来越突出，而滑动式拉格朗日插值法则保持较高的插值精度和稳定性。

对于地球同步轨道卫星来说，在取4~10阶插值时，非滑动式与滑动式拉格朗日插值法的精度差异小些，插值误差最大值处于同一量级——毫米级，都可以较好满足精密星历的精度需求。但当阶数高于10后，非滑动式拉格朗日插值法的插值精度下降的速度越来越快，到18阶时，其插值误差最大值已高达5 dm以上，而滑动式拉格朗日插值法的误差最大值仍保持在1 mm左右。

对于倾斜地球同步轨道卫星来说，在取6~10阶插值时，非滑动式与滑动式拉格朗日插值法的插值误差最大值处于毫米级，即此时精度差异较小。在取6~18阶多项式时，非滑动式拉格朗日插值法的插值误差最大值随着阶数的增加而增大，而滑动式拉格朗日插值法的插值误差最大值稳定在1.1 mm左右，故它们之间的精度差异与插值阶数在该阶数范围内属于正相关关系。

对于中地球轨道卫星来说，插值阶数当取低阶——4阶、5阶时，非滑动式与滑动式拉格朗日插值法的插值精度都差，达到分米级及以上，这远不能满足精度要求。在取7~9阶插值时，它们的插值精度差异处于较小的阶段。此后它们之间的插值精度差距越来越大，在15~18阶时，非滑动式拉格朗日插值法的插值误差最大值是滑动式拉格朗日插值法的数百倍。

2.3.2 对3类不同轨道卫星的整体分析

为方便对比分析非滑动式与滑动式拉格朗日插值法在同一插值阶数下的插值精度，图 2和图 3分别绘出了6~14阶下的非滑动式与滑动式拉格朗日插值误差最大值和均方根差。从图 2和图 3可以直观看出，在同一插值阶数下，无论插值阶数取值是低是高，滑动式拉格朗日插值法的插值精确度都要优于非滑动式拉格朗日插值法的插值精确度；在插值阶数取7~10时，它们的误差最大值和均方根差之间的差距略小，此后随着插值阶数的增加，它们的误差最大值和均方根差之间的差距呈现越来越大的趋势，滑动式拉格朗日插值法的插值精确度基本稳定不变，而非滑动式拉格朗日插值法的插值精确度越来越差。但是非滑动式拉格朗日插值法可以内插出更多数据，故插值效率略高。

3 结论

本文基于非滑动式和滑动式拉格朗日插值法，通过对北斗导航系统精密星历提供的离散北斗卫星坐标进行实例插值分析，可得到以下结论：

(1) 通过对非滑动式拉格朗日插值法在不同插值阶数下的插值精度研究可知，插值精度随插值阶数的增加先提高后降低，即一个开口向下的抛物线关系。当插值阶数取7~9时，其插值效果达到顶峰，误差最大值仅为毫米级，可以满足插值精度的需求。

(2) 通过对滑动式拉格朗日插值法在不同插值阶数下的插值精度研究可知，滑动式插值算法具有高度的稳定性和准确性，对北斗导航系统精密星历内插而言，当插值阶数取7阶时，已达到最佳插值精度，误差最大值稳定在1 mm左右。再增加插值阶数，插值精度保持稳定不变。

(3) 对于地球同步轨道、倾斜地球同步轨道、中地球轨道3类不同轨道类型卫星来说，在同一插值阶数下，滑动式拉格朗日插值法的精确度要高于非滑动式拉格朗日插值法的精确度，但是前者的插值效率略低于后者。由于目前的计算机计算性能卓越，故它们之间的效率差异可忽略不计。因此在实际应用中，建议使用滑动式拉格朗日插值法对北斗导航系统精密星历进行内插。