一、径向基函数的定义

如果 $||x_1||=||x_2||$，那么 $\varphi (x_1)=\varphi (x_2)$ 的函数 $\varphi$ 就是径向函数，即仅由 $r=||x||$ 控制的函数（径向基函数是一个取值仅仅依赖于离原点距离的实值函数，或者还可以是到任意一点 $c$ 的距离）。

给定一个一元函数 $\varphi ：R_{+}\to R$，在定义域 $x\in R^d$ 上，所有形如 $\psi (x)=\varphi (||x-c||)$ 及其线性组合张成的函数空间称为由函数 $\varphi$ 导出的径向基函数空间。

在一定的条件下，只要取 { x j } \{x_j\} {xj​} 两两不同， { ϕ ( x − x j ) } \{\phi(x-x_j)\} {ϕ(x−xj​)} 就是线性无关的，从而形成径向基函数空间中某子空间的一组基。当 { x j } \{x_j\} {xj​} 几乎充满 $R$ 时， { x j } \{x_j\} {xj​} 几乎充满 $R$ 时， { ϕ ( x − x j ) } \{\phi(x-x_j)\} {ϕ(x−xj​)} 及其线性组合可以逼近几乎任何函数。

各类文献中常用的径向基函数有：

1. Kriging 方法的 Gauss 分布函数：$\varphi (r)=e^{-c^2{r}^2}$
2. Kriging 方法的 Markoff 分布函数：$\varphi (r)=e^{-c|r|}$，及其他概率分布函数；
3. Hardy 的 Multi-Quadric 函数：$\varphi (r)=(c^2+r^2{)}^{\beta }$（其中 $\beta$ 是正的实数）；
4. Hardy 的逆 Multi-Quadric 函数：$\varphi (r)=(c^2+r^2{)}^{-\beta }$（其中 $\beta$ 是正的实数）；
5. Duchon 的薄板样条：$d$ 为偶数时，$\varphi (r)=r^{2k-d}\log r$；$d$ 为奇数时，$\varphi (r)=r^{2k-d}$；

二、径向基函数插值

定义：径向基函数插值是对于给定的多元散乱数据 { x j , f j } j = 1 n , x j ∈ R n , f j ∈ R , j = 1 , ⋯   , n \{x_j,f_j\}^n_{j=1},x_j\in R^n,f_j\in R,j=1,\cdots,n {xj​,fj​}j=1n​,xj​∈Rn,fj​∈R,j=1,⋯,n。选取径向函数 $\varphi :R_{+}\to R$ 来构造函数系 { ϕ ( ∣ ∣ x − x j ∣ ∣ ) } j = 1 n \{\phi(||x-x_j||)\}_{j=1}^n {ϕ(∣∣x−xj​∣∣)}j=1n​ 并寻找形如 $S(x)={\sum }_{j=1}^n{\lambda }_j\varphi (||x-x_j||)$ 的插值函数 $S(x)$，使其满足条件 $S(x_j)=f_j,j=1,\cdots ,n$。

