主线:Neville 递推原理 →
interppol()逐行对应 →pephpos()的 11 节点 SP3 插值 → 与拉格朗日比较 → 阶数、误差及工程注意事项。本文以 RTKLIB 官方仓库
master/src/preceph.c的所核对实现为准。
1. 澄清三个概念
Neville(内维尔)算法将相邻数据点的低次插值结果逐层合并,直接求出指定时刻的插值值,不显式计算多项式系数。RTKLIB 的interppol()使用的正是 Neville 递推,而不是 Newton(牛顿)差商插值。
| 名称 | 核心思路 | 与 RTKLIB 的关系 |
|---|---|---|
| Neville 插值 | 子节点组逐层递推,在目标时刻求值 | interppol()的实际算法 |
| Lagrange(拉格朗日)插值 | 构造基函数并加权求和 | 与 Neville 表示同一个插值多项式 |
| Newton(牛顿)差商插值 | 先求差商系数,再用嵌套形式求值 | 数学结果相同,不是该函数的实现 |
对同一组互异节点,三种方法在精确算术下得到同一个唯一插值多项式;换一种求值算法,并不能自动降低该多项式的截断误差。
RTKLIB 在preceph.c中定义NMAX=10:10 次多项式需要 11 个节点,并不意味着只读取 10 个 SP3 历元。
2. 数学原理:从一般公式推到 RTKLIB 原地递推
2.1 两节点线性插值
设节点为( x i , y i ) (x_i,y_i)(xi,yi),其中x i x_ixi是相对目标时刻的时间(s),y i y_iyi是一个待插值的卫星坐标分量(m);节点时间两两互异。目标自变量记为u uu(s)。利用相邻两个节点构造一次多项式:
P i , 1 ( u ) = y i x i + 1 − u x i + 1 − x i + y i + 1 u − x i x i + 1 − x i = ( x i + 1 − u ) y i − ( x i − u ) y i + 1 x i + 1 − x i . \begin{aligned} P_{i,1}(u) &=y_i\frac{x_{i+1}-u}{x_{i+1}-x_i} +y_{i+1}\frac{u-x_i}{x_{i+1}-x_i}\\ &=\boxed{\frac{(x_{i+1}-u)y_i-(x_i-u)y_{i+1}}{x_{i+1}-x_i}}. \end{aligned}Pi,1(u)=yixi+1−xixi+1−u+yi+1xi+1−xiu−xi=xi+1−xi(xi+1−u)yi−(xi−u)yi+1.
其中i ii为起始节点下标,P i , 1 P_{i,1}Pi,1表示由节点i , i + 1 i,i+1i,i+1构造的一次插值多项式。
2.2 从低次到高次的 Neville 递推
记P i , j ( u ) P_{i,j}(u)Pi,j(u)为使用节点i , i + 1 , … , i + j i,i+1,\ldots,i+ji,i+1,…,i+j构造的最高j jj次插值多项式在u uu处的值。起始条件为
P i , 0 ( u ) = y i . \boxed{P_{i,0}(u)=y_i}.Pi,0(u)=yi.
已有两个相互重叠的( j − 1 ) (j-1)(j−1)次插值多项式:
- P i , j − 1 ( u ) P_{i,j-1}(u)Pi,j−1(u):由节点i , … , i + j − 1 i,\ldots,i+j-1i,…,i+j−1构造;
- P i + 1 , j − 1 ( u ) P_{i+1,j-1}(u)Pi+1,j−1(u):由节点i + 1 , … , i + j i+1,\ldots,i+ji+1,…,i+j构造。
按首尾节点构造新函数:
Q ( u ) = ( x i + j − u ) P i , j − 1 ( u ) + ( u − x i ) P i + 1 , j − 1 ( u ) x i + j − x i . Q(u)=\frac{(x_{i+j}-u)P_{i,j-1}(u)+(u-x_i)P_{i+1,j-1}(u)}{x_{i+j}-x_i}.Q(u)=xi+j−xi(xi+j−u)Pi,j−1(u)+(u−xi)Pi+1,j−1(u).
在u = x i u=x_iu=xi、u = x i + j u=x_{i+j}u=xi+j和所有中间节点处,Q QQ都等于对应的观测值;同时Q QQ的次数不超过j jj。根据互异节点的插值多项式唯一性,Q = P i , j Q=P_{i,j}Q=Pi,j,得到Neville 通式:
P i , j ( u ) = ( x i + j − u ) P i , j − 1 ( u ) − ( x i − u ) P i + 1 , j − 1 ( u ) x i + j − x i , j ≥ 1. (1) \boxed{ P_{i,j}(u)=\frac{(x_{i+j}-u)P_{i,j-1}(u)-(x_i-u)P_{i+1,j-1}(u)}{x_{i+j}-x_i} },\quad j\ge 1. \tag{1}P