MQ徑向基函數(shù)的理論、方法及應(yīng)用?

喬遠陽，吳技蓮，馮新龍

(新疆大學(xué) 數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院,新疆 烏魯木齊 830046)

0 引言

徑向基函數(shù)(Radial Basis Function,簡稱RBF)是一個取值僅依賴于離原點距離(本文使用歐氏距離)的實值函數(shù),即或者還可以是到任意點的距離,其中稱為中心點,即其研究的主要工作就是研究該函數(shù)張成的空間及其性質(zhì),并如何利用這個空間來解決一般對象的描述問題.徑向基函數(shù)是處理多元函數(shù)逼近問題的一種有效方法.事實上,它是通過定義在[0,∞)上的一元函數(shù)φ與Rd上的歐幾里德范數(shù)來表示元函數(shù)其中向量因此,用徑向基函數(shù)處理多元函數(shù)逼近問題具有效率高,運算簡單,易于編程,在計算機中儲存方便以及各向同性,不依賴于網(wǎng)格,求解精度高等優(yōu)點,其在函數(shù)逼近、偏微分方程數(shù)值解、神經(jīng)網(wǎng)絡(luò)、小波、多尺度分析以及地球物理學(xué)、測繪學(xué)、遙感與信號處理、圖像處理、機器學(xué)習(xí)、天文學(xué)與采礦等領(lǐng)域得到了廣泛的研究和應(yīng)用.

目前在解決實際問題中常用到的數(shù)值方法有:有限差分法、有限元法、有限體積法、邊界元法、無網(wǎng)格法等.雖然前四種方法解決了眾多的科學(xué)和工程計算問題,但它們存在一個共同的缺點:求解過程中每次都需要剖分網(wǎng)格,計算工作量大,尤其是三維問題.所以要想徹底解決這些方法面臨的網(wǎng)格重構(gòu)問題,就應(yīng)該避免使用網(wǎng)格.于是,無網(wǎng)格思想就被提出來了,并在近年來得到了迅速發(fā)展.目前,具有代表性的無網(wǎng)格方法主要有光滑粒子動力學(xué)法、無單元Galerkin法、重構(gòu)粒子法、單位分解法、RBF法等十余種,其中無單元Galerkin法、RBF法和重構(gòu)粒子法是無網(wǎng)格法研究的主流方向.RBF方法是以徑向函數(shù)為插值基函數(shù),在求解中采用配點形式,無需背景網(wǎng)格支持的方法.此方法的優(yōu)點是不受單元形狀的限制,在場量變化劇烈的地方,可以很容易地增加中心點的密度,因此避免了網(wǎng)格剖分在求解過程中的困難.

基于RBF的諸多優(yōu)點,國內(nèi)外研究人員進行了大量的科學(xué)研究并取得了豐碩成果,在此主要介紹RBF在計算數(shù)學(xué)上的應(yīng)用.徑向基函數(shù)的研究是從徑向基函數(shù)插值開始的,以下章節(jié)中會給出徑向基函數(shù)插值的基本理論結(jié)果.20世紀60年代,由Hardy[1]提出的MQ方法是一種無網(wǎng)格配置方法,目前它是徑向基函數(shù)無網(wǎng)格法的主流.Hardy最早將其用于散點曲面插值擬合.20世紀80年代初,Franke[2]對數(shù)十種散點插值法綜合比較,得出MQ方法是具有最好效果的方法.在之后的研究結(jié)果中,Micchelli[3]給出了RBF插值的可解性證明;Schaback[4]進一步研究了插值的存在性、唯一性和收斂性問題,使得RBF插值理論日趨成熟;此外,吳瑞豐[5]等人研究并構(gòu)造了高精度RBF擬插值算子;Kansa[6]把徑向基函數(shù)引入配點法中,首次成功地用以求解雙曲、橢圓和拋物型問題;此后,RBF方法在偏微分方程求解方面引起了廣泛關(guān)注,Franke[7]和Schaback首先建立了RBF方法在求解偏微分方程上的理論基礎(chǔ);Wu[8,9]主要從理論上研究了該方法的可解性和誤差限.但是,Dubal[10]和Fornberg[11]進一步研究了這類基于RBF的方法,發(fā)現(xiàn)隨著秩的增加,系數(shù)矩陣會變得病態(tài),為了解決這個問題,Kansa[12]提出用塊分解和LU分解格式取代全局求解器,從而避免病態(tài)矩陣的產(chǎn)生.此外,在前人工作的基礎(chǔ)上,張云新[13]和張淮清[14]分別系統(tǒng)地研究了RBF方法在不可壓縮流體和電磁場計算中的應(yīng)用.

