一種面向多模態(tài)手術(shù)軌跡的快速無監(jiān)督分割方法

邵振洲，趙紅發(fā)，渠 瀛，施智平，關(guān) 永，袁慧梅

1(首都師范大學(xué) 信息工程學(xué)院，北京 100048)

2(首都師范大學(xué) 輕型工業(yè)機器人與安全驗證北京市重點實驗室，北京 100048)

3(成像技術(shù)北京市高精尖創(chuàng)新中心，北京 100048)

4(田納西大學(xué) 諾克斯維爾分校工程學(xué)院，田納西 美國 37996)

1 引 言

醫(yī)療手術(shù)軌跡分割技術(shù)是機器人輔助微創(chuàng)手術(shù)(RMIS)中的重要研究方向，在示范學(xué)習(xí)[1]、技能評估[2-4]和復(fù)雜任務(wù)自動化[5]等諸多領(lǐng)域有廣泛的應(yīng)用.以著名的“達(dá)芬奇”醫(yī)療手術(shù)機器人為例，在醫(yī)生手術(shù)過程中可以產(chǎn)生運動學(xué)和視頻記錄，每一個手術(shù)任務(wù)的過程軌跡(即運動學(xué)數(shù)據(jù)和視頻記錄)都可以分解為若干個可視為原始軌跡的子軌跡段(如圖1所示)，這些子軌跡段復(fù)雜度低、方差小、容易剔除異常值，可以通過相互組合來實現(xiàn)新的手術(shù)任務(wù).如何更加快速準(zhǔn)確地進(jìn)行手術(shù)軌跡分割在機器人輔助微創(chuàng)手術(shù)領(lǐng)域具有十分重要的研究意義.

圖1 手術(shù)軌跡分割示意圖Fig.1 Diagram of surgical trajectory segmentation

目前，解決手術(shù)軌跡分割問題的方法主要包括基于機器學(xué)習(xí)的有監(jiān)督和無監(jiān)督方法.其中有監(jiān)督的方法從人工標(biāo)注中學(xué)習(xí)子軌跡段到預(yù)定義軌跡段的匹配關(guān)系.例如，Lin等人[6]采用基于線性判別分析(LDA)的簡單貝葉斯分類器學(xué)習(xí)手術(shù)軌跡和標(biāo)注的分類匹配關(guān)系；受到隱馬爾科夫模型(HMM)在語音識別領(lǐng)域成功的啟發(fā)，文獻(xiàn)[7]提出了采用專家手術(shù)示范數(shù)據(jù)驅(qū)動的方法完成HMM的建模，對手術(shù)過程進(jìn)行分割和識別；然而，這類方法沒有考慮細(xì)粒度軌跡的臨近關(guān)系，所以捕獲手術(shù)軌跡中的長程狀態(tài)轉(zhuǎn)換[8]有更好的效果.但是這些模型在訓(xùn)練前需要預(yù)先完成人工標(biāo)注，而人工標(biāo)注是非常耗時的.于是為了擺脫對人工標(biāo)注的依賴，一些無監(jiān)督方法被提出，通常首先提取手術(shù)軌跡數(shù)據(jù)(包括運動學(xué)、視覺數(shù)據(jù))的特征，然后通過聚類的方法，估計出一系列的時間分割點或者分割子軌跡段.例如，Zhou等人[9]提出了一種無監(jiān)督分割方法，使用人為設(shè)計的特征(速度，方向，夾角)，采用基于LDA分類的分割方法，得到軌跡的轉(zhuǎn)折點，這種方法只考慮了軌跡數(shù)據(jù)中的部分信息.同樣是無監(jiān)督的分割方法，Sang[10]采用主成分分析(PCA)提取特征，然后利用高斯混合模型(GMM)對特征進(jìn)行聚類的分割方法達(dá)到了更好的效果.近年來，Krishnan等人[11]提出了一種對轉(zhuǎn)移狀態(tài)進(jìn)行聚類的方法(TSC)，通過構(gòu)造手術(shù)軌跡臨近時間序列的轉(zhuǎn)移表示，保留軌跡序列的局部時間相關(guān)性信息，提升了分割算法性能.最近，受到卷積神經(jīng)網(wǎng)絡(luò)(CNN)在計算機視覺領(lǐng)域獲得巨大成功的啟發(fā)，文獻(xiàn)[12]將上述工作擴(kuò)展到深度學(xué)習(xí)領(lǐng)域，提出一種采用VGG網(wǎng)絡(luò)提取手術(shù)視頻特征，然后結(jié)合運動學(xué)特征進(jìn)行分割的方法，提升了分割質(zhì)量.

