基于先驗(yàn)估計(jì)和相關(guān)相位的實(shí)時(shí)超聲彈性成像

崔少國(guó) 彭彩碧 胡 明 劉向陽(yáng)

1(湖北醫(yī)藥學(xué)院計(jì)算機(jī)中心，十堰 442300)

2(湖北醫(yī)藥學(xué)院附屬太和醫(yī)院，十堰 442300)

引言

超聲彈性成像是生物組織成像的新模式，是當(dāng)前超聲醫(yī)學(xué)工程學(xué)領(lǐng)域的前沿和熱點(diǎn)研究。該成像技術(shù)可以探測(cè)出人體組織與病理變化相關(guān)的硬度改變，對(duì)腫瘤(癌癥)的早期檢測(cè)和診斷提供重要依據(jù)。彈性成像研究自1991年被Ophir等提出以來(lái)的20年中取得了快速發(fā)展，在臨床中已逐漸顯示出重要診斷價(jià)值，如在乳房癌、前列腺癌、甲狀腺結(jié)節(jié)、肝纖維化、動(dòng)脈粥樣硬化等的檢測(cè)中取得了重要應(yīng)用［1－4］。

基于時(shí)間域的互相關(guān)成像技術(shù)主要分為兩步:位移估計(jì)和應(yīng)變計(jì)算。位移估計(jì)主要是通過(guò)互相關(guān)技術(shù)從組織壓縮后的回波信號(hào)里找到與壓縮前匹配的信號(hào)，匹配窗之間的時(shí)延即代表著組織的位移;應(yīng)變計(jì)算主要是通過(guò)位移的梯度操作或差分濾波器來(lái)實(shí)現(xiàn)。這兩步中位移估計(jì)是關(guān)鍵［5］。為了達(dá)到亞采樣點(diǎn)級(jí)精度位移估計(jì)，Céspedes等提出了拋物線插值、余弦插值和重構(gòu)插值等方法［6］。這些方法在一定程度上提高了位移估計(jì)精度，但插值操作計(jì)算代價(jià)比較大，降低了彈性成像的實(shí)時(shí)性。另外，位移估計(jì)的前提是從壓縮后的回波信號(hào)里去搜索壓縮前信號(hào)窗運(yùn)動(dòng)到的新位置，即運(yùn)動(dòng)追蹤。使用窮舉式搜索或啟發(fā)式窮舉搜索進(jìn)行運(yùn)動(dòng)追蹤搜索范圍都比較大，是較耗時(shí)的追蹤策略［7］。針對(duì)基于互相關(guān)成像算法中存在的這些問(wèn)題，本研究提出了使用先驗(yàn)估計(jì)在較小的范圍內(nèi)進(jìn)行運(yùn)動(dòng)追蹤，使用復(fù)互相關(guān)的相位進(jìn)行亞采樣點(diǎn)級(jí)位移精度估計(jì)的新算法。

1 算法

提出的新算法主要分為兩部分:基于先驗(yàn)估計(jì)的運(yùn)動(dòng)追蹤和基于相關(guān)相位的亞采樣點(diǎn)級(jí)形變估計(jì)?；谙闰?yàn)估計(jì)運(yùn)動(dòng)追蹤的基本思想是根據(jù)連續(xù)性原則，物理位置上相鄰點(diǎn)其運(yùn)動(dòng)的位移也是近似的，使用相鄰位置上已估計(jì)出來(lái)的位移做為它的初始種子位移，然后在其附近較小的范圍內(nèi)追蹤從而找到其真正的位移。基于相關(guān)相位的亞采樣點(diǎn)級(jí)精度位移估計(jì)是先使用互相關(guān)技術(shù)在基帶信號(hào)上進(jìn)行位移初步估計(jì)(只精確到整數(shù)個(gè)采樣點(diǎn))，然后在該點(diǎn)使用復(fù)互相關(guān)函數(shù)的相位角進(jìn)行亞采樣點(diǎn)級(jí)精度位移估計(jì)。

