是通过已知的离散数据求未知数据的方法。与拟合不同的是，要求曲线通过所有的已知数据。SciPy的interpolate模块提供了许多对数据进行插值运算的函数。

一、一维插值

一维数据的插值运算可以通过interpld( )完成。其调用形式如下，它实际上不是函数而是一个类：

interp1d(x, y, kind='linear', ...)

其中，x和y参数是一系列已知的数据点，kind参数是插值类型，可以是字符串或整数，它给出插值的B样条曲线的阶数，可以有如下候选值：

·'zero'、'nearest'：阶梯插值，相当于0阶B样条曲线。

·'slinear'、'linear'：线性插值，用一条直线连接所有的取样点，相当于一阶B样条曲线，'slinear'使用扩展库中的相关函数计算，而'linear'则直接使用Python编写的函数进行运算，其结果一样。

·'quadratic'、'cubic'：二阶和三阶B样条曲线，更高阶的曲线可以直接使用整数值指定。

interpld对象可以计算x的取值范围之内任意点的函数值。它可以像函数一样直接调用，和NumPy的ufunc函数一样能对数组中的每个元素进行计算，并返回一个新的数组。

下面的程序演示了kind参数以及与其对应的插值曲线，结果如图1所示。程序中我们使用循环对相同的数据进行4种不同阶数的插值运算。?首先使用数据点创建一个interpld对象f，通过kind参数指定其阶数。?调用f( )计算出一系列的插值结果。本例中，决定插值曲线的数据点一共有11个，插值之后的曲线数据点有101个。

图1　interpld的各阶插值

高次interp1d( )插值的运算量很大，因此对于点数较多的数据，建议使用后面介绍的UnivariateSpline( )。

from scipy import interpolate

x = np.linspace(0, 10, 11)
y = np.sin(x)

xnew = np.linspace(0, 10, 101)
pl.plot(x,y,'ro')
for kind in ['nearest', 'zero', 'slinear', 'quadratic']:
f = interpolate.interp1d(x,y,kind=kind) ?
ynew = f(xnew) ?
pl.plot(xnew, ynew, label=str(kind))

pl.legend(loc='lower right')

1．外推和Spline拟合

上节所介绍的interpld类要求其参数x是一个递增的序列，并且只能在x的取值范围之内进行内插计算，不能用它进行外推运算，即计算x的取值范围之外的数据点。UnivariateSpline类的插值运算比interpld更高级，它支持外推和拟合运算，其调用形式如下：

UnivariateSpline(x, y, w=None, bbox=[None, None], k=3, s=None)

·x、y是保存数据点的X-Y坐标的数组，其中x必须是递增序列。

·w是为每个数据点指定的权重值。

·k为样条曲线的阶数。

·s是平滑系数，它使得最终生成的样条曲线满足条件：

即当s>0时，样条曲线并不一定通过各个数据点。为了让曲线通过所有数据点，必须将s参数设置为0。此外，还可以使用InterpolatedUnivariateSpline类，它与UnivariateSpline的唯一区别就是它通过所有的数据点，相当于将s设置为0。

下面的程序演示了使用UnivariateSpline对数据进行插值、外推以及样条曲线拟合：

?如图2（上）所示，UnivariateSpline能够进行外推运算，虽然输入数据中没有X轴大于10的点，但是它能计算出X轴在0到12的插值结果。在X轴大于10的部分，样条曲线仍然呈现出和正弦波类似的形状，越远离输入数据范围，误差会越大，因此外推的范围是有限的。由于s参数为0，因此插值曲线经过所有的数据点。

?图2（下）则显示了s参数不为零时的结果，对于带噪声的输入数据，选择合适的s参数能够使得样条曲线接近无噪声时的波形，可以把它看作使用样条曲线对数据进行拟合运算。

x1 = np.linspace(0, 10, 20)
y1 = np.sin(x1)
sx1 = np.linspace(0, 12, 100)
sy1 = interpolate.UnivariateSpline(x1, y1, s=0)(sx1) ?

x2 = np.linspace(0, 20, 200)
y2 = np.sin(x2) + np.random.standard_normal(len(x2))*0.2
sx2 = np.linspace(0, 20, 2000)
spline2 = interpolate.UnivariateSpline(x2, y2, s=8) ?
sy2 = spline2(sx2)

pl.figure(figsize=(8, 5))
pl.subplot(211)
pl.plot(x1, y1, ".", label=u"数据点")
pl.plot(sx1, sy1, label=u"spline曲线")
pl.legend( )

pl.subplot(212)
pl.plot(x2, y2, ".", label=u"数据点")
pl.plot(sx2, sy2, linewidth=2, label=u"spline曲线")
pl.plot(x2, np.sin(x2), label=u"无噪声曲线")
pl.legend( )

?

图2　使用UnivariateSpline进行插值：外推（上）和数据拟合（下）

当曲线为三阶曲线时，UnivariateSpline.roots( )可以用于计算曲线与y = 0的交点横坐标。下面显示了图2（下）中的曲线与X轴的6个交点的横坐标：