然而，無監(jiān)督方法仍存在一些不足.首先，在視頻特征提取階段，由于視頻本身數(shù)據(jù)量大并且所采用的神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)復(fù)雜，提取的特征數(shù)據(jù)量很大(1Gb以上)，導(dǎo)致在視頻特征的提取與編碼階段時間耗費較長.其次，動作細(xì)致是微創(chuàng)手術(shù)的特點，其運動學(xué)數(shù)據(jù)(即末端執(zhí)行器的位置，方向，線速度，角速度，夾子的夾角)易受到噪聲干擾，由于醫(yī)生專業(yè)水平的差異性，難免存在一些手術(shù)偏差，這樣就導(dǎo)致一些噪聲的引入，即使同一個手術(shù)專家遠(yuǎn)程操作重復(fù)同一個手術(shù)任務(wù)，手術(shù)軌跡數(shù)據(jù)也有很大不同，對手術(shù)軌跡進(jìn)行平滑處理有助于分割.然而，已有的分割方法對運動學(xué)數(shù)據(jù)處理只是簡單的歸一化處理，沒有對其做進(jìn)一步的去噪處理.

針對上述問題，本文提出了一種快速手術(shù)軌跡無監(jiān)督分割方法(如圖2所示)，首先通過堆疊卷積自編碼器(SCAE)提取手術(shù)過程視頻特征，在將視頻數(shù)據(jù)映射到低維空間表示的同時保留視頻中的主要特征，其次采用小波對手術(shù)過程的運動學(xué)數(shù)據(jù)進(jìn)行多尺度濾波處理，平滑短程軌跡，得到運動學(xué)特征.為了減少由于采用單一特征進(jìn)行分割導(dǎo)致的信息缺失，將視覺特征及運動學(xué)特征進(jìn)行融合，然后對融合特征進(jìn)行時間軸加窗處理得到轉(zhuǎn)移表示，對轉(zhuǎn)移進(jìn)行聚類找出轉(zhuǎn)移狀態(tài).最后在視覺空間、運動學(xué)空間和時間上對轉(zhuǎn)移狀態(tài)進(jìn)行層次化聚類得到聚類簇，并對其中的異常簇進(jìn)行修剪，得到最終的分割結(jié)果分割子軌跡段.本文對涉及到的數(shù)學(xué)符號進(jìn)行了統(tǒng)一匯總?cè)绫?所示.

圖2 快速手術(shù)軌跡無監(jiān)督分割算法框圖Fig.2 Algorithm block diagram of fast surgical trajectory unsupervised segmentation

表1 數(shù)學(xué)符號意義Table 1 Meaning of mathematical symbols

2 基于SCAE的視覺特征提取

(1)

(2)

為了進(jìn)一步降低編碼器輸出特征的維度，在編碼器的最后一層(式3)和解碼器的第一層(式4)采用卷積核為1×1大小的卷積層，在降低維度的同時保留了空間信息：

(3)

(4)

圖3 SCAE網(wǎng)絡(luò)結(jié)構(gòu)示意圖Fig.3 Diagram of SCAE network structure

(5)

為了避免由于傳統(tǒng)轉(zhuǎn)置卷積帶來的棋盤效應(yīng)，本文采用雙線性插值的方法進(jìn)行上采樣，解碼器第l上采樣層可由式(6)計算.

(6)

網(wǎng)絡(luò)的訓(xùn)練通過最小化均方誤差L式(7)進(jìn)行.

(7)

視頻圖像特征提取網(wǎng)絡(luò)的偽代碼如下：

