基于B樣條的函數型數據表示與配準

孫 萍， 劉源昊， 王旭輝

(合肥工業大學 數學學院，合肥230601)

1 引 言

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

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

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

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

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

2 基于B樣條的函數型數據配準問題

2.1 函數型數據的近似表示

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

(1)

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

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

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

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

2.2 函數型數據的曲線配準問題

則

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

(2)

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

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

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

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

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

(3)

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

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

(i′)γ(s)的節點序列為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.

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

(4)

證由于復合函數的n階導函數公式為

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

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

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

(5)

2.3 曲線配準的算法步驟

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

3 算例研究

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

例2一組函數由下式給出：

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.

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

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

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

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

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

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

4 結 論

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

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

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

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