为了方便，我们定义
$$\{\begin{array}{l}{f}^T=(f_1,f_2,\cdots ,f_n)\\ {\varphi }^T(x)=(\varphi (||x-x_1||,\varphi (||x-x_2||,\cdots ,\varphi (||x-x_n||))\\ {\lambda }^T=({\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_n)\\ A=(\varphi (||x_j-x_k||){)}_{n\times n}\end{array}$$

上述插值方程对任意两两不同的 $x_j$ 的数据 { x j , f j } \{x_j,f_j\} {xj​,fj​} 有解的充要条件是：对任意两两不同的 $x_j$，对称矩阵 $A$ 都非奇异。

定理：函数 $\varphi :R_{+}\to R$ 是连续的，${\lim }_{r\to \infty }\varphi (r)=0$，那么对于 $n$ 元的径向基函数插值总是存在唯一解的充分条件是：矩阵 $A$ 是正定矩阵。

上面提到的径向基函数中逆 Multi-Quadric 函数和 Gauss 函数在任意维空间上都是正定函数，因此插值是唯一的。

三、用高斯函数进行散乱数据的插值

对于数据量少的情况，径向基函数（尤其是高斯函数）插值的结果较令人满意，而且计算也比较简单。

令径向基函数插值方程为
$$S(x)=\sum _{j=1}^n{\lambda }_j\varphi (||x-x_j||)$$

将已知点 $(x_j,f_j),j=1,\cdots ,n$ 代入方程，可得：
$$\left[\begin{array}{cccc}{\lambda }_1& {\lambda }_2& \cdots & {\lambda }_n\end{array}\right]\left[\begin{array}{cccc}\varphi (||x_1-x_1||)& \varphi (||x_2-x_1||)& \cdots & \varphi (||x_n-x_1||)\\ \varphi (||x_1-x_2||)& \varphi (||x_2-x_2||)& \cdots & \varphi (||x_n-x_2||)\\ ⋮& ⋮& & ⋮\\ \varphi (||x_1-x_n||)& \varphi (||x_2-x_n||)& \cdots & \varphi (||x_n-x_n||)\end{array}\right]=\left[\begin{array}{cccc}{f}_1& f_2& \cdots & f_n\end{array}\right]$$

求解上述方程，可求出 ${\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_n$ 的值，从而求出插值曲线方程。插值曲面方程类似，将 $x$ 替换成向量 $X$ 即可。

具体应用到高斯函数，设高斯函数插值方程为
$$S(x)=\sum _{j=1}^n{\lambda }_j{e}^{-\alpha ||x-x_j|{|}^2}\alpha >0$$

其中，$\alpha$ 为形状调整参数，可根据散乱数据点分布特征选取，当数据点对应的函数值变化较大时，$\alpha$ 可取的稍大些；数据点对应的函数值变化较小时，$\alpha$ 可取得稍小些。

python 代码实现

import numpy as np
import matplotlib.pyplot as plt


def gen_data(x1, x2):
    # 用于生成插值节点和总数据点 x1,x2 分别为插值节点的横坐标构成的行向量，总数据点的横坐标构成的行向量
    y_sample = np.sin(np.pi * x1 / 2) + np.cos(np.pi * x1 / 3)  # 插值节点的函数值
    y_all = np.sin(np.pi * x2 / 2) + np.cos(np.pi * x2 / 3)  # 总数据点的函数值
    return y_sample, y_all


def kernel_interpolation(y_sample, x1, sig):
    # 求解插值函数中高斯基函数的系数
    gaussian_kernel = lambda x, c, h: np.exp(-h*(x - x[c]) ** 2)  # 高斯基函数
    num = len(y_sample)
    w = np.zeros(num)
    int_matrix = np.asmatrix(np.zeros((num, num)))
    for i in range(num):
        int_matrix[i, :] = gaussian_kernel(x1, i, sig)
    w = np.asmatrix(y_sample) * int_matrix.I
    w = w.T
    return w


def kernel_interpolation_rec(w, x1, x2, sig):
    gkernel = lambda x, xc, h: np.exp(-h*(x - xc) ** 2)  # 高斯基函数
    num = len(x2)
    y_rec = np.zeros(num)
    for i in range(num):
        for k in range(len(w)):
            y_rec[i] = y_rec[i] + w[k] * gkernel(x2[i], x1[k], sig)
    return y_rec


if __name__ == '__main__':
    snum = 20  # control point数量
    ratio = 20  # 总数据点数量：snum*ratio
    sig = 0.5  # 核函数宽度
    xs = -8
    xe = 8
    x1 = np.linspace(xs, xe, snum)
    x2 = np.linspace(xs, xe, (snum - 1) * ratio + 1)
    y_sample, y_all = gen_data(x1, x2)
    plt.figure(1)
    w = kernel_interpolation(y_sample, x1, sig)
    y_rec = kernel_interpolation_rec(w, x1, x2, sig)
    plt.plot(x2, y_rec, 'k')
    plt.plot(x2, y_all, 'r:')
    plt.ylabel('y')
    plt.xlabel('x')
    for i in range(len(x1)):
        plt.plot(x1[i], y_sample[i], 'go', markerfacecolor='none')

    plt.legend(labels=['reconstruction', 'original', 'control point'], loc='lower left')
    plt.title('kernel interpolation:$y=sin(\pi x/2)+cos(\pi x/3)$')
    plt.show()

运行结果

Matlab 代码实现

clc;clear;

%% 参数
snum = 20;  % 插值节点个数
ratio = 20;  % 总数据点个数：(snum-1)*ratio + 1
sig = 0.5;  % 形状控制参数

xs = -8;  % 起点
xe = 8;  % 终点

%% 生成样本点
x1 = linspace(xs,xe,snum);  % 生成插值点坐标
x2 = linspace(xs,xe,(snum-1)*ratio + 1);
[y_sample,y_all] = gen_data(x1,x2);

figure('Name','RBF插值方法')

%% 计算基函数技术lambda
w = kernel_interpolation(y_sample,x1,sig);

%% 重构曲线
y_rec = kernel_interpolation_rec(w,x1,x2,sig);

%% 绘图
plot(x2,y_rec,'--',x2,y_all,':','LineWidth',2)
hold on
for i = 1:1:length(x1)
    scatter(x1(i),y_sample(i),70,'green')
end


title('$y=sin(\pi x/2)+cos(\pi x/3)$','Interpreter','latex')
xlabel('x')
ylabel('y')
legend('重构曲线','原始曲线','插值点','Location','southwest')

%% 生成插值节点和总数据点
function [y_sample,y_all] = gen_data(x1,x2)
y_sample = sin(pi * x1 / 2) + cos(pi * x1 / 3);
y_all = sin(pi * x2 / 2) + cos(pi * x2 / 3);
end


%% 求解插值函数中高斯基函数的系数
function w = kernel_interpolation(y_sample,x1,sig)
num = length(y_sample);
int_matrix = zeros(num,num);

for i = 1:1:num 
    int_matrix(i,:) = exp(-sig.*(x1-x1(i)).^2);
end

w = y_sample * inv(int_matrix);
end


%% 重构曲线
function y_rec = kernel_interpolation_rec(w,x1,x2,sig)
num = length(x2);
y_rec = zeros(1,num);

for i = 1:1:num
    for k = 1:1:length(w)
        y_rec(i) = y_rec(i) + w(k) .* exp(-sig.*(x2(i)-x1(k)).^2);
    end
end
end

运行结果