Scipy bisplrep双样条插值如何实现数据点区域外的外推？

问题背景 #

使用scipy.interpolate.bisplrep拟合2D点到1D轴时，整体效果优于多项式插值，但在拟合数据区域外会返回有界值，且该方法没有类似interp1d中fill_value="extrapolate"的参数。由于样条本质是分段多项式，可通过调用边界处的多项式系数实现外推。

解决方案思路 #

bisplrep返回的样条对象是元组tck = (tx, ty, c, kx, ky)，其中：

• tx/ty：x/y方向的节点数组
• c：样条系数矩阵
• kx/ky：x/y方向的样条阶数

实现外推的核心步骤：

1. 判断待计算点是否在插值区域外
2. 若在区域外，找到最近的边界样段
3. 提取该样段的系数，构造局部二元多项式计算外推值

代码实现 #

import scipy.interpolate as interpolate
import numpy as np

# 原始数据
xy = [[0, 0], [0, 1], [1, 0], [1, 1], [2, 2], [0, 2], [2, 0], [2, 1], [0.1, 0]]
x, y = tuple(list(zip(*xy)))
z =  [0, 1, 1, 1.1, 3, 2, 2, 3, 0.1]

# 生成二阶二元样条
spline = interpolate.bisplrep(x, y, z, kx=2, ky=2)
tx, ty, c, kx, ky = spline

def bisplev_extrapolate(xi, yi, tck):
    tx, ty, c, kx, ky = tck
    nx = len(tx) - kx - 1
    ny = len(ty) - ky - 1
    
    # 确定x方向的样段索引
    if xi <= tx[kx]:
        x_idx = kx
    elif xi >= tx[-kx-1]:
        x_idx = len(tx) - kx - 2
    else:
        # 区域内直接用原生bisplev计算
        return interpolate.bisplev([xi], [yi], tck)[0]
    
    # 确定y方向的样段索引
    if yi <= ty[ky]:
        y_idx = ky
    elif yi >= ty[-ky-1]:
        y_idx = len(ty) - ky - 2
    else:
        return interpolate.bisplev([xi], [yi], tck)[0]
    
    # 提取当前样段的系数矩阵
    coeff = c[y_idx - ky : y_idx + 1, x_idx - kx : x_idx + 1]
    
    # 计算点相对于样段起始节点的偏移量
    x_offset = xi - tx[x_idx]
    y_offset = yi - ty[y_idx]
    
    # 计算二元多项式值
    val = 0.0
    for i in range(kx + 1):
        for j in range(ky + 1):
            val += coeff[j, i] * (x_offset ** i) * (y_offset ** j)
    return val

# 测试用例
print(bisplev_extrapolate(0, 0, spline))  # 接近0，符合预期
print(bisplev_extrapolate(1, 1, spline))  # 输出1.1，符合预期
print(bisplev_extrapolate(2, 2, spline))  # 输出3，符合预期
print(bisplev_extrapolate(3, 3, spline))  # 输出外推值，不再返回3

说明 #

• 函数bisplev_extrapolate会自动判断点的位置：在插值区域内时直接调用原生bisplev，区域外则使用边界样段的多项式系数计算外推值
• 该方法严格遵循样条的分段多项式特性，外推逻辑与样条在边界内的插值逻辑一致

内容的提问来源于stack exchange，提问作者PJ127