平滑样条曲线#

对于插值问题，任务是构造一条通过给定数据点的曲线。如果数据存在噪声，这可能不合适：我们希望构造一条平滑曲线 g(x)，它*近似*输入数据而不精确通过每个点。为此，scipy.interpolate 允许基于 P. Dierckx 的 Fortran 库 FITPACK 构建*平滑样条曲线*。

具体来说，给定数据数组 x 和 y 以及非负*权重*数组 w，我们寻找一个样条函数 g(x)，它满足

\[\sum_j \left[ w_j (g(x_j) - y_j)\right]^2 \leqslant s\]

其中 s 是输入参数，控制结果曲线 g(x) 的平滑度与数据近似质量（即 \(g(x_j)\) 和 \(y_j\) 之间的差异）之间的权衡。

请注意，极限 s = 0 对应于插值问题，其中 \(g(x_j) = y_j\)。增加 s 会导致更平滑的拟合，而在非常大的 s 极限下，\(g(x)\) 退化为单个最佳拟合多项式。

找到 s 参数的良好值是一个试错过程。如果权重对应于输入数据标准偏差的倒数，则 s 的“良好”值预计在 \(m - \sqrt{2m}\) 和 \(m + \sqrt{2m}\) 之间，其中 \(m\) 是数据点的数量。如果所有权重等于一，一个合理的猜测可能是大约 \(s \sim m\,\sigma^2\)，其中 \(\sigma\) 是对数据标准偏差的估计。

在内部，FITPACK 库通过向样条拟合 g(x) 添加内部节点来工作，因此**结果节点不一定与输入数据重合**。

一维样条平滑#

scipy.interpolate 提供了 FITPACK 库的两个接口，一个是函数式接口，另一个是面向对象的接口。虽然它们是等价的，但这些接口有不同的默认设置。下面我们将分别讨论它们，首先从函数式接口开始——我们推荐在新代码中使用它。

过程式 (splrep)#

样条插值需要两个基本步骤：(1) 计算曲线的样条表示，(2) 在所需点处对样条进行求值。为了找到样条表示，有两种不同的方式来表示曲线并获得（平滑）样条系数：直接方式和参数化方式。直接方法使用函数 splrep 在二维平面中找到曲线的样条表示。前两个参数是唯一必需的，它们提供了曲线的 \(x\) 和 \(y\) 分量。正常输出是一个三元组 \(\left(t,c,k\right)\)，包含节点点 \(t\)、系数 \(c\) 和样条的阶数 \(k\)。默认的样条阶数是三次，但可以通过输入关键字 k 来更改。

节点数组定义了插值区间为 t[k:-k]，因此 t 数组的前 \(k+1\) 和后 \(k+1\) 个条目定义了*边界节点*。系数是一个长度至少为 len(t) - k - 1 的一维数组。一些例程会填充这个数组以使其长度为 len(c) == len(t)——这些额外的系数在样条求值时会被忽略。

tck 元组格式与 插值 B 样条 兼容：splrep 的输出可以包装成一个 BSpline 对象，例如 BSpline(*tck)，并且下面描述的求值/积分/求根例程可以互换使用 tck 元组和 BSpline 对象。

对于 N 维空间中的曲线，函数 splprep 允许定义曲线 参数化。对于此函数，只需要1个输入参数。该输入是一个包含 \(N\) 个数组的列表，表示 N-D 空间中的曲线。每个数组的长度是曲线的点数，每个数组提供 N-D 数据点的一个分量。参数变量通过关键字参数 u 给出，默认情况下是一个在 \(0\) 和 \(1\) 之间均匀分布的单调序列（即，均匀参数化）。

输出包含两个对象：一个三元组 \(\left(t,c,k\right)\)，包含样条表示和参数变量 \(u\)。

系数是一个包含 \(N\) 个数组的列表，其中每个数组对应于输入数据的一个维度。注意，节点 t 对应于曲线 u 的参数化。

