慣例聲明：本人沒有相關(guān)的工程應(yīng)用經(jīng)驗(yàn)，只是純粹對相關(guān)算法感興趣才寫此博客。所以如果有錯誤，歡迎在評論區(qū)指正，不勝感激。本文主要關(guān)注于算法的實(shí)現(xiàn)，對于實(shí)際應(yīng)用等問題本人沒有任何經(jīng)驗(yàn)，所以也不再涉及。

0 前言

插值是一個工程中非常常見的擴(kuò)展數(shù)據(jù)方法。通常數(shù)據(jù)測量數(shù)量永遠(yuǎn)是已知的，數(shù)據(jù)的儲存空間也是有限的，但工程中的數(shù)據(jù)需求卻永遠(yuǎn)是無限的。

如何用較少位置處的數(shù)據(jù)點(diǎn)來推廣到任意位置處的數(shù)據(jù)點(diǎn)，是工程中常見的問題。其中徑向基函數(shù)RBF插值具有不依賴數(shù)據(jù)網(wǎng)格的特點(diǎn)，省去了傳統(tǒng)插值的網(wǎng)格剖分，是一種基于擬合的插值。

本文主要注重于具有幾何意義上的RBF插值，所以只列舉了一維二維和三維插值，更高維插值數(shù)據(jù)可以稍加改寫代碼就可以實(shí)現(xiàn)，本文也不再涉及。

本文的參考文獻(xiàn)如下：
[1]Meshfree Approximation Methods with MATLAB.Gregory E. Fasshauer.

1 RBF思路

徑向基函數(shù)的大概原理是利用一系列函數(shù)疊加，對原函數(shù)進(jìn)行擬合。也就是:
$$F(x)=\sum w_i?f_i(x_0,x)$$
其中f(x0,x)為基函數(shù)，是中心對稱函數(shù)，函數(shù)中心點(diǎn)在x0上。w為每個函數(shù)的權(quán)重。

因此，只需要求出權(quán)重w，就可以利用基函數(shù)在各個點(diǎn)的值，計(jì)算出整個域的函數(shù)。而求解權(quán)重，也是一個簡單的線性代數(shù)問題，直接用線性方程組求逆的方式就可以得到。

因此，RBF方法也具有原理簡單，編程容易的特點(diǎn)。

下面用一個簡單的例子來解釋。

首先問題假設(shè)如下，我有下面6個點(diǎn)的數(shù)據(jù)（xi,yi），想插值出藍(lán)色的曲線結(jié)果，該如何處理？

那么第一步，我們構(gòu)造核函數(shù)。這里選用高斯函數(shù)作為核函數(shù)：

