B.6 三次样条插值

线性插值就是两个样本点之间拉一条直线。问题在哪？走到样本点，曲线方向会突然折一下。位置是连着的，但走向变了。对温度曲线、机械轨迹、传感器趋势这类数据，这种折角往往太生硬。

先看一个四点的例子：

(0,0),\quad (1,1),\quad (2,0.5),\quad (3,0)

用线性插值， $(0,0)\to(1,1)$ 是一条直线， $(1,1)\to(2,0.5)$ 又换成另一条直线。曲线穿过了所有样本点，但在 $(1,1)$ 和 $(2,0.5)$ 这两个点上方向硬拐。

右图同样四个点，三次样条穿过每个点，接得平滑。想想电梯启动时，速度不会突然从 0 跳到额定值，先缓缓加速，再缓缓收住。现实里的温度、位移、传感器读数大多是这种慢慢变的量，不会在一个采样点上瞬间改变走势。三次样条插值就是把这种假设体现在插值规则里。

连接处要对齐什么

要让曲线在样本点处拐弯不生硬，先弄清楚连接处两段要对齐到什么程度。

线性插值只让两段在位置上连上，方向各走各的，所以会折角。如果在连接处再让两段的斜率相等，方向就接上了，折角消失。如果连弯曲程度也对齐，曲线连”怎么转弯”都一致，就是一条流畅的轨迹。

图 B.6-2 从左到右就是这三层。最左只对齐位置，折角仍在；中间让斜率也相等，方向变顺；最右连弯曲都对齐，曲线完全平滑。三次样条要做到的就是最右这一档。

不过中间和最右两条几乎一个样。只对齐斜率和再对齐弯曲，差别藏在节点附近，肉眼很难分辨。把节点 $t=2$ 附近放大、再看一眼弯曲程度本身，区别就出来了：

图 B.6-2附左边把两条曲线叠在一起，放大插图里它们在节点附近确实只差一点点。差别在右边，弯曲程度就是二阶导 $S''$ 。只对齐斜率的那条（蓝线），二阶导在 $t=2$ 这个节点上跳了一个台阶。曲线方向虽然连上了，但转弯的急缓突然换挡；而三次样条（绿线）的二阶导平滑穿过节点，不存在这种突变。完全平滑指的就是连弯曲程度都连续，这一点光看曲线形状不容易察觉，看二阶导才清楚。

要同时安排位置、斜率、弯曲，每一段曲线就得有足够的自由度。一段曲线能调几个量，取决于它有几个待定系数。直线 $a+bt$ 只有两个，定好两端高度就用光了，斜率和弯曲再没法单独指定。要在一段里同时管住高度、斜率、弯曲，至少得有第三个系数；而弯曲还需要能沿段变化，于是要第四个。三次多项式 $a+bt+ct^2+dt^3$ 恰好四个系数，够分这四件事。这就是用三次的原因，二次少一个自由度，做不到弯曲连续。

把这四个系数逐个拨动一下，就能看清它们各管曲线的哪一面：

图 B.6-2c 以一条平直段为基线（虚线），每个面板只动一个系数：调 $b$ 把整段抬成斜线，管的是起点斜率；调 $c$ 让曲线凹下或凸起，管的是弯曲方向；调 $d$ 则让弯曲沿着这一段慢慢转向。还有个没画的 $a$ ，它只是整体上下平移，定起点高度。四个系数各管一件事，合起来才够拼出一段能任意调姿态的曲线。

至于”样条”这个名字，来自老式制图，绘图员把一根有弹性的细木条（英文 spline）压在几个固定钉点上，木条自己弯成一条最省力的顺滑曲线，再照着描下来。

三次样条插值就是这根木条的数学版本，把整条曲线拆成一段一段，每段用一个三次多项式接力。设样本点为 $(t_i,y_i)$ ，第 $i$ 段写成

S_i(t)=a_i+b_i(t-t_i)+c_i(t-t_i)^2+d_i(t-t_i)^3

把 $t=t_i$ 代进 $S_i$ 、 $S_i'$ 、 $S_i''$ ，各系数的含义就出来了：

S_i(t_i)=a_i,\qquad S_i'(t_i)=b_i,\qquad S_i''(t_i)=2c_i

也就是 $a_i$ 是这一段起点的高度， $b_i$ 是起点的斜率， $c_i$ 对应起点的弯曲程度，剩下的 $d_i$ 控制弯曲沿着这一段如何变化。三次曲线的二阶导是一条直线，所以弯曲程度能从左端平滑过渡到右端，这正是直线和二次曲线做不到、而样条需要的。