根據(jù)定義域的不同可將徑向基函數(shù)分為全域型和緊支型,全域型RBF具有連續(xù)性和無限可微等優(yōu)點,在近似函數(shù)方面能獲得較高的精度,相對于移動最小二乘近似,其計算量大幅減小,對散點模型適應(yīng)性強,且對空間維數(shù)不敏感.這在高維空間(二維和三維)問題分析中優(yōu)勢更加突出.雖然RBF方法基函數(shù)中心點和配點設(shè)置靈活、程序編制簡便.但由于它是一種全域型方法,離散形成的代數(shù)矩陣是滿陣,故得到系數(shù)矩陣的條件數(shù)很大,不便于求解較大規(guī)模的問題.為了解決這一問題,出現(xiàn)了多種可行的方法,如優(yōu)化調(diào)整MQ徑向基函數(shù)的形狀參數(shù)[12,15]、緊支徑向基函數(shù)[16]以及區(qū)域分解法[17].近年來,吳宗敏[18]給出了正定緊支徑向基函數(shù)的構(gòu)造法則,并構(gòu)造出一系列正定緊支徑向基函數(shù);Buhmann[19]給出了基于TPS函數(shù)的正定緊支徑向基函數(shù).Wendland[20]給出的正定緊支徑向基函數(shù)與吳宗敏的類似,但它們具有最少的自由度.這些正定緊支RBF具有系數(shù)矩陣稀疏、帶狀分布的特點,有利于解決大型問題,但它們不滿足完備性條件,甚至不能正確描述常應(yīng)變狀態(tài).

在過去的幾十年里,徑向基函數(shù)方法取得了長足的發(fā)展,但在嚴格的數(shù)學(xué)論證、計算效率、邊界條件處理等方面還不能與成熟的有限元方法相媲美.雖然徑向基函數(shù)方法的理論研究還不成熟,但因它不需要背景網(wǎng)格支持,故在許多應(yīng)用領(lǐng)域中具有廣闊的發(fā)展前景.

在本文后面的章節(jié)中,主要討論RBF的定義、分類、插值和應(yīng)用,并著重介紹了MQ-RBF的發(fā)展歷程,理論知識及迭代自適應(yīng)MQ-RBF方法在間斷檢測上的應(yīng)用.由于MQ-RBF是全域型的,離散形成的系數(shù)矩陣是滿陣且條件數(shù)很大,為解決這一難題,將采用維數(shù)分裂的方法進行處理,該方法非常適合并行計算.

本文接下來的結(jié)構(gòu)框架如下:第二節(jié)主要介紹RBF的基本理論知識和插值方法,第三節(jié)主要介紹MQRBF的相關(guān)知識及其在間斷檢測上的應(yīng)用,第四節(jié)主要對RBF方法進行了總結(jié),并探討了其研究的發(fā)展趨勢和未來可能的研究方向.

1 徑向基函數(shù)

本節(jié)主要介紹RBF的一些基本概念、分類及插值方法[14,21].

1.1 基本概念

首先給出徑向函數(shù)(RF)和徑向基函數(shù)(RBF)如下形式的基本定義:

定義1[14]對于給定的徑向函數(shù)φ滿足:若則有其中φ是僅依賴于的函數(shù),這里為Euclidean范數(shù).

定義2[21]對于給定的和非負一元函數(shù)φ:R+→R,徑向基函數(shù)是由所有形如的函數(shù)及其線性組合張成的函數(shù)空間,其中為徑向函數(shù)的中心.

簡單地說,RBF就是一個距離基函數(shù),其變量為空間距離r,故具有形式簡單和各向同性等優(yōu)點.在特定條件下,只要選取兩兩不同,則就是線性無關(guān)的,故其可形成徑向基函數(shù)空間中的一組基,當選取幾乎充滿Rd時,及其線性組合便可逼近空間內(nèi)幾乎所有的函數(shù).