算法1.基于SCAE特征提取網(wǎng)絡(luò)訓(xùn)練與預(yù)測

網(wǎng)絡(luò)訓(xùn)練輸入：手術(shù)視頻圖像幀I

1. FORlin encoder layer

//編碼器卷積層，3×3卷積核

4. END FOR

//編碼器末層，1×1卷積核

//解碼器首層，1×1卷積核

7. FORlin decoder layer

//解碼器卷積層，3×3卷積核

10. END FOR

網(wǎng)絡(luò)預(yù)測輸出：特征圖H4

1. FORlin encoder layer

4. END FOR

3 基于視頻和運動學(xué)軌跡的層次化非參聚類

3.1 基于小波變換的運動學(xué)數(shù)據(jù)尺度平滑處理

小波變換可以將手術(shù)軌跡分解為低頻部分與高頻部分，高頻部分反映了軌跡在細(xì)節(jié)上的差異，而低頻部分可以看作是原始手術(shù)軌跡的平滑逼近[14].采用多尺度濾波處理手術(shù)軌跡可以在一定程度上消除噪聲，減少短程軌跡的影響.用雙通道多采樣率濾波器組表示的小波分解過程如圖4所示.

圖4 小波分解示意圖Fig.4 Diagram of wavelet decomposition

本文采用db10小波對運動學(xué)數(shù)據(jù)進(jìn)行5級濾波處理，以去除由于專家遠(yuǎn)程操作手術(shù)所引入的小尺度噪聲.小波濾波前后的右臂末端的位置隨時間的變化如圖5所示，縱坐標(biāo)為末端所處的位置以米(m)為單位，橫坐標(biāo)為時間以幀數(shù)表示(30幀/s)，x、y、z表示末端執(zhí)行器的位置，xw、yw、zw為對應(yīng)的小波濾波結(jié)果.

圖5 手術(shù)軌跡小波濾波前后對比圖Fig.5 Comparison of surgery trajectory before and after filtering by wavelet

3.2 層次化非參混合模型

本文采用狄利克雷高斯混合模型(DPGMM)進(jìn)行聚類.高斯混合模型是生成模型，考慮一個數(shù)據(jù)集X=xi，i=1…n，GMM模型將X視為從k個高斯組分中生成的，即：X～GMM(μi，σi)，i=1…k.而如何確定參數(shù)k是比較困難的，通常通過構(gòu)建狄里克雷混合模型進(jìn)行解決，即認(rèn)為GMM的參數(shù)是從狄里克雷過程產(chǎn)生：(μi，σi，k)～DP(α，H).DPGMM采用EM算法進(jìn)行擬合.

DPGMM仍需要定義最大組分個數(shù)kmax，本文在非參聚類過程中，首先采用K-means對輸入特征進(jìn)行聚類(k=kmax)，將聚類結(jié)果作為DPGMM的初始化，然后采用EM算法對DPGMM進(jìn)行細(xì)致的迭代優(yōu)化，能夠在一定程度上減少迭代次數(shù).

結(jié)合視覺特征和運動學(xué)特征的快速手術(shù)軌跡無監(jiān)督分割算法偽代碼如下：

算法2.快速手術(shù)軌跡無監(jiān)督分割算法

輸入：手術(shù)視頻圖像幀I，手術(shù)軌跡運動學(xué)數(shù)據(jù)K

輸出：一組時間分割點Segs

//使用訓(xùn)練好的SCAE提取視覺特征，見算法1

1.V←SCAEtest(I)

2.X←wavelet(K)//使用小波變換處理運動學(xué)數(shù)據(jù)

//K-means對轉(zhuǎn)移進(jìn)行聚類，計算聚類中心cs

4.cs←Kmeans(D)

//DPGMM對轉(zhuǎn)移進(jìn)行聚類，用cs初始化

5.CD←DPGMM(D，cs)

6.S←Dt；?t，Ct≠Ct+1//轉(zhuǎn)移狀態(tài)識別

//Kmeans初始化的DPGMM在視覺空間聚類

7.CSV←KDPGMM(VS)

//Kmeans初始化的DPGMM在運動學(xué)空間聚類

