径向基函数插值

2024-01-08 17:34:45

一、径向基函数的定义

如果 ∣ ∣ x 1 ∣ ∣ = ∣ ∣ x 2 ∣ ∣ ||x_1||=||x_2|| ∣∣x1?∣∣=∣∣x2?∣∣,那么 ? ( x 1 ) = ? ( x 2 ) \phi(x_1)=\phi(x_2) ?(x1?)=?(x2?) 的函数 ? \phi ? 就是径向函数,即仅由 r = ∣ ∣ x ∣ ∣ r=||x|| r=∣∣x∣∣ 控制的函数(径向基函数是一个取值仅仅依赖于离原点距离的实值函数,或者还可以是到任意一点 c c c 的距离)。

给定一个一元函数 ? : R + → R \phi:R_+\rightarrow R ?R+?R,在定义域 x ∈ R d x\in R^d xRd 上,所有形如 ψ ( x ) = ? ( ∣ ∣ x ? c ∣ ∣ ) \psi(x)=\phi(||x-c||) ψ(x)=?(∣∣x?c∣∣) 及其线性组合张成的函数空间称为由函数 ? \phi ? 导出的径向基函数空间

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

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

  1. Kriging 方法的 Gauss 分布函数: ? ( r ) = e ? c 2 r 2 \phi(r)=e^{-c^2r^2} ?(r)=e?c2r2
  2. Kriging 方法的 Markoff 分布函数: ? ( r ) = e ? c ∣ r ∣ \phi(r)=e^{-c|r|} ?(r)=e?cr,及其他概率分布函数;
  3. Hardy 的 Multi-Quadric 函数: ? ( r ) = ( c 2 + r 2 ) β \phi(r)=(c^2+r^2)^\beta ?(r)=(c2+r2)β(其中 β \beta β 是正的实数);
  4. Hardy 的逆 Multi-Quadric 函数: ? ( r ) = ( c 2 + r 2 ) ? β \phi(r)=(c^2+r^2)^{-\beta} ?(r)=(c2+r2)?β(其中 β \beta β 是正的实数);
  5. Duchon 的薄板样条: d d d 为偶数时, ? ( r ) = r 2 k ? d log ? r \phi(r)=r^{2k-d}\log r ?(r)=r2k?dlogr d d d 为奇数时, ? ( r ) = r 2 k ? d \phi(r)=r^{2k-d} ?(r)=r2k?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。选取径向函数 ? : R + → R \phi:R_+\rightarrow R ?:R+?R 来构造函数系 { ? ( ∣ ∣ x ? x j ∣ ∣ ) } j = 1 n \{\phi(||x-x_j||)\}_{j=1}^n {?(∣∣x?xj?∣∣)}j=1n? 并寻找形如 S ( x ) = ∑ j = 1 n λ j ? ( ∣ ∣ x ? x j ∣ ∣ ) S(x)=\sum_{j=1}^n\lambda_j\phi(||x-x_j||) S(x)=j=1n?λj??(∣∣x?xj?∣∣) 的插值函数 S ( x ) S(x) S(x),使其满足条件 S ( x j ) = f j , j = 1 , ? ? , n S(x_j)=f_j,j=1,\cdots,n S(xj?)=fj?,j=1,?,n