總共構(gòu)造6個核函數(shù)f(x,xi)，每個核函數(shù)中心點(diǎn)在xi上。
$$f_i(x)=exp(?(x?x_i{)}^2/2/s^2)$$
每個函數(shù)乘以權(quán)重，可以得到當(dāng)前基函數(shù)對應(yīng)各個點(diǎn)的數(shù)值。以第一個基函數(shù)為例：
$$w_1?\left[\begin{array}{c}{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{array}\right]$$
則原函數(shù)F(x)可以計(jì)算為：
$$F(x)=\sum w_i?f_i(x)$$$$=\sum w_i?\left[\begin{array}{c}{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{array}\right]$$$$=\left[\begin{array}{cccccc}{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{array}\right]?\left[\begin{array}{c}{w}_1\\ w_2\\ w_3\\ w_4\\ w_5\\ w_6\end{array}\right]$$
然后利用線性代數(shù)方式，就可以得到系數(shù)w。
則原函數(shù)F可以利用這些個基函數(shù)疊加得到：

疊加后的函數(shù)和原函數(shù)對比見下圖：
可以看到計(jì)算結(jié)果還可以，曲線過渡也比較光滑。

如果已知的插值點(diǎn)更多，還可以擬合的更好。下圖為10個已知點(diǎn)進(jìn)行插值的結(jié)果，和預(yù)想曲線幾乎重合。

上面代碼如下：

%RBF基本原理
clear
clc
close all%插值點(diǎn)
x0=0.2:0.5:3;
y0=sin(2*pi*0.5*x0)+cos(2*pi*0.8*x0)+2;
%原函數(shù)
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構(gòu)造函數(shù)
N_k=length(x0);%構(gòu)造函數(shù)的數(shù)量
RBF_Kernel=cell(N_k,1);
for k=1:N_kx_k=x0(k);%中心位置，插值節(jié)點(diǎn)RBF_Kernel_k=@(x) exp(-(x-x_k).^2/2/0.3^2);%正態(tài)函數(shù)RBF_Kernel{k}=RBF_Kernel_k;%儲存每個函數(shù)
end%2計(jì)算出插值矩陣
InterpMat=zeros(N_k,N_k);
for k=1:N_kRBF_phi_k=RBF_Kernel{k}(x0);%計(jì)算出當(dāng)前正態(tài)函數(shù)InterpMat(:,k)=RBF_phi_k(:);%插值矩陣每一列儲存的都是對應(yīng)的正態(tài)函數(shù)
end%3利用插值矩陣求解出每個高斯函數(shù)的權(quán)重
%∑φ(k)*w(k)=InterpMat*w=y0，可以線性求逆直接得到系數(shù)w
w=InterpMat\y0';%繪制出每個正態(tài)函數(shù)
figure()
hold on
mcp=colormap('lines');
for k=1:N_kRBF_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%繪制出最終相加的結(jié)果
y_RBF=zeros(size(y2));
for k=1:N_ky_RBF=y_RBF+RBF_Kernel{k}(x2)*w(k);%疊加每一個正態(tài)函數(shù)
end
figure()
hold on
plot(x2,y2)
plot(x2,y_RBF,'--')
plot(x0,y0,'o','MarkerSize',8)
legend({'原函數(shù)','RBF插值結(jié)果'},'Location','best')

2 1維RBF函數(shù)

下面為1維RBF函數(shù)插值的代碼，基本思路和上面第一章節(jié)的一樣。但是對于計(jì)算速度和輸入輸出等方面簡單做了一些優(yōu)化。

計(jì)算結(jié)果如下：

代碼如下：

clear
clc
close all%初始已知點(diǎn)
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 offfunction y=RBF1(x0,y0,x,Method,Rs,Npoly,Error)
%一維RBF插值，輸入x0散點(diǎn)，y0值。x是輸出插值得到的點(diǎn)。
%Method方法，默認(rèn)'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大概和點(diǎn)和點(diǎn)之間的距離差不多就行。
%Npoly多項(xiàng)式擬合。默認(rèn)是1，只擬合1次項(xiàng)。
%Error誤差。在(-∞,∞)區(qū)間。默認(rèn)是0，無誤差。表示可以在部分誤差范圍內(nèi)去插值。
%Error一般在0~1之內(nèi)。如果只是為了收斂，可以比較小，在0.1以內(nèi)就可以有很好效果；如果想實(shí)現(xiàn)平滑，可適當(dāng)增大，大于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);%點(diǎn)的數(shù)量，也是RBF核的數(shù)量
%整理為列向量
x0=x0(:);
y0=y0(:);
x=x(:);%函數(shù)默認(rèn)信息
if nargin<4 || isempty(Method)Method='linear';%默認(rèn)線性核函數(shù)
elseif nargin<5 || isempty(Rs)Rs=1.1*(max(x0)-min(x0))/N;%假設(shè)空間均勻分布
elseif nargin<6 || isempty(Npoly)Npoly=1;%默認(rèn)擬合1次項(xiàng)
elseif nargin<7 || isempty(Error)Error=0;%默認(rèn)輸入數(shù)據(jù)沒有誤差
end%選擇核函數(shù)
K_method=0;
switch Methodcase 'linear'K_method=1;fun=@(RMat) Kernel_Linear(RMat,Rs);case 'gaussian'K_method=2;fun=@(RMat) Kernel_Gaussian(RMat,Rs);Error=-Error;%因?yàn)榱泓c(diǎn)不是零，遠(yuǎn)點(diǎn)是0，所以這里誤差為負(fù)case 'thin_plate'K_method=3;fun=@(RMat) Kernel_Thin_plate(RMat,Rs);case 'linearEpsR'K_method=4;fun=@(RMat) Kernel_LinearEpsR(RMat,Rs);Error=-Error;%因?yàn)榱泓c(diǎn)不是零，遠(yuǎn)點(diǎn)是0，所以這里誤差為負(fù)
end%將原始數(shù)據(jù)分離出多項(xiàng)式項(xiàng)Npoly
if Npoly==0C=ones(N,1);%常數(shù)項(xiàng)wC=mean(y0);
elseif Npoly==1C=[ones(N,1),x0];%常數(shù)項(xiàng)+一次項(xiàng)wC=C\y0;
elseif Npoly==2C=[ones(N,1),x0,x0.^2];%常數(shù)項(xiàng)+一次項(xiàng)+二次項(xiàng)wC=C\y0;
elseerror('只支持Npoly=0,1,2')
end
yC=C*wC;%多項(xiàng)式項(xiàng)
y1=y0-yC;%計(jì)算每個基函數(shù)在各個點(diǎn)的值
K2=zeros(N);%每一列對應(yīng)一個函數(shù)
%計(jì)算距離矩陣
DisMat=squareform(pdist(x0));
K3=fun(DisMat);%每一個核函數(shù)的中心點(diǎn)在節(jié)點(diǎn)上
K3=K3-Error*eye(N);%把誤差函數(shù)加入w=K3\y1;%計(jì)算權(quán)重%根據(jù)權(quán)重計(jì)算插值
y=zeros(size(x));
for k=1:N %計(jì)算每個核的貢獻(xiàn)，然后疊加R_k=abs(x-x0(k));yt=w(k)*fun(R_k );%plot(t,yt);y=y+yt;
endNOut=length(y);
%再還原回多項(xiàng)式項(xiàng)
if Npoly==0C=ones(NOut,1);%常數(shù)項(xiàng)
elseif Npoly==1C=[ones(NOut,1),x];%常數(shù)項(xiàng)+一次項(xiàng)
elseif Npoly==2C=[ones(NOut,1),x,x.^2];%常數(shù)項(xiàng)+一次項(xiàng)+二次項(xiàng)
end
y=y+C*wC;%再加上多項(xiàng)式項(xiàng)%檢查結(jié)果
if max(abs(w))/max(abs(y1))>1e3warning('結(jié)果未收斂，建議調(diào)整誤差Error');
elseif (max(y)-min(y))>5*(max(y0)-min(y0))warning('結(jié)果未收斂，建議增大區(qū)間Eps，或者調(diào)整誤差Error');
endendfunction z=Kernel_Linear(R,s)
%R距離
z=abs(R);%線性函數(shù)
endfunction z=Kernel_Gaussian(R,s)
%R距離
%s大概和采樣點(diǎn)間距差不多就行，可以略大（更胖）
z=exp(-R.^2/2/s^2);%正態(tài)函數(shù)
endfunction z=Kernel_Thin_plate(R,s)
%R距離
z=R.^2.*log(R);%thin_plate
endfunction z=Kernel_LinearEpsR(R,s)
%R距離
%s大概和采樣點(diǎn)間距差不多就行，可以略大（更胖）
s=2*s;%這里要更大
z=s-R;%線性函數(shù)
z(z<0)=0;
end