1.1 運(yùn)動(dòng)追蹤

1.1.1 種子位移選擇

先驗(yàn)估計(jì)引導(dǎo)追蹤可以使用軸向上相鄰的上一已估計(jì)點(diǎn)的位移做為該點(diǎn)的種子位移(軸向引導(dǎo))，也可以使用側(cè)向上左邊或右邊相鄰已估計(jì)點(diǎn)的位移作為該點(diǎn)的種子位移(側(cè)向引導(dǎo))。軸向引導(dǎo)算法示意圖如圖1(a)。用seed(i，j)表示第(i，j)個(gè)估計(jì)位置的初始種子位移，u(i，j)表示第(i，j)個(gè)估計(jì)位置的位移，則種子選擇方法為:

位移提煉則使用:

式中，2≤i≤M，1≤j≤N，M 和 N 是位移估計(jì)圖像的行和列維數(shù)。－S≤r≤S，S是相鄰兩個(gè)估計(jì)位置可能的最大位移差，［－S，S］是局部搜索范圍，搜索步長(zhǎng)1個(gè)采樣點(diǎn)。Coefficient＜＞是計(jì)算互相關(guān)函數(shù)的系數(shù)，argmax()是求當(dāng)互相關(guān)系數(shù)值最大時(shí)局部位移值r，round()是四舍五入函數(shù)。

側(cè)向引導(dǎo)算法如圖1(b)所示，假設(shè)最中間一列位移已被估計(jì)，則種子選擇可按如下公式進(jìn)行

位移提煉同樣使用公式(2)進(jìn)行。

1.1.2 局部搜索范圍確定

種子位移確定后，需要確定局部搜索范圍。局部搜索范圍可以依據(jù)下面公式進(jìn)行:

圖1 軸向引導(dǎo)和側(cè)向引導(dǎo)運(yùn)動(dòng)追蹤示意(～u表示種子位移，?代表以～u為種子估計(jì)出來(lái)的位移，空心箭頭代表種子選擇來(lái)源，實(shí)心箭頭代表種子傳遞方向)。(a)軸向引導(dǎo);(b)側(cè)向引導(dǎo)。Fig.1 Illustrations for various guidance tracking strategies，where～u denotes the seed displacement and ? denotes the estimated displacement using the～u as seed.(a)Axial Guidance;(b)Lateral Guidance.

式中，Rup、Rdown、Rleft、Rright分別是局部搜索范圍的向上、向下、向左、向右的維數(shù)。εmax是系統(tǒng)檢測(cè)到的最大可能應(yīng)變(屬于系統(tǒng)設(shè)定參數(shù)，一般為6%)，Axial_shift和Lateral_shift分別是估計(jì)窗軸向和側(cè)向移動(dòng)的維數(shù)。一維軸向運(yùn)動(dòng)追蹤時(shí)只采用式(4)，二維運(yùn)動(dòng)追蹤時(shí)采用式(4)和式(5)。

1.1.3 局部啟發(fā)式搜索

如果組織的全局運(yùn)動(dòng)方向已知，局部搜索范圍可進(jìn)一步縮小。全局運(yùn)動(dòng)方向的獲取可以先隨機(jī)地選取幾個(gè)估計(jì)點(diǎn)，這些點(diǎn)的位移可通過(guò)啟發(fā)式窮舉搜索獲得。然后根據(jù)這些位移的正負(fù)號(hào)確定組織全局運(yùn)動(dòng)方向。如果向上(壓縮)，則局部搜索范圍

如果向下(拉伸)，則

向左和向右搜索維數(shù)Rleft、Rright不變，使用式(5)。

1.2 位移估計(jì)

設(shè)壓縮前的回波信號(hào)為x(t)，壓縮后的回波信號(hào)為y(t)。則x(t)和y(t)的復(fù)互相關(guān)函數(shù)為

2-D離散化的歸一化復(fù)互相關(guān)函數(shù)為

