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

(1．中山大學(xué)新華學(xué)院,廣東 廣州 510520；2．中山大學(xué)科學(xué)計(jì)算與計(jì)算機(jī)應(yīng)用系,廣東 廣州 510275)

散亂數(shù)據(jù)擬合在許多領(lǐng)域中都有廣泛的應(yīng)用[1-4]。以前考慮得比較多的是平面散亂數(shù)據(jù)擬合問題,但近年來已有不少更高維的問題出現(xiàn),例如動畫設(shè)計(jì)與制作,動態(tài)醫(yī)學(xué)三維圖像演示等方面有熱烈的討論[5]。從上世紀(jì)60年代以來，眾多的科技工作者對散亂數(shù)據(jù)曲面插值問題進(jìn)行了一系列的研究[6-10]。近來，Lai、吳宗敏等有一系列總結(jié)性的工作[11-14]，但是對多元散亂數(shù)據(jù)與大規(guī)模散亂數(shù)據(jù)插值問題的解決一直不夠理想。目前在空間散亂數(shù)據(jù)插值中用得較多的三角剖分方法很難推廣到更高維去，而徑向基函數(shù)方法則在計(jì)算穩(wěn)定性與良好逼近性方面有欠缺，不像一元三次B-樣條那樣有一系列良好的性質(zhì)，能很好地解決一元散亂數(shù)據(jù)插值問題。李岳生和關(guān)履泰在1989年試圖推廣一元三次自然樣條對散亂點(diǎn)插值的方法到二元情形，提出散亂數(shù)據(jù)的二元多項(xiàng)式自然樣條插值，研究了廣義混合樣條函數(shù)空間矩形域帶連續(xù)邊界條件和離散邊界條件的多元散亂數(shù)據(jù)最優(yōu)插值問題[15]。1993年，Chui和關(guān)履泰把二元的結(jié)果全面地推廣到一般多元情形[16]。關(guān)履泰在1993年給出了一種二元散亂數(shù)據(jù)的多項(xiàng)式自然樣條光順[17]，并在1997年研究了類似B-樣條的多元局部支撐基函數(shù)[18]，并于2003年發(fā)表了關(guān)于這種局部支撐基的性質(zhì)及插值自然樣條算法[19]。但是這類樣條的目標(biāo)泛函比較復(fù)雜，帶有一系列的積分項(xiàng)，沒有明顯的物理意義，而且在很多情形下，如果區(qū)域邊界上沒有插值點(diǎn)或插值點(diǎn)太少，那么插值效果就會比較差。關(guān)履泰、許偉志和朱慶勇等研究了一類新的二元自然樣條插值方法，該方法的目標(biāo)泛函較為簡單，沒有離散邊界插值點(diǎn)，更符合實(shí)際情況,而且結(jié)果比以前的更加簡單,計(jì)算起來更加方便[20]。

