譚勇 王叢知

【摘要】位移估計是彈性成像的第一步，它決定著整個成像的時間以及質(zhì)量，目前位移估計算法有兩大類：傳統(tǒng)的歸一化互相關(guān)和相位檢測算法。相位檢測法由于計算速度快被廣泛運用，通過對相位檢測算法的深入研究提出一種新的基于原始RF信號的Doppler位移估計算法，該方法只需一次積和運算以及不用構(gòu)建IQ信號，大大縮短了計算時間。通過生物仿體研究證明，該方法比傳統(tǒng)相位檢測法和Doppler速度估計法更快，從性能上看，RF-Doppler法得到應(yīng)變圖像質(zhì)量與相位檢測法相似，但是速度卻快了接近100倍。

【關(guān)鍵詞】聲輻射力成像;位移估計;GPU;并行算法

1.引言

聲輻射力彈性成像（ARFI）是一種非侵害性的影像技術(shù)，將局部組織的彈性模量數(shù)值映射成彩色編碼的圖像信息[1]。它的原理是，首先利用聚焦超聲波束產(chǎn)生的聲輻射力在聚焦位置造成組織偏移，并形成向側(cè)向傳播的剪切波，然后由高幀頻的探測超聲波束追蹤該剪切波的傳播，從射頻（RF）數(shù)據(jù)中計算出剪切波引起的組織位移變化，進(jìn)而計算出剪切波傳播速度，最后根據(jù)剪切波速度計算出組織的彈性模量值。它可以提供一些以生物軟組織（如乳房、肝臟、腎、肌肉、眼睛、胰腺等）的彈性病變?yōu)榛A(chǔ)的診斷信息，具有重要的臨床價值和廣闊的應(yīng)用前景。

目前超聲彈性成像所使用的組織位移估計算法主要包括兩大類：傳統(tǒng)的歸一化互相關(guān)算法和基于IQ信號的Doppler相位檢測算法。傳統(tǒng)的歸一化互相關(guān)計算需要在一定范圍內(nèi)確定互相關(guān)系數(shù)的極大值位置，但直接的搜索結(jié)果只能是采樣周期的整數(shù)倍，因此，通常需要對互相關(guān)函數(shù)極大值點的前后兩點之間進(jìn)行拋物線插值，然后在插值范圍內(nèi)重新尋找更為精確的極大值點位置。插值倍數(shù)越大，搜索極大值點位置的精度就越高，但耗時也越長。基于IQ信號的Doppler相位檢測算法的原理是，對原始的超聲射頻（RF）信號進(jìn)行IQ分解后（其中一個步驟是讓信號通過低通濾波器去除載波成分），將所得到的正交IQ信號進(jìn)行互相關(guān)計算，則其互相關(guān)系數(shù)極大值點的相位，再除以相鄰兩幀數(shù)據(jù)間的時間間隔，對應(yīng)于組織的相對平均移動速度。據(jù)此便可進(jìn)一步計算出精確的組織相對位移距離。該方法不需要進(jìn)行插值，減少了計算時間并提高了計算精度，已被廣泛采用。

要解決計算量和計算速度之間的矛盾，可以通過使用GPU進(jìn)行并行計算加速實現(xiàn)。NVIDIA公司于2006年推出的“統(tǒng)一編程架構(gòu)”使得擁有NVIDIA GPU的用戶能輕松使用GPU進(jìn)行編程[2]。GPU具有比CPU更多的運算單元，更高的帶寬以及浮點運算能力，因此GPU特別適合高密集型的運算。目前，GPU已經(jīng)被廣泛用在計算機(jī)斷層掃描、核磁共振等成像上并取得了顯著效果，但在超聲成像，特別是聲輻射力超聲彈性成像上的應(yīng)用還很少見于報道。

本文提出一種應(yīng)用于ARFI中的特別適合于利用GPU進(jìn)行加速計算的快速組織位移估計算法，并通過實驗證明該算法有較快的計算速度以及不錯的應(yīng)變圖像信噪比。