根據(jù)式(9)，當(dāng)復(fù)互相關(guān)函數(shù) f(i，j，d)的值最大時(shí)所求得的d即為整數(shù)采樣點(diǎn)級(jí)位移估計(jì)。那么亞采樣點(diǎn)級(jí)精度的位移估計(jì)u為

式中，angle()是求復(fù)互相關(guān)函數(shù) f(i，j，d)的相位角，ω0是中心角頻率。

1.3 應(yīng)變計(jì)算

應(yīng)變計(jì)算使用羅建文等提出的數(shù)字低通差分濾波器SG-I DD［8］，軸向應(yīng)變計(jì)算方法為

式中，s表示應(yīng)變，u表示位移，M是濾波器參數(shù)(濾波長(zhǎng)度為2M+1)，Δk=u(n+k)－u(n－k)。

2 實(shí)驗(yàn)

在MatLab 7.0實(shí)現(xiàn)上述提出的成像算法，并在仿真數(shù)據(jù)集和彈性體模上獲取的真實(shí)回波數(shù)據(jù)集上運(yùn)行，比較提出的算法和經(jīng)典的互相關(guān)插值法的成像速度和成像質(zhì)量。插值法選用拋物線插值，比較的具體項(xiàng)目包括成像一幀應(yīng)變圖像所需要的時(shí)間、成像效果、應(yīng)變圖像的性噪比SNRe。成像所使用的計(jì)算機(jī)是一臺(tái) Pentium(R)4，CPU 2.4 GHz，RAM 1 GB，操作系統(tǒng)是 Microsoft Windows XP。

2.1 仿真

本實(shí)驗(yàn)仿真中使用的中心頻率是5 MHz，采樣頻率是20 MHz，聲速設(shè)定為常數(shù)1 540 m/s，兩種模型被仿真:均勻彈性組織模型和含有兩個(gè)硬包容物的組織模型。組織大小均為40 mm×40 mm，包容物的直徑是10 mm。散射子濃度設(shè)為45(個(gè)散射子/脈沖寬度)，硬包容物的彈性模量是背景的10倍。1D應(yīng)變模型被采用(設(shè) dy=0)。仿真的每幀數(shù)據(jù)是1039(采樣點(diǎn))×128(線)，然后加高斯白噪聲使其性噪比為50 dB。對(duì)每種模型產(chǎn)生不同整體應(yīng)變的數(shù)據(jù)，應(yīng)變分別是:0.01%，0.1%，0.5%，1%，2%，3%，4%，5%。

2.2 體模實(shí)驗(yàn)

為了驗(yàn)證提出的算法的真實(shí)成像效果，使用彈性體模上采集的超聲回波數(shù)據(jù)進(jìn)行實(shí)驗(yàn)。體模采用CIRS Inc.研發(fā)的專用于彈性成像研究的Model 049彈性體模。體模里在背景介質(zhì)中放有4個(gè)10 mm和4個(gè)20 mm直徑的球形包容物，球體和背景是等回聲的，在B模式圖像中不能區(qū)分。背景的彈性模量是30 kPa，，每種直徑的4個(gè)包容物各有不同的彈性模量，它們分別是 8 kPa(TypeⅠ，軟)，18 kPa(TypeⅡ，軟)，44 kPa(TypeⅢ，硬)和 63 kPa(TypeⅣ，硬)。實(shí)際體模中不同直徑的包容物在軸向(上下方向，即體模壓縮方向)上是錯(cuò)開(kāi)的。本實(shí)驗(yàn)中成像對(duì)象選擇直徑是10 mm、彈性模量是63 kPa的較硬包容物。

實(shí)驗(yàn)采用的探頭是型號(hào)為SA5L38B的128陣元的線陣探頭，中心頻率是5 MHz，75%分?jǐn)?shù)階帶寬。實(shí)驗(yàn)使用的超聲系統(tǒng)是聲泰特(成都)公司研發(fā)的iMago C21超聲機(jī)，系統(tǒng)采樣頻率設(shè)為 40 MHz。使用手持探頭均勻軸向壓縮體模，在壓縮過(guò)程中共采集序列回波300幀。RF信號(hào)使用積分解調(diào)器解調(diào)到基帶I/Q信號(hào)，然后下采樣到10 MHz。幀間應(yīng)變大約0.7%，幀大小是 512(采樣點(diǎn))×188(線)。

