尤鑫

（四川大學(xué)醫(yī)學(xué)影像實(shí)驗(yàn)室，成都　610065）

聲輻射力成像位移計(jì)算改進(jìn)方法

尤鑫

（四川大學(xué)醫(yī)學(xué)影像實(shí)驗(yàn)室，成都610065）

0　引言

醫(yī)學(xué)超聲振動(dòng)性彈性成像（Vibration Sonoelastography）由Fatemi和Greenleaf于1998年發(fā)明。該技術(shù)是用一個(gè)超聲波場產(chǎn)生一個(gè)低頻振動(dòng)并作用于受檢組織，組織受到激勵(lì)根據(jù)自身的彈性模量大小產(chǎn)生不同的振動(dòng)幅度，并最終通過圖像表現(xiàn)出來[1-2]。采用聲輻射力的聲輻射脈沖成像屬于振動(dòng)性彈性成像的一種，ARFI成像測量組織的反應(yīng)是通過使用由長脈沖引起的聲輻射力產(chǎn)生局部位移從而深入探察在有機(jī)體體內(nèi)和體外的軟組織局部粘彈性性能[3-4]。位移檢測是ARFI成像中最為重要的環(huán)節(jié)，是組織生物力學(xué)特性分析的基礎(chǔ)。傳統(tǒng)的方法通常采用胡相關(guān)方法，迭代相位零（Iterative phase zero）方法或者loupas算法，組織位移的估算可以通過對激勵(lì)之前（參考）和激勵(lì)之后（跟蹤）信號進(jìn)行計(jì)算獲得[5]。Loupas算法提供了一種對軸向速度推導(dǎo)得到的速度估算，它是通過利用顯示估計(jì)在每個(gè)距離選通脈沖位置的平均多普勒頻率與平均射頻頻率來轉(zhuǎn)換多普勒方程式[6]。Loupas結(jié)合了一般的較高性能的2D寬帶時(shí)域技術(shù)與有相對溫和復(fù)雜性的1D窄帶相域并滿足被提出的估算函數(shù)計(jì)算要求的速度估算函數(shù)。

現(xiàn)有醫(yī)學(xué)超聲聲輻射力成像一般采用Loupas算法得到由聲輻射力引起的微小位移。而在傳統(tǒng)Loupas算法中，計(jì)算位移時(shí)通常只選取初始時(shí)刻和之后的各個(gè)不同時(shí)刻的采樣點(diǎn)來進(jìn)行絕對位移的計(jì)算，這種計(jì)算方式在時(shí)間方向的窗口N為2，容易受噪音影響，使位移曲線不夠連續(xù)和平滑。具體而言，由于N取值小，回波相位位移計(jì)算時(shí)間短，包含信息少，計(jì)算誤差大，從3節(jié)位移曲線可看出，計(jì)算結(jié)果抖動(dòng)較大，導(dǎo)致最后成圖效果不佳；而當(dāng)N取值過大時(shí)，回波相位位移計(jì)算時(shí)間長，且由于包含過多冗余信息，而導(dǎo)致計(jì)算結(jié)果不準(zhǔn)確，成圖效果變差，因此本文通過在成像過程中比較作用于模型焦點(diǎn)位置的回波信號不同時(shí)間取樣窗口大小，研究不同時(shí)間取樣窗口的大小下，焦點(diǎn)位置的不同的回波位移-時(shí)間曲線，選擇出最適時(shí)間取樣窗口取值，將優(yōu)化后的參數(shù)反饋給系統(tǒng)從而取得更好的成像效果。

本文著重比較了成像處理過程中在焦點(diǎn)區(qū)域不同時(shí)間取樣窗口下的回波位移-時(shí)間曲線，選擇合適的時(shí)間取樣窗口大小的值，并反饋優(yōu)化后的參數(shù)到系統(tǒng)中，從而獲得更好的成圖效果。

1　研究方法

1.1Loupas算法介紹

復(fù)雜解調(diào)信號z（m，n）下2D自相關(guān)器獲得的平均速度[7]：

其中γdem是z（m，n）的2D自相關(guān)函數(shù)；

最后，依據(jù)I＆Q采樣，公式（1）和 （2）可以結(jié)合來表示平均速度，如公式 （3）：

1.2平均位移計(jì)算

在涉及使用Loupas算法計(jì)算ARFI位移的先前文獻(xiàn)中，時(shí)間取樣窗口N通常取值為2，計(jì)算位移時(shí)，由激勵(lì)脈沖產(chǎn)生的回波信號由于受到發(fā)射長波造成的干擾而沒有計(jì)算價(jià)值，因此會將其去除，且以與參考脈沖對應(yīng)的回波信號為基準(zhǔn)，分別與剩余的檢測回波信號進(jìn)行計(jì)算。

