孫 正

（華北電力大學(xué)電子與通信工程系，河北 保定 071003）

血管內(nèi)超聲(intravascular ultrasound，IVUS)成像是臨床診斷冠心病的一種基于導(dǎo)管的腔內(nèi)成像技術(shù)，具有獨(dú)特的在活體中觀察血管壁、管腔及斑塊形態(tài)性質(zhì)，甚至管壁功能狀態(tài)的特點(diǎn)，已經(jīng)成為 X射線冠狀動(dòng)脈造影(X-ray coronary angiography)的彌補(bǔ)影像技術(shù)。

由于心臟運(yùn)動(dòng)以及血壓變化所造成的運(yùn)動(dòng)偽像是影響IVUS圖像的視覺效果，進(jìn)而影響定量分析和血管三維重建精度的主要因素。目前主要采用門控的方法減少 IVUS序列中的運(yùn)動(dòng)偽像，包括3種方式：

一是心電(electrocardiogram，ECG)門控的圖像采集方法：即采用 ECG門控的馬達(dá)回撤超聲導(dǎo)管，根據(jù) ECG信號只在每個(gè)心動(dòng)周期的相同時(shí)相（一般是R波）采集圖像，改善心臟搏動(dòng)所致的偽像[1]。該方法的缺點(diǎn)在于需要專門的ECG門控采集裝置，而目前臨床采用的多數(shù)IVUS系統(tǒng)不包含此功能。并且由于每個(gè)心動(dòng)周期只采集一幀，因此與連續(xù)回撤導(dǎo)管的非門控采集方式相比，圖像采集時(shí)間延長至少3倍[2]。

二是離線 ECG門控方法：即連續(xù)回撤超聲導(dǎo)管采集圖像序列，同時(shí)記錄 ECG信號。待介入檢查過程結(jié)束后，由醫(yī)生根據(jù) ECG信號從圖像序列中選取在相同心臟時(shí)相采集的圖像。這種方式雖然不會(huì)延長介入檢查的時(shí)間，但是其不足之處在于：ECG信號是心臟的生理電特性通過骨骼、肌肉和皮膚等反映在體表的信號，而 IVUS圖像是通過插入血管腔內(nèi)的超聲探頭采集的，兩者不可能完全同步；門控結(jié)果的客觀性和可重復(fù)性差，受操作者的臨床經(jīng)驗(yàn)和主觀因素影響較大；雖然多采用 0%點(diǎn)即舒張末期作為采樣點(diǎn)，但無法確定采用該點(diǎn)是否可以獲得最穩(wěn)定的門控圖像序列；目前尚無公認(rèn)、明確的標(biāo)準(zhǔn)，來確定采用心動(dòng)周期中的哪一時(shí)相作為最佳采樣點(diǎn)，可獲得最大幀間穩(wěn)定度[3]。

三是基于圖像的離線門控方法：基于 ECG門控技術(shù)在抑制運(yùn)動(dòng)偽像方面的不足，近年來，運(yùn)用數(shù)字圖像處理技術(shù)，從連續(xù)回撤超聲導(dǎo)管采集的IVUS圖像序列中提取出隱含的心動(dòng)周期信息，并據(jù)此從每個(gè)心動(dòng)周期中選取一幀圖像，組成同步采樣序列，即基于圖像的門控(image-based gating)技術(shù)[2-8]，逐漸成為該領(lǐng)域的研究熱點(diǎn)。這種方法雖然可以避免離線 ECG門控方法的缺點(diǎn)，但是由于每個(gè)心動(dòng)周期只選擇一幀，需要拋棄大量幀，因而可能會(huì)丟失很多有診斷價(jià)值的信息。

目前，對非門控IVUS圖像序列直接進(jìn)行運(yùn)動(dòng)估計(jì)的方法包括光流法(optical flow，OF)[9]和塊匹配法(block matching，BM)[10]。此類方法既不需要門控圖像采集裝置，無需記錄ECG信號，也不需要拋棄有用幀，可保證圖像數(shù)據(jù)集合的完整性。但是光流法的估計(jì)結(jié)果極易受到血流隨機(jī)運(yùn)動(dòng)的干擾，且對噪聲十分敏感；塊匹配法需要明顯的標(biāo)志物（例如鈣化，血管分叉等），而此類標(biāo)志并非在每幀中都存在。

