王凱龍，常兆華，葉 萍

（上海理工大學(xué)健康科學(xué)與工程學(xué)院，上海 200082）

0 引言

目前在全球范圍內(nèi)，肺癌仍是癌癥中死亡率最高的疾病之一［1］，近年來以滅活腫瘤細(xì)胞為目的的冷凍消融治療手段引起了人們的廣泛關(guān)注。該治療手段需在CT 引導(dǎo)下做氬氦刀經(jīng)皮穿刺，精準(zhǔn)穿刺是治療的重要條件，但胸腹部腫瘤的位置會(huì)隨著患者的呼吸發(fā)生運(yùn)動(dòng)，這種運(yùn)動(dòng)會(huì)影響穿刺的精準(zhǔn)度，在有些情況下腫瘤運(yùn)動(dòng)會(huì)達(dá)到28mm［2］。研究表明，胸腹部外部呼吸信號(hào)與肺部腫瘤運(yùn)動(dòng)具有良好的相關(guān)性，因此利用外部信號(hào)跟蹤肺部腫瘤位置是可行的［3］。但是由于外部運(yùn)動(dòng)信號(hào)的獲取、軟件內(nèi)部算法的運(yùn)行以及機(jī)械臂延遲而導(dǎo)致的系統(tǒng)延遲會(huì)影響呼吸跟蹤的實(shí)時(shí)性［4］，通常這種系統(tǒng)延遲在200ms 左右［5-7］，這會(huì)降低機(jī)器人的穿刺精度。Ma 等［8］研究表明，利用呼吸信號(hào)進(jìn)行預(yù)測的方法能夠有效提升穿刺機(jī)器人的穿刺精度，他們利用呼吸預(yù)測模型對(duì)公共數(shù)據(jù)集中的呼吸信號(hào)進(jìn)行預(yù)測，并用運(yùn)動(dòng)體模進(jìn)行穿刺實(shí)驗(yàn)，結(jié)果表明，預(yù)測結(jié)果的平均均方根誤差（Root Mean Squared Error，RMSE）為0.58mm，穿刺誤差的平均均方根（Root Mean Square，RMS）為0.65 mm。

有關(guān)呼吸跟蹤和呼吸預(yù)測的研究在放療中展開較早，Hoogeman 等［9］的研究表明，使用呼吸跟蹤和呼吸預(yù)測技術(shù)能夠有效地減少由呼吸運(yùn)動(dòng)和系統(tǒng)延遲引起的幾何誤差。比較經(jīng)典的呼吸預(yù)測算法有：CyberKnife?的Synchrony?呼吸同步跟蹤系統(tǒng)使用的預(yù)測算法、支持向量機(jī)（Support Vector Machine，SVM）、自回歸移動(dòng)平均模型（Autoregressive Integrated Moving Average Model，ARIMA）和神經(jīng)網(wǎng)絡(luò)。

Synchrony?系統(tǒng)在臨床實(shí)踐中表現(xiàn)良好，其利用兩個(gè)正交的kV 級(jí)X 射線系統(tǒng)等中心地投射到患者的治療部位觀測內(nèi)部腫瘤運(yùn)動(dòng)并建立與外部光學(xué)標(biāo)記相關(guān)的運(yùn)動(dòng)模型，然后使用呼吸預(yù)測技術(shù)補(bǔ)償固有的系統(tǒng)延遲［10-11］。Synchrony?呼吸同步跟蹤系統(tǒng)應(yīng)用了一種混合預(yù)測算法，該混合算法結(jié)合了線性算法、模式匹配和模糊邏輯3 種方法［12-13］。支持向量回歸（Support Vector Regression，SVR）是SVM 的一個(gè)應(yīng)用分支，Riaz 等［14］利用SVR 方法對(duì)腫瘤運(yùn)動(dòng)進(jìn)行預(yù)測，在延遲為1s 的狀態(tài)下預(yù)測結(jié)果的RMSE 在2mm以內(nèi)。Khashei 等［15］研究表明，ARIMA 在時(shí)間序列預(yù)測方面表現(xiàn)良好，Ichiji 等［16］使用ARIMA 結(jié)合多層感知器（Multilayer Perceptron，MLP）對(duì)腫瘤運(yùn)動(dòng)進(jìn)行預(yù)測，其結(jié)果顯示平均絕對(duì)誤差（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ù)測呼吸運(yùn)動(dòng)信號(hào)，設(shè)定延遲為400ms 的情況下其預(yù)測結(jié)果的平均RMSE為0.097 mm、MAE為0.074 mm。

