李 仟，劉 佳，莫 英，于佳弘，張昆鵬，張 華

南方醫(yī)科大學(xué)生物醫(yī)學(xué)工程學(xué)院，廣東 廣州 510000

肺癌是最常見的致命癌癥之一［1］。胸部數(shù)字X射線平片攝影（DR）是肺部病變檢測［2］最常用的影像學(xué)方式。然而，DR產(chǎn)生的二維圖像包含重疊的解剖結(jié)構(gòu)，大大降低了DR對于肺部結(jié)節(jié)尤其是肺部微小結(jié)節(jié)或隔膜下結(jié)節(jié)檢測［3-4］的特異性。胸部數(shù)字合成體層成像（CTS）是一種新興的能夠通過去除部分重疊組織來提高解剖的可視性的成像技術(shù)，該技術(shù)利用有限角度下的掃描投影重建偽3D圖像［5-7］，因此，CTS檢測肺結(jié)節(jié)的靈敏度至少是DR的3倍［8-9］。肺癌胸部合成的檢出率與低劑量CT相當(dāng)［10］。與標(biāo)準(zhǔn)胸部CT相比，胸部數(shù)字合成體層成像分辨率更高，對患者的輻射劑量更低?？勺鳛橐环N增強(qiáng)肺實質(zhì)檢測和描述肺結(jié)節(jié)等異常的影像學(xué)工具［11］。

呼吸運(yùn)動引起的偽影在CTS檢查中較為常見。在某些情況下，由于呼吸運(yùn)動的存在，導(dǎo)致CTS對肺部病灶的特異性降至DR水平［9,12-13］。在臨床實踐中，放射技師很難直接從X線投影圖像上判斷病人是否屏住呼吸。如果能夠在掃描過程中自動識別是否屏氣良好（存在呼吸運(yùn)動），將會有非常重要的臨床意義。在CTS成像過程中，X射線管掃描在平行于探測器的平行直線路徑上［14］，投影域中的運(yùn)動物體的運(yùn)動模式與CT的不同，因此通常用于從胸錐束CT 投影中提取呼吸信號的方法［15］，例如Amsterdam Shroud（AS）［16］方法，可能不適用于CTS。

目前，對CTS中呼吸信號提取算法的研究較少，sin-quadratic模型［17］對于某些特殊狀態(tài)，如掃描過程中病人存在不規(guī)律的呼吸，可能不適合提取呼吸信號。在CTS成像中，管沿患者的頭尾方向直線移動，且每個視圖的X射線投影都可以清楚地顯示橫膈膜的邊緣線。當(dāng)前X射線球管的速度有限，故對呼吸信號提取影響不大。在此基礎(chǔ)上，本文提出了一個模型，以橫膈肌的物理位置作為替代信號來分析胸部數(shù)字合成體層成像時的呼吸運(yùn)動。該模型不僅包含了患者的呼吸運(yùn)動，而且考慮了CTS系統(tǒng)幾何結(jié)構(gòu)引起的基礎(chǔ)運(yùn)動。本文依據(jù)投影數(shù)據(jù)本身，通過對目標(biāo)對象呼吸運(yùn)動的實時分析，不需要任何其他外部呼吸檢測裝置作為參考，直接從投影圖中提取呼吸信號。

1 方法

1.1 CTS的運(yùn)動分析

1.1.1 靜止物體CTS運(yùn)動分析 在CTS掃描過程中，X射線球管與探測器是連續(xù)平行相對運(yùn)動的［18-20］，其運(yùn)動連續(xù)且均勻。當(dāng)X射線管勻速移動時，單個點的運(yùn)動軌跡將被很好地限制在一定范圍內(nèi)。如圖1所示，D是X射線管平面到探測器所在平面的距離，d是探測器與X射線管同步移動支點平面的高度。假設(shè)在某個時刻，X射線管位于坐標(biāo)c處，探測器的中心位于坐標(biāo)b處，那么位于（x，y）的物體的投影脈沖相對于探測器中點的相對位移m則有：m=b-a

圖1 CTS的幾何形狀Fig.1 Geometry of chest tomosynthesis.

由CTS系統(tǒng)幾何可以得到：

由上文可知，一旦CTS的幾何形狀固定，m就與c線性相關(guān)。對于使用固定探測器的系統(tǒng)（例如GE Volume RAD），m與c也存在線性關(guān)系。如果X射線管勻速移動，那么相鄰?fù)队皥D中單個點的移動距離幾乎是恒定的［21］，這是CTS模式所固有的。

