平滑样条

对于插值问题，任务是构造一条通过给定数据集点的曲线。如果数据有噪声，这可能不合适：我们希望构造一条平滑曲线，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\) 是数据点的数量。如果所有权重都等于 1，则一个合理的选择可能是大约 \(s \sim m\,\sigma^2\)，其中 \(\sigma\) 是对数据标准差的估计。

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

一维样条平滑

scipy.interpolate 提供了两种用于 FITPACK 库的接口，一种是函数式接口，另一种是面向对象接口。虽然等效，但这些接口具有不同的默认值。下面我们将依次讨论它们，从函数式接口开始——我们建议在新的代码中使用它。

过程式 (splrep)

样条插值需要两个基本步骤：（1）计算曲线的样条表示，以及（2）在所需点处对样条进行评估。为了找到样条表示，有两种不同的方法来表示曲线并获得（平滑）样条系数：直接法和参数法。直接法使用函数 splrep 找到二维平面中曲线的样条表示。前两个参数是唯一必需的参数，它们提供了曲线的 \(x\) 和 \(y\) 分量。正常的输出是一个 3 元组，\(\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 维空间中的曲线。每个数组的长度是曲线点的数量，每个数组提供 N 维数据点的其中一个分量。参数变量使用关键字参数 u 给出，默认情况下为 \(0\) 和 \(1\) 之间的等间距单调序列（即 均匀参数化）。

输出包含两个对象：一个 3 元组 \(\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()

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

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

操作样条对象：过程式 (splXXX)

一旦确定了数据的样条表示，就可以使用函数来评估样条 (splev) 及其导数 (splev, spalde) 在任何点以及样条在任意两点之间的积分 ( splint)。此外，对于具有 8 个或更多节点的立方样条 ( \(k=3\) )，可以估计样条的根 ( 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(['Linear', 'Cubic Spline', 'True'])
>>> plt.axis([-0.05, 6.33, -1.05, 1.05])
>>> plt.title('Cubic-spline interpolation')
>>> 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('All derivatives of a B-spline')
>>> 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])  # may vary

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

>>> x = np.linspace(-np.pi/4, 2.*np.pi + np.pi/4, 21)
>>> 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(['Linear', 'Cubic Spline', 'True'])
>>> plt.axis([-1.05, 1.05, -1.05, 1.05])
>>> plt.title('Spline of parametrically-defined curve')
>>> plt.show()

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

>>> tt, cc, k = tck
>>> cc = np.array(cc)
>>> bspl = BSpline(tt, cc.T, k)    # note the transpose
>>> xy = bspl(u)
>>> xx, yy = xy.T   # transpose to unpack into a pair of arrays
>>> np.allclose(x, xx)
True
>>> np.allclose(y, yy)
True

面向对象 (UnivariateSpline)

上述描述的样条拟合功能也可以通过面向对象的接口使用。一维样条是 UnivariateSpline 类的对象，并通过将曲线的 \(x\) 和 \(y\) 分量作为参数提供给构造函数来创建。该类定义了 __call__，允许使用 x 轴值调用该对象，在该值处应评估样条，并返回插值的 y 值。这在下面的示例中显示了子类 InterpolatedUnivariateSpline。 integral、derivatives 和 roots 方法也适用于 UnivariateSpline 对象，允许计算样条的定积分、导数和根。

UnivariateSpline 类也可以用于通过提供非零平滑参数值 s 来平滑数据，其含义与上面描述的 splrep 函数的 s 关键字相同。这将导致一个具有比数据点数量更少的节点的样条曲线，因此不再严格是插值样条曲线，而是平滑样条曲线。如果不需要这样做，可以使用 InterpolatedUnivariateSpline 类。它是 UnivariateSpline 的子类，它始终通过所有点（相当于强制平滑参数为 0）。此类在下面的示例中进行了演示。

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

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

InterpolatedUnivariateSpline

>>> 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(['Linear', 'InterpolatedUnivariateSpline', 'True'])
>>> plt.axis([-0.05, 6.33, -1.05, 1.05])
>>> plt.title('InterpolatedUnivariateSpline')
>>> plt.show()

具有非均匀节点的 LSQUnivarateSpline

>>> 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(['Linear', 'LSQUnivariateSpline', 'True'])
>>> plt.axis([-0.05, 6.33, -1.05, 1.05])
>>> plt.title('Spline with Specified Interior Knots')
>>> 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).\) 上的点。可以通过可选参数 kx 和 ky 指定 x 和 y 方向上的样条阶数。默认值为双三次样条，kx=ky=3。

bisplrep 的输出是一个列表 [tx ,ty, c, kx, ky]，其条目分别代表节点位置的各个分量、样条的系数以及每个坐标方向上的样条阶数。将此列表保存在一个名为 tck 的对象中比较方便，这样可以轻松地将其传递给 bisplev 函数。关键字 s 可用于更改在确定适当样条时对数据执行的平滑程度。推荐的 \(s\) 值取决于权重 \(w_i\)。如果将它们设为 \(1/d_i\)，其中 \(d_i\) 是 \(z_i\) 标准差的估计值，则 \(s\) 的一个良好值应该在 \(m- \sqrt{2m}, m + \sqrt{2m}\) 范围内，其中 \(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 对的功能子集，并具有不同的默认值。

它将权重数组的元素设置为等于 1，\(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('Perfect result, sampled function', 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=default', 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()

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

几个特征清晰可见。首先，s 的默认值对这些数据来说过于平滑；强制插值条件，s = 0，可以将底层函数恢复到合理的精度。其次，在插值范围之外（即白点覆盖的区域），结果使用最近邻常数进行外推。最后，我们不得不屏蔽警告（这是一种不好的做法，是的！）。

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

网格上数据的双变量样条拟合

对于网格化的 2D 数据，可以使用 RectBivariateSpline 类拟合平滑张量积样条。它的接口类似于 SmoothBivariateSpline，主要区别在于 1D 输入数组 x 和 y 被理解为定义一个 2D 网格（作为它们的外部积），而 z 数组是 2D 的，形状为 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)        # the grid is an outer product
y = np.arange(-5.01, 7.51, 0.25)        # of x and y arrays

xx, yy = np.meshgrid(x, y, indexing='ij')
z = np.sin(xx**2 + 2.*yy**2)            # z array needs to be 2-D

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()

球坐标系中数据的双变量样条曲线拟合

如果您的数据以球坐标系给出，\(r = r(\theta, \phi)\)，SmoothSphereBivariateSpline 和 RectSphereBivariateSpline 分别提供了 SmoothBivariateSpline 和 RectBivariateSpline 的便捷类比。

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