2.1 參數(shù)說明

下面說一下上面函數(shù)RBF1的后面各項(xiàng)參數(shù)結(jié)果：

y=RBF1(x0,y0,x,Method,Rs,Npoly,Error)

2.1.1 核函數(shù)選擇

RBF函數(shù)的核函數(shù)有非常多種，比如高斯函數(shù)、線性函數(shù)、薄板函數(shù)等。

下面以分別選擇高斯函數(shù)和線性函數(shù)作為核函數(shù)的RBF進(jìn)行對比：

clear
clc
close all%初始已知點(diǎn)
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

可以看到線性核函數(shù)的效果相當(dāng)于線性插值，而高斯核函數(shù)對峰值的擬合預(yù)測很好，各有側(cè)重點(diǎn)。

事實(shí)上，核函數(shù)的選擇對RBF插值有重要影響，不同核函數(shù)對應(yīng)著不同的場景。

2.1.2 作用半徑

對于某些核函數(shù)，還需要確定其作用半徑。比如對于高斯函數(shù)，半徑越大，核越寬。

對于高斯函數(shù)，作用半徑通常要大于兩個散點(diǎn)之間距離，但最好不要太大。如果太小，會導(dǎo)致每個點(diǎn)都是一個孤立的峰。如果太大，會導(dǎo)致求系數(shù)時(shí)方程求解困難。

下面代碼展示了作用半徑為0.05、0.3和0.8半徑的插值效果。

%初始節(jié)點(diǎn)
x0=0:0.3: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.01);
y3=RBF1(x0,y0,x,'gaussian',0.05,0,0.0001);
y4=RBF1(x0,y0,x,'gaussian',0.8,0,0.0001);
figure()
hold on
plot(x,y3)
plot(x,y)
plot(x,y4)
plot(x0,y0,'o')
hold off
legend({'Rs=0.05','Rs=0.3','Rs=0.7'},'location','best')
ylim([-0.5,4.5])

可以看到上面每個點(diǎn)的間距為0.2。當(dāng)Rs很小的時(shí)候，每個點(diǎn)只有一個孤立的高斯峰。當(dāng)Rs大于間距時(shí)，比如Rs=0.3和0.7，都可以求解出結(jié)果。

