基于MATLAB的径向基函数插值（RBF插值）（一维、二维、三维） |
您所在的位置：网站首页 › matlab中的插值方法 › 基于MATLAB的径向基函数插值（RBF插值）（一维、二维、三维） |
基于MATLAB的径向基函数插值（RBF插值）（一维、二维、三维）
|
基于MATLAB的径向基函数插值（RBF插值）（一维、二维、三维）
0 前言1 RBF思路2 1维RBF函数2.1 参数说明2.1.1 核函数选择2.1.2 作用半径2.1.3 多项式拟合2.1.4 误差项（光滑项）
3 2维RBF函数4 3维RBF函数
惯例声明：本人没有相关的工程应用经验，只是纯粹对相关算法感兴趣才写此博客。所以如果有错误，欢迎在评论区指正，不胜感激。本文主要关注于算法的实现，对于实际应用等问题本人没有任何经验，所以也不再涉及。 0 前言插值是一个工程中非常常见的扩展数据方法。通常数据测量数量永远是已知的，数据的储存空间也是有限的，但工程中的数据需求却永远是无限的。 如何用较少位置处的数据点来推广到任意位置处的数据点，是工程中常见的问题。其中径向基函数RBF插值具有不依赖数据网格的特点，省去了传统插值的网格剖分，是一种基于拟合的插值。 本文主要注重于具有几何意义上的RBF插值，所以只列举了一维二维和三维插值，更高维插值数据可以稍加改写代码就可以实现，本文也不再涉及。 本文的参考文献如下： [1]Meshfree Approximation Methods with MATLAB.Gregory E. Fasshauer. 1 RBF思路径向基函数的大概原理是利用一系列函数叠加，对原函数进行拟合。也就是: F ( x ) = ∑ w i ∗ f i ( x 0 , x ) F(x)=\sum w_i*f_i(x_0,x) F(x)=∑wi∗fi(x0,x) 其中f(x0,x)为基函数，是中心对称函数，函数中心点在x0上。w为每个函数的权重。 因此，只需要求出权重w，就可以利用基函数在各个点的值，计算出整个域的函数。而求解权重，也是一个简单的线性代数问题，直接用线性方程组求逆的方式就可以得到。 因此，RBF方法也具有原理简单，编程容易的特点。 下面用一个简单的例子来解释。 首先问题假设如下，我有下面6个点的数据（xi,yi），想插值出蓝色的曲线结果，该如何处理？ 如果已知的插值点更多，还可以拟合的更好。下图为10个已知点进行插值的结果，和预想曲线几乎重合。 上面代码如下： %RBF基本原理 clear clc close all %插值点 x0=0.2:0.5:3; y0=sin(2*pi*0.5*x0)+cos(2*pi*0.8*x0)+2; %原函数 x2=0:0.02:3; y2=sin(2*pi*0.5*x2)+cos(2*pi*0.8*x2)+2; figure() hold on plot(x2,y2) plot(x0,y0,'o','MarkerSize',8) hold off %1构造函数 N_k=length(x0);%构造函数的数量 RBF_Kernel=cell(N_k,1); for k=1:N_k x_k=x0(k);%中心位置，插值节点 RBF_Kernel_k=@(x) exp(-(x-x_k).^2/2/0.3^2);%正态函数 RBF_Kernel{k}=RBF_Kernel_k;%储存每个函数 end %2计算出插值矩阵 InterpMat=zeros(N_k,N_k); for k=1:N_k RBF_phi_k=RBF_Kernel{k}(x0);%计算出当前正态函数 InterpMat(:,k)=RBF_phi_k(:);%插值矩阵每一列储存的都是对应的正态函数 end %3利用插值矩阵求解出每个高斯函数的权重 %∑φ(k)*w(k)=InterpMat*w=y0，可以线性求逆直接得到系数w w=InterpMat\y0'; %绘制出每个正态函数 figure() hold on mcp=colormap('lines'); for k=1:N_k RBF_phi2_k=RBF_Kernel{k}(x2)*w(k); plot(x2,RBF_phi2_k,'Color',mcp(k,:)) plot([x0(k),x0(k)],[0,RBF_Kernel{k}(x0(k))*w(k)],'Color',mcp(k,:),'LineStyle','--') end plot(x2,y2,'Color','k','LineWidth',2) hold off %绘制出最终相加的结果 y_RBF=zeros(size(y2)); for k=1:N_k y_RBF=y_RBF+RBF_Kernel{k}(x2)*w(k);%叠加每一个正态函数 end figure() hold on plot(x2,y2) plot(x2,y_RBF,'--') plot(x0,y0,'o','MarkerSize',8) legend({'原函数','RBF插值结果'},'Location','best') 2 1维RBF函数下面为1维RBF函数插值的代码，基本思路和上面第一章节的一样。但是对于计算速度和输入输出等方面简单做了一些优化。 计算结果如下： 代码如下： clear clc close all %初始已知点 x0=0:0.1:3; y0=sin(2*pi*0.5*x0)+cos(2*pi*0.8*x0)+2+0.02*rand(size(x0)); x=(-1:0.01:4); y=RBF1(x0,y0,x,'gaussian',0.3,0,0.001); figure() hold on plot(x0,y0,'o') plot(x,y) hold off function y=RBF1(x0,y0,x,Method,Rs,Npoly,Error) %一维RBF插值，输入x0散点，y0值。x是输出插值得到的点。 %Method方法，默认'linear'。 %'linear',|R| %'gaussian',exp(-(R/Rs)^2) %'thin_plate',R^2*log(R) %'linearEpsR',|R|。分段线性插值，要求s更大一些 %'cubic',|R^3| %Rs，插值核作用半径。Rs对于'linear'，'cubic'，'thin_plate'无影响。Rs大概和点和点之间的距离差不多就行。 %Npoly多项式拟合。默认是1，只拟合1次项。 %Error误差。在(-∞,∞)区间。默认是0，无误差。表示可以在部分误差范围内去插值。 %Error一般在0~1之内。如果只是为了收敛，可以比较小，在0.1以内就可以有很好效果；如果想实现平滑，可适当增大，大于1。 %示例：x0=0:0.3:pi;y0=sin(2*x0)+0.2*x0;x=-1:0.01:4;y=RBF1(x0,y0,x,'linear',0.2,1,0.0); %示例：x0=0:0.1:3;y0=0.25*x0.^2-0.5*x0+0.1*rand(size(x0));x=-1:0.01:4;y=RBF1(x0,y0,x,'gaussian',0.3,2,0.01); N=length(x0);%点的数量，也是RBF核的数量 %整理为列向量 x0=x0(:); y0=y0(:); x=x(:); %函数默认信息 if nargin'常数项','一次项'},'location','best')
当数据采集具有一定的误差，或者计算不需要完全贴合每一个点的数据，或者拟合出的曲线太曲折需要光滑，或者计算出的系数w过大甚至不收敛，都可以适当的添加误差项来解决。 其中系数w不收敛的问题，可能是由于核半径Rs过大，或者原始数据本身不太好导致的。适当添加一点误差，1e-3甚至更小，就可以解决这个问题。 而如果想要消除误差，光滑曲线，则需要较大的数，比如0.2，甚至超过1。这个根据具体效果来看。 这个的原理是直接在矩阵K的对角线上减去一个数。这样可以略微放宽求解的约束，使得节点处的数值不一定是控制点数值。 下面展示了添加了随机项后，大误差和小误差项的对比代码和结果： x0=0:0.1:3; y0=sin(2*pi*0.5*x0)+cos(2*pi*0.8*x0)+1+0.4*randn(size(x0)); x=(-0.0:0.01:3.0); y1=RBF1(x0,y0,x,'gaussian',0.15,0,0.0001); y2=RBF1(x0,y0,x,'gaussian',0.15,0,0.3); figure() hold on plot(x,y1) plot(x,y2) plot(x0,y0,'o') hold off legend({'小误差项','大误差项'},'location','best')
二维RBF函数的原理和一维相同。 z=RBF2(x0,y0,z0,x,y,Method,Rs,Npoly,ExtrapMethod,Error)其中x0、y0是已知散点，z0是对应的值。然后要插值出x、y对应的值。 二维RBF函数额外考虑了外插的方法。通常并不建议数据外插，因为缺乏足够的信息。 本函数总共考虑了三种外插方法： 第一个是’rbf’，其实就是RBF方法算出来多少就是多少。 第二个是’nearest’，就是不外插，超出包线的点的数值直接等于相邻数据点数值。 第三个是’ploy’，是多项式外插，超出包线的点，按照多项式拟合的值计算得到，适用于波动不大的数据。 推荐’rbf’方法，因为不需要计算包线，如果不怎么考虑外插的话，非常节省时间。 下面展示了二维RBF插值计算的结果： 展示的代码如下： clear clc close all %初始节点 N=60; x0=rand(N,1)*6-3;%-3~3 y0=rand(N,1)*6-3;%-3~3 z0=peaks(x0,y0)+0.05*x0; [x,y]=meshgrid(-5:0.02:5); %二维RBF插值 z3=RBF2(x0,y0,z0,x,y,'gaussian',0.5,1,'rbf',0.001); z=zeros(size(x));z(:)=z3(:); figure() hold on surf(x,y,z,'EdgeColor','none'); scatter3(x0,y0,z0) hold off function z=RBF2(x0,y0,z0,x,y,Method,Rs,Npoly,ExtrapMethod,Error) %二维RBF插值，输入x0，y0散点，z0值。x，y是输出插值得到的点。 %Method方法，默认'linear'。 %'linear',|R| %'gaussian',exp(-(R/Rs)^2) %'thin_plate',R^2*log(R) %'linearEpsR',|R|。分段线性插值，要求s更大一些 %'cubic',|R^3| %Rs，插值核作用半径。Rs对于'linear'，'cubic'，'thin_plate'无影响。Rs大概和点和点之间的距离差不多就行。 %Npoly多项式拟合。默认是1，只拟合1次项。 %ExtrapMethod,外插方法。某个数，全部赋值为这个数。'nearest'，按照临近值外插。'ploy'，多项式外插。'rbf'，默认用RBF插值得到的值。 %Error误差。在(-∞,∞)区间。默认是0，无误差。表示可以在部分误差范围内去插值。 %Error一般在0~1之内。如果只是为了收敛，可以比较小，在0.1以内就可以有很好效果；如果想实现平滑，可适当增大，大于1。 %示例： N=numel(x0);%点的数量，也是RBF核的数量 %整理为列向量 x0=x0(:); y0=y0(:); z0=z0(:); x=x(:); y=y(:); %计算包线 [k_ch,V]=convhull(x0,y0); if V |
那么第一步，我们构造核函数。这里选用高斯函数作为核函数： 
总共构造6个核函数f(x,xi)，每个核函数中心点在xi上。
f
i
(
x
)
=
e
x
p
(
−
(
x
−
x
i
)
2
/
2
/
s
2
)
f_i(x)=exp(-(x-x_i)^2/2/s^2)
fi(x)=exp(−(x−xi)2/2/s2) 每个函数乘以权重，可以得到当前基函数对应各个点的数值。以第一个基函数为例：
w
1
∗
[
f
1
(
x
1
)
f
1
(
x
2
)
f
1
(
x
3
)
f
1
(
x
4
)
f
1
(
x
5
)
f
1
(
x
6
)
]
w_1* \begin{bmatrix} f_1(x_1)\\ f_1(x_2)\\ f_1(x_3)\\ f_1(x_4)\\ f_1(x_5)\\ f_1(x_6)\\ \end{bmatrix}
w1∗
f1(x1)f1(x2)f1(x3)f1(x4)f1(x5)f1(x6)
则原函数F(x)可以计算为：
F
(
x
)
=
∑
w
i
∗
f
i
(
x
)
F(x)=\sum w_i*f_i(x)
F(x)=∑wi∗fi(x)
=
∑
w
i
∗
[
f
i
(
x
1
)
f
i
(
x
2
)
f
i
(
x
3
)
f
i
(
x
4
)
f
i
(
x
5
)
f
i
(
x
6
)
]
=\sum w_i*\begin{bmatrix} f_i(x_1)\\ f_i(x_2)\\ f_i(x_3)\\ f_i(x_4)\\ f_i(x_5)\\ f_i(x_6)\\ \end{bmatrix}
=∑wi∗
fi(x1)fi(x2)fi(x3)fi(x4)fi(x5)fi(x6)
=
[
f
1
(
x
1
)
f
2
(
x
1
)
f
3
(
x
1
)
f
4
(
x
1
)
f
5
(
x
1
)
f
6
(
x
1
)
f
1
(
x
2
)
f
2
(
x
2
)
f
3
(
x
2
)
f
4
(
x
2
)
f
5
(
x
2
)
f
6
(
x
2
)
f
1
(
x
3
)
f
2
(
x
3
)
f
3
(
x
3
)
f
4
(
x
3
)
f
5
(
x
3
)
f
6
(
x
3
)
f
1
(
x
4
)
f
2
(
x
4
)
f
3
(
x
4
)
f
4
(
x
4
)
f
5
(
x
4
)
f
6
(
x
4
)
f
1
(
x
5
)
f
2
(
x
5
)
f
3
(
x
5
)
f
4
(
x
5
)
f
5
(
x
5
)
f
6
(
x
5
)
f
1
(
x
6
)
f
2
(
x
6
)
f
3
(
x
6
)
f
4
(
x
6
)
f
5
(
x
6
)
f
6
(
x
6
)
]
∗
[
w
1
w
2
w
3
w
4
w
5
w
6
]
= \begin{bmatrix} f_1(x_1)&f_2(x_1)&f_3(x_1)&f_4(x_1)&f_5(x_1)&f_6(x_1)\\ f_1(x_2)&f_2(x_2)&f_3(x_2)&f_4(x_2)&f_5(x_2)&f_6(x_2)\\ f_1(x_3)&f_2(x_3)&f_3(x_3)&f_4(x_3)&f_5(x_3)&f_6(x_3)\\ f_1(x_4)&f_2(x_4)&f_3(x_4)&f_4(x_4)&f_5(x_4)&f_6(x_4)\\ f_1(x_5)&f_2(x_5)&f_3(x_5)&f_4(x_5)&f_5(x_5)&f_6(x_5)\\ f_1(x_6)&f_2(x_6)&f_3(x_6)&f_4(x_6)&f_5(x_6)&f_6(x_6)\\ \end{bmatrix}* \begin{bmatrix} w_1\\ w_2\\ w_3\\ w_4\\ w_5\\ w_6\\ \end{bmatrix}
=
f1(x1)f1(x2)f1(x3)f1(x4)f1(x5)f1(x6)f2(x1)f2(x2)f2(x3)f2(x4)f2(x5)f2(x6)f3(x1)f3(x2)f3(x3)f3(x4)f3(x5)f3(x6)f4(x1)f4(x2)f4(x3)f4(x4)f4(x5)f4(x6)f5(x1)f5(x2)f5(x3)f5(x4)f5(x5)f5(x6)f6(x1)f6(x2)f6(x3)f6(x4)f6(x5)f6(x6)
∗
w1w2w3w4w5w6
然后利用线性代数方式，就可以得到系数w。 则原函数F可以利用这些个基函数叠加得到： 
叠加后的函数和原函数对比见下图： 
可以看到计算结果还可以，曲线过渡也比较光滑。
可以看到小误差项，曲线经过了每一个散点。大误差项，曲线只是大致经过了散点。
可以看到RBF插值可以较好的拟合出预期的效果。【本文地址】
今日新闻 |
推荐新闻 |