📚 课程目录第3章 RBF插值原理

3. RBF插值原理：从散乱数据点到连续曲面的映射

好，咱们进入正题。前面两章聊了RBF是什么、能干什么，这一章咱们要动真格的了——搞清楚RBF到底怎么把一堆散乱的点，变成一张光滑连续的曲面。

说白了，RBF插值就是干一件事：给你一堆已知点，让你猜未知点的值。听起来像数学版的“读心术”，对吧？

3.1 从散乱数据点到连续曲面的映射

先想想看，你手里有一堆测量数据点，比如汽车A柱表面的三维坐标 (x, y, z)。这些点分布得乱七八糟，没有规则网格。你要的是整个A柱的曲面，而不是几个孤零零的点。

RBF的做法很聪明：用一组“基函数”去拟合这些点。每个已知点都像一个“小磁铁”，对周围区域产生影响力。这个影响力的大小，由径向基函数 φ(r) 来描述，其中 r 是到该点的距离。

我记得第一次做这个的时候，脑子里总想着“这不就是加权平均吗？”后来发现，比加权平均复杂多了。每个基函数都有自己的形状参数，相互之间还会叠加影响。

核心思想：整个曲面可以表示为所有基函数的线性组合：

s(x) = Σ w_i · φ(||x - x_i||) + p(x)

其中 w_i 是权重系数，φ 是径向基函数，p(x) 是多项式项（用于保证刚体运动不变性）。

这里有个关键点：基函数只依赖于距离，不依赖于方向。这就是“径向”二字的由来。你想想看，不管从哪个方向靠近一个已知点，它的影响都是一样的——各向同性。

实际项目中，我遇到过一个问题：如果数据点分布极不均匀，有些区域点很密，有些区域点很稀。这时候基函数的形状参数就得小心调，否则稀疏区域会出现“波浪”，看着就不对劲。

3.2 插值条件与线性方程组构建

好，现在咱们有了模型形式，但权重系数 w_i 还不知道。怎么求？靠插值条件。

插值条件很简单：在已知数据点处，曲面值必须等于已知值。 也就是说：

s(x_j) = f_j，对于 j = 1, 2, ..., N

把模型代入，就得到：

Σ w_i · φ(||x_j - x_i||) + p(x_j) = f_j

这看起来像什么？没错，一个线性方程组！

写成矩阵形式就是：

[ Φ   P ] [ w ]   [ f ]
[ P^T   0 ] [ c ] = [ 0 ]

其中：

Φ 是 N×N 矩阵，元素 Φ_ji = φ(||x_j - x_i||)
P 是 N×(m+1) 矩阵，包含多项式基函数在数据点处的值
w 是权重系数向量
c 是多项式系数向量
f 是已知函数值向量

我的经验：这个矩阵通常是对称的，但不一定是正定的。如果基函数选得好（比如高斯函数、多二次函数），矩阵就是正定的，求解起来稳定得多。

这里有个坑，我踩过：矩阵条件数。当数据点很多（比如几千个）或者点之间距离很近时，矩阵会变得病态。我曾经用高斯函数做插值，形状参数设得太小，结果矩阵几乎奇异，求解出来的曲面全是毛刺。

3.3 系数求解方法

方程组建好了，怎么解？这取决于你的数据规模。

数据规模	推荐方法	特点
N < 1000	直接求解（LU分解、Cholesky分解）	精度高，速度快，内存占用小
1000 ≤ N < 10000	迭代法（共轭梯度法、GMRES）	内存可控，适合中等规模
N ≥ 10000	快速多极方法（FMM）或局部RBF	避免全局矩阵，但精度有损失

我个人习惯，对于中小规模问题（N < 5000），直接用Cholesky分解。因为矩阵是对称正定的，Cholesky比LU快一倍，而且数值稳定性好。

代码实现大概是这样的：

import numpy as np
from scipy.linalg import cho_factor, cho_solve

def rbf_interpolate(x_data, f_data, x_query, epsilon=1.0):
    """
    RBF插值求解
    x_data: 已知数据点 (N, d)
    f_data: 已知函数值 (N,)
    x_query: 查询点 (M, d)
    epsilon: 形状参数
    """
    N = x_data.shape[0]
    
    # 构建矩阵 Phi
    Phi = np.zeros((N, N))
    for i in range(N):
        for j in range(N):
            r = np.linalg.norm(x_data[i] - x_data[j])
            Phi[i, j] = np.exp(-(epsilon * r)**2)  # 高斯函数
    
    # Cholesky分解
    L, lower = cho_factor(Phi)
    
    # 求解权重系数
    w = cho_solve((L, lower), f_data)
    
    # 在查询点处插值
    f_query = np.zeros(x_query.shape[0])
    for i in range(x_query.shape[0]):
        for j in range(N):
            r = np.linalg.norm(x_query[i] - x_data[j])
            f_query[i] += w[j] * np.exp(-(epsilon * r)**2)
    
    return f_query

注意：上面这个双重循环在数据量大时效率极低。实际项目中一定要用向量化或者FMM加速。我曾经用这个代码处理5000个点，跑了整整两分钟——后来改成向量化，0.2秒就搞定了。

对于大规模问题，我建议用局部RBF方法。说白了，就是每个查询点只考虑附近K个已知点的影响，而不是全部。这样矩阵规模从 N×N 降到 K×K，求解速度飞起。

但代价是什么？精度。局部方法在数据稀疏区域可能产生不连续。我做过一个实验：用全局RBF和局部RBF分别拟合同一个自由曲面，全局方法最大误差0.01mm，局部方法（K=30）最大误差0.08mm。看你的精度要求了。

嗯，这里还要提一句：形状参数 epsilon 的选择。这个参数直接影响插值质量。太小了，矩阵病态；太大了，曲面过于平滑，细节丢失。我一般用交叉验证来选，或者用经验公式 epsilon = 0.815 * d_avg，其中 d_avg 是数据点之间的平均距离。

最后总结一下：RBF插值的本质，就是把一个看似复杂的曲面拟合问题，转化成一个线性方程组求解问题。你只要选好基函数、建好矩阵、解好方程，剩下的就交给数学了。

当然，实际工程中还有很多细节要处理——边界条件、形状参数优化、大规模加速等等。这些咱们后面章节会逐一展开。