很多相關(guān)研究實(shí)驗(yàn)所用的呼吸數(shù)據(jù)樣本數(shù)量過少，或僅截取呼吸平穩(wěn)的區(qū)段進(jìn)行實(shí)驗(yàn)，無法測試模型的魯棒性。本文使用了29 組呼吸數(shù)據(jù)進(jìn)行實(shí)驗(yàn)，在保證樣本數(shù)量的同時(shí)不特意截取呼吸平穩(wěn)區(qū)段，使實(shí)驗(yàn)結(jié)果更貼近臨床需求。本文提出的基于Bi-LSTM 的呼吸預(yù)測模型框架不但具有很高的準(zhǔn)確性，而且在面對(duì)各種復(fù)雜呼吸狀況時(shí)也表現(xiàn)出了較好的魯棒性。

1 實(shí)驗(yàn)數(shù)據(jù)

1.1 數(shù)據(jù)獲取及預(yù)處理

實(shí)驗(yàn)中，使用NDI Polaris 光學(xué)定位系統(tǒng)采集了23 組志愿者的外部呼吸數(shù)據(jù)，該系統(tǒng)通過追蹤仰臥狀態(tài)下的志愿者胸腹部體表的紅外光學(xué)定位支架以獲取呼吸數(shù)據(jù)，系統(tǒng)采樣頻率設(shè)定為20Hz，每組的采樣時(shí)間為5～7min，采集到的每組數(shù)據(jù)包括呼吸運(yùn)動(dòng)的左右、頭腳、前后3 個(gè)維度。志愿者在呼吸數(shù)據(jù)采集過程中完全處于自由呼吸狀態(tài)，因此數(shù)據(jù)中包含了基線漂移、呼吸頻率改變、深呼吸等狀態(tài)。此外，實(shí)驗(yàn)還加入了德國Lubeck 大學(xué)機(jī)器人與認(rèn)知系統(tǒng)研究所公開數(shù)據(jù)庫中的6 組呼吸數(shù)據(jù)，該數(shù)據(jù)的志愿者均在肺部有惡性腫瘤。

為了方便實(shí)驗(yàn)結(jié)果的比較并保證實(shí)驗(yàn)數(shù)據(jù)可靠，本文中統(tǒng)一選取前后方向的數(shù)據(jù)進(jìn)行預(yù)測，且并不特意截取呼吸數(shù)據(jù)的平滑段進(jìn)行試驗(yàn)，保留原數(shù)據(jù)的諸多不利于預(yù)測精度的狀況，29組數(shù)據(jù)均截取6 000個(gè)數(shù)據(jù)點(diǎn)進(jìn)行實(shí)驗(yàn)。

采集的原始呼吸數(shù)據(jù)的噪聲嚴(yán)重影響了之后數(shù)據(jù)訓(xùn)練效果，因此實(shí)驗(yàn)開始時(shí)需對(duì)每條數(shù)據(jù)進(jìn)行平滑降噪處理，在保證呼吸數(shù)據(jù)原有特性不變的前提下降低信噪比。本文選用Savitzky-Golay 濾波器對(duì)原始數(shù)據(jù)進(jìn)行平滑降噪處理。

采集到的原始呼吸數(shù)據(jù)的數(shù)值分布并不統(tǒng)一且數(shù)值較大，因此還需對(duì)呼吸數(shù)據(jù)作歸一化處理，其目的是在統(tǒng)一量綱的同時(shí)便于梯度下降的計(jì)算、加快訓(xùn)練模型的收斂。實(shí)驗(yàn)中采用Z-Score 方法對(duì)平滑降噪后的數(shù)據(jù)作標(biāo)準(zhǔn)歸一化處理，并在模型訓(xùn)練完成后將數(shù)據(jù)去標(biāo)準(zhǔn)化。

1.2 訓(xùn)練樣本構(gòu)建

因每位志愿者的呼吸運(yùn)動(dòng)幅值、頻率等諸多因素不同，甚至同一人在不同時(shí)間所測得的呼吸數(shù)據(jù)亦大不相同，在呼吸預(yù)測中應(yīng)保證數(shù)據(jù)時(shí)效性，因此每組呼吸數(shù)據(jù)都需要構(gòu)建相應(yīng)的訓(xùn)練樣本。

