基于B樣條的函數(shù)型數(shù)據(jù)表示與配準

孫 萍， 劉源昊， 王旭輝

(合肥工業(yè)大學 數(shù)學學院，合肥230601)

1 引 言

函數(shù)型數(shù)據(jù)概念最初是由J. O. Ramsay[1]和C. J. Dalzell[2-3]提出的.函數(shù)型數(shù)據(jù)是指在緊區(qū)域上實值函數(shù)的實現(xiàn)，其將觀測到的離散數(shù)據(jù)當作一個整體來考慮.從函數(shù)型數(shù)據(jù)角度，可用一個連續(xù)的函數(shù)來表示觀測數(shù)據(jù)，即離散數(shù)據(jù)來自一條光滑的曲線.

在函數(shù)型數(shù)據(jù)分析中，大多數(shù)樣本函數(shù)都呈現(xiàn)出典型的形狀特征序列(例如峰和谷)，如著名的伯克利生長數(shù)據(jù)集[4-6]，為研究生長曲線的形狀特征，Laura M. Sangalli[7]、Yu-Hsiang Cheng[8]等人引入了曲線配準方法來聚類和對齊生長曲線.該方法允許振幅和相位變化的分離，來研究生長的階段性特點與生長曲線形狀特征.振幅和相位為曲線的形狀特征，相位變化由扭曲函數(shù)(warping function)表示.

通過曲線對齊和扭曲函數(shù)復合來消除相位變化的過程稱為曲線配準[9].配準的目的是分解函數(shù)型數(shù)據(jù)曲線樣本的幅度和相位變化.相位變化通過扭曲函數(shù)捕獲，扭曲函數(shù)單調(diào)的變換函數(shù)域[10].配準后的曲線應僅顯示振幅的變化，使曲線特征突出.文獻[11]中主要講述了函數(shù)數(shù)據(jù)的配準或注冊的特定方面，討論如何處理兩個或多個函數(shù)的配準問題，提出將扭曲函數(shù)作為工具，以水平調(diào)整觀察到的函數(shù)，減少它們的方差.數(shù)學上，研究2D、3D和更高維度的曲線形狀，通過配準使函數(shù)形狀的波峰波谷對齊，為研究函數(shù)形狀帶來便利，因此配準問題對函數(shù)型數(shù)據(jù)的形狀分析至關重要.

曲線配準問題已被國內(nèi)外許多學者研究.A. Kneip(1995)等[12-13]深入研究了標志配準，將標志稱為結(jié)構(gòu)性特征，位置稱為結(jié)構(gòu)點；作為標志配準的替代方法，Silverman(1995)[14]開發(fā)了一種不需要對標志進行解釋的曲線配準技術.D. Telesca (2008)等[15]提出了貝葉斯層次曲線配準新模型.O. Collier (2015)等[16]提出了一個非參數(shù)測試框架，其中開發(fā)了廣義似然比測試以執(zhí)行曲線配準.Aishwarya Pawar (2019)等[17]提出一種基于動態(tài)水平集(使用截斷的分層B 樣條)的聯(lián)合圖像分割和配準.

本文將函數(shù)型數(shù)據(jù)用B樣條基函數(shù)擬合，充分利用了B 樣條基函數(shù)自由度高，靈活性強等特點[18-21].并構(gòu)造曲線配準能量項作為標準，從而尋找扭曲函數(shù)對擬合后曲線進行曲線配準.需要指出的是，為了將整個函數(shù)型數(shù)據(jù)配準問題統(tǒng)一在B樣條函數(shù)框架上，扭曲函數(shù)也選擇為B樣條函數(shù).

2 基于B樣條的函數(shù)型數(shù)據(jù)配準問題

2.1 函數(shù)型數(shù)據(jù)的近似表示