2.組織位移估計算法介紹

目前在商用設(shè)備和科學(xué)研究中被廣泛采用的組織位移估計算法主要包括基于RF信號的改進(jìn)互相關(guān)算法和基于IQ信號的Doppler相移檢測算法，本章將分別對現(xiàn)有算法和我們提出的新算法進(jìn)行介紹。

2.1 基于RF信號的改進(jìn)互相關(guān)算法

傳統(tǒng)的互相關(guān)算法需要進(jìn)行大量的錯位積和計算和插值計算以提高定位精度。計算量大且組織位移的結(jié)果不夠準(zhǔn)確。為解決這些問題，Cabot在1981年提出了一種通過構(gòu)造互相關(guān)系數(shù)的解析信號并對其相位進(jìn)行分析來提高定位精度的新算法[3]。首先對原始的RF數(shù)據(jù)分段并對相鄰兩幀中的對應(yīng)數(shù)據(jù)段進(jìn)行互相關(guān)計算，再對所得的互相關(guān)系數(shù)進(jìn)行希爾伯特變換，構(gòu)造其解析信號。該解析信號在互相關(guān)系數(shù)極大值附近的相位過零點位置，被證明精確地對應(yīng)于互相關(guān)系數(shù)的真實極大值點的位置。該算法具體步驟：

1）先對回波信號分段，本文每段長為100個數(shù)據(jù)點。

2）對每段進(jìn)行時域互相關(guān)：

（1）

式中：m窗口索引;s為幀序列;L為窗長;n是延遲間隔。

3）對互相關(guān)函數(shù)構(gòu)造解析信號。

（2）

（3）

h（n）是希爾伯特變換函數(shù)[4]，是的希爾伯特變換，為互相關(guān)函數(shù)的解析信號。

4）尋找相位過零點：

（4）

使用牛頓迭代法[5]尋找到τ使得，即該延時精確對應(yīng)時域互相關(guān)函數(shù)最大值的延時。

該方法的弊端除了前面提到的多次錯位運算外，運用牛頓迭代并不適合在GPU上并行運算，因為后一次的運算需要前一次運算得出的結(jié)果。

2.2 基于IQ信號的Doppler相移檢測算法

此算法由Kasai于1985年提出，最初用于Doppler血流速度估計[6]，后來又被用于組織位移估計。該算法首先對原始的RF信號進(jìn)行IQ分解，在對得到的IQ信號進(jìn)行分段，并計算相鄰兩幀中對應(yīng)數(shù)據(jù)段之間的點積和，并計算其相對于中心頻率的平均相位移動，從而計算出組織的相對位移，其公式如（5）。c是波速，f是RF信號中心頻率。

由于ARFI數(shù)據(jù)中兩幀相鄰信號間的位移幅度很小，所以運用該方法時基本不會出現(xiàn)相位卷繞的問題。但由于IQ分解所引起的系統(tǒng)復(fù)雜度上升和信號信噪比下降的問題依然存在，并最終影響到組織位移估計結(jié)果的精度。

（5）

2.3 特別適合于GPU加速的快速組織位移估計新算法

以上兩種方法雖然比傳統(tǒng)歸一化互相關(guān)算法高效和精確，但是第一種算法需要多次錯位運算和迭代運算，計算量仍舊很大，第二種構(gòu)建IQ信號需要專門的電路，若用軟件實現(xiàn)，計算量大而且繁瑣。通過分析發(fā)現(xiàn)，構(gòu)建解析信號與IQ信號雖然從構(gòu)建方法上看完全不同，但利用其相位進(jìn)行位移或者速度估計的方法從本質(zhì)上說是相似的。用軟件實現(xiàn)希爾伯特變換來構(gòu)建解析信號比用軟件構(gòu)建IQ信號簡單的多，而且可以利用傅立葉正、反變換在頻域上很方便地實現(xiàn)，不再需要完成IQ分解的硬件電路，簡化了系統(tǒng)設(shè)計，也不再需要低通濾波，保留了完整的超聲回波信息，具有很大優(yōu)勢。所以我們提出，綜合上述兩種方法的優(yōu)勢，直接利用原始的RF數(shù)據(jù)進(jìn)行計算，直接對整幀的RF數(shù)據(jù)進(jìn)行解析信號的構(gòu)建，然后再進(jìn)行分段等后續(xù)計算，這一步驟可以大大提高利用GPU進(jìn)行并行計算加速的效果。算法的大致步驟是：首先對原始RF信號構(gòu)建希爾伯特變換，分段后，對應(yīng)數(shù)據(jù)段內(nèi)的對應(yīng)數(shù)據(jù)點進(jìn)行一次求積和，公式如下：