然而，橫膈膜的運(yùn)動比單個點的運(yùn)動要復(fù)雜的多。從解剖學(xué)的角度來看，人體橫膈膜肌表現(xiàn)為丘陵狀。因此，即使是屏住呼吸的病人，每次曝光時X射線和橫膈膜的交點是不同的。式（1）所示的線性關(guān)系不足以描述橫膈膜的整體的運(yùn)動軌跡。圖2顯示了有無呼吸運(yùn)動的橫膈膜的CTS掃描示意圖。

1.1.2 CTS運(yùn)動物體的運(yùn)動分析 除了CTS系統(tǒng)成像幾何引起的基礎(chǔ)運(yùn)動之外，在CTS成像過程中無法屏住呼吸的患者，其橫膈膜投影位置也會受到呼吸運(yùn)動的影響。所以投影X線片中橫膈肌的運(yùn)動可以看作是成像系統(tǒng)引起的基礎(chǔ)運(yùn)動和患者呼吸運(yùn)動的結(jié)合。同時，本文中不再關(guān)注橫膈肌上單點的運(yùn)動狀態(tài)，而是將橫膈肌看作整體來構(gòu)建運(yùn)動模型。具體來說，我們可以將橫膈膜的運(yùn)動建模為：

其中M表示在投影中橫膈膜的總運(yùn)動，M呼吸表示患者呼吸引起的運(yùn)動，M基礎(chǔ)表示基礎(chǔ)運(yùn)動。

1.2 運(yùn)動模型

如圖2B所示，X線投影圖中橫膈膜的運(yùn)動可以概括為由系統(tǒng)幾何引起基礎(chǔ)運(yùn)動加上病人呼吸引起的呼吸運(yùn)動?？紤]到如圖2所示的，探測器、放射源異界掃描位置的關(guān)系，可以推定其線性相關(guān)。根據(jù)本文2.1分析，采用一次多項式估計橫膈膜的由CTS系統(tǒng)幾何引起的基礎(chǔ)運(yùn)動分量，并從整體運(yùn)動中減去基本運(yùn)動，呼吸信號就可以被檢測到。

圖2 CTS掃描橫膈膜分別在（A）靜息狀態(tài)（B）呼吸狀態(tài)Fig.2 CTS scan of the diaphragms in stationary state(A)and when breathing(B).

1.3 呼吸信號提取

基于以上分析，我們在圖3中給出了提取呼吸信號工作流程圖，該流程圖包括以下步驟：

步驟一：從投影中提取橫膈膜信息（圖3A～C）

通過提取膈膜的運(yùn)動來產(chǎn)生呼吸信號。按照投影順序，依次對每個投影圖像（圖像尺寸為1440*1440）進(jìn)行差分。本文中的差分圖像由每個投影圖像自身兩兩相對應(yīng)像素值相減（同一列中的第n行的像素值減去第n-2行的像素值）得到，以突出圖像中的變化部分。這種方法差分后的圖像往往能夠檢測出運(yùn)動目標(biāo)的輪廓，得到豎直方向上也就是橫膈膜運(yùn)動方向上的像素值的變化程度。一方面由于橫膈膜位于胸腔中，在投影圖像的中下部分；另一方面由于左側(cè)橫膈膜的運(yùn)動受心臟跳動的影響；所以在投影圖中選擇右邊的橫膈膜運(yùn)動來提取信號。在實現(xiàn)過程中，我們根據(jù)多幅圖像分析獲得構(gòu)造橫膈膜位置圖所需處理的圖像區(qū)域。再通過將該區(qū)域水平方向求和，獲取該幅圖像橫膈膜位置信息。

步驟二：將逐幀橫膈膜位置信息組合獲得橫膈膜位置圖（圖3D）

橫膈膜位置圖是將信息列按照其投影順序組合在一起，得到一個的大小為nb*np的矩陣。其中nb是每列的像素值的數(shù)目，np是投影圖的數(shù)目。

步驟三：從橫膈膜位置圖中提取橫膈膜運(yùn)動軌跡（圖3E）

本文中，我們利用基于動態(tài)規(guī)劃算法從橫膈膜位置圖中提取橫膈膜的運(yùn)動軌跡。限制條件是輪廓線精確地通過每列僅一次?；贒P理論［22-24］，可以建立成本函數(shù)：