8.CSVX←KDPGMM(XC)；?c∈CSV

//Kmeans初始化的DPGMM在時間軸聚類

9.CSVXT←KDPGMM(T)；T=t(c)；?c∈CSVX

//計算聚類簇中心

10.Segs←mean(T(c))；?c∈CSVXT

4 實驗結(jié)果與分析

為了評估快速手術(shù)軌跡無監(jiān)督分割算法的性能，本文進(jìn)行了2組實驗.第一組實驗主要驗證本算法自身的通用性和探究使用視頻特征的必要性.在第二組實驗中，將本文算法與同樣使用了視頻數(shù)據(jù)和運動學(xué)數(shù)據(jù)進(jìn)行分割的TSC-DL[12]，以及僅使用運動學(xué)數(shù)據(jù)進(jìn)行分割的TSC[11]，GMM[10]，HMMGMM[15]，HSMM[16]進(jìn)行了軌跡分割性能對比.

本文實驗采用約翰霍普金斯大學(xué)提供的JIGSAWS數(shù)據(jù)集[17].該數(shù)據(jù)集采集自達(dá)芬奇通用醫(yī)療機器人系統(tǒng)，共分為視頻數(shù)據(jù)和運動學(xué)數(shù)據(jù)兩部分.其中視頻數(shù)據(jù)為固定攝像機對手術(shù)區(qū)域錄制的640×480大小的視頻；運動學(xué)數(shù)據(jù)共76維，包括遠(yuǎn)程端雙臂的末端執(zhí)行器38維數(shù)據(jù)和病人端雙臂末端執(zhí)行器38維數(shù)據(jù)；38維數(shù)據(jù)(每個臂19維)包括：雙臂末端的3維笛卡爾坐標(biāo)、9維表示的旋轉(zhuǎn)矩陣、3維的速度、3維的角速度和1維的夾持器角度.視頻和運動學(xué)數(shù)據(jù)采樣頻率均為30Hz.該數(shù)據(jù)集包括不同熟練程度的8名外科醫(yī)生完成的三種手術(shù)任務(wù)：縫合，針傳遞和打結(jié).

表2 實驗環(huán)境配置

本實驗選取JIGSAWS數(shù)據(jù)集的一個子集進(jìn)行實驗，共包括2種手術(shù)任務(wù)(針傳遞和針縫合)，每種手術(shù)任務(wù)包括5個專家(E)，3個中級專家(I)，3個非專家(N)的數(shù)據(jù).運動學(xué)數(shù)據(jù)只使用病人端的38維數(shù)據(jù).對運動學(xué)數(shù)據(jù)和視覺數(shù)據(jù)進(jìn)行3倍降采樣以降低計算量.在特征提取之前分別為針傳遞和針縫合訓(xùn)練一個視頻特征提取網(wǎng)絡(luò)(SCAE)，SCAE的輸入為640×480的彩色圖像，輸出為一個70維的特征向量，采用Adam[18]進(jìn)行訓(xùn)練.本實驗的機器配置見表2.

4.1 評估方法

本文采用歸一化互信息(NMI)和整體運行時間對算法進(jìn)行評估.

歸一化互信息：給定預(yù)測聚類結(jié)果A和真實標(biāo)記B，歸一化互信息用于衡量A，B兩個聚類結(jié)果的相似度，可由式(8)計算：

(8)

其中H(A)，H(B)分別為A，B的信息熵，I(A，B)為A和B的互信息，NMI范圍在0-1之間，0表示兩個聚類結(jié)果沒有相關(guān)關(guān)系，而1表示完全相關(guān).

整體運行時間：統(tǒng)計分割所有手術(shù)示范過程中特征提取及進(jìn)行聚類分割的總運行時間，沒有使用視頻數(shù)據(jù)的只統(tǒng)計聚類分割的運行時間.

4.2 本文提出方法的性能評估

