一種基于GPU的聲輻射力成像組織位移估計算法

譚勇 王叢知

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

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

1.引言

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

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

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

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

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

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

2.1 基于RF信號的改進互相關算法

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

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

2）對每段進行時域互相關：

（1）

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

3）對互相關函數(shù)構造解析信號。

（2）

（3）

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

4）尋找相位過零點：

（4）

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

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

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

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

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

（5）

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

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

（6）

（7）

（8）

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

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

3.1 基于RF信號的改進互先關算法

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

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

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

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

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

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

4.實驗與討論

4.1 實驗

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

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

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

（1）先構造解析信號，再分段互相關。

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

4.2 結果與討論

表1 各算法的時間對比

算法 時間

NCC_1

NCC_2

New

Doppler CPU/ms

1701

2920

141

5290 GPU/ms

266

1619

2.4

1630

表1給出了各算法在不同窗長下的運行時間。其中NCC_1代表先構造解析信號，再分段互相關的實現(xiàn)方式，NCC_2代表先分段互相關，最后用互相關系數(shù)構造解析信號， Doppler代表基于IQ信號的位移計算時間，New代表基于RF信號的我們所提出的新算法的計算時間。從上表可以發(fā)現(xiàn)算法在GPU上的實現(xiàn)方式對計算速度有很大影響，先對原始RF信號進行希爾伯特變換再分段互相關比先分段再互相關最后進行希爾伯特變化速度有明顯提升，大概為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差別較明顯，這主要是由于我們程序中使用了低通濾波器，但是很難達到完全理想的濾波效果，這就造成信號波形和相位的偏移和扭曲，進而影響估算Doppler頻率和組織位移的準確程度。NCC_1、New、Doppler算法對區(qū)域2所得彈性圖像的SNRe分別為：8.67、7.5、4.26。三種算法的對比度噪聲比分別為：0.991、0.953、0.684。

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

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

5.結論

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

參考文獻

[1]彭博，諶勇，流動權.基于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]萬明習等.生物醫(yī)學超聲學[M].北京：科學出版社，2010：

483-486.

作者簡介：

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

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