（6）

（7）

（8）

是相鄰兩幀解析信號對應(yīng)窗口的積和值，m是窗口索引，L是窗長，*是共軛標(biāo)志，f是RF信號中心頻率，c代表波速。

3.基于GPU的快速組織位移估計算法實

3.1 基于RF信號的改進(jìn)互先關(guān)算法

一條A-Line信號總共有K幀每幀有N個采樣點，為了達(dá)到并行運算，我們分配K-1個線程塊，每個線程塊有M個線程（M是延遲數(shù)量），每個線程塊處理S幀和（S+1）幀信號的互相關(guān)運算。

由于做互相關(guān)時窗口覆蓋率達(dá)到98%，若使用傳統(tǒng)的互相關(guān)運算，即S和（S-1）的各個對應(yīng)窗口進(jìn)行相關(guān)運算，這樣會導(dǎo)致大量的重復(fù)運算，因此我們使用滑窗的辦法。首先每個線程計算第一個窗口的所有數(shù)據(jù)，計算完后，進(jìn)行下一個窗口的計算時只需計算新增的點積和，然后再減去被移出的點的點積和，這樣至少可以節(jié)約90%的運算量。有線程存取共享內(nèi)存比全局內(nèi)存有更高的訪問速度，所以在做互相關(guān)之前先將信號存入共享內(nèi)存，經(jīng)驗證能提高15%的運算速度。

3.2 基于IQ信號的Doppler相移檢測算法

為了要將原始信號的每一幀構(gòu)造IQ信號，需要將原始信號拷入兩個矩陣I和Q。在這兩個矩陣?yán)锓謩e對其進(jìn)行I轉(zhuǎn)換和Q轉(zhuǎn)換。在做互相關(guān)時，需要K個線程塊，每個線程塊有M個線程，K代表需要做互相關(guān)的信號，M代表窗口數(shù)量。K線程塊的序號0、1、2…，分別對應(yīng)I、Q矩陣的第0、1、2…信號幀，即線程塊0處理I信號和Q信號的第0幀信號，線程塊1處理I信號和Q信號的第1幀信號…。線程序號對應(yīng)每個窗口序號，相當(dāng)于公式（8）中的m，每個線程負(fù)責(zé)計算該窗口內(nèi)所有采樣點的點積和。為了能提高效率，同樣先將要進(jìn)行互相關(guān)運算的I、Q信號存入共享內(nèi)存，當(dāng)所有點積和運算完后，使用庫函數(shù)里的atan2f對該和求出相位值，并將其寫回GPU全局內(nèi)存。

3.3 特別適合于GPU加速的快速組織位移估計新算