本文在綜合分析IVUS圖像序列中由周期性心臟運(yùn)動(dòng)所致運(yùn)動(dòng)偽像的產(chǎn)生機(jī)制和表現(xiàn)形式的基礎(chǔ)上，對覆蓋多個(gè)心動(dòng)周期的非門控 IVUS圖像序列直接進(jìn)行剛性運(yùn)動(dòng)偽像的定量估計(jì)和補(bǔ)償。首先建立橫向視圖中的血管壁在心動(dòng)周期中的剛性運(yùn)動(dòng)模型；然后定量估計(jì)相鄰幀間的血管橫斷面剛性運(yùn)動(dòng)參數(shù)（包括位移和旋轉(zhuǎn)角），并采用譜分析的方法分離運(yùn)動(dòng)分量與幾何分量；最后通過補(bǔ)償運(yùn)動(dòng)分量，達(dá)到抑制運(yùn)動(dòng)偽像的目的。無需利用 ECG信號，也不需要拋棄任何一幀圖像，在抑制運(yùn)動(dòng)偽像的同時(shí)保證圖像數(shù)據(jù)集合的完整性。

1 方 法

1.1 運(yùn)動(dòng)偽像的產(chǎn)生機(jī)制和表現(xiàn)形式

IVUS圖像序列中由周期性心臟運(yùn)動(dòng)和搏動(dòng)血流所致的運(yùn)動(dòng)偽像的主要表現(xiàn)形式為：

1） 幀間的錯(cuò)位，即橫斷面圖像序列（圖1(a)）中相鄰幀之間血管橫截面的平移和旋轉(zhuǎn)（圖1(b)）；

2） 在沿管腔長軸方向的縱向視圖中，血管壁邊緣呈現(xiàn)鋸齒形（圖1(c)）。

圖2說明了運(yùn)動(dòng)偽像的產(chǎn)生機(jī)制。假設(shè)對一段均勻直徑的直血管段進(jìn)行IVUS成像，在回撤超聲導(dǎo)管的過程中，導(dǎo)管始終位于各幀斷層圖像的中心。若無運(yùn)動(dòng)偽像，且假設(shè)導(dǎo)管回撤路徑與管腔軸線重合，則將各幀圖像按照采集順序依次排列得到的縱向視圖中，血管壁的上下輪廓應(yīng)近似為直線（圖 2(a)）。若血管處于周期性運(yùn)動(dòng)狀態(tài)，則可能導(dǎo)致導(dǎo)管和管腔之間的相對運(yùn)動(dòng)。因而在回撤導(dǎo)管的過程中，采集到的各幀斷層圖像記錄的并非血管的真實(shí)狀態(tài)，在縱向視圖中將出現(xiàn)“鋸齒樣”的血管壁上下輪廓，如圖2(b)所示。這種偽像不僅會(huì)影響對血管腔和斑塊形態(tài)的觀察，同時(shí)也會(huì)給形態(tài)參數(shù)的準(zhǔn)確測量和血管的三 維重建帶來很大困難。

圖1 IVUS圖像序列

圖2 IVUS圖像序列示意圖

1.2 建立血管壁剛性運(yùn)動(dòng)模型

IVUS斷層圖像序列中相鄰切片之間管腔橫截面重心的位移和空間方向的改變主要由運(yùn)動(dòng)和幾何兩方面的因素造成[3,11-13]：

1） 運(yùn)動(dòng)：即外部因素，指由周期性心臟運(yùn)動(dòng)和搏動(dòng)的血流造成的超聲導(dǎo)管相對于管腔的運(yùn)動(dòng)和血管形態(tài)的變化；

2） 幾何：即內(nèi)部因素，指血管腔本身不規(guī)則的幾何形狀。

從各幀斷層圖像中分割出血管腔邊界之后，相鄰幀之間管腔橫截面的平移和旋轉(zhuǎn)可分別用管腔邊界曲線重心的位移(,)xyΔΔ和邊界曲線之間的旋轉(zhuǎn)角αΔ來表示，它們分別由兩部分組成

其中，腳標(biāo)d和g分別表示運(yùn)動(dòng)分量和幾何分量。