Laurent[21]在其著作《逼近與優(yōu)化》中推廣一元三次樣條的優(yōu)化性質(zhì)，提出Hilbert空間樣條進(jìn)行研究，Bezhaev等[22]其后在《樣條的變分理論》中對多元Hilbert空間樣條作了進(jìn)一步的討論，李岳生[4]對樣條變分法的歐拉方程有進(jìn)一步的分析。不論是一元還是多元Hilbert空間樣條，解決問題的難點(diǎn)和關(guān)鍵都在于如何提出一個合理的目標(biāo)泛函，使得用變分法得出的歐拉方程容易求解，從而能夠求得問題顯式的解。在動態(tài)醫(yī)學(xué)圖象處理，三維動畫設(shè)計(jì)等實(shí)際應(yīng)用問題中，出現(xiàn)了采樣點(diǎn)不僅要確定三維空間中位置，而且還有時間的限制，這使得采樣數(shù)據(jù)出現(xiàn)了4維散亂數(shù)據(jù)的情形，因此對于4維散亂數(shù)據(jù)的處理也要作進(jìn)一步討論，如怎樣更有效的對4維散亂數(shù)據(jù)進(jìn)行插值或光順等。本文深入進(jìn)行目標(biāo)泛函的研究，提出簡單合理的目標(biāo)泛函進(jìn)行討論，構(gòu)造出一種更符合實(shí)際的Hilbert空間三元帶自然邊界條件的樣條。這種帶自然邊界條件的樣條在簡單情況下,即不帶導(dǎo)數(shù)條件下是三奇次的多項(xiàng)式(即關(guān)于每個變元xj是2pj-1次多項(xiàng)式，j=1,2,3)，其算子T的核是〈p1,p2,p3〉階(見定義1)多項(xiàng)式空間的子空間，在一般情形下這種樣條是三元多項(xiàng)式。用這種新的三元帶自然邊界條件的樣條考慮4維空間散亂數(shù)據(jù)的光順問題，則該問題的帶自然邊界條件的樣條解與以往不同，其基函數(shù)有良好的邊界條件與簡單的表達(dá)式，且求解此問題時所用到的線性方程組的系數(shù)矩陣是對稱矩陣，這使得求解也更為方便。雖然解在表達(dá)式上與單邊基有關(guān)，但可以證明對應(yīng)的系數(shù)矩陣是正定對稱的，可以有較穩(wěn)定的迭代算法，這些我們在以后的文章中討論。

1 散亂數(shù)據(jù)三元帶自然邊界條件的光順樣條

若函數(shù)u在區(qū)域Ω的邊界xi=ai(i=1,2,3)滿足如下的條件：

u(i1,p2,0)(a1,x2,x3)=0,u(p1,i2,0)(x1,a2,x3)=0,

u(i1,i2,p3)(a1,a2,x3)=0,u(p1,p2,i3)(x1,x2,a3)=0

(1)

其中ik=0,1,…,pk-1,k=1,2,3,則稱u滿足自然邊界條件NB，簡記為u|NB=0。記Y1=L2(Ω)，令T1:X→Y1是一個從X到Y(jié)1的線性算子，定義為T1(u)=u(p)=?pu，其中

X={u|u∈X1,u|NB=0}

Au=(u(α)(x1),u(α)(x2),…,u(α)(xN))

其中u(α)(xi)=u(α)(x)|x=xi,i=1,…,N,k=1,2,3。對每個k有αk∈Ik,而Ik?Ipk={0,1,…,pk-1}都是任意指定的指標(biāo)集,總指標(biāo)數(shù)為q。

問題T 求u(x)∈X，使得如下的泛函

ρ

(2)

取Y=Y1?Q為Y1與Q的乘積空間,定義線性連續(xù)算子T:X→Y為T=T1?A,則問題T還可寫成:

問題T1 求σ(x)∈X，使

(3)

這里考慮到正數(shù)ρ的作用,定義空間Y的范數(shù)為

ρ

除以上問題中引入的記號外，為表示簡單，后文中還需要用到如下的記號

定理1 算子T1的化零子空間滿足

(4)

反之，如果u∈N(T1)且滿足邊界條件(1)，則由u(p1,p2,p3)(x1,x2,x3)=0，于是有

再由u(p1,p2,k)(x1,x2,a3)=0,k=0,1,…,p3-1，知

,

以及

μ(μ-1)…(μ-k+1),,

k=1,…,p3-1

于是可推得

再由邊界條件

u(i,p2,0)(a1,x2,x3)=0,i=0,1,…,p1-1,

u(p1,j,0)(x1,a2,x3)=0,j=0,1,…,p2-1

可知

,

i=1,…,p1-1,j=1,…,p2-1

當(dāng)a1≠0,a2≠0,時，由x3的任意性及上述各式容易解得

p1-1,ν=0,1,…,p2-1,μ=0,1,…,p3-1

于是可知

,

代入u(x1,x2,x3)的表達(dá)式，便知u可表示為

即u∈P〈p1,p2,p3〉=P〈p〉。綜上可知N(T1)=P〈p〉。

定理2 算子T的化零子空間滿足

α=(α1,α2,α3),αk∈Ik,k=1,2,3,i=1,…,N}