其中（x1，x2，...xn）是每列中路徑的大小，n是列數(shù)（也是本研究中投影的數(shù)量）。i數(shù)據(jù)項c（xi）是通過路徑xi的成本，先前項d（xi-1,xi）是xi-1到xi之間部分路徑的成本。這里，用橫膈膜位置圖的梯度大小去構(gòu)造數(shù)據(jù)項：。先前項是xi-1與xi之間的絕對距離同時滿足d（xi-1，xi）≤k。所以輪廓線也就是橫膈膜的運(yùn)動軌跡是通過求解等式（3）得到的。

步驟四：擬合橫膈膜運(yùn)動（圖3F）

對步驟三中獲取的膈肌運(yùn)動軌跡，我們采用一次函數(shù)y=ax+b來擬合橫膈膜的基本運(yùn)動分量，擬合方法為最小二乘法。最小二乘法通過等式（4）最小化誤差的平方和尋找數(shù)據(jù)的最佳函數(shù)匹配。利用最小二乘法可以簡便地求得擬合函數(shù)的參數(shù)，并使得這些求得的數(shù)據(jù)與實際數(shù)據(jù)之間誤差的平方和為最小。

步驟五：獲取呼吸信號（圖3G）

圖3 應(yīng)用一次模型提取呼吸信號流程圖Fig.3 Flowchart of applying model for respiratory signal extraction.(A)are the CTS projections.The projections of the diaphragm in the right pulmonary field of the red rectangle in(B)selected to construct the diaphragma location map.The orange circle marked curve in(F)is the diaphragm moving trajectory extracted from the diaphragm location map(E).The solid blue curve in(F)and the black circular marked curve in(G)are the model fitted diaphragm trajectory caused by the system respectively.The circle marked red curve in(G)is the estimated respiratory signal.

呼吸信號可以通過等式（5），利用橫膈膜的整體運(yùn)動減去擬合的一次函數(shù)就可以得到。

其中M表示在投影中橫膈膜的總運(yùn)動，M呼吸表示患者呼吸引起的運(yùn)動，M擬合表示擬合的基礎(chǔ)運(yùn)動。

1.4 數(shù)據(jù)采集

1.4.1 物理體模數(shù)據(jù) 物理體模實驗在LUNGMAN上進(jìn)行。該體模是一個精確的真人大小的人體軀干解剖模型，胸壁的厚度是基于臨床數(shù)據(jù)的測量。軟組織替代材料和合成骨的X射線吸收率與人體組織非常接近。體模的骨骼和血管在圖像上顯示出與真人一樣的對比度和電子管電壓由縱隔、肺血管和腹部塊組成的內(nèi)部成分很容易分離，允許插入模擬腫瘤或其他病變。物體體模數(shù)據(jù)采集過程與下文中患者數(shù)據(jù)采集過程相同。圖4給出了物理體模的圖片。

圖4 物理體模Fig.4 Multipurpose chest phantom.

1.4.2 XCAT 體模數(shù)據(jù) 采用動態(tài)數(shù)字胸部XCAT 體模［25］進(jìn)行實驗仿真。XCAT體?？梢蕴峁┱鎸嵉娜梭w解剖模型。在人體解剖模四維胸部體模中設(shè)計10個呼吸階段。按照設(shè)定的次序進(jìn)行排列10個呼吸階段。設(shè)定數(shù)字體模的體素大小為1 mm×1 mm×1 mm。仿真中使用了在下節(jié)中描述的Shimadzu系統(tǒng)的幾何結(jié)構(gòu)。我們將大約一個半呼吸周期的4D體模投影到模擬CTS系統(tǒng)上，得到74張尺寸為400×400的投影。采用快速射線追蹤算法進(jìn)行仿真。不考慮X射線散射、光束質(zhì)量和探測器特性等因素的影響，因為這些因素不會影響所提出模型的適用性。圖5為模擬投影X線照片示例。

圖5 XCAT數(shù)據(jù)仿真CTS投影圖像Fig.5 Digital XCAT phantom-based simulation.

1.4.3 病人數(shù)據(jù) 有關(guān)病人的研究得到了南方醫(yī)科大學(xué)南方醫(yī)院倫理委員會的批準(zhǔn)，在病人知情的情況下收集并使用數(shù)據(jù)。病人數(shù)據(jù)是通過Sonialvision Safire II射線照相系統(tǒng)（Shimadzu,Japan）獲取的，該系統(tǒng)配備了一個直接轉(zhuǎn)換的數(shù)字平板探測器。震源到等中心點的距離為998 mm，震源到檢測器的距離為1098毫米。數(shù)據(jù)采集時設(shè)置管電壓為110 kVp，管電流［26］為160 mA，采集管掃角為±20°、尺寸為1440×1440的X線投影圖片74張。采集3組患者數(shù)據(jù)，分別標(biāo)記為患者a,b,c。其中患者b為明顯的脊柱側(cè)彎病人。圖6顯示了3組患者數(shù)據(jù)的X線投影圖。

