一维插值:从折线、三次样条到保形曲线与批次数据

原文:SciPy 文档贡献者,1-D interpolation。本文依据 2026 年 10 月 5 日取得的 SciPy 1.18.0 手册全文翻译整理,保留主要算法、示例与适用边界;两处原文问题在相应位置明确订正。示例仅做静态审查,没有实际运行。

插值是在已知采样点之间构造函数,使它经过这些数据点。选择方法时,需要决定的并不只是“曲线够不够平滑”:还要考虑是否允许过冲、是否需要导数、数据有没有批次维度,以及区间之外应当返回什么。以下从最简单的折线开始,逐步说明这些取舍。

一维插值选择示意:线性插值强调简单,三次样条提供二阶连续导数,PCHIP强调保形;构造边界条件和区间外求值分别决定不同性质。
图:未完纪为本文绘制的方法与边界条件示意;不是数值运行结果。

1. 分段线性插值

如果只需要在相邻点之间连直线,可以使用 numpy.interp。输入 x 是采样位置,y 是对应值,xnew 是希望求值的位置:

import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(0, 10, num=11)
y = np.cos(-x**2 / 9.0)
xnew = np.linspace(0, 10, num=1001)
ynew = np.interp(xnew, x, y)

plt.plot(xnew, ynew, '-', label='linear interp')
plt.plot(x, y, 'o', label='data')
plt.legend(loc='best')
plt.show()

这种方法容易理解,但连接处通常有折角。numpy.interp 默认在区间两端之外返回端点值;它还允许指定 left、right 常量或周期参数,不过没有样条那样的边界多项式外推接口。若希望沿边缘直线继续延伸,可以使用后面介绍的 make_interp_spline(..., k=1)。

2. 三次样条及其导数

CubicSpline 把若干三次多项式拼成一条曲线,并使连接处的一阶、二阶导数连续。创建对象后,可以像调用函数一样传入标量或数组;nu 表示要求的导数阶数。

from scipy.interpolate import CubicSpline

spl = CubicSpline(x, y)
fig, ax = plt.subplots(4, 1, figsize=(5, 7))
ax[0].plot(xnew, spl(xnew), label='spline')
ax[0].plot(x, y, 'o', label='data')
ax[1].plot(xnew, spl(xnew, nu=1), '--', label='1st derivative')
ax[2].plot(xnew, spl(xnew, nu=2), '--', label='2nd derivative')
ax[3].plot(xnew, spl(xnew, nu=3), '--', label='3rd derivative')
for a in ax:
    a.legend(loc='best')
plt.tight_layout()
plt.show()

前两阶导数连续来自构造约束,三阶导数可以在采样点处跳变。这里得到的是插值函数的导数,不等于真实过程的导数;数据噪声、采样间隔和边界条件都会影响结果。

3. 需要保留形状时,检查过冲

二阶连续不一定符合数据本身。在尖峰或异常点附近,普通三次样条可能振荡,甚至超出相邻数据的范围。PchipInterpolator 和 Akima1DInterpolator 都只保证一阶连续,侧重保留局部形状。PCHIP 在单调数据区间保持单调;不能据此把 Akima 也说成对任意数据都有同样的单调性保证。

from scipy.interpolate import (
    CubicSpline, PchipInterpolator, Akima1DInterpolator,
)

x = np.array([1., 2., 3., 4., 4.5, 5., 6., 7., 8.])
y = x**2
y[4] += 101
xx = np.linspace(1, 8, 51)

plt.plot(xx, CubicSpline(x, y)(xx), '--', label='spline')
plt.plot(xx, Akima1DInterpolator(x, y)(xx), '-', label='Akima1D')
plt.plot(xx, PchipInterpolator(x, y)(xx), '-', label='pchip')
plt.plot(x, y, 'o', label='data')
plt.legend()
plt.show()

原教程用这一组带尖峰的数据对照三种方法。代码中的尖峰是有意加入的演示数据,不能把插值当成异常值检测或去噪:三种方法仍然都试图经过输入点。真正的问题若是测量误差,应先决定数据清理或平滑策略。

4. B 样条表示与非三次插值

B 样条与分段多项式在表示能力上对应,但采用不同的基函数,通常比幂基表示更利于数值计算。它除了插值,也用于回归和曲线建模。make_interp_spline 根据数据构造 BSpline 对象:

from scipy.interpolate import make_interp_spline

x = np.linspace(0, 3 / 2, 7)
y = np.sin(np.pi * x)
bspl = make_interp_spline(x, y, k=3)
der = bspl.derivative()
xx = np.linspace(0, 3 / 2, 51)
plt.plot(xx, bspl(xx), '--', label='sin(pi*x) approx')
plt.plot(x, y, 'o', label='data')
plt.plot(xx, der(xx) / np.pi, '--', label='derivative / pi')
plt.legend()
plt.show()