(5)

反之，如果u∈N(T)，則Tu=0且滿足邊界條件(1)。因此T1u=0,Au=0。

一方面由T1u=0，根據(jù)定理1可知u∈P〈p〉。

另一方面再由Au=0,知

由此得

,αk∈Ik,k=1,2,3,i=1,…,N

綜上兩方面可知定理得證。

2 散亂數(shù)據(jù)三元自然光順樣條特征及構(gòu)造

定理3 (特征定理)σ∈X是4維帶自然邊界條件的光順樣條的充分必要條件如下：

σ(p)(x)u(p)(x)dx+

對一切u∈X成立.

證明根據(jù)定義，上式即

〈T1σ,T1u〉Y1+ρ〈Aσ-z,Au〉Q=0

〈T1σ,T1u〉Y1+ρ〈Aσ-z,Au〉Q

因此有〈T1σ,T1u〉Y1+ρ〈Aσ-z,Au〉Q=0。

反之，若〈T1σ,T1u〉Y1+ρ〈Aσ-z,Au〉Q=0，則

即J(u)在σ處取極小。

為了找出問題T的解，即帶自然邊界條件的光順樣條所具有的形式，先討論在散亂點(diǎn)方型域Ω中不帶導(dǎo)數(shù)條件時的情形，此時稱問題T的解為簡單帶自然邊界條件的光順樣條。

定理4 (構(gòu)造定理) 4維散亂數(shù)據(jù)的三元簡單帶自然邊界條件的光順樣條σ(x)是三奇次多項(xiàng)式且具有如下的顯式及緊湊格式的表達(dá)式:

(6)

其中

,aj;xj)，i=1,…,N

(7)

是三奇次樣條基函數(shù)，并且有

(8)

此處

(9)

kj=0,1,…,pj-1,j=1,2,3。系數(shù)cjkl及di由方程組

(10)

決定，此處D=(di)T,C=(clmn)T,矩陣E=(eij)是N×N階方陣，且滿足

BTD=0

(11)

對任意X中的元素u，在a1點(diǎn)Taylor展開，有

類似地，由Taylor展開還有

u(i,0,0)(a1,x2,x3)=

u(i,j,0)(a1,x2,x3)=

u(p1,0,0)(t1,x2,x3)=

u(p1,p2,0)(t1,t2,x3)=

及自然邊界條件(1)。所以

再由gi(x1,x2,x3)滿足自然邊界條件

,x2,x3)=0,l=0,1,…,p1-1

可知

[G(l)(x1i,a1;a1)]G(p2)(x2i,a2;x2)G(x3i,a3;x3)=0

由此有

,

l=0,1,…,p1-1

于是可得

，k=0,1,…,p1-1

類似地可得G(x2i,a2;x2)及G(x3i,a3;x3)，其形式與G(x1i,a1;x1)類同，這里不再贅述。

再由特征定理可知在不帶導(dǎo)數(shù)的簡單情形下有

σ(p)(x)u(p)(x)dx+ρσ(xi)-zi)u(xi)=0

而由前述證明有

σ(p)(x)u(p)(x)dx=

,

i=1,…,N

這正好是方程

AD+BC=z

(12)

綜合式(11)、(12)，定理得證。

下面用類似的方法再考慮帶導(dǎo)數(shù)條件的一般情形，則有

定理5 (構(gòu)造定理)設(shè)σ(x)是4維散亂數(shù)據(jù)的三元帶自然邊界條件的光順樣條，則其具有如下的顯式及緊湊格式的表達(dá)式

(13)

其中

,aj;xj),i=1,…,N

(14)

是三維樣條基函數(shù)，并且有

(15)

此處

(16)

kj=0,…,pj-1,j=1,2,3.系數(shù)cjkl及di由方程組

(17)

BTD=0

(18)

對任意X中的元素u，在a1點(diǎn)Taylor展開，有

u(α1,α2,α3)(x1,x2,x3)=

類似地，由Taylor展開還有

u(α1+i,α2,α3)(a1,x2,x3)=