徑向基函數(shù)按定義域的不同可將其劃分為全域型和緊支型兩類,其中全域型函數(shù)在插值中具有較高的精度,但不善于求解高維問題;緊支型函數(shù)則由于求解矩陣為稀疏、帶狀分布的，故其計算量較小,適合大型問題的求解.常見的全域型RBF如表1所示.緊支型RBF可看作是全域型函數(shù)通過一個截斷多項式,在給定的連續(xù)性條件和自變量維數(shù)下截斷產(chǎn)生的,一般要求它在給定的維數(shù)空間內(nèi)正定和連續(xù).常見的緊支型RBF如表2所示.

表1 常見全域型RBF

表2 常見緊支型RBF

其中(1?r)+定義如下:

1.2 插值方法

插值法是一種非常重要且經(jīng)典的函數(shù)逼近方法,RBF的研究也是從插值開始的.插值過程是通過給定的離散數(shù)據(jù)節(jié)點,尋求一個能逼近未知函數(shù)的近似函數(shù),且要求近似函數(shù)在插值節(jié)點處的數(shù)值與原函數(shù)相等.根據(jù)給定問題條件的不同,將插值方法大體分為直接插值方法(包括拉格朗日插值和牛頓插值)和帶導(dǎo)數(shù)條件的Hermite插值方法.對于RBF插值,與之對應(yīng)的方法主要有Kansa插值、Hermite插值和擬插值等[14].RBF插值是在插值中選用徑向函數(shù)作為基函數(shù)的插值方式.它最早源于1951年Kringe在礦藏領(lǐng)域用于沉積插值的各向同性、穩(wěn)定的隨機函數(shù);隨后在1971年Hardy應(yīng)用多重二次插值方法,實現(xiàn)飛機外形曲面擬合問題并給出了MQ函數(shù);另外,1975年Duchon從樣條彎曲能最小的理論出發(fā),提出了多元問題的薄板樣條函數(shù)逼近方法.以上應(yīng)用于不同領(lǐng)域的插值方法都為形成RBF插值奠定了基礎(chǔ).RBF插值的數(shù)學(xué)語言表述如下:

給定個已知節(jié)點和未知函數(shù)且其在節(jié)點上的函數(shù)值為若想對未知節(jié)點進行插值,首先需找到簡單函數(shù)使得在已知節(jié)點滿足插值條件:

其中aj是待定系數(shù),則方程(1)的矩陣形式如下:

2 基于MQ-RBF插值的間斷檢測方法

圖像邊緣是以圖像局部特性的不連續(xù)性形式出現(xiàn)的.圖像邊緣檢測在圖像處理中是眾多科學(xué)家感興趣的課題,經(jīng)典的檢測方法是給出一個滿足初始條件的微分方程模型,使得初始曲線向著圖像邊緣移動.當?shù)竭_圖像邊緣時,曲線停止移動,這樣就找到了圖像的邊緣.數(shù)值擬合該微分方程常用的方法是水平集方法,該方法的優(yōu)點是在不知道檢測邊緣的前提下可以很自然地檢測到不同拓撲結(jié)構(gòu)的圖像邊緣.但其缺點是,求解高維問題會耗費大量的計算時間.而MQ-RBF邊緣檢測方法在檢測高維問題時,使用維數(shù)分裂的方法，這將在某種程度上降低求解的復(fù)雜度.

2.1 MQ函數(shù)發(fā)展歷程

MQ函數(shù)是Hardy于1968年提出的一種徑向基函數(shù).Hardy[22]總結(jié)了MQ函數(shù)自發(fā)現(xiàn)以來20年間所取得的研究成果,并列舉了MQ函數(shù)在眾多工程方面的應(yīng)用.他將MQ函數(shù)的發(fā)展歷程分為以下三個階段:

?1968-1972年是研究MQ函數(shù)基本方法的時期,在這一時期得到的比較實用的方法包括配置模型、最小二乘模型、密切模型、區(qū)域分解等;

?1972-1979年主要表現(xiàn)在MQ函數(shù)在各方面的廣泛應(yīng)用,如在大地測量學(xué)、地球物理學(xué)、測繪學(xué)、攝影測量學(xué)、遙感與信號處理、數(shù)字地形模型、水文學(xué)、地質(zhì)學(xué)與采礦等領(lǐng)域的應(yīng)用;