由于要改變時(shí)間取樣窗口，且需保證計(jì)算中的每個(gè)回波信號在時(shí)間方向的間隔保持一致，因此在位移計(jì)算前，需要對回波信號進(jìn)行重采樣。本文采用的是線性插值法，分別采用相鄰檢測脈沖對應(yīng)的回波信號進(jìn)行線性插值替換3組激勵(lì)脈沖對應(yīng)的回波信號；并利用公式（3）并改變式中N的取值來分別計(jì)算相鄰的回波在時(shí)間取樣窗口N內(nèi)的平均位移速度。之后，運(yùn)用公式（5），（6）計(jì)算出回波在時(shí)間方向的具體位移：

其中ui為組織時(shí)間方向?qū)?yīng)的局部位移，fc是信號中心頻率，fs是采樣頻率，di為被測組織在時(shí)間放向?qū)?yīng)的局部位移，TPRF為脈沖重復(fù)周期，i的取值范圍為0到取樣容積數(shù)之間，d0為初始時(shí)刻，因?yàn)闆]有聲輻射力的激勵(lì)，所以d0=0。

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

二維的ARFI圖像由多個(gè)掃描位置計(jì)算得到，它們之間間隔0.5mm，且具有相同的焦點(diǎn)深度。接收到的檢測信號與參考信號通過使用Loupas算法進(jìn)行相關(guān)性計(jì)算得到由聲輻射力引起的微小位移。

在聲輻射力檢測過程中，首先發(fā)射的脈沖作為位移計(jì)算的參考信號，信號使用較高的電壓和較短的脈沖長度，接下來發(fā)射了3組激勵(lì)脈沖和檢測脈沖。單脈沖序列中包含有參考檢測脈沖、激勵(lì)脈沖（長脈沖）、檢測脈沖（短脈沖），回波信號中具有與上述脈沖一一對應(yīng)的回波信號。掃描包的大小為24，脈沖序列的示意圖見圖1。

圖1　脈沖序列的示意圖

實(shí)驗(yàn)數(shù)據(jù)采集時(shí)，探頭中心頻率為5MHz，采樣頻率為40MHz，重復(fù)頻率為8.88kHz，取樣容積數(shù)為24。實(shí)驗(yàn)采用CIRS公司的Model049模型，模型由4種力學(xué)特征不同的球形組織和彈性背景組織構(gòu)成，類型1-4的彈性系數(shù)依次為：8kPa，14kPa，45 kPa，80 kPa，背景區(qū)域?yàn)?5 kPa。本實(shí)驗(yàn)采用類型1彈性系數(shù)為8 kPa的小球進(jìn)行聲輻射力激勵(lì)與回波檢測。

２.1位移曲線

本實(shí)驗(yàn)的取樣容積數(shù)為24，所以自相關(guān)算法時(shí)間取樣窗口N的取值范圍為2≤N≤24，由于N取值過大，會加入多余信息干擾實(shí)驗(yàn)結(jié)果并影響計(jì)算效率，所以本實(shí)驗(yàn)只選取了2≤N≤10，9組數(shù)值進(jìn)行研究考察，由于N取值相近時(shí)結(jié)果相差不大，因此本文只選取了原Loupas算法、N取值為4，6，8的結(jié)果進(jìn)行展示。通過第二節(jié)中的算法分析，利用本文提出的自相關(guān)Loupas算法的改進(jìn)，由公式（5）可計(jì)算得到在實(shí)驗(yàn)環(huán)境下所選ROI內(nèi)位置由聲輻射力引起的微小位移。由于ARFI成像技術(shù)只在焦點(diǎn)區(qū)域能量最強(qiáng)、信息最為準(zhǔn)確，且為了方便研究聲輻射力所引起的位移規(guī)律，本文根據(jù)公式（7）選擇了焦點(diǎn)有效區(qū)域位置內(nèi)的點(diǎn)進(jìn)行線性平均計(jì)算并顯示其回波位移-時(shí)間曲線。

其中rD為焦點(diǎn)有效區(qū)域邊長，f#為震源深度與孔徑寬度之比，λ為波長，本實(shí)驗(yàn)中f#=1.5，因此rD為0.14cm。

圖2　焦點(diǎn)有效區(qū)域回波位移-時(shí)間曲線

圖2為原Loupas算法與改進(jìn)算法不同時(shí)間取樣窗口大小（N=4，6，8）下焦點(diǎn)有效區(qū)域的平均回波位移-時(shí)間曲線圖。由于原算法中，激勵(lì)脈沖所對應(yīng)的回波無效被刪除，所以位移曲線中本文使用線性插值的方法將3個(gè)回波所對應(yīng)的位移值插出。由圖可知，在相應(yīng)的取樣容積數(shù)目（Ensemble size）下，小球的運(yùn)動(dòng)情況均為從靜止運(yùn)到到一定位移，最后歸于平衡。而隨著N的值地不斷增加，位移曲線也更加平滑。