為了探究使用視頻數(shù)據(jù)對本文方法的影響，本文在同一數(shù)據(jù)集上對僅使用運動學(xué)數(shù)據(jù)、僅使用視頻數(shù)據(jù)和使用運動學(xué)數(shù)據(jù)及視頻數(shù)據(jù)三種軌跡分割方法進(jìn)行了實驗，NMI見表3.使用的數(shù)據(jù)集是JIGSAWS數(shù)據(jù)集的子集，包括全專家E(5個專家示范)，專家及中級專家E+I(3個專家示范和3個中等專家示范)，所有級別專家E+I+N(3個專家示范、3個中等專家示范及3個非專家示范).

表3 采用不同特征進(jìn)行分割的NMI對比Table 3 NMI of segmentation using different features

對于縫合任務(wù)，從表3中，我們可以看出，結(jié)合運動學(xué)和視頻數(shù)據(jù)進(jìn)行分割的NMI比僅使用運動學(xué)數(shù)據(jù)提高了至少4.7%；對于僅使用運動學(xué)數(shù)據(jù)來說，NMI隨著非專家示范的比例增長呈現(xiàn)逐漸下降的趨勢，是因為專家示范通常比非專家完成的更加順利和迅速.然而在使用運動學(xué)數(shù)據(jù)和視頻數(shù)據(jù)進(jìn)行分割時，卻沒有一致的趨勢.我們認(rèn)為這與縫合任務(wù)的特殊性和新手的一些不一致的錯誤有關(guān).

對于針傳遞任務(wù)，從表3中，我們可以看出，結(jié)合運動學(xué)和視頻數(shù)據(jù)進(jìn)行分割的NMI比僅使用運動學(xué)數(shù)據(jù)提高了至少37.9%，這比縫合任務(wù)要顯著很多.與縫合任務(wù)相似，對于僅使用運動學(xué)數(shù)據(jù)的分割方法來說，NMI隨著非專家示范的比例增長呈現(xiàn)逐漸下降的趨勢.然而在使用運動學(xué)數(shù)據(jù)和視頻數(shù)據(jù)進(jìn)行分割時，與縫合任務(wù)不同，針傳遞任務(wù)NMI表現(xiàn)出與僅使用運動學(xué)數(shù)據(jù)進(jìn)行分割相同的趨勢.我們猜測這是因為針傳遞任務(wù)的軌跡相比于縫合任務(wù)的軌跡更加規(guī)則.

圖6對無監(jiān)督聚類的分割結(jié)果進(jìn)行了可視化，從圖中我們發(fā)現(xiàn)使用視覺和運動學(xué)數(shù)據(jù)進(jìn)行分割(圖6(b)，圖6(d))的效果要好于僅使用運動學(xué)數(shù)據(jù)進(jìn)行分割(圖6(a)，圖6(c))；僅使用運動學(xué)數(shù)的分割方法通常會導(dǎo)致過度分割，從而導(dǎo)致一些分割片段無法被正確分割；而結(jié)合視覺和運動學(xué)數(shù)據(jù)的分割方法結(jié)合了兩種互補的信息，能夠從整體進(jìn)行分割，而減少過度分割的發(fā)生.

4.3 不同軌跡分割方法的對比

實驗中分別采用HSMM，HMMGMM，GMM，TSC，TSC-DL以及本文算法對同一手術(shù)示范數(shù)據(jù)集進(jìn)行分割；手術(shù)示范數(shù)據(jù)集包括全專家E(5個專家示范)，專家及中級專家E+I(3個專家示范和3個中等專家示范)，所有級別專家E+I+N(3個專家示范、3個中等專家示范及3個非專家示范).不同方法的分割精度和運行時間如表4和表5所示.表5對不同分割方法的總運行時間進(jìn)行了統(tǒng)計，對于僅使用運動學(xué)數(shù)據(jù)的分割方法(HSMM、HMMGMM、GMM、TSC)，運行時間為其聚類分割的運行時間.對于使用視覺數(shù)據(jù)和運動學(xué)數(shù)據(jù)的分割方法(TSC-DL、本文方法)，運行時間為視頻特征提取和聚類分割運行時間的總和.

表4 不同分割方法的NMI對比Table 4 NMI of segmentation using different approachs

