數值分析：資料插值方法

阿新 • • 發佈：2019-01-08

插值、擬合和逼近的區別

據維基百科，科學和工程問題可以通過諸如取樣、實驗等方法獲得若干離散的資料，根據這些資料，我們往往希望得到一個連續的函式（也就是曲線）或者更加密集的離散方程與已知資料相吻合，這過程就叫做擬合。通過擬合得到的函式獲得未知點的資料的方法，叫做插值。其中，擬合函式經過所有已知點的插值方法，叫做內插。
擬合是已知點列，從整體上靠近它們；插值是已知點列並且完全經過點列；逼近是已知曲線，或者點列，通過逼近使得構造的函式無限靠近它們。
最小二乘意義下的擬合，是要求擬合函式與原始資料的均方誤差達到極小，是一種整體意義的逼近，對區域性性質沒有要求。而所謂“插值”，就是要在原有離散資料之間“插入”一些值，這就要求插值函式必須通過所有的離散點，插值函式在離散點之外的那些點都相當於“插入”的值。插值有全域性插值，也有區域性插值（比如分段線性插值），插值誤差通常考慮的是逐點誤差或最大模誤差，插值的好壞往往通過某些區域性的性質來體現，比如龍格現象或吉布斯振盪。

插值方法

多項式插值

對於大部分多項式插值函式，插值點的高度值可以視為所有（或某些）節點高度值的線性組合，而線性組合的係數一般是x座標的多項式函式，稱作基函式。對於一個節點的基函式，它在x等於該節點的x時等於1，在x等於其他節點的x時等於0。這就保證曲線必定經過所有節點，所以屬於內插方法。

在本小節，均以一組隨機數作為已知的高度值，使它們對應於間隔固定的x座標，使用不同的插值函式獲得各已知點（稱為插值函式的節點）之外其它x座標所對應的高度值，畫出這些點所對應的曲線。再把所有高度值轉換成灰度值，以顏色的變化比較各插值函式。

原點列如圖：（假定橫向為x，縱向為y。各點x座標的間隔是固定的，但y座標是隨機的）

線性插值

線性插值是用一系列首尾相連的線段依次連線相鄰各點，每條線段內的點的高度作為插值獲得的高度值。

以(xi,yi)表示某條線段的前一個端點，(x(i+1),y(i+1))表示該線段的後一個端點，則對於在[xi,x(i+1)]範圍內的橫座標為x的點，其高度y為：

p(x) = f(x_0) + \frac{f(x_1)-f(x_0)}{x_1-x_0}(x-x_0). \,\!

為便於與後面各函式比較，寫成比較對稱的形式：

其中，yi和y(i+1)的兩個引數稱為基函式，二者之和為1，分別代表yi和y(i+1)對插值點高度的權值。

插值影象如下：

將高度轉化為灰度，得到如下條帶：

線性插值的特點是計算簡便，但光滑性很差。如果用線性插值擬合一條光滑曲線，對每一段線段，原曲線在該段內二階導數絕對值的最大值越大，擬合的誤差越大。

二次插值

如果按照線性插值的形式，以每3個相鄰點做插值，就得到了二次插值：

OpenGL實現程式碼如下：

void quadratic(float p[20][2])
{
float x,y;
int i;
float x01,x02,x12;
glColor3f(0.0,0.0,1.0);
glBegin(GL_LINE_STRIP);
for(i=0;i<20;i+=2)
{
x01=p[i][0]-p[i+1][0];
x02=p[i][0]-p[i+2][0];
x12=p[i+1][0]-p[i+2][0];
for(x=p[i][0];x<=p[i+2][0];x+=1.0)
{
y=(x-p[i+1][0])*(x-p[i+2][0])/x01/x02*p[i][1]-(x-p[i][0])*(x-p[i+2][0])/x01/x12*p[i+1][1]+(x-p[i][0])*(x-p[i+1][0])/x02/x12*p[i+2][1];
glVertex2f(x,y);
}
}
glEnd();
}

二次（分段）插值影象如下：

轉換成灰度值如圖：

二次插值在每段二次曲線內是光滑的，但在每條曲線的連線處其光滑性可能甚至比線性插值還差。二次插值只適合3個節點的情形，當節點數超過3個時，就需要分段插值了。

Cubic插值

如果想要保證各段曲線連線處光滑（一階導數相同），並且不想使用除法運算，可以考慮Cubic插值函式：

其中，v代表插值點，v0、v1、v2、v3代表4個連續的節點。t取值為[0,1]，將會產生一段連線v1和v2的曲線。也就是說，如果有n個節點，Cubic插值函式將會產生(n-2)段曲線，位於首尾兩端的節點不會納入曲線。

實現程式碼如下：