真正把曲线缝顺的，是连接点处的三个条件。位置对上：

S_{i-1}(t_i)=S_i(t_i)=y_i

斜率对上：

S'_{i-1}(t_i)=S'_i(t_i)

弯曲对上：

S''_{i-1}(t_i)=S''_i(t_i)

如果一共有 $n$ 个小区间，每段四个系数，就是 $4n$ 个未知数。上面三组连续条件加上“曲线必须过样本点”，给出的方程比未知数少两个。差的两个由端点条件补上，比如自然边界让两端二阶导为 0，钳制边界指定两端斜率。

MATLAB 里 spline 是一行调用，背后就是解这组系数。

下面这段代码可以跑出图 B.6-1 的对比：

t = [0 1 2 3];                              % 样本点横坐标
y = [0 1 0.5 0];                            % 样本点纵坐标
tq = linspace(0, 3, 100);                   % 加密的查询位置

y_linear = interp1(t, y, tq, 'linear');     % 线性插值：两点间直线
y_spline = spline(t, y, tq);                % 三次样条：解系数后逐段求值

plot(tq, y_linear, '--', 'LineWidth', 1.4); hold on;
plot(tq, y_spline, 'LineWidth', 1.8);
plot(t, y, 'ko', 'MarkerFaceColor', 'k');   % 黑点标出原样本
legend('线性插值', '三次样条', '样本点');
grid on;

spline(t, y, tq) 这一行，就是把上面那组连接条件交给求解器，再在 tq 上把每段三次多项式求值出来。跑完就是图 B.6-1。

样条在测量曲线、轨迹规划、动画关键帧过渡里都常用，适合点稀又希望中间走得自然的场景，能把离散点串成一条顺滑的趋势。

动手算一个中间点

挑 $t=1.5$ 手算一次，看样条给出的数和线性插值差在哪里。

这套方法叫三弯矩法，给每个节点的二阶导各列一个方程联立求解（节点处的二阶导在力学里叫弯矩，方法因此得名）。

还用开头那四个样本点，算 $t=1.5$ 处的值，它落在 $[1,2]$ 段。线性插值好算， $t=1.5$ 正好在中点，取两端平均：

y_{\text{linear}}(1.5)=\frac{1+0.5}{2}=0.75

样条要分两步。

第一步，求各节点的二阶导。样条的未知量是各节点处的二阶导 $S''(t_0),\dots,S''(t_3)$ 。两端由自然边界直接给出： $S''(t_0)=S''(t_3)=0$ 。中间两个要解方程，方程本身由三个连接条件推出来。推导分三小步：

步骤①:写出一段内的 $S''$ 。

令 $u=t-t_i$ , $u$ 从 $0$ 跑到 $h$ (等间隔时 $h=1$ )。对第 $i$ 段逐项求导：

S_i(u)=a_i+b_i u+c_i u^2+d_i u^3 \quad\Longrightarrow\quad S_i''(u)=2c_i+6d_i u

这是一条关于 $u$ 的直线。它在 $u=0$ (节点 $t_i$ )处等于 $S''(t_i)$ ，在 $u=h$ (节点 $t_{i+1}$ )处等于 $S''(t_{i+1})$ 。

下图把这条直线画了出来，同时标出了段内任意一点到两端的距离 $u$ 和 $h-u$ :

权重正比于到对面端点的距离，归一化后：

S_i''(u)=S''(t_i)\cdot\frac{h-u}{h}+S''(t_{i+1})\cdot\frac{u}{h}

验证一下： $u=0$ 代入得 $S''(t_i)$ , $u=h$ 代入得 $S''(t_{i+1})$ ，两端对上了。

步骤②:积分一次，得到斜率。

对上式关于 $u$ 积分：

S_i'(u)=S''(t_i)\!\left(u-\frac{u^2}{2h}\right)+S''(t_{i+1})\cdot\frac{u^2}{2h}+C_1

再积分一次，得到曲线本身：

S_i(u)=S''(t_i)\!\left(\frac{u^2}{2}-\frac{u^3}{6h}\right)+S''(t_{i+1})\cdot\frac{u^3}{6h}+C_1 u+C_2

出现两个积分常数 $C_1,C_2$ ，由”两端必须过样本点”各定一个：

