王凱龍,常兆華,葉 萍
(上海理工大學(xué)健康科學(xué)與工程學(xué)院,上海 200082)
目前在全球范圍內(nèi),肺癌仍是癌癥中死亡率最高的疾病之一[1],近年來以滅活腫瘤細胞為目的的冷凍消融治療手段引起了人們的廣泛關(guān)注。該治療手段需在CT 引導(dǎo)下做氬氦刀經(jīng)皮穿刺,精準(zhǔn)穿刺是治療的重要條件,但胸腹部腫瘤的位置會隨著患者的呼吸發(fā)生運動,這種運動會影響穿刺的精準(zhǔn)度,在有些情況下腫瘤運動會達到28mm[2]。研究表明,胸腹部外部呼吸信號與肺部腫瘤運動具有良好的相關(guān)性,因此利用外部信號跟蹤肺部腫瘤位置是可行的[3]。但是由于外部運動信號的獲取、軟件內(nèi)部算法的運行以及機械臂延遲而導(dǎo)致的系統(tǒng)延遲會影響呼吸跟蹤的實時性[4],通常這種系統(tǒng)延遲在200ms 左右[5-7],這會降低機器人的穿刺精度。Ma 等[8]研究表明,利用呼吸信號進行預(yù)測的方法能夠有效提升穿刺機器人的穿刺精度,他們利用呼吸預(yù)測模型對公共數(shù)據(jù)集中的呼吸信號進行預(yù)測,并用運動體模進行穿刺實驗,結(jié)果表明,預(yù)測結(jié)果的平均均方根誤差(Root Mean Squared Error,RMSE)為0.58mm,穿刺誤差的平均均方根(Root Mean Square,RMS)為0.65 mm。
有關(guān)呼吸跟蹤和呼吸預(yù)測的研究在放療中展開較早,Hoogeman 等[9]的研究表明,使用呼吸跟蹤和呼吸預(yù)測技術(shù)能夠有效地減少由呼吸運動和系統(tǒng)延遲引起的幾何誤差。比較經(jīng)典的呼吸預(yù)測算法有:CyberKnife?的Synchrony?呼吸同步跟蹤系統(tǒng)使用的預(yù)測算法、支持向量機(Support Vector Machine,SVM)、自回歸移動平均模型(Autoregressive Integrated Moving Average Model,ARIMA)和神經(jīng)網(wǎng)絡(luò)。
Synchrony?系統(tǒng)在臨床實踐中表現(xiàn)良好,其利用兩個正交的kV 級X 射線系統(tǒng)等中心地投射到患者的治療部位觀測內(nèi)部腫瘤運動并建立與外部光學(xué)標(biāo)記相關(guān)的運動模型,然后使用呼吸預(yù)測技術(shù)補償固有的系統(tǒng)延遲[10-11]。Synchrony?呼吸同步跟蹤系統(tǒng)應(yīng)用了一種混合預(yù)測算法,該混合算法結(jié)合了線性算法、模式匹配和模糊邏輯3 種方法[12-13]。支持向量回歸(Support Vector Regression,SVR)是SVM 的一個應(yīng)用分支,Riaz 等[14]利用SVR 方法對腫瘤運動進行預(yù)測,在延遲為1s 的狀態(tài)下預(yù)測結(jié)果的RMSE 在2mm以內(nèi)。Khashei 等[15]研究表明,ARIMA 在時間序列預(yù)測方面表現(xiàn)良好,Ichiji 等[16]使用ARIMA 結(jié)合多層感知器(Multilayer Perceptron,MLP)對腫瘤運動進行預(yù)測,其結(jié)果顯示平均絕對誤差(Mean Absolute Error,MAE)在500ms 的延遲下為0.795 3±0.024 3 mm。Wang等[17]應(yīng)用七層雙向長短期記憶神經(jīng)網(wǎng)絡(luò)(Bidirectional Long Short-term Memory,Bi-LSTM)預(yù)測呼吸運動信號,設(shè)定延遲為400ms 的情況下其預(yù)測結(jié)果的平均RMSE為0.097 mm、MAE為0.074 mm。
很多相關(guān)研究實驗所用的呼吸數(shù)據(jù)樣本數(shù)量過少,或僅截取呼吸平穩(wěn)的區(qū)段進行實驗,無法測試模型的魯棒性。本文使用了29 組呼吸數(shù)據(jù)進行實驗,在保證樣本數(shù)量的同時不特意截取呼吸平穩(wěn)區(qū)段,使實驗結(jié)果更貼近臨床需求。本文提出的基于Bi-LSTM 的呼吸預(yù)測模型框架不但具有很高的準(zhǔn)確性,而且在面對各種復(fù)雜呼吸狀況時也表現(xiàn)出了較好的魯棒性。
實驗中,使用NDI Polaris 光學(xué)定位系統(tǒng)采集了23 組志愿者的外部呼吸數(shù)據(jù),該系統(tǒng)通過追蹤仰臥狀態(tài)下的志愿者胸腹部體表的紅外光學(xué)定位支架以獲取呼吸數(shù)據(jù),系統(tǒng)采樣頻率設(shè)定為20Hz,每組的采樣時間為5~7min,采集到的每組數(shù)據(jù)包括呼吸運動的左右、頭腳、前后3 個維度。志愿者在呼吸數(shù)據(jù)采集過程中完全處于自由呼吸狀態(tài),因此數(shù)據(jù)中包含了基線漂移、呼吸頻率改變、深呼吸等狀態(tài)。此外,實驗還加入了德國Lubeck 大學(xué)機器人與認知系統(tǒng)研究所公開數(shù)據(jù)庫中的6 組呼吸數(shù)據(jù),該數(shù)據(jù)的志愿者均在肺部有惡性腫瘤。
為了方便實驗結(jié)果的比較并保證實驗數(shù)據(jù)可靠,本文中統(tǒng)一選取前后方向的數(shù)據(jù)進行預(yù)測,且并不特意截取呼吸數(shù)據(jù)的平滑段進行試驗,保留原數(shù)據(jù)的諸多不利于預(yù)測精度的狀況,29組數(shù)據(jù)均截取6 000個數(shù)據(jù)點進行實驗。
采集的原始呼吸數(shù)據(jù)的噪聲嚴重影響了之后數(shù)據(jù)訓(xùn)練效果,因此實驗開始時需對每條數(shù)據(jù)進行平滑降噪處理,在保證呼吸數(shù)據(jù)原有特性不變的前提下降低信噪比。本文選用Savitzky-Golay 濾波器對原始數(shù)據(jù)進行平滑降噪處理。
采集到的原始呼吸數(shù)據(jù)的數(shù)值分布并不統(tǒng)一且數(shù)值較大,因此還需對呼吸數(shù)據(jù)作歸一化處理,其目的是在統(tǒng)一量綱的同時便于梯度下降的計算、加快訓(xùn)練模型的收斂。實驗中采用Z-Score 方法對平滑降噪后的數(shù)據(jù)作標(biāo)準(zhǔn)歸一化處理,并在模型訓(xùn)練完成后將數(shù)據(jù)去標(biāo)準(zhǔn)化。
因每位志愿者的呼吸運動幅值、頻率等諸多因素不同,甚至同一人在不同時間所測得的呼吸數(shù)據(jù)亦大不相同,在呼吸預(yù)測中應(yīng)保證數(shù)據(jù)時效性,因此每組呼吸數(shù)據(jù)都需要構(gòu)建相應(yīng)的訓(xùn)練樣本。
將每條經(jīng)過標(biāo)準(zhǔn)歸一化后數(shù)據(jù)的前80%作為訓(xùn)練集,后20%作為驗證集。將預(yù)測時長設(shè)定為250ms,能夠補償大部分系統(tǒng)的延遲時間。自測呼吸數(shù)據(jù)的采樣頻率為20Hz,德國Lubeck 大學(xué)公開數(shù)據(jù)集中的呼吸數(shù)據(jù)采樣頻率也為20Hz 左右(17.5~21.3Hz),即預(yù)測步長為5。實驗設(shè)定通過前50 個數(shù)據(jù)點的位置變化預(yù)測后5 個數(shù)據(jù)點的位置信息。
利用滑動窗口將訓(xùn)練集X=[x1,x2,x3,…,xend]構(gòu)成有效的訓(xùn)練樣本,即輸入矩陣Input和輸出矩陣Output?;瑒哟翱跇?gòu)建訓(xùn)練樣本的過程如圖1 所示。輸入步長為50,預(yù)測步長為5,從訓(xùn)練集X中的x1開始將50 個數(shù)據(jù)點設(shè)定為輸入矩陣Input的第一列,從x51開始將5 個數(shù)據(jù)點設(shè)定為輸出矩陣Output的第一列,隨后滑動窗口右移繼續(xù)采樣。輸入矩陣Input和輸出矩陣Output如下:
Fig.1 Schematic diagram of sliding window圖1 滑動窗口示意圖
式(1)、式(2)中,n代表滑動窗口右移的步長,每一個訓(xùn)練集X可構(gòu)建k+1組訓(xùn)練樣本。
循環(huán)神經(jīng)網(wǎng)絡(luò)(Recurrent Neural Network,RNN)在時間序列數(shù)據(jù)預(yù)測方面有著較好的表現(xiàn)并且廣泛應(yīng)用于多個領(lǐng)域,但因梯度消失或爆炸問題,RNN 在處理長距離依賴問題上的能力有限。Hochreiter 等[18]提出的長短期記憶網(wǎng)絡(luò)(Long Short-Term Memory,LSTM)引入記憶單元和遺忘門,克服了RNN 的局限性。如圖2 所示,LSTM 有4 個門,即輸入門i、遺忘門f、控制門c和輸出門o[19]。
輸入門決定哪些信息可以傳遞給單元格,定義為:
遺忘門決定了輸入中的哪些信息應(yīng)該從先前的記憶中忽略,定義為:
控制門控制單元狀態(tài)從Ct-1到Ct的更新,如式(5)、式(6):
輸出門負責(zé)生成輸出和更新隱藏向量ht-1,定義為:
其中,σ是Sigmoid 激活函數(shù),W是相應(yīng)的權(quán)重矩陣,tanh用于將值縮放到-1~1的范圍內(nèi)。
Fig.2 LSTM unit structure圖2 LSTM 單元結(jié)構(gòu)
Bi-LSTM 屬于雙向循環(huán)神經(jīng)網(wǎng)絡(luò)(Bi-RNN)的一種,Bi-RNN 于1997 年由Schuster[20]首次提出。Bi-RNN 結(jié)構(gòu)如圖3 所示,由兩個正負時間方向的RNN 疊加而成,前向狀態(tài)和后向狀態(tài)連接在一個輸出層,在訓(xùn)練過程中Bi-RNN 在前向和后向兩個方向同時進行訓(xùn)練。Bi-RNN 的結(jié)構(gòu)決定了每一時刻的輸出都能獲取整個輸入序列的完整信息而不僅僅是當(dāng)前時間之前的信息。這種雙向結(jié)構(gòu)可以應(yīng)用于RNN 的變種網(wǎng)絡(luò),本實驗中就使用了Bi-LSTM。
Fig.3 The general structure of BRNN圖3 BRNN的一般結(jié)構(gòu)
本實驗基于MATLAB 平臺利用29 組呼吸數(shù)據(jù)對LSTM、門控循環(huán)單元(gated recurrent unit,GRU)、Bi-LSTM 3 種算法進行對照比較。實驗用計算機裝有Intel(R)Core(TM)i7-9750H CPU 和NVIDIA GeForce GTX 1650 GPU。
利用RMSE 和MAE 作為主要的實驗結(jié)果評價標(biāo)準(zhǔn),RMSE 表示預(yù)測值和真實值之間殘差的樣本標(biāo)準(zhǔn)差,MAE表示預(yù)測值和觀測值之間絕對誤差的平均值。
式(9)、式(10)中,N代表呼吸樣本數(shù)據(jù)點個數(shù),xi代表呼吸樣本數(shù)據(jù)點的實際值,ˉxi代表該點的算法預(yù)測值。此外,本實驗還統(tǒng)計了每一例樣本預(yù)測結(jié)果中誤差大于0.5mm 的數(shù)據(jù)點個數(shù)。
首先,對滑動窗口右移的步長n 取何值做了實驗比較,結(jié)果如圖4 所示。橫坐標(biāo)表示制作訓(xùn)練集時滑動窗口右移的步長,縱坐標(biāo)表示不同n 值時的預(yù)測精度。選用樣本1 進行實驗,預(yù)測算法采用LSTM。實驗結(jié)果顯示,其RMSE 與MAE 均呈現(xiàn)先減小后增大的趨勢,可觀察到當(dāng)滑動窗口右移步長為5 時,RMSE 與MAE 最小,預(yù)測效果最佳。
Fig.4 Experiment of sliding window right step size n圖4 滑動窗口右移步長n實驗測試
Fig.5 Experiment of network layer’s number圖5 網(wǎng)絡(luò)層數(shù)實驗測試
其次,對預(yù)測模型中的網(wǎng)絡(luò)層數(shù)設(shè)置進行實驗,實驗結(jié)果如圖5 所示。實驗選用樣本5 進行,預(yù)測算法采用Bi-LSTM。橫坐標(biāo)表示從1 至5 不同的Bi-LSTM 網(wǎng)絡(luò)層數(shù),縱坐標(biāo)表示不同網(wǎng)絡(luò)層數(shù)的預(yù)測精度??梢悦黠@看出,隨著網(wǎng)絡(luò)層數(shù)的增加,雖然網(wǎng)絡(luò)結(jié)構(gòu)逐漸變得復(fù)雜,但RMSE與MAE 卻呈現(xiàn)上升趨勢,因此當(dāng)預(yù)測算法的網(wǎng)絡(luò)層數(shù)為1時,預(yù)測效果最佳。
構(gòu)建訓(xùn)練集時,n 越小所得到的訓(xùn)練樣本更多,但卻導(dǎo)致樣本特征不明顯,在網(wǎng)絡(luò)訓(xùn)練中很容易導(dǎo)致網(wǎng)絡(luò)過擬合。因此,設(shè)置合適的n 值不但能減少訓(xùn)練時長,還可提高預(yù)測準(zhǔn)確性。設(shè)置網(wǎng)絡(luò)層數(shù)時發(fā)現(xiàn)更深的網(wǎng)絡(luò)層使預(yù)測精度下降且徒增訓(xùn)練時長,得出復(fù)雜的網(wǎng)絡(luò)結(jié)構(gòu)并非適用所有情況,本實驗中隱藏層設(shè)為單層預(yù)測效果最佳。此外,本實驗中網(wǎng)絡(luò)模型的超參數(shù)如下:訓(xùn)練次數(shù)、神經(jīng)元個數(shù)、批處理大小、初始學(xué)習(xí)率、學(xué)習(xí)率下降因子、下降迭代次數(shù)、求解器分別為350、200、all group、0.002、0.1、250、adam。
表1 對29 組呼吸數(shù)據(jù)預(yù)測結(jié)果的RMSE 進行統(tǒng)計并在結(jié)尾作平均值計算。相較于LSTM、GRU,Bi-LSTM 算法的平均值更小。29 組數(shù)據(jù)中有19 組,Bi-LSTM 算法預(yù)測結(jié)果的RMSE 小于LSTM、GRU。
Table 1 RMSE statistics of prediction results(mm)表1 預(yù)測結(jié)果的RMSE統(tǒng)計(mm)
表2 對預(yù)測結(jié)果的MAE 進行統(tǒng)計并在結(jié)尾作平均值計算。29 組數(shù)據(jù)中,Bi-LSTM 算法預(yù)測結(jié)果的MAE 均小于LSTM、GRU。并且,Bi-LSTM 算法預(yù)測結(jié)果MAE 的平均值明顯小于LSTM、GRU。
Table 2 MAE statistics of prediction results(mm)表2 預(yù)測結(jié)果的MAE統(tǒng)計(mm)
此外,針對臨床關(guān)注的預(yù)測誤差略大情況出現(xiàn)的次數(shù)也做了統(tǒng)計,結(jié)果如表3 所示。表3 中統(tǒng)計了預(yù)測結(jié)果中誤差大于0.5mm 的個數(shù),29 組數(shù)據(jù)中的18 組,Bi-LSTM 算法預(yù)測結(jié)果中誤差大于0.5mm 的個數(shù)小于LSTM、GRU,且Bi-LSTM 算法在該統(tǒng)計中平均值亦最小。
圖4 中展示的呼吸數(shù)據(jù)其呼吸頻率與幅值均不穩(wěn)定,且具有較為明顯的基線漂移;圖5 中展示的呼吸數(shù)據(jù)其呼吸狀態(tài)較為穩(wěn)定,頻率與幅值均無較大變化。此外,圖6-圖7 分別直觀展示了兩組呼吸數(shù)據(jù)應(yīng)用不同算法的預(yù)測結(jié)果。橫坐標(biāo)表示呼吸數(shù)據(jù)的點數(shù),縱坐標(biāo)表示呼吸位置數(shù)據(jù),藍色實線代表真實值,橙色實線代表預(yù)測值(彩圖掃OSID 碼可見)??梢灾庇^看出,LSTM 與GRU 算法的預(yù)測結(jié)果中,在呼吸的波峰和波谷處均有小幅度的搖擺,Bi-LSTM 的預(yù)測結(jié)果更擬合于真實值且更為平滑,明顯優(yōu)于LSTM 與GRU。
對實驗結(jié)果進行分析,從RMSE、MAE 和誤差大于0.5mm 的數(shù)據(jù)點個數(shù)這3 個參考標(biāo)準(zhǔn)可知,各項指標(biāo)Bi-LSTM 最佳。并且,通過直觀觀察發(fā)現(xiàn),在預(yù)測較不穩(wěn)定的呼吸數(shù)據(jù)時Bi-LSTM 算法依然結(jié)果較好。與很多相關(guān)研究相比,本實驗并沒有特意截取呼吸數(shù)據(jù)中較平穩(wěn)的區(qū)域進行實驗,而是將諸多不利因素(如呼吸頻率、幅值明顯變化,基線漂移等)都考慮在內(nèi),在比較準(zhǔn)確性的同時也驗證了算法魯棒性,使結(jié)果更貼近臨床需求。
Table 3 Statistics of number of data points with error>0.5mm表3 誤差>0.5mm的數(shù)據(jù)點個數(shù)統(tǒng)計
Fig.6 Examples of unstable breathing data prediction results圖6 不穩(wěn)定呼吸數(shù)據(jù)預(yù)測結(jié)果示例
Fig.7 Examples of relatively stable respiratory data prediction results圖7 較穩(wěn)定呼吸數(shù)據(jù)預(yù)測結(jié)果示例
本文提出一種基于Bi-LSTM 的呼吸預(yù)測模型框架,調(diào)整優(yōu)化了模型各參數(shù)后用29 組數(shù)據(jù)進行實驗。實驗結(jié)果表明,其在250ms 的延遲下平均RMSE、MAE 分別為0.091mm、0.042mm,在提高預(yù)測準(zhǔn)確性的同時也具有較強的魯棒性,能夠應(yīng)對復(fù)雜的呼吸狀況。并且,其單層隱藏層的網(wǎng)絡(luò)結(jié)構(gòu)使得訓(xùn)練速度很快,具有一定的臨床應(yīng)用可行性。后續(xù)將從如下方面進行深入研究:①建立內(nèi)外呼吸關(guān)聯(lián)模型;②搭建仿真穿刺實驗平臺;③進行仿真穿刺實驗。