但當(dāng)Rs比較大，比如Rs=0.7時(shí)，會出現(xiàn)計(jì)算很奇怪的解，比如系數(shù)w特別大。這種如果適當(dāng)?shù)募右恍↑c(diǎn)誤差Error（有一點(diǎn)就行，1e-3量級的已經(jīng)足夠了），放寬求解精度，可以避免這種系數(shù)w很大不收斂的情況。

2.1.3 多項(xiàng)式擬合

本RBF插值函數(shù)默認(rèn)帶一個常數(shù)項(xiàng)。因?yàn)橄窀咚购瘮?shù)等函數(shù)無窮大處是0，對于本身數(shù)值在0附近的函數(shù)插值效果還可以，但是函數(shù)距離x軸很遠(yuǎn)時(shí)，再用這些基函數(shù)去擬合效果就會非常差。

因此在RBF插值之前，需要先求出常數(shù)項(xiàng)C，然后將所有散點(diǎn)的yi減去這個常數(shù)項(xiàng)C。之后計(jì)算完系數(shù)再插值后，再把這個常數(shù)項(xiàng)C加上。

一次項(xiàng)也是同樣的原理，就是先擬合出y’=kx+b，然后把所有散點(diǎn)yi-y’，得到計(jì)算系數(shù)用的y值。之后再把這個一次項(xiàng)加上即可。

通常常數(shù)項(xiàng)一般就夠，除非數(shù)據(jù)有很明顯的線性度。

下面代碼展示了常數(shù)項(xiàng)和一次項(xiàng)對比插值的效果：

x0=0:0.3:3;
x=(-1:0.01:4);
y0=sin(2*pi*0.5*x0)+cos(2*pi*0.8*x0)+1+2*x0+0.02*rand(size(x0));
y1=RBF1(x0,y0,x,'gaussian',0.3,0,0.0001);
y2=RBF1(x0,y0,x,'gaussian',0.3,1,0.0001);
figure()
hold on
plot(x,y1)
plot(x,y2)
plot(x0,y0,'o')
hold off
legend({'常數(shù)項(xiàng)','一次項(xiàng)'},'location','best')

2.1.4 誤差項(xiàng)（光滑項(xiàng)）