將每條經(jīng)過標(biāo)準(zhǔn)歸一化后數(shù)據(jù)的前80%作為訓(xùn)練集，后20%作為驗(yàn)證集。將預(yù)測時(shí)長設(shè)定為250ms，能夠補(bǔ)償大部分系統(tǒng)的延遲時(shí)間。自測呼吸數(shù)據(jù)的采樣頻率為20Hz，德國Lubeck 大學(xué)公開數(shù)據(jù)集中的呼吸數(shù)據(jù)采樣頻率也為20Hz 左右（17.5～21.3Hz），即預(yù)測步長為5。實(shí)驗(yàn)設(shè)定通過前50 個(gè)數(shù)據(jù)點(diǎn)的位置變化預(yù)測后5 個(gè)數(shù)據(jù)點(diǎn)的位置信息。

利用滑動(dòng)窗口將訓(xùn)練集X=[x1，x2，x3，…，xend]構(gòu)成有效的訓(xùn)練樣本，即輸入矩陣Input和輸出矩陣Output?；瑒?dòng)窗口構(gòu)建訓(xùn)練樣本的過程如圖1 所示。輸入步長為50，預(yù)測步長為5，從訓(xùn)練集X中的x1開始將50 個(gè)數(shù)據(jù)點(diǎn)設(shè)定為輸入矩陣Input的第一列，從x51開始將5 個(gè)數(shù)據(jù)點(diǎn)設(shè)定為輸出矩陣Output的第一列，隨后滑動(dòng)窗口右移繼續(xù)采樣。輸入矩陣Input和輸出矩陣Output如下：

Fig.1 Schematic diagram of sliding window圖1 滑動(dòng)窗口示意圖

式（1）、式（2）中，n代表滑動(dòng)窗口右移的步長，每一個(gè)訓(xùn)練集X可構(gòu)建k+1組訓(xùn)練樣本。

2 預(yù)測模型

2.1 LSTM 模型

循環(huán)神經(jīng)網(wǎng)絡(luò)（Recurrent Neural Network，RNN）在時(shí)間序列數(shù)據(jù)預(yù)測方面有著較好的表現(xiàn)并且廣泛應(yīng)用于多個(gè)領(lǐng)域，但因梯度消失或爆炸問題，RNN 在處理長距離依賴問題上的能力有限。Hochreiter 等［18］提出的長短期記憶網(wǎng)絡(luò)（Long Short-Term Memory，LSTM）引入記憶單元和遺忘門，克服了RNN 的局限性。如圖2 所示，LSTM 有4 個(gè)門，即輸入門i、遺忘門f、控制門c和輸出門o［19］。

輸入門決定哪些信息可以傳遞給單元格，定義為：

遺忘門決定了輸入中的哪些信息應(yīng)該從先前的記憶中忽略，定義為：

控制門控制單元狀態(tài)從Ct-1到Ct的更新，如式（5）、式（6）：

輸出門負(fù)責(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)

2.2 Bi-LSTM 模型

Bi-LSTM 屬于雙向循環(huán)神經(jīng)網(wǎng)絡(luò)（Bi-RNN）的一種，Bi-RNN 于1997 年由Schuster［20］首次提出。Bi-RNN 結(jié)構(gòu)如圖3 所示，由兩個(gè)正負(fù)時(shí)間方向的RNN 疊加而成，前向狀態(tài)和后向狀態(tài)連接在一個(gè)輸出層，在訓(xùn)練過程中Bi-RNN 在前向和后向兩個(gè)方向同時(shí)進(jìn)行訓(xùn)練。Bi-RNN 的結(jié)構(gòu)決定了每一時(shí)刻的輸出都能獲取整個(gè)輸入序列的完整信息而不僅僅是當(dāng)前時(shí)間之前的信息。這種雙向結(jié)構(gòu)可以應(yīng)用于RNN 的變種網(wǎng)絡(luò)，本實(shí)驗(yàn)中就使用了Bi-LSTM。

Fig.3 The general structure of BRNN圖3 BRNN的一般結(jié)構(gòu)

3 實(shí)驗(yàn)結(jié)果與分析

