平滑样条曲线#

对于插值问题，任务是构造一条穿过给定数据集点的曲线。如果数据有噪声，这可能不合适：然后我们希望构造一条平滑曲线 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)。此外，对于三次样条 ( \(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(['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 = \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(['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 对的功能子集，并且具有不同的默认值。

它将权重数组的元素取为单位，\(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 文档字符串）。

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

对于网格化二维数据，可以使用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)        # 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]\)的周期性，并提供对极点连续性的控制。有关详细信息，请参阅这些类的文档字符串。