为了方便,我们定义
{ f T = ( f 1 , f 2 , ? ? , f n ) ? T ( x ) = ( ? ( ∣ ∣ x ? x 1 ∣ ∣ , ? ( ∣ ∣ x ? x 2 ∣ ∣ , ? ? , ? ( ∣ ∣ x ? x n ∣ ∣ ) ) λ T = ( λ 1 , λ 2 , ? ? , λ n ) A = ( ? ( ∣ ∣ x j ? x k ∣ ∣ ) ) n × n \begin{cases} \pmb{f^T}=(f_1,f_2,\cdots,f_n)\\[2ex] \pmb{\phi^T}(x)=\Big(\phi(||x-x_1||,\phi(||x-x_2||,\cdots,\phi(||x-x_n||)\Big)\\[2ex] \pmb{\lambda^T}=(\lambda_1,\lambda_2,\cdots,\lambda_n)\\[2ex] \pmb{A}=\Big(\phi(||x_j-x_k||)\Big)_{n\times n} \end{cases} ? ? ??fT=(f1?,f2?,?,fn?)?T(x)=(?(∣∣x?x1?∣∣,?(∣∣x?x2?∣∣,?,?(∣∣x?xn?∣∣))λT=(λ1?,λ2?,?,λn?)A=(?(∣∣xj??xk?∣∣))n×n??

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

定理:函数 ? : R + → R \phi:R_+\rightarrow R ?:R+?R 是连续的, lim ? r → ∞ ? ( r ) = 0 \lim_{r\rightarrow\infty}\phi(r)=0 limr??(r)=0,那么对于 n n n 元的径向基函数插值总是存在唯一解的充分条件是:矩阵 A \pmb A A 是正定矩阵。

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

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

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

令径向基函数插值方程为
S ( x ) = ∑ j = 1 n λ j ? ( ∣ ∣ x ? x j ∣ ∣ ) S(x)=\sum_{j=1}^n\lambda_j\phi(||x-x_j||) S(x)=j=1n?λj??(∣∣x?xj?∣∣)

将已知点 ( x j , f j ) , j = 1 , ? ? , n (x_j,f_j),j=1,\cdots,n (xj?,fj?),j=1,?,n 代入方程,可得:
[ λ 1 λ 2 ? λ n ] [ ? ( ∣ ∣ x 1 ? x 1 ∣ ∣ ) ? ( ∣ ∣ x 2 ? x 1 ∣ ∣ ) ? ? ( ∣ ∣ x n ? x 1 ∣ ∣ ) ? ( ∣ ∣ x 1 ? x 2 ∣ ∣ ) ? ( ∣ ∣ x 2 ? x 2 ∣ ∣ ) ? ? ( ∣ ∣ x n ? x 2 ∣ ∣ ) ? ? ? ? ( ∣ ∣ x 1 ? x n ∣ ∣ ) ? ( ∣ ∣ x 2 ? x n ∣ ∣ ) ? ? ( ∣ ∣ x n ? x n ∣ ∣ ) ] = [ f 1 f 2 ? f n ] \left[ \begin{matrix} \lambda_1 & \lambda_2 & \cdots & \lambda_n\\ \end{matrix} \right] \left[ \begin{matrix} \phi(||x_1-x_1||) & \phi(||x_2-x_1||) & \cdots & \phi(||x_n-x_1||)\\ \phi(||x_1-x_2||) & \phi(||x_2-x_2||) & \cdots & \phi(||x_n-x_2||)\\ \vdots & \vdots & & \vdots\\ \phi(||x_1-x_n||) & \phi(||x_2-x_n||) & \cdots & \phi(||x_n-x_n||)\\ \end{matrix} \right]= \left[ \begin{matrix} f_1& f_2& \cdots& f_n\\ \end{matrix} \right] [λ1??λ2????λn??] ??(∣∣x1??x1?∣∣)?(∣∣x1??x2?∣∣)??(∣∣x1??xn?∣∣)??(∣∣x2??x1?∣∣)?(∣∣x2??x2?∣∣)??(∣∣x2??xn?∣∣)??????(∣∣xn??x1?∣∣)?(∣∣xn??x2?∣∣)??(∣∣xn??xn?∣∣)? ?=[f1??f2????fn??]

求解上述方程,可求出 λ 1 , λ 2 , ? ? , λ n \lambda_1,\lambda_2,\cdots,\lambda_n λ1?,λ2?,?,λn? 的值,从而求出插值曲线方程。插值曲面方程类似,将 x x x 替换成向量 X X X 即可。

具体应用到高斯函数,设高斯函数插值方程为
S ( x ) = ∑ j = 1 n λ j e ? α ∣ ∣ x ? x j ∣ ∣ 2 α > 0 S(x)=\sum_{j=1}^n\lambda_je^{-\alpha||x-x_j||^2}\quad\alpha>0 S(x)=j=1n?λj?e?α∣∣x?xj?2α>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 代码实现

文章来源:https://blog.csdn.net/m0_72748751/article/details/135453652
本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。