3.1 實(shí)驗(yàn)平臺(tái)及評(píng)價(jià)標(biāo)準(zhǔn)

本實(shí)驗(yàn)基于MATLAB 平臺(tái)利用29 組呼吸數(shù)據(jù)對(duì)LSTM、門控循環(huán)單元（gated recurrent unit，GRU）、Bi-LSTM 3 種算法進(jìn)行對(duì)照比較。實(shí)驗(yàn)用計(jì)算機(jī)裝有Intel（R）Core（TM）i7-9750H CPU 和NVIDIA GeForce GTX 1650 GPU。

利用RMSE 和MAE 作為主要的實(shí)驗(yàn)結(jié)果評(píng)價(jià)標(biāo)準(zhǔn)，RMSE 表示預(yù)測值和真實(shí)值之間殘差的樣本標(biāo)準(zhǔn)差，MAE表示預(yù)測值和觀測值之間絕對(duì)誤差的平均值。

式（9）、式（10）中，N代表呼吸樣本數(shù)據(jù)點(diǎn)個(gè)數(shù)，xi代表呼吸樣本數(shù)據(jù)點(diǎn)的實(shí)際值，ˉxi代表該點(diǎn)的算法預(yù)測值。此外，本實(shí)驗(yàn)還統(tǒng)計(jì)了每一例樣本預(yù)測結(jié)果中誤差大于0.5mm 的數(shù)據(jù)點(diǎn)個(gè)數(shù)。

3.2 模型參數(shù)實(shí)驗(yàn)

首先，對(duì)滑動(dòng)窗口右移的步長n 取何值做了實(shí)驗(yàn)比較，結(jié)果如圖4 所示。橫坐標(biāo)表示制作訓(xùn)練集時(shí)滑動(dòng)窗口右移的步長，縱坐標(biāo)表示不同n 值時(shí)的預(yù)測精度。選用樣本1 進(jìn)行實(shí)驗(yàn)，預(yù)測算法采用LSTM。實(shí)驗(yàn)結(jié)果顯示，其RMSE 與MAE 均呈現(xiàn)先減小后增大的趨勢，可觀察到當(dāng)滑動(dòng)窗口右移步長為5 時(shí)，RMSE 與MAE 最小，預(yù)測效果最佳。

Fig.4 Experiment of sliding window right step size n圖4 滑動(dòng)窗口右移步長n實(shí)驗(yàn)測試

Fig.5 Experiment of network layer’s number圖5 網(wǎng)絡(luò)層數(shù)實(shí)驗(yàn)測試

其次，對(duì)預(yù)測模型中的網(wǎng)絡(luò)層數(shù)設(shè)置進(jìn)行實(shí)驗(yàn)，實(shí)驗(yàn)結(jié)果如圖5 所示。實(shí)驗(yàn)選用樣本5 進(jìn)行，預(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時(shí)，預(yù)測效果最佳。

構(gòu)建訓(xùn)練集時(shí)，n 越小所得到的訓(xùn)練樣本更多，但卻導(dǎo)致樣本特征不明顯，在網(wǎng)絡(luò)訓(xùn)練中很容易導(dǎo)致網(wǎng)絡(luò)過擬合。因此，設(shè)置合適的n 值不但能減少訓(xùn)練時(shí)長，還可提高預(yù)測準(zhǔn)確性。設(shè)置網(wǎng)絡(luò)層數(shù)時(shí)發(fā)現(xiàn)更深的網(wǎng)絡(luò)層使預(yù)測精度下降且徒增訓(xùn)練時(shí)長，得出復(fù)雜的網(wǎng)絡(luò)結(jié)構(gòu)并非適用所有情況，本實(shí)驗(yàn)中隱藏層設(shè)為單層預(yù)測效果最佳。此外，本實(shí)驗(yàn)中網(wǎng)絡(luò)模型的超參數(shù)如下：訓(xùn)練次數(shù)、神經(jīng)元個(gè)數(shù)、批處理大小、初始學(xué)習(xí)率、學(xué)習(xí)率下降因子、下降迭代次數(shù)、求解器分別為350、200、all group、0.002、0.1、250、adam。

3.3 呼吸數(shù)據(jù)預(yù)測實(shí)驗(yàn)

表1 對(duì)29 組呼吸數(shù)據(jù)預(yù)測結(jié)果的RMSE 進(jìn)行統(tǒng)計(jì)并在結(jié)尾作平均值計(jì)算。相較于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)計(jì)（mm）

