孫 萍， 劉源昊， 王旭輝

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

1 引 言

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

在函數型數據分析中，大多數樣本函數都呈現出典型的形狀特征序列(例如峰和谷)，如著名的伯克利生長數據集[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]開發(fā)了一種不需要對標志進行解釋的曲線配準技術.D. Telesca (2008)等[15]提出了貝葉斯層次曲線配準新模型.O. Collier (2015)等[16]提出了一個非參數測試框架，其中開發(fā)了廣義似然比測試以執(zhí)行曲線配準.Aishwarya Pawar (2019)等[17]提出一種基于動態(tài)水平集(使用截斷的分層B 樣條)的聯(lián)合圖像分割和配準.

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

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

2.1 函數型數據的近似表示

一個隨機變量在無窮維空間或函數空間取值，則稱其為一個函數型變量[22]，稱其觀測值為函數型數據，函數型數據一般以離散數據的形式存在.對一個連續(xù)的函數過程，將其表示為{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]提出一種動態(tài)規(guī)劃算法用于解決式(3)中所述的最小化問題，該動態(tài)編程算法通過動態(tài)調整離散參數格點方法得到扭曲函數.

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

(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.

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

(4)

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

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

例1給定節(jié)點向量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正態(tài)分布參數，圖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)中的參數越多，配準后函數階數也就越高，函數過于平滑，影響函數的配準效果.

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