3 結(jié)果及分析

3.1 時(shí)間比較

新算法與經(jīng)典的互相關(guān)插值算法的時(shí)間性能比較分別如表1和表2。從表中可看出，不論是仿真還是體模實(shí)驗(yàn)，新算法比經(jīng)典算法成像速度都大大提高，這主要是因?yàn)槭褂昧讼闰?yàn)位移進(jìn)行位移估計(jì)，使運(yùn)動(dòng)追蹤的范圍大大減小。經(jīng)典互相關(guān)插值算法實(shí)時(shí)性較差，每秒只能達(dá)到約2～3幀，而新算法實(shí)時(shí)性較好，每秒可以達(dá)30～50幀，可滿足臨床超聲成像的實(shí)時(shí)性要求。

3.2 成像效果

仿真成像效果如圖2所示。使用含兩個(gè)硬包容物仿真模型，成像中估計(jì)窗大小為40個(gè)采樣點(diǎn)，窗軸向重疊率0.75。從圖2可以看出，新算法所得到的應(yīng)變圖像(圖2(b))比較細(xì)膩平滑，噪聲顯著減少。這主要是因?yàn)椴蓸酉嚓P(guān)相位進(jìn)行亞采樣點(diǎn)級(jí)精度位移估計(jì)比拋物線插值法更準(zhǔn)確。

表1 仿真實(shí)驗(yàn)中經(jīng)典算法和新算法的時(shí)間性能比較Tab.1 Comparison of computation time for conventional and presented algorithms in simulations

表2 體模實(shí)驗(yàn)中經(jīng)典算法和新算法的時(shí)間性能比較Tab.2 Comparison of computation time for conventional and presented algorithms in phantom

圖2 仿真彈性成像。(a)拋物線插值法;(b)復(fù)相關(guān)相位法Fig.2 Simulation strain images. (a)Parabola interpolation;(b)Cross－correlation phase

體模成像效果如圖3所示，成像中估計(jì)窗大小20個(gè)采樣點(diǎn)，窗軸向重疊率0.75。從圖3可以看出，體模實(shí)驗(yàn)結(jié)果與仿真實(shí)驗(yàn)結(jié)果一致，新算法得到的應(yīng)變圖像(圖3(b))有更好的性能，彈性圖像中的噪聲明顯減少。

圖3 體模彈性成像。(a)拋物線插值法;(b)復(fù)相關(guān)相位法Fig.3 Phantom strain images. (a)Parabola interpolation;(b)Cross－correlation phase

3.3 應(yīng)變圖像的性噪比

圖4顯示了在仿真的均勻彈性組織模型數(shù)據(jù)集上應(yīng)用提出的新算法所得應(yīng)變圖像的性噪比與拋物線插值法的比較。性噪比計(jì)算公式為