float cubic(float v0,float v1,float v2,float v3,float x)
{
float v32=v3-v2;
float v01=v0-v1;
float v20=v2-v0;
float temp=(v32-v01)*x;
temp=(temp+v01+v01-v32)*x;
temp=(temp+v20)*x+v1;
return temp;
}
void drawCubic(float p[20])
{
float x,y;
int step;
float delta;
glColor3f(0.0,0.0,1.0);
glBegin(GL_LINE_STRIP);
for(step=0;step<17;step++)
{
for(delta=0.0;delta<=1.0;delta+=0.05)
{
x=delta*(p[step+1][0]-p[step][0])+p[step][0];
y=cubic(p[step],p[step+1],p[step+2],p[step+3],delta);
glVertex2f(x,y);
}
}
glEnd();
}

Cuibc插值影象如下：

轉化成灰度如圖：

Cubic插值可以產生整體上光滑的曲線，但容易產生較劇烈的波動，使得曲線的最高點比最高的節點還高、曲線的最低點比最低的節點還低。所以對顏色等取值有嚴格的上下界的資料進行插值時，會造成曲線的擷取，破壞其光滑性。比如顏色（RGB三個分量取值範圍都是[0,255]），如果最高的節點是255，最低的節點是0，那麼插值後負數會被擷取成0，大於255的數會被擷取成255。

拉格朗日多項式插值

解決插值函式的直接構造問題以及插值誤差的估計問題，但是同時也使得插值函式值的計算變得複雜。

依照線性插值和二次插值的思路，可以增加基函式分子和分母的階數，構造拉格朗日插值多項式：

一個n次的拉格朗日插值函式可以繪製經過(n+1)個節點的曲線，但運算量非常大。而且在次數比較高時，容易產生劇烈的震盪（龍格現象）。所以要選擇位置特殊的節點（比如切比雪夫多項式的零點）進行插值，或使用多個次數較低的拉格朗日函式分段插值。（關於拉格朗日多項式和龍格現象，詳見維基百科連結）

分段插值實現程式碼如下：

//n為次數，nmax為節點的總數。n不大於nmax
void lagrange(float p[][2],int n,int nmax)
{
float x,y;
int i,j,t;
float temp;
glColor3f(0.0,0.0,1.0);
for(i=0;i<=(nmax-1);i+=(n-1))
{
glBegin(GL_LINE_STRIP);
for(x=p[i][0];x<=p[i+n-1][0];x+=1.0)
{
y=0.0;
for(j=0;j<n;j++)
{
temp=1.0;
for(t=0;t<n;t++)
{
if(t==j)
continue;
temp*=(x-p[i+t][0])/(p[i+j][0]-p[i+t][0]);
}
y+=temp*p[i+j][1];
}
glVertex2f(x,y);
}
glEnd();
}
}

使用4次拉格朗日多項式分段插值：

轉化為灰度：

拉格朗日多項式插值也存在連線處不光滑的問題。

如果直接使用20次的拉格朗日插值，得到的影象如下：

這樣的插值曲線顯然是不能容忍的~

牛頓插值

從演算法角度考慮，提出便於計算的插值方法，也就是牛頓插值公式。拉格朗日插值每增加一個節點，整個函式要重新計算，計算量巨大。而牛頓插值每增加一個點只需要在多項式的最後增加一項，而且各基函式的係數可以遞迴計算，減少了很多計算量。

均差（Divided differences）

是遞迴除法過程。在數值分析中，也稱差商（Difference quotient），可用於計算牛頓多項式形式的多項式插值的係數。

定義

[數值分析鍾爾傑電子科大]

差商表（高階差商是兩個低一階差商的差商）
$0$ 階差商	$1$ 階差商	$2$ 階差商	$3$ 階差商	$\ldots$	$k-1$ 階差商
$x_{0}$	$f[x_{{0}}]$
$x_{{1}}$	$f[x_{{1}}]$	$f[x_{{0}},x_{{1}}]$
$x_{{2}}$	$f[x_{{2}}]$	$f[x_{{1}},x_{{2}}]$	$f[x_{{0}},x_{{1}},x_{{2}}]$
$x_{{3}}$	$f[x_{{3}}]$	$f[x_{{2}},x_{{3}}]$	$f[x_{{1}},x_{{2}},x_{{3}}]$	$f[x_{{0}},x_{{1}},x_{{2}},x_{{3}}]$
$\ldots$	$\ldots$	$\ldots$	$\ldots$	$\ldots$	$\ldots$
$x_{{k}}$	$f[x_{{k}}]$	$f[x_{{k-1}},x_{{k}}]$	$f[x_{{k-2}},x_{{k-1}},x_{{k}}]$	$f[x_{{k-3}},x_{{k-2}},x_{{k-1}},x_{{k}}]$	$\ldots$	$f[x_{{0}},\ldots ,x_{{k}}]$

由函式表出發，從上往下從左往右依次計算出一階，二階，。。。，各階均差。