由第2節(jié)可知利用Loupas算法計(jì)算出的結(jié)果為回波在時(shí)間取樣窗口內(nèi)的平均位移速度，而當(dāng)位移出現(xiàn)負(fù)值時(shí)，說明小球往反向運(yùn)動(dòng)從而產(chǎn)生了位移；然而小球在取樣容積數(shù)目范圍內(nèi)，延時(shí)間方向，應(yīng)趨于平衡，而非向負(fù)向發(fā)生位移，從原算法所計(jì)算出的時(shí)間-位移曲線圖可以看出，在時(shí)間最末小球位移趨于0，為非產(chǎn)生負(fù)向位移。因此，當(dāng)位移出現(xiàn)負(fù)值時(shí)，說明時(shí)間取樣窗口過大，計(jì)算位移速度時(shí)包含過多信息，且延時(shí)間方向向后的信息信噪比低，導(dǎo)致計(jì)算出的反向平均速度過大且不準(zhǔn)確，從而產(chǎn)生負(fù)向位移，使得最終計(jì)算結(jié)果不準(zhǔn)確。

所以，當(dāng)位移曲線中出現(xiàn)負(fù)值時(shí)，說明N選取過大，使得計(jì)算結(jié)果不準(zhǔn)確，最終降低了成圖效果。由圖2可知，當(dāng)N=4時(shí)曲線平滑度較高，且在其相應(yīng)的最后時(shí)刻，小球位移趨于0（小球趨于平衡），未產(chǎn)生負(fù)向位移，為最合適的取值，成圖效果為最佳。

２.2MSATLAB仿真成圖與性能參數(shù)比較

本文對不同的時(shí)間取樣窗口大小，選擇其相應(yīng)的取樣容積范圍內(nèi)，小球位移達(dá)到最大值的時(shí)刻成像，而非像現(xiàn)有的醫(yī)學(xué)超聲聲輻射力成像一般提供的固定時(shí)刻位移成像，并與原Loupas算法的成像結(jié)果進(jìn)行比較，仿真成像結(jié)果如圖3所示。

由圖3可以看出，當(dāng)N取值較小時(shí)（原算法），成圖效果不佳；當(dāng)N取值過大后，成圖效果沒有隨著N的取值增大而線性改善，反而成像效果逐漸變差;可看出當(dāng)N取值為4時(shí)，成圖效果最佳。

為進(jìn)一步確認(rèn)是否Loupas時(shí)間取樣窗口為4時(shí)，成圖效果最佳，本文對所成圖像利用公式（8）進(jìn)行CNR[8]計(jì)算與比較：

圖3　Loupas算法仿真成像與改進(jìn)后不同時(shí)間取樣窗口對應(yīng)仿真成像

其中μt與σt分別表示目標(biāo)區(qū)域應(yīng)變的均值與標(biāo)準(zhǔn)差，μb與σb分別表示背景區(qū)域應(yīng)變的均值和標(biāo)準(zhǔn)差。

本文以球內(nèi)的區(qū)域?yàn)槟繕?biāo)區(qū)域，球外的4個(gè)區(qū)域分別為背景區(qū)域，分別計(jì)算出相應(yīng)的CNR，再將4個(gè)CNR進(jìn)行平均（5個(gè)區(qū)域范圍大小幾乎相等），如圖4:

圖4　CNR計(jì)算區(qū)域選擇示意圖

結(jié)果如表1所示。

表1　各仿真圖像的CNR值

從表1可以看出，時(shí)間取樣窗口為4時(shí)，CNR值最高，再次證明4為最佳值。

3　結(jié)語

本文研究聲輻射力的成像算法問題，對聲輻射力成像現(xiàn)使用的Loupas算法進(jìn)行了改進(jìn)，利用焦點(diǎn)位置區(qū)域中心的時(shí)間-位移曲線選取合適的算法時(shí)間取樣窗口的值，達(dá)到更好的成像效果。

由于現(xiàn)使用的Loupas算法中，時(shí)間取樣窗口N取值為最小值2，N取值小，回波相位位移計(jì)算時(shí)間短，計(jì)算誤差大，因此本文對此進(jìn)行了改進(jìn)，并研究N的最佳取值。在觀察計(jì)算得出的不同時(shí)間取樣窗口大小下焦點(diǎn)有效區(qū)域內(nèi)的平均回波位移-時(shí)間曲線圖后分析發(fā)現(xiàn)，當(dāng)N取值增加到一定值時(shí)，位移出現(xiàn)負(fù)值時(shí)，說明小球往反向運(yùn)動(dòng)從而產(chǎn)生了位移，與小球在取樣容積數(shù)目范圍內(nèi)，沿時(shí)間方向，最后時(shí)刻應(yīng)趨于平衡，而非向負(fù)向發(fā)生位移相悖，說明此時(shí)N取值過大，受過多信息干擾，使得計(jì)算結(jié)果不準(zhǔn)確，最終降低成圖效果。最終根據(jù)不同時(shí)間取樣窗口大小下焦點(diǎn)有效區(qū)域內(nèi)的平均回波位移-時(shí)間曲線圖分析，MATLAB仿真成圖效果以及CNR計(jì)算比較，得出時(shí)間取樣窗口N取值為4時(shí)，成圖效果最佳。