关键字参数 s 用于指定在样条拟合过程中执行的平滑量。\(s\) 的默认值是 \(s=m-\sqrt{2m}\)，其中 \(m\) 是拟合的数据点数。因此，如果不希望进行平滑，应将值 \(\mathbf{s}=0\) 传递给例程。

一旦确定了数据的样条表示，可以使用 splev 函数或通过将 tck 元组包装成 BSpline 对象来对其进行评估，如下所示。

我们首先通过在某些合成噪声数据上展示 s 参数对平滑效果的影响来开始说明。

>>> import numpy as np
>>> from scipy.interpolate import splrep, BSpline

生成一些噪声数据

>>> x = np.arange(0, 2*np.pi+np.pi/4, 2*np.pi/16)
>>> rng = np.random.default_rng()
>>> y =  np.sin(x) + 0.4*rng.standard_normal(size=len(x))

构建两个具有不同 s 值的样条。

>>> tck = splrep(x, y, s=0)
>>> tck_s = splrep(x, y, s=len(x))

并绘制它们

>>> import matplotlib.pyplot as plt
>>> xnew = np.arange(0, 9/4, 1/50) * np.pi

>>> plt.plot(xnew, np.sin(xnew), '-.', label='sin(x)')
>>> plt.plot(xnew, BSpline(*tck)(xnew), '-', label='s=0')
>>> plt.plot(xnew, BSpline(*tck_s)(xnew), '-', label=f's={len(x)}')
>>> plt.plot(x, y, 'o')
>>> plt.legend()
>>> plt.show()
../../_images/smoothing_splines-1.png

我们看到，s=0 曲线跟随数据点的（随机）波动， 而 s > 0 曲线接近于基础的正弦函数。 还要注意，外推值根据 s 的值变化很大。

s 的默认值取决于是否提供了权重， 并且对于 splrep 和 splprep 也不同。因此，我们建议始终 显式提供 s 的值。

操作样条对象：过程式（splXXX）#

一旦确定了数据的样条表示， 就可以使用函数在任何点评估样条 （splev）及其导数 （splev，spalde）， 以及样条在任意两点之间的积分（ splint）。此外，对于三次样条（\(k=3\)） 且有8个或更多节点，可以估计样条的根（ sproot）。这些函数在下面的示例中进行了演示。

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from scipy import interpolate

三次样条

>>> x = np.arange(0, 2*np.pi+np.pi/4, 2*np.pi/8)
>>> y = np.sin(x)
>>> tck = interpolate.splrep(x, y, s=0)
>>> xnew = np.arange(0, 2*np.pi, np.pi/50)
>>> ynew = interpolate.splev(xnew, tck, der=0)

请注意，最后一行等价于 BSpline(*tck)(xnew)。

>>> plt.figure()
>>> plt.plot(x, y, 'x', xnew, ynew, xnew, np.sin(xnew), x, y, 'b')
>>> plt.legend(['线性', '三次样条', '真实'])
>>> plt.axis([-0.05, 6.33, -1.05, 1.05])
>>> plt.title('三次样条插值')
>>> plt.show()

样条的导数

>>> yder = interpolate.splev(xnew, tck, der=1)   # or BSpline(*tck)(xnew, 1)
>>> plt.figure()
>>> plt.plot(xnew, yder, xnew, np.cos(xnew),'--')
>>> plt.legend(['Cubic Spline', 'True'])
>>> plt.axis([-0.05, 6.33, -1.05, 1.05])
>>> plt.title('Derivative estimation from spline')
>>> plt.show()

样条的所有导数

>>> yders = interpolate.spalde(xnew, tck)
>>> plt.figure()
>>> for i in range(len(yders[0])):
...    plt.plot(xnew, [d[i] for d in yders], '--', label=f"{i} derivative")
>>> plt.legend()
>>> plt.axis([-0.05, 6.33, -1.05, 1.05])
>>> plt.title('B样条的所有导数')
>>> plt.show()

样条的积分