print np.array_str( spline2.roots( ), precision=3 )
[  3.288   6.329   9.296  12.578  15.75   18.805]

如果需要计算曲线与任意横线的交点，可以事先将曲线的拟合数据在Y轴方向上进行平移。但若要计算与多条y=c横线的交点，则需要对同样的数据进行平移和拟合多次。

?下面的root_at( )通过直接修改拟合曲线的参数，实现曲线的平移，从而可以计算与任意横线的交点。?将roots_at( )动态添加为UnivariateSpline类的方法。?对多条横线求交点，并进行绘图，其结果如图3所示。

def roots_at(self, v): ?
coeff = self.get_coeffs( )
coeff -= v
try:
root = self.roots( )
return root
finally:
coeff += v

interpolate.UnivariateSpline.roots_at = roots_at ?

pl.plot(sx2, sy2, linewidth=2, label=u"spline曲线")

ax = pl.gca( )
for level in [0.5, 0.75, -0.5, -0.75]:
ax.axhline(level, ls=":", color="k")
xr = spline2.roots_at(level) ?
pl.plot(xr, spline2(xr), "ro")

?

图3　计算Spline与水平线的交点

2．参数插值

前面介绍的插值函数都需要X轴的数据是按照递增顺序排列的，就像一般的y=f(x)函数曲线一样。数学上还有一种参数曲线，它使用参数t和两个函数x=f(t)、y=g(t)来定义二维平面上的一条曲线，例如圆形、心形等曲线都是参数曲线。参数曲线的插值可以通过splprep( )和splev( )实现，这组函数支持高维空间的曲线的插值，这里以二维曲线为例介绍其用法。

?首先调用splprep( )，其第一个参数为一组一维数组，每个数组是各点在对应轴上的坐标。s参数为平滑系数，与UnivariateSpline的含义相同。splprep( )返回两个对象，其中tck是一个元组，它包含了插值曲线的所有信息。t是自动计算出的参数曲线的参数数组。

?调用splev( )进行插值运算，其第一个参数为一个新的参数数组，这里将t的取值范围等分200份，第二个参数为splprep( )返回的第一个对象。实际上，参数数组t是正规化之后的各个线段长度的累计，因此t的范围为0到1。

其结果如图4所示，图中比较了平滑系数为0和1e-4时的插值曲线。当平滑系数为0时，插值曲线通过所有的数据点。

x = [ 4.913,  4.913,  4.918,  4.938,  4.955,  4.949,  4.911,
4.848,  4.864,  4.893,  4.935,  4.981,  5.01 ,  5.021]

y = [ 5.2785,  5.2875,  5.291 ,  5.289 ,  5.28  ,  5.26  ,  5.245 ,
5.245 ,  5.2615,  5.278 ,  5.2775,  5.261 ,  5.245 ,  5.241]

pl.plot(x, y, "o")

for s in (0, 1e-4):
tck, t = interpolate.splprep([x, y], s=s) ?
xi, yi = interpolate.splev(np.linspace(t[0], t[-1], 200), tck) ?
pl.plot(xi, yi, lw=2, label=u"s=%g" % s)

pl.legend( )

图4　使用参数插值连接二维平面上的点

3．单调插值

前面介绍的几种插值方法不能保证数据点的单调性，即曲线的最值可能出现在数据点之外的地方。PchipInterpolator类（别名pchip）使用单调三次插值，能够保证曲线的所有最值都出现在数据点之上。下面的程序用pchip( )对数据点进行插值，并绘制其一阶导数曲线，由图5的导数曲线可知，所有最值点处的导数都为0。

x = [0, 1, 2, 3, 4, 5]
y = [1, 2, 1.5, 2.5, 3, 2.5]
xs = np.linspace(x[0], x[-1], 100)
curve = interpolate.pchip(x, y)
ys = curve(xs)
dys = curve.derivative(xs)
pl.plot(xs, ys, label=u"pchip")
pl.plot(xs, dys, label=u"一阶导数")
pl.plot(x, y, "o")
pl.legend(loc="best")
pl.grid( )
pl.margins(0.1, 0.1)

图5　单调插值能保证两个点之间的曲线为单调递增或递减

二、多维插值

使用interp2d( )可以进行二维插值运算。它的调用形式如下：

interp2d(x, y, z, kind='linear', ...)

其中x、y、z都是一维数组，如果传入的是多维数组，则先将其转为一维数组。kind参数指定插值运算的阶数，可以为'linear'、'cubic'、'quintic'。

下面的例子对某个函数曲面上的网格点进行二维插值，效果如图6所示。其中左图显示插值之前的数据，而右图显示插值运算后得到的结果。

def func(x, y): ?
return (x+y)*np.exp(-5.0*(x**2 + y**2))

# X-Y轴分为15*15的网格
y, x = np.mgrid[-1:1:15j, -1:1:15j] ?
fvals = func(x, y) # 计算每个网格点上的函数值