當(dāng)數(shù)據(jù)采集具有一定的誤差，或者計(jì)算不需要完全貼合每一個點(diǎn)的數(shù)據(jù)，或者擬合出的曲線太曲折需要光滑，或者計(jì)算出的系數(shù)w過大甚至不收斂，都可以適當(dāng)?shù)奶砑诱`差項(xiàng)來解決。

其中系數(shù)w不收斂的問題，可能是由于核半徑Rs過大，或者原始數(shù)據(jù)本身不太好導(dǎo)致的。適當(dāng)添加一點(diǎn)誤差，1e-3甚至更小，就可以解決這個問題。

而如果想要消除誤差，光滑曲線，則需要較大的數(shù)，比如0.2，甚至超過1。這個根據(jù)具體效果來看。

這個的原理是直接在矩陣K的對角線上減去一個數(shù)。這樣可以略微放寬求解的約束，使得節(jié)點(diǎn)處的數(shù)值不一定是控制點(diǎn)數(shù)值。

下面展示了添加了隨機(jī)項(xiàng)后，大誤差和小誤差項(xiàng)的對比代碼和結(jié)果：

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({'小誤差項(xiàng)','大誤差項(xiàng)'},'location','best')

可以看到小誤差項(xiàng)，曲線經(jīng)過了每一個散點(diǎn)。大誤差項(xiàng)，曲線只是大致經(jīng)過了散點(diǎn)。

3 2維RBF函數(shù)

二維RBF函數(shù)的原理和一維相同。

z=RBF2(x0,y0,z0,x,y,Method,Rs,Npoly,ExtrapMethod,Error)

其中x0、y0是已知散點(diǎn)，z0是對應(yīng)的值。然后要插值出x、y對應(yīng)的值。

二維RBF函數(shù)額外考慮了外插的方法。通常并不建議數(shù)據(jù)外插，因?yàn)槿狈ψ銐虻男畔ⅰ?/p>

本函數(shù)總共考慮了三種外插方法：
第一個是’rbf’，其實(shí)就是RBF方法算出來多少就是多少。
第二個是’nearest’，就是不外插，超出包線的點(diǎn)的數(shù)值直接等于相鄰數(shù)據(jù)點(diǎn)數(shù)值。
第三個是’ploy’，是多項(xiàng)式外插，超出包線的點(diǎn)，按照多項(xiàng)式擬合的值計(jì)算得到，適用于波動不大的數(shù)據(jù)。
推薦’rbf’方法，因?yàn)椴恍枰?jì)算包線，如果不怎么考慮外插的話，非常節(jié)省時(shí)間。

下面展示了二維RBF插值計(jì)算的結(jié)果：

可以看到RBF插值可以較好的擬合出預(yù)期的效果。

展示的代碼如下：

clear
clc
close all%初始節(jié)點(diǎn)
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 offfunction z=RBF2(x0,y0,z0,x,y,Method,Rs,Npoly,ExtrapMethod,Error)
%二維RBF插值，輸入x0，y0散點(diǎn)，z0值。x，y是輸出插值得到的點(diǎn)。
%Method方法，默認(rèn)'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大概和點(diǎn)和點(diǎn)之間的距離差不多就行。
%Npoly多項(xiàng)式擬合。默認(rèn)是1，只擬合1次項(xiàng)。
%ExtrapMethod,外插方法。某個數(shù)，全部賦值為這個數(shù)。'nearest'，按照臨近值外插。'ploy'，多項(xiàng)式外插。'rbf'，默認(rèn)用RBF插值得到的值。
%Error誤差。在(-∞,∞)區(qū)間。默認(rèn)是0，無誤差。表示可以在部分誤差范圍內(nèi)去插值。
%Error一般在0~1之內(nèi)。如果只是為了收斂，可以比較小，在0.1以內(nèi)就可以有很好效果；如果想實(shí)現(xiàn)平滑，可適當(dāng)增大，大于1。%示例：N=numel(x0);%點(diǎn)的數(shù)量，也是RBF核的數(shù)量
%整理為列向量
x0=x0(:);
y0=y0(:);
z0=z0(:);
x=x(:);
y=y(:);
%計(jì)算包線
[k_ch,V]=convhull(x0,y0);
if V<=0warning('輸入散點(diǎn)過于集中');
end
%函數(shù)默認(rèn)信息
narginchk(5,10)
if nargin<7 || isempty(Method)Method='linear';%默認(rèn)線性核函數(shù)
elseif nargin<8 || isempty(Rs)Rs=1.1*(max(x0)-min(x0))/N;%假設(shè)空間均勻分布
elseif nargin<9 || isempty(Npoly)Npoly=1;%默認(rèn)擬合1次項(xiàng)
elseif nargin<10 || isempty(Npoly)ExtrapMethod='rbf';%默認(rèn)不外插
elseif nargin<11 || isempty(Error)Error=0;%默認(rèn)輸入數(shù)據(jù)沒有誤差
end%選擇核函數(shù)
switch Methodcase 'linear'fun=@(RMat) Kernel_Linear(RMat,Rs);case 'gaussian'fun=@(RMat) Kernel_Gaussian(RMat,Rs);Error=-Error;%因?yàn)榱泓c(diǎn)不是零，遠(yuǎn)點(diǎn)是0，所以這里誤差為負(fù)case 'thin_plate'fun=@(RMat) Kernel_Thin_plate(RMat,Rs);case 'linearEpsR'fun=@(RMat) Kernel_LinearEpsR(RMat,Rs);Error=-Error;%因?yàn)榱泓c(diǎn)不是零，遠(yuǎn)點(diǎn)是0，所以這里誤差為負(fù)
end%將原始數(shù)據(jù)分離出多項(xiàng)式項(xiàng)Npoly
if Npoly==0C=ones(N,1);%常數(shù)項(xiàng)wC=mean(z0);
elseif Npoly==1if N<2;warning('點(diǎn)數(shù)太少，建議Npoly等于0');endC=[ones(N,1),x0,y0];%常數(shù)項(xiàng)+一次項(xiàng)wC=C\z0;
elseif Npoly==2if N<5;warning('點(diǎn)數(shù)太少，建議Npoly等于1');endC=[ones(N,1),x0,y0,x0.^2,y0.^2,x0.*y0];%常數(shù)項(xiàng)+一次項(xiàng)+二次項(xiàng)wC=C\z0;
elseerror('只支持Npoly=0,1,2')
end
zC=C*wC;%多項(xiàng)式項(xiàng)
z1=z0-zC;%計(jì)算距離矩陣
DisMat=squareform(pdist([x0,y0]));
K3=fun(DisMat);%每一個核函數(shù)的中心點(diǎn)在節(jié)點(diǎn)上
K3=K3-Error*eye(N);%把誤差函數(shù)加入w=K3\z1;%計(jì)算權(quán)重%根據(jù)權(quán)重計(jì)算插值
z=zeros(size(x));
for k=1:N %計(jì)算每個核的貢獻(xiàn)，然后疊加R_k=sqrt((x-x0(k)).^2+(y-y0(k)).^2);zt=w(k)*fun(R_k );z=z+zt;
endNOut=length(z);
%再還原回多項(xiàng)式項(xiàng)
if Npoly==0C=ones(NOut,1);%常數(shù)項(xiàng)
elseif Npoly==1C=[ones(NOut,1),x,y];%常數(shù)項(xiàng)+一次項(xiàng)
elseif Npoly==2C=[ones(NOut,1),x,y,x.^2,y.^2,x.*y];%常數(shù)項(xiàng)+一次項(xiàng)+二次項(xiàng)
end
zC2=C*wC;%再加上多項(xiàng)式項(xiàng)%判斷外插方法
Inpoly=inpolygon(x,y,x0(k_ch),y0(k_ch));%多邊形內(nèi)
indxInpoly=find(Inpoly);
if isnumeric(ExtrapMethod)z(Inpoly)=z(Inpoly)+zC1(Inpoly);%內(nèi)部的正常z(~Inpoly)=ExtrapMethod;%外部的固定值
elseif ischar(ExtrapMethod)switch ExtrapMethodcase 'rbf'z=z+zC2;%再加上多項(xiàng)式項(xiàng)case 'ploy'z(Inpoly)=z(Inpoly)+zC2(Inpoly);%內(nèi)部的正常z(~Inpoly)=zC2(~Inpoly);%外部的直接用多項(xiàng)式值case 'nearest'z(Inpoly)=z(Inpoly)+zC2(Inpoly);%內(nèi)部的正常%找到最近的點(diǎn)F = scatteredInterpolant(x(Inpoly),y(Inpoly),z(Inpoly),'nearest','nearest');%外部的直接插值z(~Inpoly)=F(x(~Inpoly),y(~Inpoly));end
elseerror('未識別ExtrapMethod')
end%檢查結(jié)果
if max(abs(w))/max(abs(z1))>1e3warning('結(jié)果未收斂，建議調(diào)整誤差Error');
elseif (max(z)-min(z))>5*(max(z0)-min(z0))warning('結(jié)果未收斂，建議增大區(qū)間Eps，或者調(diào)整誤差Error');
endendfunction z=Kernel_Linear(R,s)
%R距離
z=abs(R);%線性函數(shù)
endfunction z=Kernel_Gaussian(R,s)
%R距離
%s大概和采樣點(diǎn)間距差不多就行，可以略大（更胖）
z=exp(-R.^2/2/s^2);%正態(tài)函數(shù)
endfunction z=Kernel_Thin_plate(R,s)
%R距離
z=R.^2.*log(R);%thin_plate
endfunction z=Kernel_LinearEpsR(R,s)
%R距離
%s大概和采樣點(diǎn)間距差不多就行，可以略大（更胖）
s=2*s;%這里要更大
z=s-R;%線性函數(shù)
z(z<0)=0;
end

4 3維RBF函數(shù)

計(jì)算方法和二維的方法一致。其實(shí)更高維的數(shù)據(jù)計(jì)算方法也相差不多。
對于實(shí)際工程的幾何插值使用，最高維度也就是三維。對于其它應(yīng)用可能還是需要高維的插值，本文暫不涉及，由于RBF代碼較為簡單，因此稍加改動即可更改為高維RBF插值結(jié)果。

三維RBF函數(shù)的原理和二維也相同。

P=RBF3(x0,y0,z0,P0,x,y,z,Method,Rs,Npoly,ExtrapMethod,Error)

其中x0、y0、z0是已知散點(diǎn)，P0是對應(yīng)的值。然后要插值出x、y、z對應(yīng)的值。

下面展示了三維RBF插值計(jì)算的結(jié)果和對應(yīng)代碼：

clear
clc
close all%初始節(jié)點(diǎn)
N=1000;
x0=rand(N,1)*6-3;%-3~3
y0=rand(N,1)*6-3;%-3~3
z0=rand(N,1)*6-3;%-3~3
P0=sin(1*pi*x0)+cos(1.5*pi*y0)+sin(1*pi*z0)+1;[x,y,z]=meshgrid(-5:0.2:5);P3=RBF3(x0,y0,z0,P0,x,y,z,'linear',0.5,1,'nearest',0.01);%RBF插值P=zeros(size(x));P(:)=P3(:);%實(shí)際插值結(jié)果
figure()
hold on
s=slice(x,y,z,P,[0],[0],[0]);
set(s,'LineStyle','none');
scatter3(x0,y0,z0)
hold off
xlabel('x');ylabel('y');zlabel('z');
view(3)%理論目標(biāo)結(jié)果
figure()
P3=sin(1*pi*x)+cos(1.5*pi*y)+sin(1*pi*z)+1;
P3(x>3)=0;P3(y>3)=0;P3(z>3)=0;P3(x<-3)=0;P3(y<-3)=0;P3(z<-3)=0;
s=slice(x,y,z,P3,[0],[0],[0]);
xlabel('x');ylabel('y');zlabel('z');
view(3)function P=RBF3(x0,y0,z0,P0,x,y,z,Method,Rs,Npoly,ExtrapMethod,Error)
%二維RBF插值，輸入x0，y0散點(diǎn)，z0值。x，y是輸出插值得到的點(diǎn)。
%Method方法，默認(rèn)'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大概和點(diǎn)和點(diǎn)之間的距離差不多就行。
%Npoly多項(xiàng)式擬合。默認(rèn)是1，只擬合1次項(xiàng)。
%ExtrapMethod,外插方法。某個數(shù)，全部賦值為這個數(shù)。'nearest'，按照臨近值外插。'rbf'，默認(rèn)用RBF插值得到的值。三維外插判據(jù)復(fù)雜，不建議用。
%Error誤差。在(-∞,∞)區(qū)間。默認(rèn)是0，無誤差。表示可以在部分誤差范圍內(nèi)去插值。
%Error一般在0~1之內(nèi)。如果只是為了收斂，可以比較小，在0.1以內(nèi)就可以有很好效果；如果想實(shí)現(xiàn)平滑，可適當(dāng)增大，大于1。
N=numel(x0);%點(diǎn)的數(shù)量，也是RBF核的數(shù)量
%整理為列向量
x0=x0(:);
y0=y0(:);
z0=z0(:);
P0=P0(:);
x=x(:);
y=y(:);
z=z(:);
%計(jì)算包線
[k_ch,V]=convhull(x0,y0,z0);
%函數(shù)默認(rèn)信息
narginchk(7,12)
if nargin<7 || isempty(Method)Method='linear';%默認(rèn)線性核函數(shù)
elseif nargin<8 || isempty(Rs)Rs=1.1*(max(x0)-min(x0))/N;%假設(shè)空間均勻分布
elseif nargin<9 || isempty(Npoly)Npoly=1;%默認(rèn)擬合1次項(xiàng)
elseif nargin<10 || isempty(Npoly)ExtrapMethod='rbf';%默認(rèn)不外插
elseif nargin<11 || isempty(Error)Error=0;%默認(rèn)輸入數(shù)據(jù)沒有誤差
end%選擇外插包線
if strcmp(ExtrapMethod,'rbf')Inpoly=[];
elsek_ch=unique(k_ch);Inpoly=IsPointInShape3([x0(k_ch),y0(k_ch),z0(k_ch)],[x,y,z]);
end
if V<=0warning('輸入散點(diǎn)過于集中');
end
%選擇核函數(shù)
switch Methodcase 'linear'fun=@(RMat) Kernel_Linear(RMat,Rs);case 'gaussian'fun=@(RMat) Kernel_Gaussian(RMat,Rs);Error=-Error;%因?yàn)榱泓c(diǎn)不是零，遠(yuǎn)點(diǎn)是0，所以這里誤差為負(fù)case 'thin_plate'fun=@(RMat) Kernel_Thin_plate(RMat,Rs);case 'linearEpsR'fun=@(RMat) Kernel_LinearEpsR(RMat,Rs);Error=-Error;%因?yàn)榱泓c(diǎn)不是零，遠(yuǎn)點(diǎn)是0，所以這里誤差為負(fù)
end
%將原始數(shù)據(jù)分離出多項(xiàng)式項(xiàng)Npoly
if Npoly==0C=ones(N,1);%常數(shù)項(xiàng)wC=mean(P0);
elseif Npoly==1if N<3;warning('點(diǎn)數(shù)太少，建議Npoly等于0');endC=[ones(N,1),x0,y0,z0];%常數(shù)項(xiàng)+一次項(xiàng)wC=C\P0;
elseif Npoly==2if N<9;warning('點(diǎn)數(shù)太少，建議Npoly等于1');endC=[ones(N,1),x0,y0,z0,x0.^2,y0.^2,z0.^2,x0.*y0,x0.*z0,z0.*y0];%常數(shù)項(xiàng)+一次項(xiàng)+二次項(xiàng)wC=C\P0;
elseerror('只支持Npoly=0,1,2')
end
PC=C*wC;%多項(xiàng)式項(xiàng)
P1=P0-PC;
%計(jì)算距離矩陣
DisMat=squareform(pdist([x0,y0,z0]));
K3=fun(DisMat);%每一個核函數(shù)的中心點(diǎn)在節(jié)點(diǎn)上
K3=K3-Error*eye(N);%把誤差函數(shù)加入w=K3\P1;%計(jì)算權(quán)重
%根據(jù)權(quán)重計(jì)算插值
P=zeros(size(x));
for k=1:N %計(jì)算每個核的貢獻(xiàn)，然后疊加R_k=sqrt((x-x0(k)).^2+(y-y0(k)).^2+(z-z0(k)).^2);Pt=w(k)*fun(R_k );P=P+Pt;
end
NOut=length(P);
%再還原回多項(xiàng)式項(xiàng)
if Npoly==0C=ones(NOut,1);%常數(shù)項(xiàng)
elseif Npoly==1C=[ones(NOut,1),x,y,z];%常數(shù)項(xiàng)+一次項(xiàng)
elseif Npoly==2C=[ones(NOut,1),x,y,z,x.^2,y.^2,z.^2,x.*y,x.*z,z.*y];%常數(shù)項(xiàng)+一次項(xiàng)+二次項(xiàng)
end
PC2=C*wC;%再加上多項(xiàng)式項(xiàng)%判斷外插方法
if isnumeric(ExtrapMethod)P(Inpoly)=P(Inpoly)+zC1(Inpoly);%內(nèi)部的正常P(~Inpoly)=ExtrapMethod;%外部的固定值
elseif ischar(ExtrapMethod)switch ExtrapMethodcase 'rbf'P=P+PC2;%再加上多項(xiàng)式項(xiàng)case 'ploy'P(Inpoly)=P(Inpoly)+PC2(Inpoly);%內(nèi)部的正常P(~Inpoly)=PC2(~Inpoly);%外部的直接用多項(xiàng)式值case 'nearest'P(Inpoly)=P(Inpoly)+PC2(Inpoly);%內(nèi)部的正常%找到最近的點(diǎn)F = scatteredInterpolant(x(Inpoly),y(Inpoly),z(Inpoly),P(Inpoly),'nearest','nearest');%外部的直接插值P(~Inpoly)=F(x(~Inpoly),y(~Inpoly),z(~Inpoly));end
elseerror('未識別ExtrapMethod')
end
%檢查結(jié)果
if max(abs(w))/max(abs(P1))>1e3warning('結(jié)果未收斂，建議調(diào)整誤差Error');
elseif (max(P)-min(P))>5*(max(P0)-min(P0))warning('結(jié)果未收斂，建議增大區(qū)間Eps，或者調(diào)整誤差Error');
end
endfunction z=Kernel_Linear(R,s)
%R距離
z=abs(R);%線性函數(shù)
endfunction z=Kernel_Gaussian(R,s)
%R距離
%s大概和采樣點(diǎn)間距差不多就行，可以略大（更胖）
z=exp(-R.^2/2/s^2);%正態(tài)函數(shù)
endfunction z=Kernel_Thin_plate(R,s)
%R距離
z=R.^2.*log(R);%thin_plate
endfunction z=Kernel_LinearEpsR(R,s)
%R距離
%s大概和采樣點(diǎn)間距差不多就行，可以略大（更胖）
s=2*s;%這里要更大
z=s-R;%線性函數(shù)
z(z<0)=0;
endfunction IsInShape=IsPointInShape3(A,points)
%判斷點(diǎn)是否在某個三維凸多面體內(nèi)
%A，凸多面體的頂點(diǎn)。points，N行3列。
%例子：A=[0,0,0;0,0,1;0,1,0;0,1,1;1,0,0;1,0,1;1,1,0;1,1,1];points=rand(1000,3)*3-1;
%IsInShape=IsPointInShape3(A,points)NPoints=size(points,1);%計(jì)算有多少個點(diǎn)
IsInShape=false(NPoints,1);
NConvhull=size(A,1);%計(jì)算凸包上有多少個點(diǎn)
BD_A=[min(A,[],1);max(A,[],1)];
for kp=1:NPointsP_k=points(kp,:);if any(P_k<BD_A(1,:)) || any(P_k>BD_A(2,:)) %如果超過A的矩形邊界，說明肯定不在A內(nèi)continueendNewA=[A;P_k];%把這個點(diǎn)加入新凸包計(jì)算NewP=max(unique(convhull(NewA,'Simplify',false)));%查看新凸包這個點(diǎn)會不會涉及到新點(diǎn)IsInShape(kp)=(NewP==NConvhull);
end
end