>>> def integ(x, tck, constant=-1):
...     x = np.atleast_1d(x)
...     out = np.zeros(x.shape, dtype=x.dtype)
...     for n in range(len(out)):
...         out[n] = interpolate.splint(0, x[n], tck)
...     out += constant
...     return out

>>> yint = integ(xnew, tck)
>>> plt.figure()
>>> plt.plot(xnew, yint, xnew, -np.cos(xnew), '--')
>>> plt.legend(['Cubic Spline', 'True'])
>>> plt.axis([-0.05, 6.33, -1.05, 1.05])
>>> plt.title('Integral estimation from spline')
>>> plt.show()

样条的根

>>> interpolate.sproot(tck)
array([3.1416])  # 可能会有所不同

注意，sproot 可能无法在逼近区间的边缘找到明显的解，即 \(x = 0\)。如果我们在稍大的区间上定义样条，我们将恢复两个根 \(x = 0\) 和 \(x = \pi\)：

>>> x = np.linspace(-np.pi/4, np.pi + np.pi/4, 51)
>>> y = np.sin(x)
>>> tck = interpolate.splrep(x, y, s=0)
>>> interpolate.sproot(tck)
array([0., 3.1416])

参数样条

>>> t = np.arange(0, 1.1, .1)
>>> x = np.sin(2*np.pi*t)
>>> y = np.cos(2*np.pi*t)
>>> tck, u = interpolate.splprep([x, y], s=0)
>>> unew = np.arange(0, 1.01, 0.01)
>>> out = interpolate.splev(unew, tck)
>>> plt.figure()

>>> plt.plot(x, y, 'x', out[0], out[1], np.sin(2*np.pi*unew), np.cos(2*np.pi*unew), x, y, 'b')
>>> plt.legend(['线性', '三次样条', '真实'])
>>> plt.axis([-1.05, 1.05, -1.05, 1.05])
>>> plt.title('参数化定义曲线的样条')
>>> plt.show()

注意在最后一个例子中，splprep 返回的样条系数是一个数组列表，其中每个数组对应输入数据的一个维度。因此，要将其输出包装为 BSpline，我们需要对系数进行转置（或使用 BSpline(..., axis=1)）：

>>> tt, cc, k = tck
>>> cc = np.array(cc)
>>> bspl = BSpline(tt, cc.T, k)    # 注意转置
>>> xy = bspl(u)
>>> xx, yy = xy.T   # 转置以解包成一对数组
>>> np.allclose(x, xx)
True
>>> np.allclose(y, yy)
True

比数据点的数量更多，因此不再严格是一个插值样条，而是一个平滑样条。如果不希望如此，可以使用 InterpolatedUnivariateSpline 类。它是 UnivariateSpline 的子类，总是通过所有点（相当于强制平滑参数为0）。这个类在下面的示例中进行了演示。

LSQUnivariateSpline 类是 UnivariateSpline 的另一个子类。它允许用户通过参数 t 明确指定内部节点的数量和位置。这使得可以创建具有非线性间距的自定义样条，以在某些域中进行插值，在其他域中进行平滑，或改变样条的特性。

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from scipy import interpolate

插值一元样条

>>> x = np.arange(0, 2*np.pi+np.pi/4, 2*np.pi/8)
>>> y = np.sin(x)
>>> s = interpolate.InterpolatedUnivariateSpline(x, y)
>>> xnew = np.arange(0, 2*np.pi, np.pi/50)
>>> ynew = s(xnew)

>>> plt.figure()
>>> plt.plot(x, y, 'x', xnew, ynew, xnew, np.sin(xnew), x, y, 'b')
>>> plt.legend(['线性', '插值一元样条', '真实'])
>>> plt.axis([-0.05, 6.33, -1.05, 1.05])
>>> plt.title('插值一元样条')
>>> plt.show()

具有非均匀节点的 LSQ 一元样条

>>> t = [np.pi/2-.1, np.pi/2+.1, 3*np.pi/2-.1, 3*np.pi/2+.1]
>>> s = interpolate.LSQUnivariateSpline(x, y, t, k=2)
>>> ynew = s(xnew)