式中μs和 σs分別是均勻彈性組織應(yīng)變圖像的應(yīng)變平均值和應(yīng)變標(biāo)準(zhǔn)方差。圖中曲線上每一個(gè)點(diǎn)是5次獨(dú)立仿真實(shí)驗(yàn)數(shù)據(jù)的平均值。圖中最上面的一條曲線是使用相關(guān)相位(新算法)法在各種應(yīng)變數(shù)據(jù)上成像的性噪比，中間一條是拋物線插值法對(duì)應(yīng)的曲線，最下面一條是沒(méi)有進(jìn)行亞采樣點(diǎn)級(jí)位移估計(jì)(估計(jì)精度是整數(shù)個(gè)采樣點(diǎn))的產(chǎn)生的曲線。從圖中可以看出當(dāng)應(yīng)變s＜1%時(shí)，性噪比隨應(yīng)變減小而減小，這主要是因?yàn)閼?yīng)變較小時(shí)，回波信號(hào)的噪聲是應(yīng)變估計(jì)的主要噪聲源，回波噪聲使較小應(yīng)變的數(shù)據(jù)窗匹配產(chǎn)生較大誤差;當(dāng)應(yīng)變 s＞3%時(shí)，性噪比隨應(yīng)變?cè)龃蠖陆?，這主要是因?yàn)閼?yīng)變較大時(shí)，解相關(guān)是應(yīng)變估計(jì)的主要噪聲源，高的解相關(guān)使位移(時(shí)延)估計(jì)產(chǎn)生較大的方差。另外，從圖中可以明顯看出，在任何應(yīng)變下，新算法的性噪比SNRe最高，拋物線插值法次之，整數(shù)個(gè)采樣點(diǎn)級(jí)算法所得的性噪比最低。尤其在1%和2%的應(yīng)變時(shí)，新算法的性噪比較拋物插值法顯著提高。

圖4 拋物線插值法、相關(guān)相位法在仿真的各應(yīng)變數(shù)據(jù)上應(yīng)變成像的性噪比Fig.4 SNRe of strain images from the parabola interpolation and correlation phase algorithms in simulations

4 結(jié)論

本研究提出了一種新的實(shí)時(shí)高精度彈性成像算法，該算法根據(jù)連續(xù)性原則使用先驗(yàn)估計(jì)進(jìn)行運(yùn)動(dòng)追蹤，大大縮小運(yùn)動(dòng)追蹤范圍，提高成像速度;同時(shí)使用復(fù)相關(guān)函數(shù)相位進(jìn)行亞采樣點(diǎn)級(jí)精度位移估計(jì)，使應(yīng)變估計(jì)精度大大提高，圖像品質(zhì)明顯改善。仿真實(shí)驗(yàn)和體模實(shí)驗(yàn)顯示，新算法與傳統(tǒng)的基于互相關(guān)的拋物線插值法相比在成像速度、應(yīng)變圖像的性噪比等方面均有較大提高。實(shí)驗(yàn)證明了提出算法的正確性和高效性。該算法的提出將對(duì)組織彈性實(shí)時(shí)成像系統(tǒng)研發(fā)和臨床推廣應(yīng)用具有重要意義。

［1］Ophir J，Céspedes I，Ponnekanti H，et al.Elastography:a quantitative method for imaging the elasticity of biological tissue［J］.Ultrasonic Imaging，1991，13(3):111 －134.

［2］Céspedes I，Ophir J，Ponnekanti H，et al.Elasticity imaging using ultrasound with application to muscle and breast in vivo［J］.Ultrasonic Imaging，1993，15(2):73－88.

［3］RighettiR.Elastographic characterization ofHIFU-induced lesions in canine livers［J］. Ultrasound in Medicine and Biology，1999，25(7):1099－1113.

［4］Konofagou E，Ophir J.Myocardial elastography—A feasibility study［J］.Ultrasound in Medicine and Biology，2002，28(4):1113－1124.

［5］Lindop JE，Treece GE，Gee A，et al.Phase-based ultrasonic deformation estimation［J］.IEEE Transactions on Ultrasonics，F(xiàn)erroelectrics，and Frequency Control，2008，55(1):94 －111.

［6］Céspedes I，Huang Y， Ophir J.Methods for estimation of subsample time delays of digitized echo signals［J］.Ultrasonic Imaging，1995，17:142 －171.

［7］Chen Lujie，Treece GM， Lindop JE.A quality-guided displacement tracking algorithm for ultrasonic elasticity imaging［J］.Medical Image Analysis，2009，13:286 －296.

［8］Luo Jianwen，Bai Jing，He Ping，et al.Axial Strain Calculation Using a Low-Pass Digital Differentiator in Ultrasound Elastography ［J］. IEEE Transactions on Ultrasonics，F(xiàn)erroelectrics，and Frequency Control，2004，51(9):1119 －1127.