如圖3所示，設(shè)相鄰兩幀中的管腔邊界曲線分別為1γ和2γ，其重心分別為，在運(yùn)動(dòng)和幾何兩方面因素的共同作用下，1γ和2γ之間的位移為，1γ和2γ之間的旋轉(zhuǎn)角為Δα，旋轉(zhuǎn)中心為導(dǎo)管(x0,y0)，即

其中，Δx、Δy和Δα如式(1)所示分別由運(yùn)動(dòng)和幾何分量兩部分組成。對于圖像序列中的各相鄰幀，本文通過計(jì)算出位移(Δx,Δy)和旋轉(zhuǎn)角Δα，并從中分離出Δxd、Δyd和Δαd，完成對運(yùn)動(dòng)偽像的補(bǔ)償。

圖3 相鄰幀IVUS圖像運(yùn)動(dòng)參數(shù)示意圖

1.3 運(yùn)動(dòng)偽像的估計(jì)和補(bǔ)償

1.3.1 提取IVUS圖像中的血管壁邊界

對IVUS圖像進(jìn)行分割，提取血管內(nèi)腔邊界和血管壁中-外膜邊界，不僅是本文抑制運(yùn)動(dòng)偽像方法的重要步驟，而且也是量化分析[14]、三維重建[15]以及 IVUS與其它冠狀動(dòng)脈圖像融合[16]的基礎(chǔ)和保證其精度的關(guān)鍵。近年來，提出了很多計(jì)算機(jī)輔助的自動(dòng)分割方法[17-18]。考慮到處理效率等因素，本文采用一種三維分割方法[19]，與二維分割方法是對圖像序列進(jìn)行逐幀串行處理相比，三維分割技術(shù)可實(shí)現(xiàn)對圖像序列的并行處理，提高處理效率。

1.3.2 計(jì)算運(yùn)動(dòng)參數(shù)

對于圖像序列中的各相鄰幀，例如第k-1和k幀（k=2,3,…,M，M是IVUS圖像序列的總幀數(shù)），計(jì)算血管內(nèi)腔邊界曲線的幾何中心作為對其重心的近似，并如§1.2中所述方法計(jì)算出重心之間的位移量和旋轉(zhuǎn)角和kαΔ包含運(yùn)動(dòng)分量和幾何分量)兩部分，需從中分離出運(yùn)動(dòng)分量,dkxΔ、,dkyΔ和,dkαΔ。

1.3.3 分離運(yùn)動(dòng)分量和幾何分量

本文對計(jì)算出的kxΔ、kyΔ和kαΔ（k=2,3,…,M）分別進(jìn)行傅立葉變換，得到其幅度譜ΔXk、ΔYk和ΔAk。由于由血管本身的不規(guī)則幾何形狀所引起的相鄰幀之間管腔橫截面空間方向和重心位置的變化速度，遠(yuǎn)小于由周期性心臟運(yùn)動(dòng)所致的變化速度，因此在ΔXk、ΔYk和ΔAk中，高頻分量對應(yīng)于剛性運(yùn)動(dòng)分量，且其變化頻率應(yīng)等于心率，而低頻分量對應(yīng)于幾何分量 ( Δxk,g, Δyk,g, Δαk,g)。分別對ΔXk、ΔYk和ΔAk進(jìn)行高通濾波，則濾波器的輸出即是剛性運(yùn)動(dòng)分量的幅度譜(ΔXk,d, ΔYk,d, Δ Ak,d)，對其進(jìn)行逆傅立葉變換即可得到 ( Δxk,d, Δyk,d, Δαk,d)的估計(jì)值。

由于運(yùn)動(dòng)分量 ( Δxk,d,Δyk,d, Δαk,d)主要由周期性心臟運(yùn)動(dòng)引起，因此本文將高通濾波器的通帶截止頻率設(shè)定為病人的心率值（單位：次/秒，即Hz）。采用文獻(xiàn)[6]中的方法從IVUS圖像序列中估計(jì)平均心率的近似值，原理如下：

首先，對IVUS圖像序列進(jìn)行逐幀比較，通過計(jì)算兩幀圖像之間灰度值的歸一化互相關(guān)(NCC)

計(jì)算各幀圖像之間的差異值