u(α1+i,α2+j)(a1,x2,x3)=

u(p1,α2,α3)(t1,x2,x3)=

u(p1,p2,α3)(t1,t2,x3)=

及自然邊界條件(1)。所以

余下討論與定理4類同，這里不再贅述。

3 數(shù)值例子

例1 取三元函數(shù)

u(x1,x2,x3)=cosx1sinx2sinx3

邊界取為a1=a2=a3=-0.5,b1=b2=b3=4.5,即區(qū)域Ω為[-0.5,4.5]×[-0.5,4.5]×[-0.5,4.5]。散亂數(shù)據(jù)位于[0.5,3.5]×[0.5,3.5]×[0.5,3.5]，由隨機(jī)函數(shù)產(chǎn)生。用〈2,2,2〉階三奇次簡單帶自然邊界條件光順樣條擬合函數(shù)u。在此僅列出當(dāng)參數(shù)ρ=0.5，x3=2時，取500個散亂點(diǎn)和1 500個散亂點(diǎn)的光順曲面圖形，最大誤差和平均誤差(圖1，圖2及表1)。

表1 x3=2時的誤差

圖1 x3=2時500個散亂點(diǎn),光順與誤差曲面

圖2 x3=2時1 500個散亂點(diǎn),光順與誤差曲面

例2 取三元函數(shù)

u(x1,x2,x3)=

邊界取為a1=a2=a3=-0.1,b1=b2=b3=2.2,即區(qū)域Ω為[-0.1,2.2]×[-0.1,2.2]×[-0.1,2.2]。散亂點(diǎn)位于[0.9,1.2]×[0.9,1.2]×[0.9,1.2]，由隨機(jī)函數(shù)產(chǎn)生。分別用〈4,4,4〉階，參數(shù)ρ=100的三奇次簡單帶自然邊界條件的光順樣條與三奇次簡單帶自然邊界條件的插值樣條分別擬合函數(shù)u。在此僅列出當(dāng)x3=1時，取500個散亂點(diǎn)(其中10個點(diǎn)上有值為1的噪聲)和1 500個散亂點(diǎn)(其中15個點(diǎn)上有值為1的噪聲)時光順與插值的最大誤差和平均誤差(表2)。由表2可以看出，當(dāng)測得的數(shù)據(jù)有噪聲時，光順的效果要比插值好很多,這說明光順具有一定的去噪效果。

表2 x3=1時的誤差

例3 取三元函數(shù)

相應(yīng)的邊界取為a1=a2=a3=-3,b1=b2=b3=3,即區(qū)域Ω為[-3,3]×[-3,3]×[-3,3]。散亂點(diǎn)位于[-2,2]×[-2,2]×[-2,2]，由隨機(jī)函數(shù)產(chǎn)生。用參數(shù)ρ=100為〈4,4,6〉階三奇次簡單帶自然邊界條件的光順樣條與〈4,4,6〉階三奇次簡單帶自然邊界條件的插值樣條擬合函數(shù)u。在此僅列出當(dāng)x3=1時，取1 500個散亂點(diǎn)的光順曲面，插值曲面及曲面真實(shí)的等高線圖(圖3，圖4，圖5)。由等高線圖可知，光順曲面的等高線圖更接近曲面真實(shí)等高線圖。

圖3 光順曲面等高線

圖4 插值曲面等高線

圖5 曲面真實(shí)等高線

4 小 結(jié)