>>> plt.figure()
>>> plt.plot(x, y, 'x', xnew, ynew, xnew, np.sin(xnew), x, y, 'b')
>>> plt.legend(['线性', 'LSQ 一元样条', '真实'])
>>> plt.axis([-0.05, 6.33, -1.05, 1.05])
>>> plt.title('具有指定内部节点的样条')
>>> plt.show()

二维平滑样条#

除了平滑一维样条外，FITPACK 库还提供了平滑二维样条的方法。 拟合二维数据到二维*曲面*。这些曲面可以被视为两个参数的函数，\(z = g(x, y)\)，通过一维样条的张量积构造而成。

假设数据保存在三个数组中，x、y 和 z，这些数据数组有两种解释方式。第一种——*散点*插值问题——假设数据是成对的，即 x[i] 和 y[i] 的值对表示点 i 的坐标，对应于 z[i]。

构造曲面 \(g(x, y)\) 以满足

\[\sum_i \left[ w_i (g(x_i, y_i) - z_i)\right]^2 \leqslant s\]

其中 \(w_i\) 是非负权重，s 是输入参数，称为*平滑因子*，它控制了所得函数 g(x, y) 的平滑度与数据逼近质量（即 \(g(x_i, y_i)\) 与 \(z_i\) 之间的差异）之间的平衡。\(s = 0\) 的极限形式上对应于插值，此时曲面通过输入数据，\(g(x_i, y_i) = z_i\)。但请参见下面的注释。

第二种情况——*矩形网格*插值问题——假设数据点位于由 x 和 y 数组的所有元素对定义的矩形网格上。对于这个问题，假设 z 数组是二维的，z[i, j] 对应于 (x[i], y[j])。构造二元样条函数 \(g(x, y)\) 以满足

\[\sum_i \sum_j \left[ (g(x_i, y_j) - z_{i,j})\right]^2 \leqslant s\]

其中 s 是平滑因子。这里 \(s=0\) 的极限形式上也对应于插值，\(g(x_i, y_j) = z_{i, j}\)。

备注

在内部，平滑曲面 \(g(x, y)\) 是通过在由数据数组定义的边界框中放置样条节点来构造的。节点通过 FITPACK 算法自动放置，直到达到所需的平滑度。 已达到。

结点可以放置在数据点之外。

虽然 \(s=0\) 形式上对应于双变量样条插值， FITPACK 算法并不适用于插值，可能会导致意外结果。

对于散点数据插值，建议使用 griddata；对于规则网格上的数据，建议使用 RegularGridInterpolator。

备注

如果输入数据 x 和 y 使得输入维度 具有不一致的单位并且相差多个数量级，插值函数 \(g(x, y)\) 可能会有数值伪影。考虑在插值前对数据进行缩放。

我们现在依次考虑两个样条拟合问题。

散点数据的双变量样条拟合#

FITPACK 库有两个接口，一个是过程式接口，另一个是面向对象的接口。

过程式接口 (bisplrep)#

对于拟合二维曲面的（平滑）样条，可以使用函数 bisplrep。该函数需要输入 一维 数组 x、y 和 z，它们表示曲面 \(z=f(x, y)\) 上的点。 x 和 y 方向上的样条阶数可以通过可选参数 kx 和 ky 指定。默认是双三次样条，kx=ky=3。

bisplrep 的输出是一个列表 [tx, ty, c, kx, ky]，其条目分别表示结点位置的分量、样条的系数以及每个坐标中样条的阶数。 将这个列表保存在一个单一对象 tck 中很方便，这样它可以很容易地传递给函数 bisplev。 关键字 s 可以用来改变在确定适当样条时对数据进行的平滑量。推荐的 \(s\) 值取决于权重 \(w_i\)。如果这些权重取为 \(1/d_i\)， 其中 \(d_i\) 是 \(z_i\) 标准差的估计值，应在范围 \(m- \sqrt{2m}, m + \sqrt{2m}\) 内找到一个好的 \(s\) 值，其中 \(m\) 是 x、y 和 z 向量中的数据点数量。