圖6 3個病人的X線投影圖Fig.6 Projection radiographs of 3 patients.A:Patient a.B:Patient b.C:Patient c.

2 結(jié)果

2.1 物理體模數(shù)據(jù)的結(jié)果

從圖7中可以看出，橫隔膜的運(yùn)動軌跡是線性的。

圖7 物理體模數(shù)據(jù)結(jié)果Fig.7 Results of multipurpose chest phantom data.A:Projections of multipurpose chest phantom.B:Extracted diaphragm trajectory in stationary state.C:The trajectory of the diaphragm in stationary state is linear.

2.2 XCAT數(shù)字體模數(shù)據(jù)的結(jié)果

圖8顯示了XCAT數(shù)字體模橫隔膜位置圖和提取的橫隔膜運(yùn)動軌跡。從結(jié)果可以看出，橫膈膜的邊緣可以準(zhǔn)確的代表運(yùn)動?；趧討B(tài)規(guī)劃的技術(shù)可以從橫膈膜位置圖很好的提取橫膈膜運(yùn)動軌跡信息。圖9C紅色帶圈標(biāo)記曲線為本文模型提取出的呼吸運(yùn)動軌跡。圖10為設(shè)計的呼吸信號（圓形標(biāo)記紅色實線）與本文方法提取的呼吸信號（圓形標(biāo)記黑色實線）以及陶等人方法提取的呼吸信號（‘+’標(biāo)記綠色實線）三者之間的視覺對比。通過計算得到陶等人提出的方法提取的呼吸信號曲線與設(shè)計的呼吸信號曲線之間的相關(guān)系數(shù)為0.9683，而本文方法提取的呼吸信號曲線與設(shè)計的呼吸信號曲線之間的相關(guān)系數(shù)為0.9797。

圖8 數(shù)字體模數(shù)據(jù)的膜片位置圖和提取的膜片運(yùn)動軌跡（紅色曲線）Fig.8 Diaphragm location map of the phantom data and the extracted diaphragm motion trajectory(red curve).

圖9 數(shù)字體模數(shù)據(jù)結(jié)果Fig.9 Results of digital phantom data.A:Diaphragm motion trajectory.B:Model fitted trajectory (blue curve).C:Diaphragm trajectory(circle-marked black curve)and respiratory motion trajectory(circle-marked red curve).

圖10 本文方法提取的呼吸信號曲線Fig.10 Comparison of the respiratory signal curve extracted by the method in this study (black curve marked by diamond),the respiratory signal curve extracted by Tao et al.'s method (green curve marked by '+') and the respiratory signal curve designed by volume model(red curve marked by circle).

2.3 患者數(shù)據(jù)結(jié)果

3例病患利用本文所提出的一次模型提取呼吸信號的結(jié)果如圖11～12所示?？梢钥闯?，三位患者都存在不同程度的呼吸。患者a和患者b都是有規(guī)律的呼吸，而患者c掃描前期屏息良好，后期突然劇烈呼吸。

圖11 3位患者數(shù)據(jù)利用一次模型擬合基本運(yùn)動分量結(jié)果Fig.11 Model fitting of basic motion components based on data from 3 patients.A:Patient a.B:Patient b.C:Patient c.

3 討論

本文在充分考慮到CTS系統(tǒng)的掃描幾何形狀引起的基礎(chǔ)運(yùn)動以及患者自身的呼吸運(yùn)動后，提出了一種新的呼吸運(yùn)動分析方法。在實驗中，用一次多項式模擬了由CTS系統(tǒng)掃描幾何形狀引起的線性基礎(chǔ)運(yùn)動。由于在實際應(yīng)用過程中病人呼吸狀態(tài)的不確定性，我們并未用某一種特定的模型去擬合呼吸運(yùn)動。呼吸運(yùn)動信息是由實際的橫膈膜的整體運(yùn)動信號減去擬合的橫膈膜基礎(chǔ)運(yùn)動得到。