?1980-1988年是研究MQ-B函數(shù)及其數(shù)學(xué)論證的時期.

在這期間里,雖然MQ函數(shù)在微分方程和積分方程中有一定的應(yīng)用,但主要還是將MQ函數(shù)用于插值方面.Franke在文[2]中指出:就精度、穩(wěn)定性、有效性、內(nèi)存需要及易于實現(xiàn)而言,MQ函數(shù)在所有的29種散亂數(shù)據(jù)插值格式中是最好的.對于MQ函數(shù)及其逆MQ函數(shù)[23]研究了奇Buhmann數(shù)維Euclidean空間上插值的收斂階,得到如下結(jié)果:當n>3為奇數(shù)時,它們的收斂階分別為n+1和n?1階.此外他進一步研究了MQ函數(shù)單變量擬插值的收斂性問題.與此同時,他和Micchem[24]也研究了多變量整數(shù)網(wǎng)格上MQ函數(shù)的局部化性質(zhì)的改進問題.這在某種程度上豐富了MQ函數(shù)的理論.

2.2 迭代自適應(yīng)間斷檢測方法

圖像的邊緣檢測是圖像處理中一項重要的預(yù)處理技術(shù),是圖像分析與識別的重要環(huán)節(jié),也是處理許多復(fù)雜問題的關(guān)鍵,被廣泛應(yīng)用于圖像輪廓、特征的提取和紋理分析等方面.邊緣檢測法的種類有很多,如微分算子法、樣板匹配法、小波檢測法、神經(jīng)網(wǎng)絡(luò)法等.目前常用的高效邊緣檢測方法包括傅里葉共軛和增強共軛邊緣檢測方法,它們分別由Gelb[25]和Tadmor[26]提出,這兩種邊緣檢測方法都是基于傅里葉方法,并已經(jīng)成功地應(yīng)用于多種領(lǐng)域.以下將給出一種更加精確的邊緣檢測方法:基于MQ-RBF的ε迭代自適應(yīng)方法[26]

.雖然RBF方法在數(shù)值逼近間斷的地方會產(chǎn)生吉布斯現(xiàn)象[27],但可以利用該現(xiàn)象構(gòu)造一種更加精確的邊緣檢測方法.迭代自適應(yīng)MQ-RBF方法主要利用擴散系數(shù)不同的增長率,來確定間斷點的所在位置.即對于中心點總數(shù)N,如果形參固定,那么擴散系數(shù)的絕對值在間斷或間斷附近會達到局部最大并且會隨著N的增大呈指數(shù)增長.同時,在遠離間斷的地方,擴散系數(shù)的量級成指數(shù)級衰減;但是如果在間斷或間斷附近把形式參數(shù)賦空,擴散系數(shù)的量級將成線性衰減.這種迭代自適應(yīng)MQ-RBF方法應(yīng)用于邊緣檢測的成功之處在于擴散系數(shù)的增長率隨著形式參數(shù)自適應(yīng)的變化而顯著變化,并且無論所考慮的區(qū)域連續(xù)與否,擴散系數(shù)的變化都能容易的檢測到.在Jung[26,28]等人的研究工作基礎(chǔ)上,以一維和二維間斷問題為例來研究自適應(yīng)邊緣檢測方法.對于二維問題,采用維數(shù)分裂的方法進行求解,此方法的優(yōu)點是程序簡單易實現(xiàn),計算效率高.

2.2.1 一維迭代自適應(yīng)RBF方法

ε迭代自適應(yīng)方法的關(guān)鍵思想是:固定網(wǎng)格形式和N,MQ-RBF逼近的精度僅依賴于εj.理論上,當被逼近函數(shù)f(x)足夠光滑時,隨著εj的增加會得到好的收斂結(jié)果.但事實上,εj并不能隨意增大.當εj=0時,局部MQ-RBF就轉(zhuǎn)變?yōu)榉制€性函數(shù)φj=|x?xj|.但當所有的εj=0時,由于所有的基函數(shù)都是線性的,所以整體的逼近也呈線性收斂.由此,局部ε自適應(yīng)方法可簡單概括為:若中心點xj為間斷點,則εj=0;若中心點xj為非間斷點,則εj0.