該方法在GPU上的實現(xiàn)相當(dāng)簡單。首先我們利用CUDA自帶的CUFFT庫中的傅立葉變換函數(shù)對整幀信號做希爾伯特變換。如果依然采用傳統(tǒng)的先分段再構(gòu)造解析信號的方式，希爾伯特變換會反復(fù)多次調(diào)用cufftExecC2C函數(shù)，且很多數(shù)據(jù)點會被反復(fù)多次計算（因為步長僅為2%，相鄰窗口之間有98%的數(shù)據(jù)點重疊），效率低、耗時長。根據(jù)Nvidia的CUFFT庫函數(shù)使用指南，由于其采用了高效的FFTW快速傅立葉變換算法，若需要進(jìn)行變換的數(shù)據(jù)長度固定，可以用cufftPlan1d函數(shù)進(jìn)行預(yù)優(yōu)化處理，并且信號長度越長，計算速度提升的幅度就越大。因此，采用對整幀信號進(jìn)行解析信號構(gòu)造的方式，會極大地提高計算效率。同時，經(jīng)過公式推導(dǎo)和實驗驗證，我們確認(rèn)，對于ARFI中所使用的超聲RF回波信號來說，這種調(diào)換計算順序的方式對最終的組織位移估計結(jié)果的影響極小，可以忽略不計。

4.實驗與討論

4.1 實驗

本文使用NIVIDIA GT630M，它有96個CUDA核和1GB顯存，編程平臺為Visual Studio2010。實驗數(shù)據(jù)采集自我們自制的彈性仿體，三個算法分別對兩個區(qū)域進(jìn)行檢測。區(qū)域1是一個較軟與較硬區(qū)域位置的交界面的位置，區(qū)域（2）是一個硬度分布均勻且較軟的區(qū)域。以上兩個區(qū)域大小為10*10cm，區(qū)域內(nèi)每個推動位置的大小為1*1mm，即每一排每一列的推動位置數(shù)量為10。為了比較各個算法成像的差異，選擇聲輻射力彈性成像后的圖像質(zhì)量進(jìn)行評價。評價彈性圖質(zhì)量的方式主要有信噪比（SNRe）[7]和對比度噪聲比（CNRe）：

S是所測區(qū)域彈性系數(shù)的平均值，δe是所測區(qū)域彈性系數(shù)的標(biāo)準(zhǔn)差。st、sb分別是硬度區(qū)和背景區(qū)的均值，δt、δb分別為硬度區(qū)和背景區(qū)的標(biāo)準(zhǔn)差。信噪比的計算方式是選取硬度均勻的區(qū)域，如圖2-2，對比度噪聲比的計算需選擇包含不同硬度值的區(qū)域，如圖2-1。信噪比值越大說明算法抑制噪聲的能力越強(qiáng)，對比度噪聲比值越大說明算法彈性檢測能力越強(qiáng)。

在編程平臺上實現(xiàn)算法時我們采用了在CPU上單獨編程和基于CPU和GPU的混合編程，并且對基于RF信號的改進(jìn)互相關(guān)算法采取了兩種實現(xiàn)手段：

（1）先構(gòu)造解析信號，再分段互相關(guān)。

（2）先分段互相關(guān)，再由互相關(guān)系數(shù)構(gòu)造解析信號。

4.2 結(jié)果與討論

表1 各算法的時間對比

算法 時間

NCC_1

NCC_2

New

Doppler CPU/ms

1701

2920

141

5290 GPU/ms

266

1619

2.4

1630

表1給出了各算法在不同窗長下的運行時間。其中NCC_1代表先構(gòu)造解析信號，再分段互相關(guān)的實現(xiàn)方式，NCC_2代表先分段互相關(guān)，最后用互相關(guān)系數(shù)構(gòu)造解析信號， Doppler代表基于IQ信號的位移計算時間，New代表基于RF信號的我們所提出的新算法的計算時間。從上表可以發(fā)現(xiàn)算法在GPU上的實現(xiàn)方式對計算速度有很大影響，先對原始RF信號進(jìn)行希爾伯特變換再分段互相關(guān)比先分段再互相關(guān)最后進(jìn)行希爾伯特變化速度有明顯提升，大概為6倍。Doppler無論在CPU上還是在GPU上，速度都比其它算法慢很多。本文提出的新方法，相對CPU速度提升了近58倍，并且無論在CPU還是在GPU上速度均是最快的。