牛頓多項式

將n次插值多項式寫作如下形式：

N(x)=[y_{0}]+[y_{0},y_{1}](x-x_{0})+\cdots +[y_{0},\ldots ,y_{k}](x-x_{0})(x-x_{1})\cdots (x-x_{{k-1}})

也就是牛頓插值公式中的各項係數就是均差表中計算出的各階均差的第一個資料！

即牛頓插值公式為：

實現程式碼如下：

void newton(float p[][2],int n,int nmax)
{
float x,y;
int i,j,t;
float temp;
float f[20][20];
glColor3f(0.0,0.0,1.0);
for(i=0;i<=(nmax-1);i+=(n-1))
{
for(t=0;t<n;t++)
{
f[t][0]=(p[i+t][1]-p[i+t+1][1])/(p[i+t][0]-p[i+t+1][0]);
for(j=1;j<=t;j++)
f[t][j]=(f[t-1][j-1]-f[t][j-1])/(p[i+t-j][0]-p[i+t+1][0]);
}
glBegin(GL_LINE_STRIP);
for(x=p[i][0];x<=p[i+n-1][0];x+=1.0)
{
y=p[i][1];
temp=1.0;
for(j=0;j<n;j++)
{
temp*=x-p[i+j][0];
y+=temp*f[j][j];
}
glVertex2f(x,y);
}
glEnd();
}
}

可能由於計算精度的原因，牛頓插值繪製的曲線與拉格朗日插值的曲線略有不同。但次數較高時，牛頓插值也會產生劇烈的震盪。分段4次牛頓插值影象如下：

轉化成灰度如下：

牛頓插值也存在連線處不光滑的缺陷。

埃爾米特插值

以上各多項式插值方法的插值條件都是各節點的座標，在以低階函式分段插值時雖然可以保持曲線的穩定（比較平緩），但在各分段曲線的連線處無法保持光滑（一階導數相等）。埃爾米特插值方法不但規定了各節點的座標值，還規定了曲線在每個節點的各階導數，這樣就可以既保持曲線的穩定，又保證在連線處足夠光滑。

以3次二重（"m重"就是規定座標和曲線在所有節點處1到m-1階導數的值）埃爾米特插值為例：

4個基函式滿足分別只在y0，y1，y0的導數，y1的導數處等於1，而在其他3個條件下等於0。可以把埃爾米特插值看作對座標和導數的插值的組合。

曲線在每個節點的導數可以根據相鄰節點和它的相對位置確定，也可以完全隨機生成。只要位置比較高和比較低的節點的導數絕對值不是很大，就可以使整條曲線落在最高與最低節點定義的帶狀區域內，這樣就可以避免對插值的擷取。

對於本節的示例節點，一種可能的導數值如下：grad[20]={-4.0,4.0,0.5,-2.0,1.0,-2.0,-2.0,2.0,1.0,1.0,0.0,-1.0,-4.0,3.0,0.0,-3.0,3.0,0.0,-4.0,-4.0};

實現程式碼如下：

//p[i][0],p[i][1]為點i的座標，p[i][2]為曲線在點i的導數
//nmax為節點的總數
void hermite(float p[][3],int nmax)
{
float x,y;
float a1,a0,b1,b0;
float x00,x11;
int i;
float temp;
glColor3f(0.0,0.0,1.0);
for(i=0;i<nmax;i++)
{
glBegin(GL_LINE_STRIP);
x00=p[i][0];
x11=p[i+1][0];
for(x=p[i][0];x<=p[i+1][0];x+=1.0)
{
temp=(x11-x00)*(x11-x00);
a0=(x11-3*x00+2*x)*(x11-x)*(x11-x)/temp/(x11-x00);
a1=(-x00+3*x11-2*x)*(x-x00)*(x-x00)/temp/(x11-x00);
b0=(x-x00)*(x-x11)*(x-x11)/temp;
b1=(x-x00)*(x-x00)*(x-x11)/temp;
y=a0*p[i][1]+a1*p[i+1][1]+b0*p[i][2]+b1*p[i+1][2];
glVertex2f(x,y);
}
glEnd();
}
}

分段3次埃爾米特插值影象如下：

轉化為灰度如下：

可見雖然n節點的埃爾米特插值是由(n-1)段曲線構成，但在每一個連線處都是光滑的。

樣條插值

B-樣條

B-樣條是Bezier樣條的一般化，也可推廣為NURBS樣條。先來介紹一下B-樣條：

該式定義了一個n次的B-樣條，它有(m+1)個控制點（樣條曲線的“節點”稱作控制點，因為這些點控制曲線的彎曲方向和程度，但曲線不一定經過這些點），也就有(m+1)個基函式。一般繪製的是完整的曲線，u(min)取0，u(max)取1。當u取值均勻時，該樣條稱作均勻B-樣條。當m=n時，B-樣條退化為Bezier樣條。