k=3 是默认的三次样条,其一阶导数为二次样条,即 bspl.k 与 der.k 分别为 3 和 2。默认情况下,make_interp_spline(x, y) 与 CubicSpline(x, y) 产生等价的插值函数;前者还允许用 k 选择次数、用 t 指定节点。两者默认边界条件都是 not-a-knot。

把次数设为 1,就能得到带线性外推的折线:

x = np.linspace(0, 5, 11)
y = 2 * x
spl = make_interp_spline(x, y, k=1)
outside_linear = spl([-1, 6])
outside_constant = np.interp([-1, 6], x, y)

按公式,前者应得到 [-2, 12],后者默认返回 [0, 10]。这是代码语义和原教程所给结果的说明,本次没有运行验证。线性外推也不意味着区间之外的预测有精度保证。

5. 一次处理多组 y:弄清插值轴

一元插值器不一定只接受一维 y。若在同一组位置 x[i] 上测量多条曲线,可以让 y[i, j] 保存第 j 条曲线的值。默认沿 axis=0 插值,其余维度作为批次维度:

n = 11
x = 2 * np.pi * np.arange(n) / n
y = np.stack((np.sin(x)**2, np.cos(x)), axis=1)
# x.shape 为 (11,),y.shape 为 (11, 2)
spl = make_interp_spline(x, y)
xv = np.linspace(0, 2 * np.pi, 51)
values = spl(xv)  # 输出形状为 (51, 2)
plt.plot(x, y, 'o')
plt.plot(xv, values, '-')
plt.show()

这看起来像 NumPy 广播,但有两个区别:x 仍必须是一维数组,不能把它与 y 自动广播;默认使用尾部维度存放批次,也不同于许多 NumPy 接口把前部维度作为批次的约定。指定其他插值轴时,必须满足 y.shape[axis] == x.size。

原文订正:教程给出的输出形状公式包含了不正确的切片顺序。依据 BSpline.__call__ API 文档,输出应当用求值位置的形状替换插值轴,而不是重新排列整个数组。将 axis 规范化为非负索引后,公式为:

y.shape[:axis] + xv.shape + y.shape[axis + 1:]

例如 y.shape == (2, 11, 3)、axis=1、xv.shape == (51,) 时,输出为 (2, 51, 3)。PCHIP、Akima、CubicSpline,底层的 PPoly、BPoly、BSpline,以及多种最小二乘和平滑样条也支持类似批次用法,具体参数仍要以各自 API 为准。

6. 把平面坐标变成参数曲线

前面的例子认为 y 是 x 的函数。参数曲线则把 (x, y) 都视为坐标,用另外的参数 u 描述曲线经过这些点的顺序。这个方法也适用于三维或更高维坐标。困难在于:同一组点采用不同的参数化方式,得到的曲线可能明显不同。

x = [0, 1, 2, 3, 4, 5, 6]
y = [0, 0, 0, 9, 0, 0, 0]
p = np.stack((x, y))  # shape: (2, 7),每一列是一个点

u_uniform = np.arange(p.shape[1], dtype=float)
dp = p[:, 1:] - p[:, :-1]
squared_lengths = (dp**2).sum(axis=0)
u_chord = np.r_[0, np.sqrt(squared_lengths).cumsum()]
u_centripetal = np.r_[0, (squared_lengths**0.25).cumsum()]

fig, ax = plt.subplots(1, 3, figsize=(8, 3))
for a, u, name in zip(
    ax,
    [u_uniform, u_chord, u_centripetal],
    ['uniform', 'chord length', 'centripetal'],
):
    spl = make_interp_spline(u, p, axis=1)
    uu = np.linspace(u[0], u[-1], 51)
    xx, yy = spl(uu)
    a.plot(xx, yy, '--')
    a.plot(p[0], p[1], 'o')
    a.set_title(name)
plt.show()

均匀参数化令 u[j]=j;弦长参数化逐段累加相邻点的欧氏距离;向心参数化则累加距离的平方根。上面第四次方根是因为 squared_lengths 已经是距离的平方。对其他数据,重复的相邻点可能让弦长或向心参数出现重复值,应先按业务含义处理,不能直接交给要求严格递增参数的构造器。

7. 缺失值不是插值器自动替你处理的空位

scipy.interpolate 并不提供统一的缺失数据插值支持。无论用 numpy.ma 的掩码数组,还是用 NaN 标记缺失,都不能假定插值器会自动忽略。库总体遵循 IEEE 754 的数值语义,NaN 表示不是有效数字,而不是带有“请补齐”含义的特殊记录。

构造前应检查数值是否有限、横坐标是否有序且满足接口的严格递增要求、是否存在重复采样位置,以及删去缺失点后是否还够构造所选次数的样条。这是数据建模步骤;删除、聚合或填补不同,得到的曲线也会不同。

8. 边界条件与区间外求值是两件事

bc_type 决定构造时的边界约束,extrapolate 决定求值超出数据范围时的行为。两者有关联,但不是同一个参数。