圖2-1、2-2、2-3分別是NCC_1、New、Doppler算法對圖2中區(qū)域1成像（硬度分布不均勻）的比較;圖2-1、2-2、2-3是對區(qū)域2（彈性分布穩(wěn)定且較軟）的成像圖。圖3給出了各算法在圖1中區(qū)域各個推動位置測得的具體彈性模量值。從分布來看NCC_1與New比較吻合，而Doppler差別較明顯，這主要是由于我們程序中使用了低通濾波器，但是很難達(dá)到完全理想的濾波效果，這就造成信號波形和相位的偏移和扭曲，進(jìn)而影響估算Doppler頻率和組織位移的準(zhǔn)確程度。NCC_1、New、Doppler算法對區(qū)域2所得彈性圖像的SNRe分別為：8.67、7.5、4.26。三種算法的對比度噪聲比分別為：0.991、0.953、0.684。

從信噪比的結(jié)果來看，NCC和New算法所得應(yīng)變圖像都有較高的圖像信噪比，而Doppler算法所得應(yīng)變圖像信噪比略低。并且前兩種算法均比Doppler算法有更高的對比度噪聲比值，即NCC_1和New有更好的彈性檢測能力。

圖3 各算法對區(qū)域1各推動位置的彈性值

5.結(jié)論

本文通過對傳統(tǒng)超聲彈性成像位移算法的研究以及結(jié)合聲輻射力位移的特點，提出了一種特別適合于GPU加速的快速組織位移估計新算法。并介紹了基于解析信號的互相關(guān)相位改進(jìn)算法、Kasais的Doppler運動估計算法和本文提出的算法在GPU上實現(xiàn)的一些技巧。最后通過實驗表明本文提出的新算法能取得與常用算法相似的應(yīng)變圖像信噪比以及對比度噪聲比，但是速度卻提高了近100倍。除此以外本文提出的新算法并不局限于被應(yīng)用在基于聲輻射力的二維超聲彈性成像系統(tǒng)，也可應(yīng)用于任何需要利用超聲回波射頻信號進(jìn)行組織位移估計的設(shè)備或方法中，比如準(zhǔn)靜態(tài)超聲彈性成像（quasi-static elastography）、瞬態(tài)超聲彈性成像（transient elastography）等其他超聲彈性成像方法，或者超聲溫度成像（temperature imaging）、超聲熱應(yīng)變成像（thermal strain imaging）等方法。

參考文獻(xiàn)

[1]彭博，諶勇，流動權(quán).基于GPU的超聲彈性成像并行研究[J]，光電工程，2013，40（5）：97-105.

[2]CUDA Tookit Document.NVIDIA.

[3]RICHARD C.CABOT.A Note on the Application of the Hilbert Transform to time delay estimation[J].IEEE TRANSICTION ACOUSTICS，SPEECH，AND SIGNAL PROCESSINF，1981，29（3）：607-609.

[4]SCHWARTZ M，BENNRTT W R，STEIN S.Communicaion Systems and Techniques[M].New-York：McGraw-Hill，1965.

[5]PESAVENTO A，PERREY C，KRUEGER M，et at.A Time-Efficient and Accurate Strain Estimation Concept for Ultrasonic Elastography Using Iterative Phase Zero Estimation[J].IEEE TRANSAXTIONS ON ULTRASONICS，F(xiàn)ERROELECTICS，AND FREQUENCY CONTROL，1999，46（5）：1057-1067.

[6]DE C，NAMEKAWA K，KOYANO K，et al.Real-Time two-dimensionalblood flow imaging using autocorrelation technique[J].IEEETRANSACTION SONICS ULTRASON，1985，42（4）：458-463.

[7]萬明習(xí)等.生物醫(yī)學(xué)超聲學(xué)[M].北京：科學(xué)出版社，2010：

483-486.

作者簡介：

譚勇（1985—），男，成都理工大學(xué)碩士研究生在讀，研究方向：嵌入式系統(tǒng)研究。

通訊作者：王叢知（1977—），男，博士，中國科學(xué)院深圳先進(jìn)技術(shù)研究院助理研究員，研究方向：生物醫(yī)學(xué)信號處理。