式中，i,j= 1 ,2,… ,M；Ii和Ij分別表示IVUS圖像序列中的第i幀和第j幀，其大小均為N1×N2像素，平均灰度值分別為μi和μj；M是圖像序列的總幀數(shù)。

然后，計(jì)算間隔為k幀的兩幀圖像Im和Im+k（m= 1,2,… ,M-k）的平均差異值

得到曲線～k，k=0,1,2,… ,M-1，且= 0 。該曲線具有近似周期的形狀，這是由心臟的周期性運(yùn)動(dòng)所造成的，該曲線的重復(fù)頻率，即其幅度譜曲線峰值所對應(yīng)的頻率即為病人心率的近似值R（單位：次/分鐘，即Hz）。

1.3.4 補(bǔ)償運(yùn)動(dòng)偽像

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

為了驗(yàn)證方法的可行性并估計(jì)精度，本文分別利用計(jì)算機(jī)模擬的血管腔輪廓序列和臨床采集的IVUS圖像序列進(jìn)行實(shí)驗(yàn)，并對實(shí)驗(yàn)結(jié)果進(jìn)行分析。

2.1 模擬圖像實(shí)驗(yàn)

圖4 7幀模擬管腔輪廓序列

為了評價(jià)§1.3中所述對血管腔輪廓序列剛性運(yùn)動(dòng)參數(shù)的估計(jì)精度，本文將用封閉曲線表示的管腔輪廓序列按照已知的位移和旋轉(zhuǎn)角進(jìn)行平移和旋轉(zhuǎn)，模擬生成運(yùn)動(dòng)參數(shù)已知的管腔輪廓序列。圖4(a)和圖4 (b)分別為7幀模擬管腔輪廓的橫向視圖和完成重心計(jì)算的縱向排列圖，圖4(c) ～ (e)為各偏移量的柱狀圖。

圖5和圖6分別為對20幀和35幀模擬管腔輪廓曲線序列進(jìn)行剛性運(yùn)動(dòng)補(bǔ)償前后的橫截面視圖和縱向視圖?？梢钥闯?，完成運(yùn)動(dòng)補(bǔ)償后的橫截面視圖中相鄰幀之間的扭動(dòng)大大減小，在縱向視圖中體現(xiàn)為管腔邊界的鋸齒形尖峰與水平基線之間的距離大大縮小。

圖6 35幀模擬圖像實(shí)驗(yàn)結(jié)果

2.2 臨床圖像實(shí)驗(yàn)

臨床圖像序列是采用Jomed Endosonic超聲成像儀采集的，采用2.9F 30MHz單軌機(jī)械超聲探頭，由馬達(dá)驅(qū)動(dòng)自動(dòng)回撤探頭導(dǎo)管，回撤速度為0.5mm/s，當(dāng)探頭以1800轉(zhuǎn)/分作360°旋轉(zhuǎn)時(shí)，可以30幀/秒的速率連續(xù)獲得血管橫軸實(shí)時(shí)切面圖像。

圖7 62幀臨床圖像實(shí)驗(yàn)結(jié)果

圖7和圖8分別為62幀和80幀臨床圖像序列的運(yùn)動(dòng)偽像補(bǔ)償前后的橫向視圖（血管腔邊界曲線）及縱向視圖?？梢钥闯觯瓿蛇\(yùn)動(dòng)偽像抑制后的橫截面視圖中相鄰幀管腔輪廓之間的扭動(dòng)大大減小，在縱向視圖中體現(xiàn)為外膜輪廓的鋸齒形邊緣變的相對平滑，因而很大程度上抑制了圖像序列中存在的運(yùn)動(dòng)偽像。

圖8 80幀臨床圖像實(shí)驗(yàn)結(jié)果

由于目前還沒有一種客觀、定量評價(jià)抑制IVUS運(yùn)動(dòng)偽像方法精度的體系[6]，因此，我們將采用本文方法進(jìn)行運(yùn)動(dòng)偽像抑制的結(jié)果和采用基于圖像的離線門控[8]結(jié)果進(jìn)行比較，評價(jià)本文方法的精度。圖9是分別對120幀和200幀的臨床圖像序列進(jìn)行實(shí)驗(yàn)的結(jié)果，可見雖然基于圖像的離線門控法可實(shí)現(xiàn)對運(yùn)動(dòng)偽像的抑制，但是門控序列的長度相對于原序列而言明顯變短，即丟棄了大量的幀數(shù)據(jù)。而本文方法在完成運(yùn)動(dòng)偽像抑制的同時(shí)，不丟棄任何一幀數(shù)據(jù)，保證了原始圖像數(shù)據(jù)集合的完整性。