$S_i(0)=y_i$ ，代入 $u=0$ ，立刻得 $C_2=y_i$ 。
$S_i(h)=y_{i+1}$ ，代入 $u=h$ ，整理得

C_1=\frac{y_{i+1}-y_i}{h}-\frac{h}{6}\bigl[2S''(t_i)+S''(t_{i+1})\bigr]

$C_1$ 就是 $u=0$ 处的斜率，即这一段左端的斜率：

\boxed{S_i'(t_i)=\frac{y_{i+1}-y_i}{h}-\frac{h}{6}\bigl[2S''(t_i)+S''(t_{i+1})\bigr]}

同理，把 $u=h$ 代入斜率表达式，得右端斜率：

\boxed{S_i'(t_{i+1})=\frac{y_{i+1}-y_i}{h}+\frac{h}{6}\bigl[S''(t_i)+2S''(t_{i+1})\bigr]}

框里这两个式子，未知量只有 $S''(t_i)$ 和 $S''(t_{i+1})$ , $y$ 和 $h$ 都是已知的。

步骤③:令内节点处两段斜率相等。

节点 $t_i$ 是段 $i-1$ 的右端，也是段 $i$ 的左端。斜率必须接上：

S_{i-1}'(t_i)=S_i'(t_i)

左边套段 $i-1$ 的右端斜率（把步骤②右端公式里的下标 $i$ 换成 $i-1$ ）:

S_{i-1}'(t_i)=\frac{y_i-y_{i-1}}{h}+\frac{h}{6}\bigl[S''(t_{i-1})+2S''(t_i)\bigr]

右边套段 $i$ 的左端斜率：

S_i'(t_i)=\frac{y_{i+1}-y_i}{h}-\frac{h}{6}\bigl[2S''(t_i)+S''(t_{i+1})\bigr]

令二者相等，两边乘以 $6/h$ ，移项整理， $y$ 差分合并后得到：

S''(t_{i-1})+4\,S''(t_i)+S''(t_{i+1})=6\,(y_{i-1}-2y_i+y_{i+1})

右边只是样本值的二阶差分，直接代数即可。 $i=1$ 和 $i=2$ 各写一个方程，并用上 $S''(t_0)=S''(t_3)=0$ :

4\,S''(t_1)+S''(t_2)=6(0-2\cdot 1+0.5)=-9

S''(t_1)+4\,S''(t_2)=6(1-2\cdot 0.5+0)=0

这是一个简单的二元一次方程组，留给你动手解一下，由第二式得 $S''(t_1)=-4\,S''(t_2)$ ，代回第一式即可。解出来：

S''(t_1)=-2.4,\qquad S''(t_2)=0.6

$S''(t_1)<0$ 说明曲线在 $t=1$ 附近向下弯， $S''(t_2)>0$ 说明到 $t=2$ 附近弯曲方向已缓过来。这两个数决定 $[1,2]$ 段怎么弯。

第二步，由节点二阶导定出这一段的四个系数。把 $[1,2]$ 段写成

S(t)=a+b(t-1)+c(t-1)^2+d(t-1)^3

令 $u=t-1$ ，则 $S(u)=a+bu+cu^2+du^3$ ，逐项求导：

S'(u)=b+2cu+3du^2,\qquad S''(u)=2c+6du

先求 $a$ 。 $u=0$ 代入 $S$ ，得 $S(0)=a$ ，而曲线起点必须过 $y_1=1$ ，所以

a=y_1=1

再求 $c$ 。 $u=0$ 代入 $S''$ ，得 $S''(0)=2c$ ;这个值等于节点二阶导 $S''(t_1)=-2.4$ ，所以

2c=S''(t_1)=-2.4\quad\Longrightarrow\quad c=-1.2

再求 $d$ 。 $u=h=1$ 代入 $S''$ ，得 $S''(1)=2c+6d$ ，这等于 $S''(t_2)=0.6$ :

2(-1.2)+6d=0.6\quad\Longrightarrow\quad 6d=3\quad\Longrightarrow\quad d=0.5

也可以直接用步骤②的结论：二阶导在段内是直线，从 $S''(t_1)$ 线性变到 $S''(t_2)$ ，所以斜率就是 $\dfrac{S''(t_2)-S''(t_1)}{6h}=\dfrac{0.6-(-2.4)}{6}=0.5$ ，结果一致。

最后求 $b$ 。现在 $a,c,d$ 都有了， $b$ 由右端必须过 $y_2$ 定出。令 $u=h=1$ , $S(1)=y_2$ :

a+b\cdot 1+c\cdot 1^2+d\cdot 1^3=y_2

1+b+(-1.2)+0.5=0.5

b=0.5-1+1.2-0.5=0.2

四个系数齐了， $[1,2]$ 段的多项式就是一条确定的曲线：

S(t)=1+0.2(t-1)-1.2(t-1)^2+0.5(t-1)^3

顺手验一下两端： $t=1$ 代进去得 $1$ , $t=2$ 代进去得 $1+0.2-1.2+0.5=0.5$ ，正好是两个样本值，系数没解错。

现在代入 $t=1.5$ ，即 $(t-1)=0.5$ ，逐项算：

S(1.5)=\underbrace{1}_{a}+\underbrace{0.2(0.5)}_{0.1}-\underbrace{1.2(0.5)^2}_{0.3}+\underbrace{0.5(0.5)^3}_{0.0625}=0.8625

同一个中间点，线性插值给 $0.75$ ，自然三次样条给 $0.8625$ 。

橙色虚线是线性结果 $0.75$ ，绿色曲线是样条结果 $0.8625$ 。样条把前后几段的弯曲趋势算了进去，把这段曲线整体往上抬了一点，原因就在刚才解出的那组节点二阶导里。

另一种做法是把所有未知系数全部列出来，再把所有条件全部写成方程：四个样本点分三段，每段四个系数，共 12 个未知数；位置连续、斜率连续、弯曲连续加上端点条件，正好凑出 12 个方程，解出来就是所有系数。哪条方程对应哪个条件，清楚直接。但 12 个方程联立是一个 $12\times12$ 的线性方程组，手算几乎不现实；节点数一多，规模按 $4n$ 增长，很快就不好处理了。

三弯矩法的规模要小得多。系数 $a,c,d$ 都能从节点二阶导直接读出， $b$ 由右端条件补一个方程定死；需要联立求解的只有内节点处的二阶导。四个节点只有 2 个内节点，方程组缩成了 $2\times2$ ; $n$ 个节点时也只是 $(n-2)\times(n-2)$ 的三对角系统，稀疏、对称、正定，解起来又快又稳。所以实际的样条求解器几乎都走三弯矩这条路，而不是直接列 $4n$ 个方程。

三次样条的局限

样条的好处是时域上顺，它的限制也几乎都从这里来。

端点要单独处理。中间节点左右都有邻居，连续条件自动管住；端点只有一边数据，所以前面才需要额外补两个边界条件。常见的有自然边界、钳制边界和 not-a-knot：自然边界让两端少受力（二阶导为 0），钳制边界告诉曲线从端点出发的方向，not-a-knot 让端点附近两段共用一个多项式。

图 B.6-5 把三种边界画在同一组样本上。中间区域三条曲线几乎重合，但在阴影标出的两端附近明显分开。用样条做外推或关心端点行为时，这一点要留神。

样条追求平滑，反而是个问题。如果数据本身有尖峰或阶跃，三次样条为了保持平滑，会在突变附近冲出样本值。

图 B.6-6 用一段阶跃数据演示：所有样本都是 0 或 1，线性插值老老实实停在 $[0,1]$ 之间，三次样条却在跳变前后冲到 1 以上、又跌到 0 以下。它把你不想抹掉的边缘，先替你抹圆了。数据里有真实突变时，这种过冲可能带来假信号。

样条不是奔着频域设计的。B.4、B.5 的低通插值在设计上很明确，保留基带、压制镜像，能回答阻带压了多少、通带掉了多少。三次样条是一种时域形状控制方法，标准是连接处顺不顺，不是频谱干不干净。所以样条适合曲线应该顺的场景，比如稀疏测量点的平滑、轨迹和动画插帧；要做高保真的采样率转换，主链还是回到低通插值那条路。

也可以反过来，把样条放进频域里观察。比如固定一组采样间隔，对一段输入算出样条输出，再比对两者的频谱，凑出一个等效的频率响应。这条路能走，但凑出来的滤波器会随输入信号变化，平滑信号和含尖峰信号的等效响应不同，不是定常线性滤波器。拿它做频域设计工具并不合适。

配套资源的 MATLAB 代码走完整流程，从三对角方程到各段系数再到任意位置求值。改动样本点或边界条件，直接观察曲线变化。