一個隨機變量在無窮維空間或函數(shù)空間取值，則稱其為一個函數(shù)型變量[22]，稱其觀測值為函數(shù)型數(shù)據(jù)，函數(shù)型數(shù)據(jù)一般以離散數(shù)據(jù)的形式存在.對一個連續(xù)的函數(shù)過程，將其表示為{f(t),t∈I}，一維情形時，I?1，二維情形時，I?2.容量為N的樣本記為fi(t),i=1,2,…,N，它們獨立且與f(t)同分布[23]. 觀測到數(shù)據(jù)對(tij,fi(tij))，其中tij,j=1,2,…,ni為觀測點，且對函數(shù)fi(t)觀測ni次.獲得的數(shù)據(jù)對為離散數(shù)據(jù)，隨著函數(shù)個數(shù)和觀察次數(shù)的增加，函數(shù)數(shù)據(jù)趨于無窮維.故在函數(shù)型數(shù)據(jù)分析過程中，函數(shù)具有高維的特征，需要對函數(shù)型數(shù)據(jù)進行降維.常見的方法是利用基函數(shù)表示函數(shù)型數(shù)據(jù)[24].設有一組基函數(shù)φi(t),i=0,…,l，則函數(shù)f(t)在函數(shù)空間span{φ0(t),φ1(t),…,φl(t)}上的估計可表示為

(1)

其中k=1,2,…,p. 特別地，當p=3時，U=[u0,u1,…,un+4].為了保證端點插值性質(zhì)，取

u0=u1=u2=u3,un+1=un+2=un+3=un+4，

系數(shù)ci,j,j=0,…,n可由最小二乘法計算，即

(ci0,…,cin)T=(ATA)-1(fi(t0),…,fi(tn0))T，

2.2 函數(shù)型數(shù)據(jù)的曲線配準問題

則

故本文中，不妨設給定函數(shù)f1(s),f2(t),s,t∈[0,1]. 成對配準問題是指尋找恰定的扭曲函數(shù)t=γ(s)，使f1,f2°γ在某種意義下誤差極小.即

(2)

其中E為某種能量函數(shù)，Γ為滿足一定條件的函數(shù)集合.函數(shù)γ(s)須滿足：

(i)γ∶[0,1]→[0,1];

(ii)γ(0)=0,γ(1)=1;

(iii)γ可逆，且γ和γ-1為光滑函數(shù).

設ΓI表示所有滿足條件(i)-(iii)的函數(shù)集合，I為γ的定義域.可知ΓI是一個李群[11]，即對任何γ1,γ2∈ΓI，有群算子(γ1°γ2)(s)=γ1(γ2(s)).一般情況下，可將能量函數(shù)取為L2范數(shù)，即問題(2)轉(zhuǎn)換為

(3)

文獻[11]提出一種動態(tài)規(guī)劃算法用于解決式(3)中所述的最小化問題，該動態(tài)編程算法通過動態(tài)調(diào)整離散參數(shù)格點方法得到扭曲函數(shù).

由2.1節(jié)的函數(shù)型數(shù)據(jù)表示方法，可設f1(s),f2(t)的B樣條函數(shù)近似為

(i′)γ(s)的節(jié)點序列為V=[0,0,0,s3,…,sm,1,1,1],其中s3>0,sm<1,si+l>si,i=3,…,m-1;

(ii′)r0=0,r1=1,ri≤ri+1,i=0,…,m-1.

根據(jù)樣條函數(shù)的性質(zhì)[27]，條件(i′) (ii′)保證了γ(s)滿足條件(i),且在區(qū)間[0,1]上單調(diào)遞增,進而為可逆函數(shù)，此外γ(0)=0,γ(1)=1. 記滿足條件(i′) (ii′)的2次B樣條函數(shù)為Π. 則在B樣條的表示框架下，成對函數(shù)型數(shù)據(jù)的匹配問題(3)可轉(zhuǎn)換為極小化問題：

(4)

證由于復合函數(shù)的n階導函數(shù)公式為

下面給出定理1的一個算例.

例1給定節(jié)點向量U=[0,0,0,0,1/3,2/3,1,1,1,1]，定義在U上的三次B樣條函數(shù)為

其中c0=0,c1=2,c2=3,c3=1,c4=2,c5=0. 扭曲函數(shù)γ(s)為

(5)

2.3 曲線配準的算法步驟

根據(jù)上述的分析，給出基于B樣條的曲線配準的具體步驟如下：

3 算例研究

為驗證本文算法的有效性，分別用文獻[30]中使用的數(shù)據(jù)和fdasrvf庫函數(shù)中的數(shù)據(jù)運行算法.通過調(diào)整函數(shù)參數(shù)獲取多組數(shù)據(jù)進行驗證本文算法，例如算例2-4來自同一函數(shù)擬合數(shù)據(jù).