函式繪製的曲線始終落在其控制點的凸包（包含一個點集所有點的凸多邊形，而且該凸多邊形所有頂點都來自這個點集）中。對於一個n次的B-樣條，改變一個控制點的位置，只改變它所在的n段曲線（由n+1個控制點定義，且以該點起始）的形狀，而不對其餘的(m-n)段曲線產生影響。

Bezier樣條

Bezier（貝塞爾）樣條是法國工程師皮埃爾·貝塞爾（Pierre Bézier）在1962年為了設計汽車主體外形曲線而提出的。其一般式表示式為：

其中，u取值為[0,1]，pk(k=0,...,n)為(n+1)個節點，n稱為階數。

Bezier樣條還可以遞迴定義為：

含義是n階Bezier樣條是兩條(n-1)階Bezier樣條的插值。

當階數降為1時，Bezier插值退化成線性插值。改變任意一個控制點的位置，整條曲線的形狀都會發生變化。

比較常用的Bezier樣條是3次Bezier：

Beizer樣條在首尾端的切線是前兩個點和最後兩個點的連線。除了第一個點和最後一個點，其他控制點一般都不在曲線上。

3階（4控制點）的Bezier函式影象如下：（黑色曲線為Bezier曲線，藍色折線為控制點的連線）

Bezier樣條可以實現非常快速的運算。為了方便說明，將3次Bezier樣條表示成如下形式：

對於任意給定的u，可以通過合併的方式將原來的7次乘法、4次加法減少為3次乘法、3次加法：

一般情況下，應用Bezier樣條繪製曲線時，都是先給定一個很小的步長t（步長足夠小才能保證Bezier曲線的精確），讓t從0取到1，從頭到尾繪製整條曲線。在t不變的條件下，可以使用更快速的差分方法，利用前一個點計算出下一個點的值，將每步的計算量減小到只有3次加法：

只需要在繪製曲線之前計算4個常數，就可以很快地計算出曲線上的所有點：

採用這種方式，Bezier樣條的運算量只隨階數線性增長。

實現程式碼如下：

void bezier3(float xx0,float yy0,float xx1,float yy1,float xx2,float yy2,float xx3,float yy3)
{
GLint i;
GLfloat x,y;
GLdouble ax,bx,cx,ay,by,cy;
GLdouble t;
GLdouble dx,d2x,dy,d2y,delta;
GLdouble d3x,d3y;
ax=-xx0+3.0*xx1-3.0*xx2+xx3;
bx=3.0*xx0-6.0*xx1+3.0*xx2;
cx=-3.0*xx0+3.0*xx1;
ay=-yy0+3.0*yy1-3.0*yy2+yy3;
by=3.0*yy0-6.0*yy1+3.0*yy2;
cy=-3.0*yy0+3.0*yy1;
delta=0.0005;
d3x=6.0*ax*delta*delta*delta;
d3y=6.0*ay*delta*delta*delta;
x=xx0;
y=yy0;
glBegin(GL_LINE_STRIP);
glVertex2f(x,y+200);
d2x=6.0*ax*delta*delta*delta+2.0*bx*delta*delta;
dx=ax*delta*delta*delta+bx*delta*delta+cx*delta;
d2y=6.0*ay*delta*delta*delta+2.0*by*delta*delta;
dy=ay*delta*delta*delta+by*delta*delta+cy*delta;
for(t=0.0;t<=1.0;t+=delta)
{
d2x+=d3x;
dx+=d2x;
x+=(float)dx;
d2y+=d3y;
dy+=d2y;
y+=(float)dy;
glVertex2f(x,y);
}
glEnd();
}

NURBS樣條

NURBS（Non-Uniform Rational B-Splines 非均勻有理B-樣條），是貝塞爾樣條的推廣。“非均勻”的意思是控制點的間隔可以是不均勻的，“有理”的意思是各控制點帶有不同的權值。和Bezier樣條相比，它對曲線形狀的控制更自由：

其基函式B(k,d)與B-樣條的基函式相同，w(k)為各點的權因子。和B-樣條一樣，改變一個控制點的位置，只改變它所在的n段曲線的形狀，而不對其餘的(m-n)段曲線產生影響。

插值的一些應用

噪聲

噪聲影象一般需要二維或三維插值函式，不過了解了各一維插值函式，就很容易擴充套件到二維和三維了。關於二維噪聲影象的產生，請參見《二維噪聲影象》。

水波

下面這個影象就是用插值函式繪製的模擬水波的三維網格（靜態圖）：

雖然這張網格看起來似乎需要二維插值，但這張網格是由6個不同波長和頻率的Gestner合成的，每個獨立的波形只需要一維插值生成。

ref:

數值分析：資料插值方法

插值、擬合和逼近的區別

插值方法

多項式插值

線性插值