对 CubicSpline 和 make_interp_spline,bc_type='periodic' 会对边界施加周期约束,需要首尾数据满足周期要求,并使曲线在接缝处满足相应的光滑条件。相反,只设 extrapolate='periodic' 会把区间外位置绕回一个周期,不能替你修复端点的值或导数不匹配。

x = np.array([0, 1, 2.5, 3, 4])
y = np.array([0.0, 1.0, 0.8, 1.2, 0.0])
x_eval = np.linspace(-4, 8, 300)

wrapped = CubicSpline(x, y, extrapolate='periodic')
periodic = CubicSpline(x, y, bc_type='periodic')
plt.plot(x_eval, wrapped(x_eval), ':', label='wrap evaluation')
plt.plot(x_eval, periodic(x_eval), '-.', label='periodic construction')
plt.plot(x, y, 'ko', label='data')
plt.legend()
plt.show()

多数相关插值器接受以下区间外策略:

  • extrapolate=True:延伸第一段或最后一段多项式。
  • extrapolate=False:区间之外返回 NaN。
  • extrapolate='periodic':在支持此值的接口中,周期性地映射求值位置。

多数接口默认允许外推;教程特别指出 Akima 默认不外推,而周期边界构造通常默认采用周期求值。PCHIP 没有 bc_type,给它周期求值参数并不能保证接缝处光滑:

for mode in [True, False, 'periodic']:
    pchip = PchipInterpolator(x, y, extrapolate=mode)
    plt.plot(x_eval, pchip(x_eval), label=f'extrapolate={mode}')
plt.plot(x, y, 'ko')
plt.legend()
plt.show()

求值时还可以临时覆盖构造时的策略,而不修改对象的默认值:

spl = CubicSpline([0, 1, 2], [0, 1, 0], extrapolate=False)
x_eval = np.array([-1, 0.5, 3])
within_only = spl(x_eval)
allow_outside = spl(x_eval, extrapolate=True)

原教程给出的结果分别为 [nan, 0.75, nan] 和 [-3, 0.75, -3]。传入 None 表示沿用对象默认值。尤其要记住:“自然”等端点导数约束也不能简单理解成“区间外必然是直线”;远离采样域的外推应单独论证。

9. 旧接口 interp1d 与替代方式

interp1d 已标为 legacy,新代码建议选择更明确的插值器;这不等于接口已经移除。官方仍说明会支持既有用法,没有移除计划。旧代码常这样选择方法:

from scipy.interpolate import interp1d
x = np.linspace(0, 10, 11)
y = np.cos(-x**2 / 9.0)
linear = interp1d(x, y)
cubic = interp1d(x, y, kind='cubic')
nearest = interp1d(x, y, kind='nearest')
previous = interp1d(x, y, kind='previous')
next_value = interp1d(x, y, kind='next')

线性情况通常改用 numpy.interp,有多维 y 或特殊外推要求时可考虑一次样条;二次、三次模式可直接使用 make_interp_spline。替换接口时需要核对默认插值轴和边界行为,不能只替换函数名。

分段常量中的 previous 对应使用前一个采样值,在区间内可用 make_interp_spline(x, y, k=0) 表达;nearest 取最近点,next 取下一个点。原文订正:若横坐标代表时间,nearest 和 next 都可能用到未来采样,不能笼统称为因果滤波;previous 才符合“只使用当前或过去样本”的常见保持方式,仍需明确边界约定。

这些模式也可以通过 numpy.searchsorted 构造。下面保留教程的最近邻思路,并以中点平分相邻点;恰好落在中点时,side='left' 选择左侧样本:

x = np.arange(8)
y = x**2
x_new = np.linspace(0, 7, 101)
x_bds = x[:-1] / 2.0 + x[1:] / 2.0
idx = np.searchsorted(x_bds, x_new, side='left')
idx = np.clip(idx, 0, len(x) - 1)
y_new = y[idx]

这段手工版本把超出两端的位置夹到端点值,与旧接口可能报错或填充值的行为不同;这种差异必须在迁移时写清楚。先选择需要保证的性质,再选择接口,通常比单纯追求更高次数稳妥。

来源、修订与许可

原文版权归 SciPy community;页面署名为 SciPy 文档贡献者。SciPy 采用 BSD 3-Clause 许可,完整项目许可文本保留在本文下方。本文补充输入检查、区间外语义和静态审查说明,订正原文的批次形状公式及 nearest/next 的因果性表述;图为未完纪原创。参考:原教程、BSpline 求值形状、SciPy 许可。

SciPy BSD 3-Clause 许可证

Copyright (c) 2001-2002 Enthought, Inc. 2003, SciPy Developers.
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions
are met:

1. Redistributions of source code must retain the above copyright
   notice, this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above
   copyright notice, this list of conditions and the following
   disclaimer in the documentation and/or other materials provided
   with the distribution.
3. Neither the name of the copyright holder nor the names of its
   contributors may be used to endorse or promote products derived
   from this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
© 版权声明
THE END
喜欢就支持一下吧
点赞0 分享
评论 抢沙发

请登录后发表评论

    暂无评论内容