例2一組函數(shù)由下式給出：

yi(t)=zi,1e-(t-1.0-λ1)2/2+zi,2e-(t+2.0-λ2)2/2,t∈[-3,3],i=1,2,…,N.

對于上述有明顯極值點的函數(shù)，可以比較極值點的方差，計算原始函數(shù)配準后極值點方差分別為0.00786645，0.0315444，0.00883703，配準前極值點方差分別為0.0250637，0.084134，0.0140217.極值點的方差配準后明顯減小，配準后函數(shù)在波峰和波谷一定程度集中，即函數(shù)配準有明顯效果.

例3對于例2，可以通過改變zi,1和zi,2的均值和方差或改變λ1和λ2的均值和方差來改變形狀，在保持其他不變的情況下，將λ1和λ2的方差分別設置為0.5和0.2.其與圖2相對應的圖形如圖3.通過計算，配準后的L2距離小于配準前的L2距離，配準效果顯著.圖3d中藍色實線帶標記X表示升階后f1(s)的控制系數(shù)的折線圖，圖3c中配準函數(shù)的控制點與f1(s)升階后的控制點距離達到最小.

例4調(diào)整例2中λ1，λ2正態(tài)分布參數(shù)，圖4a表示調(diào)整后的函數(shù)，圖4b是其配準后的函數(shù)，配準后的極值點方差分別為0.00322853、0.000136298、0.022294，配準前的極值點方差分別為0.0294151、0.0187822、0.0891169，極值點的方差明顯減小，可以判斷配準后函數(shù)在波峰明顯較配準前集中.計算函數(shù)配準前后的L2距離，配準后整體呈減小趨勢.

例5從fdasrvf庫函數(shù)中獲取模擬數(shù)據(jù)集，圖5a為原始數(shù)據(jù)B樣條擬合后的形狀，圖5b中虛線是3階B樣條基擬合γ(s)與圖5a函數(shù)配準后結(jié)果.調(diào)整擬合γ(s)的B樣條階數(shù)，圖5c、圖5d分別為4階、5階B樣條基擬合的γ(s)與原始函數(shù)配準結(jié)果.

通過計算配準前后的極值點方差和配準前后函數(shù)與f1(s)的L2距離，可以得到，函數(shù)配準后的極值點方差都是減小的.當γ(s)是3階的，其配準前后的極值點方差分別為32.3751，7.50638；當γ(s)是4階的，其配準前后的極值點方差分別為32.3751，4.63976.當γ(s)是5階的，配準后的函數(shù)形狀趨于離散，極值點的范圍更廣，無法確定極值點.隨著γ(s)階數(shù)的增加，配準后函數(shù)在一定程度較配準前更集中，但γ(s)的階數(shù)并不是越大越好.

通過算例驗證分析，得出以下幾點結(jié)論：

4 結(jié) 論

針對基于B樣條的函數(shù)型數(shù)據(jù)配準問題，本文結(jié)合B 樣條基函數(shù)擬合函數(shù)型數(shù)據(jù)，扭曲函數(shù)γ(s)同樣基于B樣條函數(shù)，配準后的函數(shù)相當于對原函數(shù)升階，將配準函數(shù)與原始直接升階的函數(shù)f1(s)進行比較.

(1) 基于B樣條函數(shù)擬合的函數(shù)型數(shù)據(jù)，對其進行函數(shù)配準，扭曲函數(shù)γ(s)滿足B樣條基函數(shù)結(jié)構(gòu)，函數(shù)配準后使其在波峰波谷更集中；

(2) 在函數(shù)型數(shù)據(jù)的配準算例中，扭曲函數(shù)γ(s)的B樣條基的次數(shù)對配準函數(shù)有一定影響.一般情況是越大越好，但到特定階之后，γ(s)的階數(shù)越高，使得γ(s)中的參數(shù)越多，配準后函數(shù)階數(shù)也就越高，函數(shù)過于平滑，影響函數(shù)的配準效果.

致謝作者非常感謝相關文獻對本文的啟發(fā)以及審稿專家提出的寶貴意見.