采用本文提供的基于聲輻射力的回波位移檢測方法，可以有效選取位移檢測過程中的時(shí)間取樣窗口，使位移檢測可靠性更高，成圖效果更好，對改進(jìn)現(xiàn)有醫(yī)學(xué)超聲聲輻射力成像系統(tǒng)，為醫(yī)生診斷提供更準(zhǔn)確結(jié)果顯示具有一定意義。本文研究只針對激勵(lì)脈沖數(shù)為3，較為局限，后續(xù)工作中要進(jìn)一步對不同的激勵(lì)脈沖數(shù)對應(yīng)的時(shí)間取樣窗口大小最佳取值進(jìn)行研究，進(jìn)行自適應(yīng)地調(diào)整位移檢測過程中時(shí)間取樣窗口。

[1]M F，JF.G.Vibro-Acoustography:an Imaging Modality Based on Ultrasound-Stimulated Acoustic Emission.[J].Proc Natl Acad Sci U S A，1999，96（12）:6603-6608.

[2]Fatemi M，Greenleaf J F.Ultrasound-Stimulated Vibro-Acoustic Spectrography[J].Science，1998，280（5360）:82-85.

[3]K N，R B，G.T.Observations of Tissue Response to Acoustic Radiation Force:Opportunities for Imaging[J].Ultrasonic Imaging，2002，24（3）:129-138.

[4]Nightingale K，Soo M.Acoustic Radiation Force Impulse Imaging:in Vivo Demonstration of Clinical Feasibility[J].Ultrasound in Medicine&Biology，2002，28（2）:227-235（9）.

[5]JR D，GE T，KR N，et al.Acoustic Radiation Force Elasticity Imaging in Diagnostic Ultrasound.[J].IEEE Trans Ultrason Ferroelectr Freq Control，2013，60（4）:685-701.

[6]Loupas T，Powers J T，Gill R W，et al.An Axial Velocity Estimator for Ultrasound Blood Flow Imaging，Based on a Full Evaluation of the Doppler Equation by Means of a Two-Dimensional Autocorrelation Approach[J].Ultrasonics Ferroelectrics and Frequency Control，IEEE Transactions on，1995，42（4）:672-688.

[7]Loupas T，Gill R W.Multifrequency Doppler:Improving the Quality of Spectral Estimation by Making full use of the Information Present in the Backscattered RF Echoes[J].Ultrasonics，F(xiàn)erroelectrics and Frequency Control，IEEE Transactions on，1994，41（4）: 522-531.

[8]Shaoguo Cui and Dong C.Liu，Noise Reduction for Ultrasonic Elastography Using Transmit-Side Frequency Compounding:A Preliminary Study，IEEE Trans Ultrason Ferroelec Freq Control，2011，58，3.

ARFI;Elasticity Imaging;2D Autocorrelator Algorithm;Sampling Time Window

An Improved and Estimated Displacement Algorithm for Acoustic Radiation Force Impulse Imaging

YOU Xin
（Medical Imaging Lab，Sichuan University，Chengdu 610065）

1007-1423（2016）02-0068-05

10.3969/j.issn.1007-1423.2016.28.016

針對超聲聲輻射力成像中使用的傳統(tǒng)Loupas算法計(jì)算誤差較大，成像效果不佳的問題，提出一種改進(jìn)方法。該方法在Loupas算法基礎(chǔ)上，提出通過改變時(shí)間取樣窗口大小，并利用二維自相關(guān)計(jì)算平均速度，進(jìn)而得到更為準(zhǔn)確的位移計(jì)算結(jié)果。實(shí)驗(yàn)結(jié)果與傳統(tǒng)Loupas算法相比，提高彈性圖像對比度噪聲比，有更好的成像效果。

超聲輻射力；彈性成像；二維自相關(guān)算法；時(shí)間取樣窗口

With the purpose of solving the issue that there are errors in the traditional Loupas algorithm used in ARFI and the effects of imaging are not as good as expected,proposes an improved method to solve this proplem.Based on Loupas algorithm,it is presented that we can get more accurate computing result of the displacement by changing the size of the sampling time windows and calculating the average speed according to the 2D autocorrelator.As a result,not only improves Elastographic Contrast-to-Noise Ratio,but also the image quality.