許偉志，關(guān)履泰，徐應(yīng)祥

(中山大學(xué)科學(xué)計算與計算機應(yīng)用系，廣東 廣州 510275)

多元散亂數(shù)據(jù)擬合廣泛用于數(shù)據(jù)壓縮、汽車外形設(shè)計、船體放樣、飛機機翼機身設(shè)計、服裝設(shè)計、地質(zhì)探礦、醫(yī)學(xué)圖像處理等諸多領(lǐng)域，從上世紀(jì)60年代開始，眾多研究工作者對散亂數(shù)據(jù)曲面擬合問題進(jìn)行了一系列的研究[1-6]。但是都不如一元B樣條處理一元散亂數(shù)據(jù)那樣理想。李岳生和關(guān)履泰在1989年試圖推廣一元三次自然樣條對散亂點插值的方法到二元情形，提出了散亂數(shù)據(jù)的二元多項式自然樣條插值，研究了廣義混合樣條函數(shù)空間矩形域帶連續(xù)邊界條件和離散邊界條件的多元散亂數(shù)據(jù)最優(yōu)插值問題[7]。胡日章研究了相應(yīng)的光順逼近問題[8]，韓國強從計算的角度提出了散亂數(shù)據(jù)光順逼近二步法[9]。但是這類樣條的目標(biāo)泛函比較復(fù)雜，帶有一系列的積分項[7- 8, 10]，而且如果離散邊界插值點不夠多的話插值效果會比較差。為了克服這些缺點，關(guān)履泰、許偉志等研究了一類新的二元自然樣條插值方法，該方法的目標(biāo)泛函較為簡單，沒有離散邊界插值點，更符合實際情況[11]。然而當(dāng)測量數(shù)據(jù)受誤差和噪聲影響時，插值方法不如光順方法那樣能夠很好地反映數(shù)據(jù)的本質(zhì)。為此，本文在文[11]的基礎(chǔ)上提出了一種多項式自然樣條光順方法，并研究了解的特征性質(zhì)和解的結(jié)構(gòu)。

1 問題的提出

給定矩形區(qū)域R=[a1,a2]×[c1,c2],a

Y=L2(R)，T:X→Y是一個從X到Y(jié)的線性微分算子，在矩形區(qū)域R內(nèi)定義為

在邊界x=a,y=c上，它帶邊界條件

記Z=RN是N維歐幾里得空間。假設(shè)A:X→Z是一個插值算子，定義為

Au=(u(x1,y1),…,u(xN,yN))

為了方便，令z=(z1,…,zN)。考慮以下在希爾伯特空間Hm,n(R)中的樣條光順問題，稱之為散亂點二元多項式自然樣條光順問題。

問題T尋找一個函數(shù)σ(x,y)∈X，滿足

(1)

其中z=(z1,z2,…,zN)，ρ>0是事先給定的常數(shù)。這里‖Tu‖2=?R(u(m,n)(x,y))2dxdy，帶邊界條件

(2)

即尋找一個函數(shù)σ(x,y)∈X，滿足

其中ρ>0，是事先給定的常數(shù)。帶邊界條件

2 帶邊界條件的自然樣條基函數(shù)的構(gòu)造

容易驗證下列引理。

引理1 算子T的化零子空間滿足