# 二维插值
newfunc = interpolate.interp2d(x, y, fvals, kind='cubic') ?

# 计算100*100的网格上的插值
xnew = np.linspace(-1,1,100)
ynew = np.linspace(-1,1,100)
fnew = newfunc(xnew, ynew) ?

图6　使用interp2d类进行二维插值

?func是计算曲面上各点高度的函数。?计算X、Y轴在-1到1范围之内，大小为15×15的等间距网格上各点的高度。注意所得到的二维数组fvals的第0轴与Y轴对应，第一轴与X轴对应。?使用网格上各点的X、Y和Z轴的坐标创建interp2d对象，这里我们使用二阶插值曲面。?interp2d对象可以像函数一样调用，我们用它计算插值曲面在一个更密网格中的高度值。注意这里的参数是两个一维数组，分别指定网格的X-Y轴坐标，而不需要通过mgrid创建网格坐标数组。

1．griddata

interp2d类只能对网格形状的取样值进行插值运算，如果需要对随机散列的取样点进行插值，则可以使用griddata( )。其调用形式如下：

griddata(points, values, xi, method='linear', fill_value=nan)

其中points表示K维空间中的坐标，它可以是形状为(N, k)的数组，也可以是一个有k个数组的序列，N为数据的点数。values是points中每个点对应的值。xi是需要进行插值运算的坐标，其形状为(M, k)。method有三个选项——'nearest'、'linear'、'cubic'，分别对应0阶、1阶以及3阶插值。

下面是griddata( )的演示程序，其输出如图7所示。左图与'nearest'算法对应，平面上每个点都被填充为与它最近的采样点的数据，因此图中由许多相同颜色的色块构成。'linear'和'cubic'算法只对采样点构成的凸包区域进行插值，区域之外采用fill_value进行填充。中图和右图中的白色区域为插值的凸包区域之外。

griddata( )使用欧几里得距离计算插值。如果K维空间中每个维度的取值范围相差较大，则应先将数据正规化，然后使用griddata( )进行插值运算。

# 计算随机N个点的坐标，以及这些点对应的函数值
N = 200
np.random.seed(42)
x = np.random.uniform(-1, 1, N)
y = np.random.uniform(-1, 1, N)
z = func(x, y)

yg, xg = np.mgrid[-1:1:100j, -1:1:100j]
xi = np.c_[xg.ravel( ), yg.ravel( )]

methods = 'nearest', 'linear', 'cubic'

zgs = [interpolate.griddata((x, y), z, xi, method=method).reshape(100, 100)
for method in methods]

?

图7　使用gridata进行二维插值

2．径向基函数插值

径向基函数(radial basis function)插值算法也可以用于高维随机散布点的插值。所谓径向基函数，是指函数值只与某特定点的距离相关的一类函数，其中xi是某个给定取样点的坐标。使用这些?函数，可以近似表示N维空间中的函数：

为了方便读者理解RBF，下面先看一个一维插值的例子，结果如图8所示。图中显示了三种?函数对应的插值曲线：multiquadric、gaussian和linear。?

from scipy.interpolate import Rbf

x1 = np.array([-1, 0, 2.0, 1.0])
y1 = np.array([1.0, 0.3, -0.5, 0.8])

funcs = ['multiquadric', 'gaussian', 'linear']
nx = np.linspace(-3, 4, 100)
rbfs = [Rbf(x1, y1, function=fname) for fname in funcs] ?
rbf_ys = [rbf(nx) for rbf in rbfs] ?

图8　一维RBF插值

?使用表示取样点的x和y数组创建一个rbf对象，并通过function参数指定所使用的径向基函数。?rbf对象可以像函数一样被调用，我们用它计算以nx为横坐标对应的插值曲线的值。

rbf对象的nodes属性保存wi系数：

for fname, rbf in zip(funcs, rbfs):
print fname, rbf.nodes
multiquadric [ -3.79570791   9.82703701   5.08190777 -11.13103777]
gaussian [ 1.78016841 -1.83986382 -1.69565607  2.5266374 ]
linear [-0.26666667  0.6         0.73333333 -0.9       ]

下面的程序演示二维径向基函数插值，效果如图9所示。

rbfs = [Rbf(x, y, z, function=fname) for fname in funcs]
rbf_zg = [rbf(xg, yg).reshape(xg.shape) for rbf in rbfs]

图9　二维径向基函数插值

某些径向基函数可以通过epsilon参数指定其作用范围，该值越大每个插值点的作用范围越广，所得到的曲面也就越平滑。下面的代码演示gaussian径向基函数的epsilon参数与插值结果的关系，效果如图10所示。

epsilons = 0.1, 0.15, 0.3
rbfs = [Rbf(x, y, z, function="gaussian", epsilon=eps) for eps in epsilons]
zgs = [rbf(xg, yg).reshape(xg.shape) for rbf in rbfs]

?

图10　epsilon参数指定径向基函数中数据点的作用范围?