從表4中我們可以看出，使用了視頻數(shù)據(jù)的軌跡分割方法(本文算法和TSC-DL)比沒有使用視頻數(shù)據(jù)的軌跡分割方法(HSMM，HMMGMM，GMM，TSC)的NMI有明顯的提升.在所有的軌跡分割任務(wù)中，本文的方法達(dá)到了最好的NMI.在針傳遞任務(wù)中，本文方法比同樣使用了視頻數(shù)據(jù)的TSC-DL算法的NMI提高了8.9%-18.7%；在針縫合任務(wù)中，本文方法比TSC-DL算法的NMI提高了2.0%-7.5%.因為本文采用小波對運動學(xué)數(shù)據(jù)進(jìn)行多尺度平滑處理，濾除小尺度噪聲而降低了短程軌跡的干擾，從而間接提升了整體分割準(zhǔn)確性；基于卷積自編碼的特征提取器能夠提取圖像中的主要特征，而不受細(xì)節(jié)干擾.

我們也注意到在針傳遞任務(wù)中，大多數(shù)算法(包括本文算法)的NMI隨著非專家手術(shù)示范比例的增加都有相同的變化趨勢，即NMI隨著非專家手術(shù)示范比例的增加而越來越低，其中E+I的比僅有E的要好一些的原因可能是選取的數(shù)據(jù)集有一定特殊性，而TSC-DL因為采用了基于專家等級進(jìn)行加權(quán)的聚類簇裁剪機制，能夠在非專家手術(shù)示范比例的增加的情況下表現(xiàn)的越來越好.

表5 不同分割方法運行時間對比(單位：秒)Table 5 Comparison of running time using different segmentation methods(unit：s)

圖7 不同分割方法運行時間柱狀圖Fig.7 Bar graph of the running time of segmentation using different methods

從圖7中可以明顯看出，結(jié)合視頻數(shù)據(jù)和運動學(xué)數(shù)據(jù)的分割方法都會比僅使用運動學(xué)數(shù)據(jù)的分割方法慢10倍左右，因為結(jié)合視頻數(shù)據(jù)和運動學(xué)數(shù)據(jù)的軌跡分割方法往往需要從視頻中提取特征，然而這一步是比較耗時的.基于層次化聚類的分割方法(TSC-DL、TSC、本文方法)比基于單次聚類的分割方法(GMM、HMMGMM、HSMM)要慢.

在結(jié)合視頻數(shù)據(jù)和運動學(xué)數(shù)據(jù)的分割方法中，本文提出的方法幾乎是TSC-DL速度的10倍，甚至優(yōu)于僅使用運動學(xué)數(shù)據(jù)進(jìn)行分割的TSC算法，因為本文方法采用了一個更簡單的無監(jiān)督深度模型對視頻數(shù)據(jù)進(jìn)行特征提取，大大降低了視頻特征提取階段的時間消耗.盡管采用K-means對混合模型進(jìn)行初始化有助于速度的提升，但是提升有限.

本實驗表明，相較于TSC-DL方法，本文方法能夠?qū)⒔Y(jié)合視覺特征和運動學(xué)特征進(jìn)行分割的速度提升10倍，并且分割結(jié)果得到較高的NMI.

5 總 結(jié)

本文針對手術(shù)軌跡分割中視頻空間的高維特征和分割效果容易受到短程軌跡干擾的問題，提出了一種快速手術(shù)軌跡無監(jiān)督分割方法.該算法通過堆疊卷積自編碼器提取手術(shù)軌跡視頻特征，在對視頻特征進(jìn)行降維的同時學(xué)習(xí)到視頻的主要特征，然后結(jié)合小波多尺度濾波處理后的運動學(xué)軌跡進(jìn)行聚類分割，減少了短程軌跡對分割的干擾.在JIGSAWS數(shù)據(jù)集上實驗表明，該算法能夠在保證分割精度的情況下，能夠有效縮短分割時間.該分割算法可以大幅縮短外科醫(yī)生的技能評估用時、推進(jìn)機器人手術(shù)自動化，在機器人輔助醫(yī)療領(lǐng)域具有廣泛的應(yīng)用場景.