默认值为 \(s=m-\sqrt{2m}\)。因此，如果不希望进行平滑处理，则应将 ``s=0`` 传递给 `bisplrep`。（但请参阅上面的注释）。

要评估二维样条及其偏导数（最高到样条的阶数），需要使用函数 bisplev。该函数的前两个参数是 两个一维数组，它们的叉积指定了评估样条的域。第三个参数是 bisplrep 返回的 tck 列表。如果需要，第四和第五个参数分别提供 \(x\) 和 \(y\) 方向上的偏导数的阶数。

备注

需要注意的是，二维插值不应用于寻找图像的样条表示。所使用的算法不适用于大量输入点。scipy.signal 和 scipy.ndimage 包含更适合寻找图像样条表示的算法。

二维插值命令旨在用于插值二维函数，如以下示例所示。此示例使用了 NumPy 中的 mgrid 命令，该命令在多维中定义“网格”非常有用。（如果不需要完整的网格，也可以使用 ogrid 命令）。输出参数的数量和每个参数的维度数量由传递给 mgrid 的索引对象的数量决定。

>>> import numpy as np
>>> from scipy import interpolate
>>> import matplotlib.pyplot as plt

在稀疏的 20x20 网格上定义函数

>>> x_edges, y_edges = np.mgrid[-1:1:21j, -1:1:21j]

>>> x = x_edges[:-1, :-1] + np.diff(x_edges[:2, 0])[0] / 2.
>>> y = y_edges[:-1, :-1] + np.diff(y_edges[0, :2])[0] / 2.
>>> z = (x+y) * np.exp(-6.0*(x*x+y*y))

>>> plt.figure()
>>> lims = dict(cmap='RdBu_r', vmin=-0.25, vmax=0.25)
>>> plt.pcolormesh(x_edges, y_edges, z, shading='flat', **lims)
>>> plt.colorbar()
>>> plt.title("Sparsely sampled function.")
>>> plt.show()
" "

在新 70x70 网格上插值函数

>>> xnew_edges, ynew_edges = np.mgrid[-1:1:71j, -1:1:71j]
>>> xnew = xnew_edges[:-1, :-1] + np.diff(xnew_edges[:2, 0])[0] / 2.
>>> ynew = ynew_edges[:-1, :-1] + np.diff(ynew_edges[0, :2])[0] / 2.
>>> tck = interpolate.bisplrep(x, y, z, s=0)
>>> znew = interpolate.bisplev(xnew[:,0], ynew[0,:], tck)

>>> plt.figure()
>>> plt.pcolormesh(xnew_edges, ynew_edges, znew, shading='flat', **lims)
>>> plt.colorbar()
>>> plt.title("Interpolated function.")
>>> plt.show()
" "

面向对象接口 (SmoothBivariateSpline)#

面向对象接口用于散点数据的二元样条平滑，SmoothBivariateSpline 类，实现了 bisplrep / bisplev 对功能的一个子集，并且有不同的默认值。

它将权重数组的元素设为统一值，\(w_i = 1\)，并根据平滑因子 s 的输入值自动构建节点向量——默认值是 \(m\)，即数据点的数量。

样条在 x 和 y 方向的阶数由可选参数 kx 和 ky 控制，默认值为 kx=ky=3。

我们使用以下示例来说明平滑因子的效果：

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import SmoothBivariateSpline

import warnings
warnings.simplefilter('ignore')

train_x, train_y = np.meshgrid(np.arange(-5, 5, 0.5), np.arange(-5, 5, 0.5))
train_x = train_x.flatten()
train_y = train_y.flatten()

def z_func(x, y):
    return np.cos(x) + np.sin(y) ** 2 + 0.05 * x + 0.1 * y

train_z = z_func(train_x, train_y)
interp_func = SmoothBivariateSpline(train_x, train_y, train_z, s=0.0)
smth_func = SmoothBivariateSpline(train_x, train_y, train_z)