特別地,若εj=0,在xi=xj,i=1,···,N時,令

現(xiàn)定義如下集合:

這里η是給定的已知數(shù).

給定:{εi},1>η>0,δ>0,aold

步驟1:由(2)和(5)式計算a,

步驟2:求得間斷點集合

步驟3:更新ε,即：若則εi=0;若為初始值,重新計算C;

步驟4:如果更新aold,重復(fù)步驟1→步驟3.

需要提及的是,對于集合C,可以只用一階導(dǎo)數(shù)或擴散系數(shù)(Ci=|ai|,i=1,···,N)來作為步驟2的判斷標準,但是,在某些特殊情況下，這二者都沒有步驟1里定義的判斷標準效果好.

現(xiàn)給出一個含有多個間斷點的一維算例,在區(qū)間[?1,1]上其間斷函數(shù)f(x)和間斷點[f]如下:

圖1 第一行是重構(gòu)函數(shù)(y),第二行是判斷標準C,第三行是形式參數(shù)ε,第四行是error =log10|(y)?(x)|

2.2.2 二維迭代自適應(yīng)RBF方法

由于整體方法需要求解的系數(shù)矩陣規(guī)模較大且為滿陣的方程組,因此,對于二維問題將采用維數(shù)分裂的方法來降低求解的規(guī)模,并提高計算效率.由文獻[26]得到,對于固定的ε=0.1,迭代自適應(yīng)方法穩(wěn)定條件就是求解的最小區(qū)域xmin和中心點個數(shù)N需滿足關(guān)系式:xmin≈0.006N?0.1.

以下將給出一些二維的算例,并且所有例子中取δ=1e-10,形式參數(shù)的初值都為ε=0.1.圖2中,從左到右依次給出一個大小為500×500的鐘表表盤的圖片以及x方向檢測的邊緣圖,y方向檢測的邊緣圖和兩個方向合起來檢測的結(jié)果,其中η=0.1,ε=0.1.圖3中原圖片大小為600×600,η=0.1.圖4中原圖片大小為465×500,η=0.03.圖5中原圖片大小為540×540,η=0.1.圖6中原圖片大小為554×435,η=0.1.圖7中原圖片大小為906×700,η=0.02.

需要說明的是,由于原始圖片本身有一些噪點,就是較弱的間斷點,所以在這些二維的例子中,其檢測出來的邊緣輪廓會有一些噪點,這是很正常的情況而不是說迭代自適應(yīng)方法不好;相反地說明用這個方法來檢測間斷點是非常精確的.

圖2 左圖是原圖,往右依次為x,y方向檢測邊緣圖和兩個方向合起來檢測的輪廓圖

圖3 左圖是一個阿貍的原圖,右圖是檢測的輪廓圖

圖4 左圖是一個蘋果的原圖,右圖是檢測的輪廓圖

圖5 左圖是一個二維碼的原圖,右圖是檢測的輪廓圖

圖6 左圖是一個指紋的原圖,右圖是檢測的輪廓圖

圖7 左圖是一個鬧鐘的原圖,右圖是檢測的輪廓圖

3 結(jié)論

通過本文的研究發(fā)現(xiàn),RBF方法具有不依賴于網(wǎng)格剖分、表達形式簡單、易于編程、計算效率和精度較高等優(yōu)點.但是,所要求解方程組的系數(shù)矩陣是滿陣,故對其進行求解時需花費較長時間,并且該方法在理論分析上不及有限差分方法、有限元方法等那么完善,但這反而為今后的研究工作提供了一個方向.在今后的工作中，將會嘗試研究高維情況下不同類型的RBF方法(如高斯型RBF,逆MQ-RBF等)在間斷檢測上的應(yīng)用.由于MQ-RBF產(chǎn)生的插值矩陣是滿陣,因此尋求一些高效的數(shù)值計算方法是很有必要的,如區(qū)域分解、并行計算等方法.此外,RBF方法是一種無網(wǎng)格法,故用它求解偏微分方程和反問題也將是今后研究的重點.

致謝本文撰寫的順利完成,作者要感謝美國Brown大學(xué)的曾維新教授和中國海洋大學(xué)的高振教授在相關(guān)題材方面進行的討論和合作,并寄予了許多寶貴的建議和精心的指導(dǎo),謹在此表示衷心的感謝.