本文主要針對4維散亂數(shù)據(jù)，提出了一種帶自然邊界條件的樣條光順方法。在使得即定的目標(biāo)泛函達(dá)到極小的要求下，給出了解的構(gòu)造過程，得到了光順問題的解——稱為4維散亂數(shù)據(jù)帶自然邊界條件的光順樣條。這種光順樣條是三元的多項(xiàng)式，對每一個變量都奇次多項(xiàng)式，但并不要求關(guān)于每個變量的次數(shù)都是一樣的，這使得在具體問題中對帶自然邊界條件的光順樣條可以根據(jù)需要選擇其對每個變量的次數(shù)，如可選擇帶自然邊界條件的光順樣條是〈4,4,2〉階或者〈2,2,4〉階的等等。從構(gòu)造過程可知，這種三元帶自然邊界條件的光順樣條并不是由一元奇次多項(xiàng)式通過張量積的方法得到，是一種新的非張量積多元樣條。由于是多項(xiàng)式，所以在計(jì)算上來說是較為簡單的。但是由于篇幅的限制，本文只討論和研究了這種帶自然邊界條件的光順樣條的特征與構(gòu)造方法，而未能將如下的問題一一作出詳細(xì)討論：

(i)光順問題解的存在唯一性；

(ii)光順問題與插值問題的關(guān)系；

(iii)光順問題解的收斂性與誤差估計(jì)；

(iv)求光順問題解時得到的線性方程組的特點(diǎn)及如何求解等。

對于以上問題將另文討論和研究。

參考文獻(xiàn)：

[1]ANDREW R W.統(tǒng)計(jì)模式識別[M].王萍，楊培龍，羅穎昕,譯.2版.北京：電子工業(yè)出版社,2004.

[2]TONY F C,SHEN J H.圖像處理與分析：變分、PDE、小波及隨機(jī)方法[M].北京：科學(xué)出版社,2009.

[3]唐澤圣.三維數(shù)據(jù)場可視化[M].北京：清華大學(xué)出版社,1999.

[4]李岳生.分布?xì)W拉方程與分片函數(shù)的表示[J].計(jì)算數(shù)學(xué),2006,28(3): 225-236.

[5]KAZHDAN M,BOLITHO M,HOPPE H.Poisson surface reconstruction[C]∥Proceeding of Eurogrouphics Synposium on Geometry Processing,Cagliari,Italy,2006: 61-70.

[6]王仁宏.多元樣條及其應(yīng)用[M].北京:科學(xué)出版社,1992.

[7]崔錦泰.多元樣條理論用其應(yīng)用[M].程正興,譯.西安:西安交通大學(xué)出版社,1991.

[8]GUAN L T,LIU B.Surface design by natural splines over refined grid points[J].Journal of Computational and Applied Mathematics,2004,163(1): 107-115.

[9]關(guān)履泰，羅笑南，黎羅羅.計(jì)算機(jī)輔助幾何圖形設(shè)計(jì)[M].北京：高等教育出版社∥海德爾堡:施普林格出版社,1999.

[10]關(guān)履泰，覃廉，張健.用參數(shù)樣條插值挖補(bǔ)方法進(jìn)行大規(guī)模散亂數(shù)據(jù)曲面造型[J].計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報，2006,18(3): 372-377.

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

[12]LAI M J.Multivariate splines for data fitting and approximation,approximation theory[M].Brentwood: Nashboro Press,2008: 210-228

[13]ZHOU T H,HAN D F,LAI M J.Energy minimization method for scattered data hermit interpolation[J].Applied Numerical Mathematics,2008,58: 646-659.

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

[15]李岳生，關(guān)履泰.散亂數(shù)據(jù)的二元多項(xiàng)式自然樣條插值[J].計(jì)算數(shù)學(xué)，1990,23(1):135-146.

[16]CHUI C K,GUAN L T.Multivariate polynomial natural spline for interpolation of scattered data and other applications[C]∥ Workship on Computational Geometry,World Sciedtific,1993:77-98.

[17]關(guān)履泰.散亂數(shù)據(jù)的多項(xiàng)式自然樣條光順與廣義插值[J].計(jì)算數(shù)學(xué)，1993,26(4): 383-401.

[18]GUAN L T.A local basis for bivariate polynomial natural splines of scattered data[C]∥Guangzhou International Symposium of Computational Mathematics,1997.

[19]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.

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

[21]LAURENT P J.Approximation et optimization[M].Paris: Hermann,1972.

[22]BEZHAEV A Y,VASILENKO V A.Variational theory of splines[M].New York: Kluwer Academic/Plenum Publishers,2001.