2.3 誤差源分析

本文的運(yùn)動(dòng)偽像抑制方法包括圖像的采集和預(yù)處理、血管腔輪廓提取、管腔輪廓重心及偏移量的計(jì)算、運(yùn)動(dòng)分量的濾出、對運(yùn)動(dòng)分量的補(bǔ)償?shù)炔襟E。除去計(jì)算誤差之外，誤差主要來自兩方面：

一是圖像的分割：由于本文方法是直接對分割出的管腔輪廓進(jìn)行運(yùn)動(dòng)估計(jì)的，因而IVUS圖像的分割精度直接影響對運(yùn)動(dòng)偽像的補(bǔ)償效果。因此，進(jìn)一步提高二維輪廓提取的準(zhǔn)確性有助于提高運(yùn)動(dòng)偽像的補(bǔ)償精度。

二是濾波器的設(shè)計(jì)：由于本文采用高通濾波的方法分離由心臟的周期性運(yùn)動(dòng)引起的偏移量，因而高通濾波器的設(shè)計(jì)（包括原型的選擇和參數(shù)的設(shè)定）直接影響對運(yùn)動(dòng)偽像的抑制效果?？紤]到運(yùn)算復(fù)雜度等因素，本文選擇了 Butterworth濾波器，并將通帶截止頻率設(shè)定為從IVUS圖像序列中估計(jì)出來的平均心率值。對于濾波器原型的選擇，文獻(xiàn)[3]首次探討了濾波器的通帶形狀對采樣精度的影響，并提出采用統(tǒng)計(jì)的方法確定哪種濾波器更適合于心臟時(shí)相信息的提取。在后續(xù)工作中，我們將著重分析濾波器的性能對于運(yùn)動(dòng)偽像抑制效果的影響。對于平均心率值的估算，若圖像采集過程中同時(shí)記錄了 ECG信號，則可直接利用ECG記錄的心率值；若無ECG信號可供利用，則采用文獻(xiàn)[6]的方法，在病人心率保持恒定的情況下，可得到精確的估算結(jié)果。但是若病人的心率發(fā)生波動(dòng)，則其估計(jì)結(jié)果存在一定誤差，詳細(xì)討論參見文獻(xiàn)[6]。

圖9 臨床IVUS圖像序列運(yùn)動(dòng)偽像抑制實(shí)驗(yàn)結(jié)果

3 結(jié) 論

本文提出了一種抑制冠狀動(dòng)脈內(nèi)超聲圖像序列中由心臟運(yùn)動(dòng)所致運(yùn)動(dòng)偽像的方法，與現(xiàn)有的心電門控方法（包括聯(lián)機(jī)和脫機(jī)兩種方式）相比較，既不需要專門的心電門控圖像采集裝置，也不需要同步記錄心電信號，而是運(yùn)用圖像分析技術(shù)，對連續(xù)回撤超聲導(dǎo)管采集到的IVUS圖像序列數(shù)據(jù)進(jìn)行處理，因而適用于沒有記錄 ECG信號的IVUS記錄序列。與基于圖像的門控方法相比，本文方法不需要拋棄圖像序列中的有用幀，可在抑制運(yùn)動(dòng)偽像的同時(shí)，保證圖像數(shù)據(jù)集合的完整性。

[1]胡廣祿, 陳興新. 心電門控技術(shù)及其應(yīng)用[J]. 醫(yī)療裝備, 2000, 13(1): 17.

[2]de Winter S A, Hamers R, Degertekin M, et al. A novel retrospective gating method for intracoronary ultrasound images based on image properties [C]//Proceedings of IEEE International Conference on Computers in Cardiology, Kassandra Chalkidiki,Greece: IEEE Press, 2003, 30: 13-16.