N(T)=P〈m,n〉={u|u(x,y)=

定理1(特征定理)σ(x,y)∈X是二元多項式自然樣條光順問題解的充要條件為

ρ?Rσ(m,n)(x,y)u(m,n)(x,y)dxdy+

(3)

對一切u(x,y)∈X成立。

證明(必要性) 設(shè)σ(x,y)∈X是二元多項式自然樣條光順問題的解，任意的u(x,y)∈X,ε為微小參數(shù)。記uε(x,y)=σ(x,y)+εu(x,y)代入問題(1)式，得

F(ε)=ρ?R(σ(m,n)(x,y)+εu(m,n)(x,y))2dxdy+

由于σ(x,y)∈X為問題T的解，所以F′(0)=0，即

F′(ε)=2ρ?Ru(m,n)(x,y)·

(σ(m,n)(x,y)+εu(m,n)(x,y))dxdy+

即

ρ?Rσ(m,n)(x,y)u(m,n)(x,y)dxdy+

J(u)=J(σ+u-σ)=

J(σ)+J1(u-σ)+a(σ,u-σ)

由充分性條件，有

ρ?Rσ(m,n)(x,y)u(m,n)(x,y)dxdy+

對一切u(x,y)∈X都成立。那么對于u-σ∈X, 有

ρ?Rσ(m,n)(x,y)(u(m,n)(x,y)-σ(m,n)(x,y))dxdy+

引理2 設(shè)σ(x,y)∈X是二元多項式自然樣條光順的解，那么必存在系數(shù)λi與ki(i=1,…,N)使得

證明由特征定理可知

ρ?Rσ(m,n)(x,y)u(m,n)(x,y)dxdy+

即

Au(x,y)=(u(x1,y1),…,u(xN,yN))=

((k1,u),…,(kN,u))

定理2 二元多項式自然樣條光順函數(shù)σ(x,y)具有如下的顯式及緊湊格式的表達(dá)式

其中g(shù)i(x,y)=G(xi,m,a;x)G(yi,n,c;y),i=1,…,N是(2m-1,2n-1)次自然樣條基函數(shù)，并且有

對任意X中的元素u，在a點Taylor展開

u(x,y)=u(a,y)+(x-a)u(1,0)(a,y)+…+

然后在c點Taylor展開u(a,y),u(1,0)(a,y),…,u(m-1,0)(a,y),u(m,0)(τ,y)得

u(xi,yi)=u(a,c)+(yi-c)u(0,1)(a,c)+…+

(xi-a)u(1,0)(a,c)+…+

利用上面結(jié)論，并注意邊界條件，對任意X中的元素u有

f(x,y)∈N(T)

其中g(shù)i(x,y)=G(xi,m,a;x)G(yi,n,c;y),i=1,…,N。滿足

所以

其中g(shù)i(x,y)=G(xi,m,a;x)G(yi,n,c;y),i=1,…,N是(2m-1,2n-1)次自然樣條基函數(shù)，并且有

且

v=0,…,n-1

從求n-1階導(dǎo)數(shù)開始，由上式推出

然后，利用這個結(jié)果代入求n-2階導(dǎo)數(shù)的式子中得到

(n-2)!cn-2(yi,n,c;y)+

cn-2(yi,n,c;y)=

再用上述這些結(jié)果，代入求n-3階導(dǎo)數(shù)的式子中，有

(n-3)!cn-3(yi,n,c;y)+

(n-3)!(cn-2(yi,n,c;y)(y-c)+

cn-3(yi,n,c;y)=

繼續(xù)下去……，得到

3 散亂點自然樣條光順問題求解

定理3 如果相應(yīng)的齊次多項式光順問題只有恒零解，則對任意散亂數(shù)據(jù)(xi,yi,zi),i=1,…,N，(2m-1,2n-1)次多項式自然樣條光順問題有唯一解，可表示為

其中g(shù)i(x,y)=G(xi,m,a;x)G(yi,n,c;y),i=1,…,N是(2m-1,2n-1)次自然樣條基函數(shù)，并且有

且

cij及λi由以下的方程組來確定

其中

Λ=(λ1,…,λN)T，z=(z1,…,zN)T，

C=(c00,c01,c02,…,c0(n-1),

c10,c11,…,c1(n-1),…,c(m-1)(n-1))T，

A=(aij)N×N，

證明由特征定理可知，σ(x,y)∈X是二元多項式自然樣條光順問題解的充要條件為

ρ?Rσ(m,n)(x,y)u(m,n)(x,y)dxdy+

(5)

對一切u(x,y)∈X成立。又由引理2知

代入(5)式得

即

再由u(x,y)的任意性可知ρλi+σ(xi,yi)-zi=0,i=1,2,…,N。即

得到前N個方程。注意在定理2的證明過程中，我們已經(jīng)得出

這就是后面m×n個方程。寫成矩陣形式后則證得本定理。由二項式定理可證以下結(jié)論。

4 求光順函數(shù)算法

由定理3和定理4可得求光順函數(shù)算法如下：

(ⅰ)輸入散亂數(shù)據(jù)點(xi,yi),i=1,…,N與要在其上擬合的數(shù)據(jù)zi,i=1,…,N及合適的光順參數(shù)ρ；

(ⅱ)計算矩陣A；

(ⅴ)求出光順函數(shù)

5 算 例

例1 擬合函數(shù)：z=1/(1+x2+y2)，區(qū)間：[0,1]×[0,1]，邊界：a=-1,c=-1。對301個散亂點(由隨機數(shù)產(chǎn)生)相應(yīng)的函數(shù)值，分別運用雙三次自然樣條光順方法(ρ=0.005)與雙三次自然樣條插值方法去擬合函數(shù)， 比較它們在所在區(qū)域30×30個格子點上的函數(shù)值Fi,導(dǎo)函數(shù)值的平均誤差與最大誤差(見表1)。

表1 301個散亂數(shù)據(jù)插值和光順的函數(shù)值、導(dǎo)函數(shù)值誤差

一般地，選取合適參數(shù)的光順樣條方法得出的曲面較用插值方法得到的曲面光滑。

例2 某地區(qū)重金屬釷的3 580個測量值[6]，分布如圖1，分別用雙三次自然樣條插值和光順(ρ=0.012)，應(yīng)用極小殘量法(Generalized Minimal Residual Algorithm)迭代求解對稱不定的線性方程組[12]，得到的結(jié)果(如圖2)。

圖1 數(shù)據(jù)點的分布

圖2 (a)插值圖像 (b)光順圖像(ρ=0.012)

由于測量數(shù)據(jù)存在誤差和噪聲，插值圖像并不能很好的反映重金屬釷的分布情況，而光順圖像很好地反映了重金屬釷的實際分布情況。

例3 擬合函數(shù)：z=cos(x)×sin(y)區(qū)間：[0,4]×[0,4]，對50、100、1 000個散亂點(由隨機數(shù)產(chǎn)生)相應(yīng)的函數(shù)值，分別用文[11]和文[8]中的雙三次自然樣條插值方法去擬合函數(shù)，相應(yīng)的系數(shù)矩陣的條件數(shù)比較如表2。

表2 方法[8]和[11]中的系數(shù)矩陣條件數(shù)比較

對相同散亂數(shù)據(jù)進(jìn)行雙三次樣條插值擬合，本文中方程組的系數(shù)矩陣條件數(shù)比文[8]中的系數(shù)矩陣條件數(shù)小,更適合大規(guī)模散亂數(shù)據(jù)擬合。

結(jié)論：選取合適參數(shù)的光順樣條方法得出的曲面較平滑，外形比用插值方法好。在實際應(yīng)用中，當(dāng)數(shù)據(jù)有誤差和噪聲時，建議選取適當(dāng)參數(shù)應(yīng)用光順樣條方法獲得擬合曲面。另外，本文的方法容易推廣到高維情形。本文中方程組的系數(shù)矩陣條件數(shù)比文[8]中的系數(shù)矩陣條件數(shù)小,更適合大規(guī)模散亂數(shù)據(jù)擬合。

參考文獻(xiàn)：

[1] WAHBA G. Spline models for observational data[J]. Soc Ind Appl Math, 1990.

[2] GUAN L T, LIU B. Surface design by natural splines over refined grid points[J]. J Comp Anal and Appl, 2004, 163(1):107-115.

[3] LAI M J. Multivariate splines for data fitting and approximation, Approximation Theory XII[M]. San Antonio, 2007, edited by Neamtu M and Schumaker L L, Brentwood :Nashboro Press, 2008: 210-228.

[4] LAI M J, SCHUMAKER L L. Spline functions over triangulations[M]. London: Cambridge University Press, 2007.

[5] 吳宗敏. 散亂數(shù)據(jù)擬合的模型、方法和理論[M]. 北京：科學(xué)出版社, 2008.

[6] BILLINGS S D, NEWSAM G N, BEATSON R K. Smooth fitting of geophysical data using continuous global surfaces[J]. Geophysics, 2002, 67 (6): 1810-1822.

[7] 李岳生, 胡日章. 多元散亂數(shù)據(jù)的樣條插值法[J]. 高等學(xué)校計算數(shù)學(xué)學(xué)報, 1990, 3: 215-226.

[8] 胡日章. 矩形域上散亂數(shù)據(jù)的離散邊界條件光順逼近[J]. 中山大學(xué)學(xué)報：自然科學(xué)版, 1990, 29(4): 17-22.

[9] 韓國強. 矩形域上散亂數(shù)據(jù)樣條光順法[J]. 華南理工大學(xué)學(xué)報：自然科學(xué)版, 1993, 21(4): 102-109.

[10] GUAN L T. Bivariate polynomial natural spline interpolation algorithms with local basis for scattered data[J]. J Comp Anal and Appl, 2003, 1: 77-101.

[11] 關(guān)履泰，許偉志，朱慶勇. 一種散亂點雙三次多項式自然樣條插值[J]. 中山大學(xué)學(xué)報：自然科學(xué)版, 2008, 47(5): 1-4.

[12] SAAD Y, SCHULTZ M H. GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems[J]. SIAM J Sci Stat Comput, 1986, 7(3): 856-869.