具體地說，呼吸信號是通過從每個二維投影中提取橫膈膜運(yùn)動信息估計得到的。這是通過將投影圖進(jìn)行差分后手動選擇大致位于右隔膜中央?yún)^(qū)域的像素水平求和得到信息列，作為替代呼吸信號來完成的。因為在這個位置的運(yùn)動比其他位置更明顯。提取的信息列按照投影順序一起排列，以產(chǎn)生一個隔膜位置圖。采用基于動態(tài)規(guī)劃邊緣提取方法提取光闌運(yùn)動軌跡。然后利用一次模型擬合運(yùn)動軌跡，估計橫膈膜由系統(tǒng)幾何引起的基礎(chǔ)運(yùn)動分量，并從整體運(yùn)動中減去基礎(chǔ)運(yùn)動，可以很容易地估計呼吸信號。整個過程的計算時間遠(yuǎn)小于1 s，因此本文方法可以看作是一種實時的方法。

本研究中，物理體模數(shù)據(jù)提取結(jié)果可以看出，利用物理體模無呼吸運(yùn)動的情況下提取橫膈膜運(yùn)動軌跡來測定用一次線性模型標(biāo)定基礎(chǔ)運(yùn)動的可行性，本文模型能夠很好地估計橫膈膜總體的由系統(tǒng)幾何引起的線性基礎(chǔ)運(yùn)動分量。定量分析表明，XCAT仿真實驗中本文提出的方法和2018年Tao等［17］提出的sin-quadratic模型都能有效地提取設(shè)計的呼吸信號，本文模型所提取的呼吸信號與原設(shè)計的呼吸信號之間的相關(guān)系數(shù)為0.9797，陶等人提出的模型所提取的呼吸信號與原設(shè)計的呼吸信號之間的相關(guān)系數(shù)為0.9683，與Tao等［17］的方法相比本文所提出的模型精確度更高。定性和定量結(jié)果均可表明，該模型能有效地從CTS投影中提取呼吸信號。3組患者數(shù)據(jù)結(jié)果也表明，在不同的呼吸狀態(tài)下，本文模型都可以有效地對呼吸信號進(jìn)行估計，表現(xiàn)得較為穩(wěn)健。在病人有規(guī)律的呼吸模式下，本文提出的方法和2018年Tao等［17］提出的sin-quadratic模型都能有效地估計患者的呼吸信號。但有一點值得注意：對于不規(guī)則呼吸模式的患者，相比Tao等［17］利用正弦函數(shù)去擬合呼吸運(yùn)動，本文并未用某一種特定的模型去擬合呼吸運(yùn)動。

圖12 3位患者數(shù)據(jù)提取呼吸信號結(jié)果Fig.12 Results of respiratory signal extraction from the datasets of 3 patients.A:Patient a.B:Patient b.C:Patient c.

因為CTS成像的掃描時間為6～12 s，也就是大約有1～1.5個呼吸周期。在真實的臨床病例中，部分患者在整個數(shù)據(jù)采集過程中可能無法屏住呼吸如同患者數(shù)據(jù)中的患者c，在采集過程中突然劇烈呼吸，很明顯這種不規(guī)則的呼吸狀態(tài)并不能用正弦函數(shù)來表達(dá)［27-28］。在這種情況下，Tao等［17］提出的sin-quadratic模型表現(xiàn)的不夠穩(wěn)健不適合提取呼吸信號。然而本文并未用某一種特定的模型去擬合呼吸運(yùn)動，不規(guī)則呼吸信號仍然可以檢測到，本文模型表現(xiàn)得更為穩(wěn)健，適合大多數(shù)真實的臨床病例掃描狀態(tài)下。

總的來說，我們開發(fā)了一種實時的呼吸信號運(yùn)動分析方法。放射技師很難直接從X線投影圖像上判斷病人是否屏住呼吸，本文的方法可以在掃描過程中自動識別是否屏氣良好（存在呼吸運(yùn)動），直接從CTS投影中提取橫膈肌運(yùn)動軌跡，而不需要任何額外的呼吸檢測設(shè)備或其他參考圖像，具有重要的臨床意義。據(jù)所獲得的呼吸信號，醫(yī)生可立即對患者進(jìn)行重新掃描或指導(dǎo)無運(yùn)動CTS圖像重建。一旦檢測到強(qiáng)烈紊亂的呼吸信號，醫(yī)生可立即對存在呼吸運(yùn)動的病人進(jìn)行重新掃描。在呼吸運(yùn)動存在的情況下模糊的肋骨和脊柱在無呼吸運(yùn)動4D-CTS重建［29-30］中可以清晰地呈現(xiàn)出來，提高醫(yī)生的診斷精確度。