表2 對(duì)預(yù)測結(jié)果的MAE 進(jìn)行統(tǒng)計(jì)并在結(jié)尾作平均值計(jì)算。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)計(jì)（mm）

此外，針對(duì)臨床關(guān)注的預(yù)測誤差略大情況出現(xiàn)的次數(shù)也做了統(tǒng)計(jì)，結(jié)果如表3 所示。表3 中統(tǒng)計(jì)了預(yù)測結(jié)果中誤差大于0.5mm 的個(gè)數(shù)，29 組數(shù)據(jù)中的18 組，Bi-LSTM 算法預(yù)測結(jié)果中誤差大于0.5mm 的個(gè)數(shù)小于LSTM、GRU，且Bi-LSTM 算法在該統(tǒng)計(jì)中平均值亦最小。

圖4 中展示的呼吸數(shù)據(jù)其呼吸頻率與幅值均不穩(wěn)定，且具有較為明顯的基線漂移；圖5 中展示的呼吸數(shù)據(jù)其呼吸狀態(tài)較為穩(wěn)定，頻率與幅值均無較大變化。此外，圖6-圖7 分別直觀展示了兩組呼吸數(shù)據(jù)應(yīng)用不同算法的預(yù)測結(jié)果。橫坐標(biāo)表示呼吸數(shù)據(jù)的點(diǎn)數(shù)，縱坐標(biāo)表示呼吸位置數(shù)據(jù)，藍(lán)色實(shí)線代表真實(shí)值，橙色實(shí)線代表預(yù)測值（彩圖掃OSID 碼可見）。可以直觀看出，LSTM 與GRU 算法的預(yù)測結(jié)果中，在呼吸的波峰和波谷處均有小幅度的搖擺，Bi-LSTM 的預(yù)測結(jié)果更擬合于真實(shí)值且更為平滑，明顯優(yōu)于LSTM 與GRU。

對(duì)實(shí)驗(yàn)結(jié)果進(jìn)行分析，從RMSE、MAE 和誤差大于0.5mm 的數(shù)據(jù)點(diǎn)個(gè)數(shù)這3 個(gè)參考標(biāo)準(zhǔn)可知，各項(xiàng)指標(biāo)Bi-LSTM 最佳。并且，通過直觀觀察發(fā)現(xiàn)，在預(yù)測較不穩(wěn)定的呼吸數(shù)據(jù)時(shí)Bi-LSTM 算法依然結(jié)果較好。與很多相關(guān)研究相比，本實(shí)驗(yàn)并沒有特意截取呼吸數(shù)據(jù)中較平穩(wěn)的區(qū)域進(jìn)行實(shí)驗(yàn)，而是將諸多不利因素（如呼吸頻率、幅值明顯變化，基線漂移等）都考慮在內(nèi)，在比較準(zhǔn)確性的同時(shí)也驗(yàn)證了算法魯棒性，使結(jié)果更貼近臨床需求。

Table 3 Statistics of number of data points with error＞0.5mm表3 誤差＞0.5mm的數(shù)據(jù)點(diǎn)個(gè)數(shù)統(tǒng)計(jì)

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é)果示例

4 結(jié)語

本文提出一種基于Bi-LSTM 的呼吸預(yù)測模型框架，調(diào)整優(yōu)化了模型各參數(shù)后用29 組數(shù)據(jù)進(jìn)行實(shí)驗(yàn)。實(shí)驗(yàn)結(jié)果表明，其在250ms 的延遲下平均RMSE、MAE 分別為0.091mm、0.042mm，在提高預(yù)測準(zhǔn)確性的同時(shí)也具有較強(qiáng)的魯棒性，能夠應(yīng)對(duì)復(fù)雜的呼吸狀況。并且，其單層隱藏層的網(wǎng)絡(luò)結(jié)構(gòu)使得訓(xùn)練速度很快，具有一定的臨床應(yīng)用可行性。后續(xù)將從如下方面進(jìn)行深入研究：①建立內(nèi)外呼吸關(guān)聯(lián)模型；②搭建仿真穿刺實(shí)驗(yàn)平臺(tái)；③進(jìn)行仿真穿刺實(shí)驗(yàn)。