test_x = np.arange(-9, 9, 0.01)
test_y = np.arange(-9, 9, 0.01)
grid_x, grid_y = np.meshgrid(test_x, test_y)

interp_result = interp_func(test_x, test_y).T
smth_result = smth_func(test_x, test_y).T
perfect_result = z_func(grid_x, grid_y)

fig, axes = plt.subplots(1, 3, figsize=(16, 8))
extent = [test_x[0], test_x[-1], test_y[0], test_y[-1]]
opts = dict(aspect='equal', cmap='nipy_spectral', extent=extent, vmin=-1.5, vmax=2.5)

im = axes[0].imshow(perfect_result, **opts)
fig.colorbar(im, ax=axes[0], orientation='horizontal')
axes[0].plot(train_x, train_y, 'w.')
axes[0].set_title('完美结果，采样函数', fontsize=21)

im = axes[1].imshow(smth_result, **opts)
axes[1].plot(train_x, train_y, 'w.')
fig.colorbar(im, ax=axes[1], orientation='horizontal')
axes[1].set_title('s=默认值', fontsize=21)

im = axes[2].imshow(interp_result, **opts)
fig.colorbar(im, ax=axes[2], orientation='horizontal')
axes[2].plot(train_x, train_y, 'w.')
axes[2].set_title('s=0', fontsize=21)

plt.tight_layout()
plt.show()
../../_images/smoothing_splines-5.png

这里我们取一个已知函数（显示在最左边的面板），在网格点上对其进行采样（由白点表示），并使用默认平滑（中间面板）和强制插值（最右边面板）构建样条拟合。

几个特征清晰可见。首先，默认的 s 值对这些数据提供了过多的平滑；强制插值条件，s = 0，可以以合理的精度恢复底层函数。其次， 超出插值范围（即白色点覆盖的区域）的结果使用最近邻常数进行外推。最后，我们不得不静默警告（是的，这是不好的形式！）。

这里的警告在 s=0 的情况下发出，表示当我们强制插值条件时，FITPACK 遇到了内部困难。如果在代码中看到此警告，请考虑切换到 bisplrep 并增加其 nxest、nyest 参数（详见 bisplrep 文档字符串）。

二维网格数据的二元样条拟合#

对于网格化的二维数据，可以使用 RectBivariateSpline 类进行平滑张量积样条拟合。它的接口与 SmoothBivariateSpline 类似，主要区别在于一维输入数组 x 和 y 被理解为定义了一个二维网格（作为它们的外积），而 z 数组是二维的，形状为 len(x) 乘以 len(y)。

x 和 y 方向上的样条阶数由可选参数 kx 和 ky 控制，默认值为 kx=ky=3，即双三次样条。

平滑因子的默认值为 s=0。但我们仍然建议始终显式指定 s。

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import RectBivariateSpline

x = np.arange(-5.01, 5.01, 0.25)        # 网格是 x 和 y 数组的外积
y = np.arange(-5.01, 7.51, 0.25)        #

xx, yy = np.meshgrid(x, y, indexing='ij')
z = np.sin(xx**2 + 2.*yy**2)            # z 数组需要是二维的

func = RectBivariateSpline(x, y, z, s=0)

xnew = np.arange(-5.01, 5.01, 1e-2)
ynew = np.arange(-5.01, 7.51, 1e-2)
znew = func(xnew, ynew)

plt.imshow(znew)
plt.colorbar()
plt.show()
../../_images/smoothing_splines-6.png

球坐标系中数据的二元样条拟合#

如果你的数据是以球坐标形式给出的，即 \(r = r(\theta, \phi)\)， SmoothSphereBivariateSpline 和 RectSphereBivariateSpline 分别提供了与 SmoothBivariateSpline 和 RectBivariateSpline 类似的便捷方法。

这些类确保了在 \(\theta \in [0, \pi]\) 和 \(\phi \in [0, 2\pi]\) 范围内的样条拟合的周期性，并且提供了一些对极点处连续性的控制。有关详细信息，请参阅这些类的文档字符串。

On this page