[3]Hernàndez-Sabaté A, Gil D, Garcia-Barnés J, et al.Image-based cardiac phase retrieval in intravascular ultrasound sequences [J]. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control,2011, 58(1): 60-72.

[4]Zhu H, Oakeson K D, Friedman M H. Retrieval of cardiac phase from IVUS sequences [C]//Proceedings of SPIE Conference on Medical Imaging 2003:Ultrasonic Imaging and Signal Processing, Bellingham WA: SPIE Press, 2003, 5035: 135-146.

[5]Nadkarni S K, Boughner D R, Fenster A. Image-based cardiac gating for three-dimensional intravascular ultrasound [J]. Ultrasound in Medicine and Biology,2005, 31(1): 53- 63.

[6]O’Malley S M, Granada J F, Carlier S, et al.Image-based gating of intravascular ultrasound pullback sequences [J]. IEEE Transactions on Information Technology in Biomedicine, 2008, 12(3):299-306.

[7]Carlo G, Oriol P, Oriol R L, et al. Robust image-based IVUS pullbacks gating [C]//Proceedings of 11th International Conference on Medical Image Computing and Computer-Assisted Intervention Part II(MICCAI 2008), New York City, USA: MICCAI Society Press, 2008, 5242: 518-525.

[8]孫 正. 抑制冠狀動(dòng)脈內(nèi)超聲圖像序列運(yùn)動(dòng)偽影的離線門控方法[J]. 光電子·激光. 2010, 21(4):632-638.

[9]Danilouchkine M G, Mastik F, Steen A F W. Accuracy in prediction of catheter rotation in IVUS with feature-based optical flow——a phantom study [J].IEEE Transactions on Information Technology in Biomedicine, 2008, 12(3): 356-365.

[10]Danilouchkine M G, Mastik F, Steen A F W. Two methods for catheter motion correction in IVUS palpography [C]// Proceedings of 2007 IEEE International Ultrasonic Symposium (IUS), New York NY: IEEE Press, 2007: 1724-1727.

[11]Nadkarni S K, Mills G, Boughner D R, et al. A pulsating coronary vessel phantom for two-and threedimensional intravascular ultrasound studies [C]//Proceedings of the 4th International Conference on Medical Image Computing and Computer-Assisted Intervention (MICCAI 2001), Utrecht, Netherlands:MICCAI Society Press, 2001, 2208: 1164-1165.

[12]Holzapfel G, Gasser T, Stadler M. A structural model for the viscoelastic behavior of arterial walls:continuum formulation and finite element analysis [J].European Journal of Mechanics A-Solids, 2002, 21:441-463.

[13]Rosales M, Radeva P, Rodriguez-Leor O, et al.Modeling of image-catheter motion for 3-D IVUS [J].Medical Image Analysis, 2009, 13(1): 91-104.

[14]Baldewsing R A, Danilouchkine M G, Mastik F, et al.An inverse method for imaging the local elasticity of atherosclerotic coronary plaques [J]. IEEE Transactions on Information Technology in Biomedicine, 2008, 12(3): 277-289.

[15]Paul S. 3-D Intravascular ultrasound (IVUS) and IVUS-Palpography: insights into the mechanical behavior of the coronary vessel wall [J]. The International Journal of Cardiovascular Imaging,2006, 22: 153-155.

[16]Sun Zheng. Three-dimensional reconstruction of vessels from intravascular ultrasound sequence and X-ray angiograms [C]//Proceedings of 2nd International Conference on BioMedical Engineering and Informatics (BMEI 2009), Tianjin,China: IEEE Press, 2009, 1: 444-448.

[17]Bovenkamp E G P, Dijkstra J, Bosch J G,et al.User-agent cooperation in multiagent IVUS image segmentation [J]. IEEE Transactions on Medical Imaging, 2009, 28(1): 94-105.

[18]Zhang Qi, Wang Yuanyuan, Wang Weiqi, et al.Automatic segmentation of calcifications in intravascular ultrasound images using snakes and the contourlet transform [J]. Ultrasound in Medicine &Biology, 2010, 36(1): 111-129.

[19]孫 正, 楊 宇. 一種基于snake模型的IVUS圖像序列三維分割方法[J]. 工程圖學(xué)學(xué)報(bào), 2011, 22(6):25-32.