胸部數(shù)字合成體層成像時(shí)的呼吸信號(hào)提取

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

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

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

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

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

1 方法

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

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

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

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

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

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

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

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

1.2 運(yùn)動(dòng)模型

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

圖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 呼吸信號(hào)提取

基于以上分析，我們?cè)趫D3中給出了提取呼吸信號(hào)工作流程圖，該流程圖包括以下步驟：

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

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

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

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

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

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

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

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

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

步驟五：獲取呼吸信號(hào)（圖3G）

圖3 應(yīng)用一次模型提取呼吸信號(hào)流程圖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.

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

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

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

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

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

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

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

1.4.3 病人數(shù)據(jù) 有關(guān)病人的研究得到了南方醫(yī)科大學(xué)南方醫(yī)院倫理委員會(huì)的批準(zhǔn)，在病人知情的情況下收集并使用數(shù)據(jù)。病人數(shù)據(jù)是通過Sonialvision Safire II射線照相系統(tǒng)（Shimadzu,Japan）獲取的，該系統(tǒng)配備了一個(gè)直接轉(zhuǎn)換的數(shù)字平板探測(cè)器。震源到等中心點(diǎn)的距離為998 mm，震源到檢測(cè)器的距離為1098毫米。數(shù)據(jù)采集時(shí)設(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個(gè)病人的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)動(dòng)軌跡是線性的。

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

圖8 數(shù)字體模數(shù)據(jù)的膜片位置圖和提取的膜片運(yùn)動(dòng)軌跡（紅色曲線）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 本文方法提取的呼吸信號(hào)曲線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例病患利用本文所提出的一次模型提取呼吸信號(hào)的結(jié)果如圖11～12所示?？梢钥闯?，三位患者都存在不同程度的呼吸?；颊遖和患者b都是有規(guī)律的呼吸，而患者c掃描前期屏息良好，后期突然劇烈呼吸。

圖11 3位患者數(shù)據(jù)利用一次模型擬合基本運(yùn)動(dòng)分量結(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)動(dòng)以及患者自身的呼吸運(yùn)動(dòng)后，提出了一種新的呼吸運(yùn)動(dòng)分析方法。在實(shí)驗(yàn)中，用一次多項(xiàng)式模擬了由CTS系統(tǒng)掃描幾何形狀引起的線性基礎(chǔ)運(yùn)動(dòng)。由于在實(shí)際應(yīng)用過程中病人呼吸狀態(tài)的不確定性，我們并未用某一種特定的模型去擬合呼吸運(yùn)動(dòng)。呼吸運(yùn)動(dòng)信息是由實(shí)際的橫膈膜的整體運(yùn)動(dòng)信號(hào)減去擬合的橫膈膜基礎(chǔ)運(yùn)動(dòng)得到。

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

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

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

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

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